/****************** Inverse CDFN *****************************
Author: Pat Mackin sagcorp@sagcorp.com
Date: 11/11/96
Provided without guarantees for public, non-commercial use.

Purpose: computes the inverse function of CDFN using
the approximation method based on Abramowitz and
Stegun (1970), p. 933.
************************************************************/

proc findz(rate);
    local probnum,c,d,t,x,x1,x2,probtog,badrate,pmat,zerotog;

    c = { 2.515517, 0.802853, 0.010328 };
    d = { 1.432788, 0.189269, 0.001308 };

    badrate = ((rate .< 0) .or (rate .> 1));
    rate = (rate .* (badrate .== 0)) + (0 .* badrate);
    badrate = miss(badrate,1);
    probtog = (rate .<= 0.5).*(-1) + (rate .> 0.5);
    probnum = 1 - rate;
    pmat = minc(rate'|probnum');
    zerotog = (pmat .== 0);
    pmat = pmat + .01 .* zerotog;
    t = sqrt(ln(1./(pmat^2)));
    x1 = (c[1]+c[2].*t+c[3].*(t^2));
    x2 = (1+d[1].*t+d[2].*(t^2)+d[3].*(t^3));
    x = probtog .* ((t-(x1./x2)) .* (1 - zerotog) + (7 .* zerotog));
    x = x + badrate;
    retp(x);
endp;


/****************** Inverse CDFN *****************************
Author: Joachim Inkmann          Joachim.Inkmann@uni-konstanz.de
	University of Konstanz
	SFB 178
	P.O. Box 5560 D124
	78434 Konstanz
	Germany
	Phone: +49.7531.88.3753
	Fax:   +49.7531.88.4450
Date: 11/11/96
Provided without guarantees for public, non-commercial use.

Purpose: computes the inverse function of CDFN using
the the more precise 4th order approximation by Odeh and Evans (1974)
instead of the 3rd order polynomial approximation to the inverse
standard normal cdf introduced by Hastings (1955) and reprinted in
the Abramowitz/ Stegun Handbook
The maximum approximation error is less than 1.5*10^-8
(compared to 4.5*10-4 in case of Hastings's approximation).
************************************************************/

proc(1) = incdfn_(p);
@ inverse of the standard normal cdf -------------------------------------- @
@ Odeh and Evans (1974), Applied Statistics, 23, 96-97 / max. error 1.5*E-8 @
local x,t,n,d,tlow,thgh;
  tlow = 1e-307;
  thgh = 1 - 1e-16;
  p = tlow.*(p .< tlow) + p.*(p .ge tlow);
  p = thgh.*(p .> thgh) + p.*(p .le thgh);
  t = sqrt(-2*ln(abs((p.>0.5)-p)));
  n = -0.322232431088-t-0.342242088547*t^2-0.0204231210245*t^3
      -0.0000453642210148*t^4;
  d = 0.0993484626060+0.588581570495*t+0.531103462366*t^2+
      0.103537752850*t^3+0.0038560700634*t^4;
  x = t+(n./d);
  x = (p.>0.5).*x-(p.<=0.5).*x;
  retp(x);
endp;




/*****************************************************************
Note that GAUSS now includes the proc cdfni.
From the help files
Purpose:    Computes the inverse of the cdf of the normal distribution.
Format:     x = cdfni(p);
Input:      p       NxK real matrix, normal probability levels, 0 <= p <= 1
Output:     x       NxK real matrix, normal deviates, such that cdfn(x)
                    equals p
cdfn(cdfni(p)) equals p to within the errors given below:

            p <= 4.6..e-308             -37.5 is returned
            4.6..e-308 < p < 5e-24      accurate to +/-5 in 12th digit
            5e-24 < p < 0.5             accurate to +/-1 in 13th digit
            0.5 < p < 1-2.22045e-16     accurate to +/-5 in 15th digit
            p >= 1-2.22045e-16          8.12589... is returned
*************************************************************************/





