Evaluating differential entropies with Matlab: NaN issue

Question

With Matlab I am trying to evaluate differential entropies. These are integrals like

$$\int_\mathbb{R} p(x) \log (p(x)) \mathrm{d}x$$

where $p(x)$ is a probability density function. My $p(x)$ is computed by an external function. I tried to define

g = @(x) x.*log(x);

and compute the entropy using

integral(@(x) g(p(x)),0,inf)

(my $p(x)$ is defined on $[0;+\infty)$). But this fails due to NaN issues. However, even if I take the precaution to define g as an external function and set NaN outputs to zero by a conditional statement, the computation fails. Why?

Consider the following simple example:

I define g as an external function to avoid the $0 \log 0$ issue as follows:

function y = g(x)

if isnan(x.*log(x))
        y = 0;
else
        y = x.*log(x);
end

Now let's try to compute the entropy integral for $p(x) = e^{-x}$, which should obviously return $-1$. But it doesn't:

>> p = @(x) exp(-x);
>> integral(@(x) g(p(x)),0,inf)
Warning: Infinite or Not-a-Number value encountered. 
> In funfun/private/integralCalc>iterateScalarValued at 349
  In funfun/private/integralCalc>vadapt at 133
  In funfun/private/integralCalc at 84
  In integral at 89

What is Matlab's problem?

Answer 1

The reason is because you've not defined your function g correctly in Matlab. The way it's written will not trap NaN cases. As per the documentation for the integral function:

For scalar-valued problems, the function y = fun(x) must accept a vector argument, x, and return a vector result, y. This generally means that fun must use array operators instead of matrix operators. For example, use .* (times) rather than * (mtimes). If you set the 'ArrayValued' option to true, then fun must accept a scalar and return an array of fixed size.

The integral function is vectorized and passes multiple values to your integrand function on each iteration. You've correctly vectorized x.*log(x), but your if statement that checks for NaN will only be triggered if the first element happens to be NaN. Using so-called logical indexing, you can rewrite g concisely as

function y = g(x)
y = x.*log(x);
y(isnan(y)) = 0;

Then the following code

p = @(x)exp(-x);
integral(@(x)g(p(x)),0,Inf)

now returns -1.000000000000000.

Also, though it may not be numerically scaled as well, you could use the mathematically equivalent g = @(x)log(x.^x); instead, which cannot evaluate to NaN.

Answer 2

MATLAB does not simplify your integrand - it treats the integrand exactly as specified. So in your first example, MATLAB tries to build a table of the function

$$ \exp(-x)\exp(-x)\log(\exp(-x)) $$

rather than $$ -x \exp(-x)\exp(-x) $$

as you would expect.

The first integrand will cause trouble as $x$ becomes very large due to floating point issues. For instance, if you get rid of all the function handles and just try

integral(@(x) exp(-2.*x).*log(exp(-x)),0, inf)

you get a Nan and the same message.

However, considering how rapidly your integral dies off,

integral(@(x) exp(-2.*x).*log(exp(-x)),0, 20)

does the job just fine. In fact, using an upper limit of 15 returns a double precision integral of -0.249999999999275, which is far more accurate than anything you might need.

Edited to add: If you want more technical details, for values of x more negative than a certain number, you cause an underflow - a number whose exponent in base 2 is smaller than -1022 cannot be represented in the usual exponent-mantissa format. Underflows can be rounded off to 0 automatically (this is probably what MATLAB does), so numbers like log(exp(-1400)) evaluate to log(0), throwing the NaN.