Evaluating normal inverse when probability is close to zero
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (2)
John D'Errico
am 17 Nov. 2018
Bearbeitet: John D'Errico
am 17 Nov. 2018
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 Kommentare
David Goodmanson
am 18 Nov. 2018
Bearbeitet: David Goodmanson
am 18 Nov. 2018
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 Kommentare
John D'Errico
am 18 Nov. 2018
Bearbeitet: John D'Errico
am 18 Nov. 2018
+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.
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!