Evaluating normal inverse when probability is close to zero

Question

Mike Pennell 2018-11-17

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/430466-evaluating-normal-inverse-when-probability-is-close-to-zero

评论： Mike Pennell 2018-11-19

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?

2 个评论
显示无隐藏无

Matt J 2018-11-19

编辑：Matt J 2018-11-19

在 MATLAB Online 中打开

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 2018-11-19

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 2018-11-17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/430466-evaluating-normal-inverse-when-probability-is-close-to-zero#answer_347574

编辑：John D'Errico 2018-11-17

在 MATLAB Online 中打开

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...

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

David Goodmanson 2018-11-18

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/430466-evaluating-normal-inverse-when-probability-is-close-to-zero#answer_347585

编辑：David Goodmanson 2018-11-18

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

John D'Errico 2018-11-18

编辑：John D'Errico 2018-11-18

+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 2018-11-19