/*
The following code was taken from the one for "Simulation of Multivariate
Normal Orthant Probabilities" written by V.A. Hajivassiliou,
D.L. McFadden and P.A. Ruud.


Author: Yoshihisa Suzuki, Univ. of Rochester
Date: May 1995
Code provided without guarantee for public non-commercial use.
*/

/*
** CDFINVN --   computes the inverse normal cdf. The algorithm used is taken
**              from Kennedy and Gentle, STATISTICAL COMPUTING, p. 95.
**
** vah 22:57:35.78, Sun, Jul 23, 1989
** Taken intact from V149b MAXPROC3.ARC
**
**              xp=cdfinvn(p);
**
**             INPUT
**              p -- nx1 vector of probabilities (that is, 0 < p < 1)
**             OUTPUT
**             xp -- the corresponding vector of values such that: cdfn(xp)=p;
**
**              This algorithm seems to be very accurate.  It is
**              rated to be accurate to 1e-7.  However, it appears
**              to be substantially more accurate -- approximately
**              1e-10.
**
*/

proc 1=cdfinvn(p);
local lim, p0, p1, p2, p3, p4, q0, q1, q2, q3, q4, maskgt, maskeq, sgn,
y, xp, pn;

@ Constants @

lim =  1e-20;

p0  = - 0.322232431088;             q0  =   0.0993484626060;
p1  = - 1.0;                        q1  =   0.588581570495;
p2  = - 0.342242088547;             q2  =   0.531103462366;
p3  = - 0.0204231210245;            q3  =   0.103537752850;
p4  = - 0.453642210148*1e-4;        q4  =   0.38560700634*1e-2;


@ Main body of code @

 @ Create masks for handling p > .5 and p >= .5 @

maskgt=(p .> 0.5);
maskeq=(p .ne 0.5);
sgn=missrv(miss(maskgt,0),-1);

 @ Convert p > .5 --> 1-p @

pn=(maskgt-p).*sgn;   clear maskgt;

 @ Computation of function for p < 0.5 @

y=sqrt(-2*ln(pn));     clear pn;
xp=y + ((((y*p4 + p3).*y + p2).*y + p1).*y + p0)./
       ((((y*q4 + q3).*y + q2).*y + q1).*y + q0); clear y;

 @ Convert results for p > .5 and p = .5 @

xp=(xp.*sgn).*maskeq;
clear sgn, maskeq;
retp(xp);
endp;
