Hi All,
Thanks to those who have responded to my question about how to get the
inverse of the standard normal cdf, i.e. the inverse of the phi() function.
Below is the summary of the reponses.
Larry Gould suggests a Hastings approximation, e.g.,
FUNCTION InvNormal(p : DOUBLE) : DOUBLE;
CONST zero = 0.0; one = 1.0; half = 0.5; split = 0.42;
a : ARRAY[1..4] OF DOUBLE
= (2.50662823884, -18.61500062529, 41.39119773534,
-25.44106049637);
b : ARRAY[1..4] OF DOUBLE
= (-8.47351093090, 23.08336743743, -21.06224101826, 3.13082909833);
c : ARRAY[1..4] OF DOUBLE
= (-2.78718931138, -2.29796479134, 4.85014127135, 2.32121276858);
d1 : DOUBLE = 3.54388924762; d2 : DOUBLE = 1.63706781897;
VAR q,r : DOUBLE;
BEGIN
q := p - half;
IF (ABS(q)<=split) THEN
BEGIN
r := SQR(q);
InvNormal := q*(((a[4]*r + a[3])*r + a[2])*r + a[1])/
((((b[4]*r + b[3])*r + b[2])*r + b[1])*r + one);
END
ELSE BEGIN
IF q > zero THEN r := one - p ELSE r := p;
IF r<= zero THEN InvNormal := -99.0
ELSE BEGIN
r := SQRT(-LN(r));
r := (((c[4]*r + c[3])*r + c[2])*r + c[1])/((d2*r + d1)*r + one);
IF q<zero THEN InvNormal := -r ELSE InvNormal := r;
END;
END;
END;
George Woodworth suggests a pretty good rational-polynomial approximation in
Abramowitz and Stegun (formula 26.2.23, Page 933) which is accurate
to 4.5 x 10(-4) over the entire domain of the function. Part of the
handbook is available on the following website
http://jove.prohosting.com/~skripty/frameindex.htm
Best,
Yuanping
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|