how to improve the precision when using the cumulative density function

Dear fellows,

I am calulating the log likelihood of bivariate normal distribution using the the buillt in function mvncdf. But I always get a lot of -Inf. Is there any solution to get the precise answer? I will show you an example to make my question clear:)

Here is the code

mu = [181.665319288082,145.332255430466]; SIGMA = [2 1; 1 2]; hatmu=[25,23]; X1=[hatmu(1)-0.5 hatmu(2)-0.5]; X2=[hatmu(1)+0.5 hatmu(2)+0.5]; X = [X1, X2]; p = log(mvncdf(X1,X2,mu,SIGMA))

When I run the above code, the result I get is p=-Inf.(the answer for mvncdf(X1,X2,mu,SIGMA)) is 0 so it is not surpriising to get p=-Inf)

Appreciate any help:)

Answer 1

[Tom Lane]

This may be too elaborate an answer, but it's all I have.

You are computing the probability of falling into a rectangle. If the lowest value of the pdf in this rectangle is L, and the highest value is H, then the answer will be bounded by A*L < answer < A*H where A is the area of the rectangle. To get a rough approximation, you could compute A*T where T is a typical value of the pdf in the rectangle.

So it remains to compute the pdf. If you try using mvnpdf you'll get zero. Edit that file, and at the end change the "y = exp(...)" stuff to remove the exponential and return the result on the log scale. Save the file under your own name and run it. When I try that, I get these log pdf values at the four corners of your rectangle:

   1.0e+03 *
   -6.8303
   -6.7373
   -6.8011
   -6.7668

The rectangle has area 4, but that doesn't matter much -- your answer is around exp(-6800). That's really small. So small that I suspect it's not going to do you much good.

Just to be clear, I don't recommend this as a general procedure, but it's one way to try to get a rough idea of the value.

[xueqi]

Hi sorry for late reply. How can you let matlab to give you the result in the form e+03? Besides how can you get the result exp(-6800)? I know it is really small number, but I need call this function in the optimiazation (fmincon), so when the result is so small then fmincon will keeps searching to get the optimal answer.

[Tom Lane]

Sorry, I don't understand. In my answer I described how I arrived at the number exp(-6800) as a rough estimate of the order of magnitude. The 1.0e+03 in my display is just the way MATLAB displays the vector when all elements are small.