Evaluating normal inverse when probability is close to zero

I am trying to find the inverse of the normal cdf for a very small probability: exp(-750). In R there is an option to input probabilities on the log scale but that option doesn't exist in the norminv function in Matlab and I get "-Inf" when I try to use p=exp(-750). Is there another function that allows this or is there a work around?

[Matt J]

Out of curiousity, why do you think you need a better solution than -Inf? It gives excellent numerical agreement with your target when run through normcdf:

>> isequal(  exp(-750), normcdf(-inf) )

ans =

  logical

   1

[Mike Pennell]

I need to sample from a normal distribution truncated over positive values in an MCMC procedure. Every once in a while, you can get a negative mean and small sd which can cause the program to crash.

Answer 1

[John D'Errico]

You need to understand that this probability is impossible to express in double precision. However, there are some things you can do. Lets see...

First, we have a very tiny value there.

vpa(exp(sym(-750)))
ans =
1.901684963475006439995456236734e-326

The simple answer is to just use syms. For example, since the inverse normal is directly related to the error function, we see that syms has no problems for very small arguments in erfcinv.

So, we know that the normal CDF, phi(X) is related to erf as:

phi(X) = (1 + erf(X/sqrt(2)))/2

Or, to erfc, the relation is

phi(X) = erfc(-X/sqrt(2)))/2

So, now we can do the necessary computation as:

double(-erfcinv(2*exp(sym(-750)))*sqrt(2))
ans =
          -38.611574423848

Ok, do you want a non-symbolic toolbox solution? Hmm. I'd probably need to dive into some asymptotics for that, so probably a little light reading in Abramowitz & Stegun...

Answer 2

[David Goodmanson]

Hi Mike,

Of course there is always the question if the normal pdf could still apply that far out in the tails, but that is assumed as true.

While John was doing things the sensible easy way, I was off into Abramowitz and Stegun.

This calculation is for a standard normal. The cdf is an integral from –inf to -x for large positive x. To simplify things, look at the equivalent calculation, the integral from +x to infinity. An asymptotic approximation to that is

P = 1/(sqrt(2*pi))*exp(-x^2/2)*A
where
A = 1/x - 1/x^3 + 3/x^5 -3*5/x^7 +3*5*7/x^9 ... 

taking logs gives

log(P) = -x^2/2 - log(sqrt(2*pi)) + log(A). 
Let
mlogP = -log(P) = 750
then
x = sqrt(2*(mlogP - logC + log(A)));

Have to solve for x where A is a function of x, but simple iteration works really quickly. In the 750 case it takes six terms in A for the answer to settle out in double precision:

mlogP = 750;
n = 6;
c = cumprod(-(1:2:2*n-1));
c = [1 c(1:end-1)]                  % asymptotic expansion coefficients 
logC = log(sqrt(2*pi));
x = sqrt(2*(mlogP - logC));          % initial guess
for k = 1:10
  A = sum(c.*(x.^-(1:2:(2*n-1))));           
  x = sqrt(2*(mlogP - logC + log(A)));
end
format long
x = -x;                             % back to the original cdf
x  
double(-erfcinv(2*exp(sym(-750)))*sqrt(2))
ans = -38.611574423848019
ans = -38.611574423848019

[John D'Errico]

+1. Yes. as I said, this is where I would have gone if the symbolic solution was inadequate. Thanks to David for doing the work that I was too lazy to do at the time.

[Mike Pennell]

David and John,

Thank you for your help. I will try David's code.

Mike