#include #include #include #include #define maxit 100 #define r2pi 0.3989422804014326e0 #define nhalf (-0.5e0) #define dennor(x) (r2pi*exp(nhalf*(x)*(x))) /* taken from DCDFLIB, Compiled and Written by: Barry W. Brown James Lovato Kathy Russell Department of Biomathematics, Box 237 The University of Texas, M.D. Anderson Cancer Center 1515 Holcombe Boulevard Houston, TX 77030 This work was supported by grant CA-16672 from the National Cancer Institute. DCDFLIB.C can be obtained by anonymous ftp to odin.mda.uth.tmc.edu (129.106.3.17) where it is available as /pub/unix/dcdflib.c.tar.Z The Fortran version of DCDFLIB is available as /pub/unix/dcdflib.f.tar.Z on the same machine. */ /* Modified by Allen Downey, October 11, 2000 tried to make it less obvious that this was translated from Fortran. Checked the results against a table for values 0.0, 0.2, ... 3.0. All agree to four decimal places. */ /* ********************************************************************** double cumnor (double x) Function Computes the cumulative of the normal distribution, i.e., the integral from -infinity to x of (1/sqrt(2*pi)) exp(-u*u/2) du X --> Upper limit of integration. X is DOUBLE PRECISION RESULT <-- Cumulative normal distribution. RESULT is DOUBLE PRECISION Renaming of function ANORM from: Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN Package of Special Function Routines and Test Drivers" acm Transactions on Mathematical Software. 19, 22-32. with slight modifications to return ccum and to deal with machine constants. ********************************************************************** Original Comments: ------------------------------------------------------------------ This function evaluates the normal distribution function: / x 1 | -t*t/2 P(x) = ----------- | e dt sqrt(2 pi) | /-oo The main computation evaluates near-minimax approximations derived from those in "Rational Chebyshev approximations for the error function" by W. J. Cody, Math. Comp., 1969, 631-637. This transportable program uses rational functions that theoretically approximate the normal distribution function to at least 18 significant decimal digits. The accuracy achieved depends on the arithmetic system, the compiler, the intrinsic functions, and proper selection of the machine-dependent constants. ******************************************************************* ******************************************************************* Explanation of machine-dependent constants. MIN = smallest machine representable number. EPS = argument below which anorm(x) may be represented by 0.5 and above which x*x will not underflow. A conservative value is the largest machine number X such that 1.0 + X = 1.0 to machine precision. ******************************************************************* ******************************************************************* Error returns The program returns ANORM = 0 for ARG .LE. XLOW. Intrinsic functions required are: ABS, AINT, EXP Author: W. J. Cody Mathematics and Computer Science Division Argonne National Laboratory Argonne, IL 60439 Latest modification: March 15, 1992 ------------------------------------------------------------------ */ double cumnor (double x) { static double a[5] = { 2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03, 1.8154981253343561249e04,6.5682337918207449113e-2 }; static double b[4] = { 4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04, 4.5507789335026729956e04 }; static double c[9] = { 3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01, 5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03, 1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8 }; static double d[8] = { 2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03, 6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04, 3.8912003286093271411e04,1.9685429676859990727e04 }; static double p[6] = { 2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2, 1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2 }; static double q[5] = { 1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2, 3.78239633202758244e-3,7.29751555083966205e-5 }; static double sqrpi = 0.39894228040143267794; static double thrsh = 0.66291; static double root32 = 5.656854248; static int i; static double del,eps,temp,xden,xnum,y,xsq,min; double result; // Machine dependent constants eps = DBL_EPSILON; min = DBL_MIN; y = fabs(x); // Evaluate anorm for |X| <= 0.66291 if(y <= thrsh) { xsq = 0.0; if(y > eps) xsq = x*x; xnum = a[4]*xsq; xden = xsq; for(i=0; i<3; i++) { xnum = (xnum+a[i])*xsq; xden = (xden+b[i])*xsq; } result = x*(xnum+a[3])/(xden+b[3]); temp = result; result = 0.5+temp; } // Evaluate anorm for 0.66291 <= |X| <= sqrt(32) else if(y <= root32) { xnum = c[8]*y; xden = y; for(i=0; i<7; i++) { xnum = (xnum+c[i])*y; xden = (xden+d[i])*y; } result = (xnum+c[7])/(xden+d[7]); xsq = floor(y*1.6)/1.6; del = (y-xsq)*(y+xsq); result = exp(-(xsq*xsq*0.5)) * exp(-(del*0.5)) * result; if(x > 0.0) { result = 1-result; } } // Evaluate anorm for |X| > sqrt(32) else { result = 0.0; xsq = 1.0/(x*x); xnum = p[5]*xsq; xden = xsq; for(i=0; i<4; i++) { xnum = (xnum+p[i])*xsq; xden = (xden+q[i])*xsq; } result = xsq*(xnum+p[4])/(xden+q[4]); result = (sqrpi-result)/y; xsq = floor(x*1.6)/1.6; del = (x-xsq)*(x+xsq); result = exp(-(xsq*xsq*0.5))*exp(-(del*0.5))*result; if(x > 0.0) { result = 1-result; } } if (result < min) result = 0.0; return result; } /* ********************************************************************** Double precision EVALuate a PoLynomial at X Function returns A(1) + A(2)*X + ... + A(N)*X**(N-1) Arguments A --> Array of coefficients of the polynomial. A is DOUBLE PRECISION(N) N --> Length of A, also degree of polynomial - 1. N is INTEGER X --> Point at which the polynomial is to be evaluated. X is DOUBLE PRECISION ***********************************************************************/ double devlpl (double a[], int n, double x) { static double devlpl,term; static int i; term = a[n-1]; for(i= n-1-1; i>=0; i--) { term = a[i]+term*x; } devlpl = term; return devlpl; } /* ********************************************************************** STarting VALue for Newton-Raphon calculation of Normal distribution Inverse Function Returns X such that CUMNOR(X) = P, i.e., the integral from - infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P Arguments P --> The probability whose normal deviate is sought. P is DOUBLE PRECISION Method The rational function on page 95 of Kennedy and Gentle, Statistical Computing, Marcel Dekker, NY , 1980. ********************************************************************** */ double stvaln (double p) { static double xden[5] = { 0.993484626060e-1,0.588581570495e0,0.531103462366e0,0.103537752850e0, 0.38560700634e-2 }; static double xnum[5] = { -0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1, -0.453642210148e-4 }; int i = 5; double stvaln, sign, y, z; if (p < 0.5) { sign = -1.0; z = p; } else { sign = 1.0; z = 1.0-p; } y = sqrt(-(2.0e0*log(z))); stvaln = y+devlpl(xnum, i, y) / devlpl (xden, i, y); stvaln = sign*stvaln; return stvaln; } /* ********************************************************************** Double precision NoRmal distribution INVerse Function Returns X such that CUMNOR(X) = P, i.e., the integral from - infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P Arguments P --> The probability whose normal deviate is sought. P is DOUBLE PRECISION Q --> 1-P P is DOUBLE PRECISION Method The rational function on page 95 of Kennedy and Gentle, Statistical Computing, Marcel Dekker, NY , 1980 is used as a start value for the Newton method of finding roots. Note If P or Q .lt. machine EPS returns +/- DINVNR(EPS) ********************************************************************** */ double invnor (double p) { double dinvnr,strtx,xcur,cum,pp,dx; int i, flag; double q = 1.0 - p; // find min of p and q flag = p <= q; if (flag) { pp = p; } else { pp = q; } // initialize strtx = stvaln(pp); xcur = strtx; // iterate for(i=1; i<=maxit; i++) { cum = cumnor (xcur); dx = (cum-pp) / dennor(xcur); xcur -= dx; if (fabs(dx/xcur) < DBL_EPSILON) { dinvnr = xcur; if (!flag) dinvnr = -dinvnr; return dinvnr; } } dinvnr = strtx; //IF WE GET HERE, NEWTON HAS FAILED if(!flag) dinvnr = -dinvnr; return dinvnr; }