Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_39.c

00001 #include "cdflib.h"
00002 double dinvnr(double *p,double *q)
00003 /*
00004 **********************************************************************
00005  
00006      double dinvnr(double *p,double *q)
00007      Double precision NoRmal distribution INVerse
00008  
00009  
00010                               Function
00011  
00012  
00013      Returns X  such that CUMNOR(X)  =   P,  i.e., the  integral from -
00014      infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
00015  
00016  
00017                               Arguments
00018  
00019  
00020      P --> The probability whose normal deviate is sought.
00021                     P is DOUBLE PRECISION
00022  
00023      Q --> 1-P
00024                     P is DOUBLE PRECISION
00025  
00026  
00027                               Method
00028  
00029  
00030      The  rational   function   on  page 95    of Kennedy  and  Gentle,
00031      Statistical Computing, Marcel Dekker, NY , 1980 is used as a start
00032      value for the Newton method of finding roots.
00033  
00034  
00035                               Note
00036  
00037  
00038      If P or Q .lt. machine EPS returns +/- DINVNR(EPS)
00039  
00040 **********************************************************************
00041 */
00042 {
00043 #define maxit 100
00044 #define eps (1.0e-13)
00045 #define r2pi 0.3989422804014326e0
00046 #define nhalf (-0.5e0)
00047 #define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
00048 static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
00049 static int i;
00050 static unsigned long qporq;
00051 /*
00052      ..
00053      .. Executable Statements ..
00054 */
00055 /*
00056      FIND MINIMUM OF P AND Q
00057 */
00058     qporq = *p <= *q;
00059     if(!qporq) goto S10;
00060     pp = *p;
00061     goto S20;
00062 S10:
00063     pp = *q;
00064 S20:
00065 /*
00066      INITIALIZATION STEP
00067 */
00068     strtx = stvaln(&pp);
00069     xcur = strtx;
00070 /*
00071      NEWTON INTERATIONS
00072 */
00073     for(i=1; i<=maxit; i++) {
00074         cumnor(&xcur,&cum,&ccum);
00075         dx = (cum-pp)/dennor(xcur);
00076         xcur -= dx;
00077         if(fabs(dx/xcur) < eps) goto S40;
00078     }
00079     dinvnr = strtx;
00080 /*
00081      IF WE GET HERE, NEWTON HAS FAILED
00082 */
00083     if(!qporq) dinvnr = -dinvnr;
00084     return dinvnr;
00085 S40:
00086 /*
00087      IF WE GET HERE, NEWTON HAS SUCCEDED
00088 */
00089     dinvnr = xcur;
00090     if(!qporq) dinvnr = -dinvnr;
00091     return dinvnr;
00092 #undef maxit
00093 #undef eps
00094 #undef r2pi
00095 #undef nhalf
00096 #undef dennor
00097 } /* END */
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

cdf_39.c
