00001 #include "cdflib.h"
00002 void cdfnor(int *which,double *p,double *q,double *x,double *mean,
00003             double *sd,int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdfnor(int *which,double *p,double *q,double *x,double *mean,
00007             double *sd,int *status,double *bound)
00008 
00009                Cumulative Distribution Function
00010                NORmal distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the normal
00017      distribution given values for the others.
00018 
00019 
00020                               Arguments
00021 
00022 
00023      WHICH  --> Integer indicating  which of the  next  parameter
00024      values is to be calculated using values  of the others.
00025      Legal range: 1..4
00026                iwhich = 1 : Calculate P and Q from X,MEAN and SD
00027                iwhich = 2 : Calculate X from P,Q,MEAN and SD
00028                iwhich = 3 : Calculate MEAN from P,Q,X and SD
00029                iwhich = 4 : Calculate SD from P,Q,X and MEAN
00030 
00031      P <--> The integral from -infinity to X of the normal density.
00032             Input range: (0,1].
00033 
00034      Q <--> 1-P.
00035             Input range: (0, 1].
00036             P + Q = 1.0.
00037 
00038      X < --> Upper limit of integration of the normal-density.
00039              Input range: ( -infinity, +infinity)
00040 
00041      MEAN <--> The mean of the normal density.
00042                Input range: (-infinity, +infinity)
00043 
00044      SD <--> Standard Deviation of the normal density.
00045              Input range: (0, +infinity).
00046 
00047      STATUS <-- 0 if calculation completed correctly
00048                -I if input parameter number I is out of range
00049                 1 if answer appears to be lower than lowest
00050                   search bound
00051                 2 if answer appears to be higher than greatest
00052                   search bound
00053                 3 if P + Q .ne. 1
00054 
00055      BOUND <-- Undefined if STATUS is 0
00056 
00057                Bound exceeded by parameter number I if STATUS
00058                is negative.
00059 
00060                Lower search bound if STATUS is 1.
00061 
00062                Upper search bound if STATUS is 2.
00063 
00064 
00065                               Method
00066 
00067 
00068 
00069 
00070      A slightly modified version of ANORM from
00071 
00072      Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
00073      Package of Special Function Routines and Test Drivers"
00074      acm Transactions on Mathematical Software. 19, 22-32.
00075 
00076      is used to calulate the  cumulative standard normal distribution.
00077 
00078      The rational functions from pages  90-95  of Kennedy and Gentle,
00079      Statistical  Computing,  Marcel  Dekker, NY,  1980 are  used  as
00080      starting values to Newton's Iterations which compute the inverse
00081      standard normal.  Therefore no  searches  are necessary for  any
00082      parameter.
00083 
00084      For X < -15, the asymptotic expansion for the normal is used  as
00085      the starting value in finding the inverse standard normal.
00086      This is formula 26.2.12 of Abramowitz and Stegun.
00087 
00088 
00089                               Note
00090 
00091 
00092       The normal density is proportional to
00093       exp( - 0.5 * (( X - MEAN)/SD)**2)
00094 
00095 **********************************************************************/
00096 {
00097 static int K1 = 1;
00098 static double z,pq;
00099 /*
00100      ..
00101      .. Executable Statements ..
00102 */
00103 /*
00104      Check arguments
00105 */
00106     *status = 0;
00107     if(!(*which < 1 || *which > 4)) goto S30;
00108     if(!(*which < 1)) goto S10;
00109     *bound = 1.0e0;
00110     goto S20;
00111 S10:
00112     *bound = 4.0e0;
00113 S20:
00114     *status = -1;
00115     return;
00116 S30:
00117     if(*which == 1) goto S70;
00118 /*
00119      P
00120 */
00121     if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
00122     if(!(*p <= 0.0e0)) goto S40;
00123     *bound = 0.0e0;
00124     goto S50;
00125 S40:
00126     *bound = 1.0e0;
00127 S50:
00128     *status = -2;
00129     return;
00130 S70:
00131 S60:
00132     if(*which == 1) goto S110;
00133 /*
00134      Q
00135 */
00136     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00137     if(!(*q <= 0.0e0)) goto S80;
00138     *bound = 0.0e0;
00139     goto S90;
00140 S80:
00141     *bound = 1.0e0;
00142 S90:
00143     *status = -3;
00144     return;
00145 S110:
00146 S100:
00147     if(*which == 1) goto S150;
00148 /*
00149      P + Q
00150 */
00151     pq = *p+*q;
00152     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S140;
00153     if(!(pq < 0.0e0)) goto S120;
00154     *bound = 0.0e0;
00155     goto S130;
00156 S120:
00157     *bound = 1.0e0;
00158 S130:
00159     *status = 3;
00160     return;
00161 S150:
00162 S140:
00163     if(*which == 4) goto S170;
00164 /*
00165      SD
00166 */
00167     if(!(*sd <= 0.0e0)) goto S160;
00168     *bound = 0.0e0;
00169     *status = -6;
00170     return;
00171 S170:
00172 S160:
00173 /*
00174      Calculate ANSWERS
00175 */
00176     if(1 == *which) {
00177 /*
00178      Computing P
00179 */
00180         z = (*x-*mean)/ *sd;
00181         cumnor(&z,p,q);
00182     }
00183     else if(2 == *which) {
00184 /*
00185      Computing X
00186 */
00187         z = dinvnr(p,q);
00188         *x = *sd*z+*mean;
00189     }
00190     else if(3 == *which) {
00191 /*
00192      Computing the MEAN
00193 */
00194         z = dinvnr(p,q);
00195         *mean = *x-*sd*z;
00196     }
00197     else if(4 == *which) {
00198 /*
00199      Computing SD
00200 */
00201         z = dinvnr(p,q);
00202         *sd = (*x-*mean)/z;
00203     }
00204     return;
00205 } /* END */