00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"   
00008 
00009 #undef USE_OLD_METHOD
00010 
00011 #ifdef USE_OLD_METHOD
00012    double stas4( double , double ) ;
00013    double stinv( double , double ) ;
00014 #else
00015    double lnbeta( double , double ) ;
00016    double incbeta( double,double,double,double ) ;
00017    double incbeta_inverse( double,double,double,double ) ;
00018 #endif
00019 
00020 double qginv( double ) ;  /* prototype */
00021 
00022 int main( int argc , char * argv[] )
00023 {
00024    double pp , dof , tt , bb , binv ;
00025    double nn,ll,mm ;
00026    int quiet = 0 ;
00027 
00028    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ){
00029       printf("\n*** NOTE: This program has been superseded by program 'cdf' ***\n\n") ;
00030 
00031       printf("Usage #1: p2t p dof\n"
00032              "  where p   = double sided tail probability for t-distribution\n"
00033              "        dof = number of degrees of freedom to use\n"
00034              "  OUTPUT = t value that matches the input p\n"
00035              "\n"
00036              "Usage #2: p2t p N L M\n"
00037              "  where p   = double sided tail probability of beta distribution\n"
00038              "        N   = number of measured data points\n"
00039              "        L   = number of nuisance parameters (orts)\n"
00040              "        M   = number of fit parameters\n"
00041              "  OUTPUT = threshold for correlation coefficient\n"
00042              "\n"
00043              "Usage #3: p2t p\n"
00044              "  where p   = one sided tail probability of Gaussian distribution\n"
00045              "  OUTPUT = z value for which P(x>z) = p\n"
00046              "\n"
00047              "Usage #4: p2t p dof N\n"
00048              "  where p   = double sided tail probability for distribution of\n"
00049              "                the mean of N  iid zero-mean t-variables\n"
00050              "        dof = number of degrees of freedom of each t-variable\n"
00051              "        N   = number of t variables averaged\n"
00052              "  OUTPUT = threshold for the t average statistic\n"
00053              "  N.B.: The method used for this calculation is the Cornish-\n"
00054              "        Fisher expansion in N, and is only an approximation.\n"
00055              "        This also requires dof > 6, and the results will be\n"
00056              "        less accurate as dof approaches 6 from above!\n"
00057             ) ;
00058       exit(0) ;
00059    }
00060 
00061    if( strcmp(argv[1],"-q") == 0 ){
00062       argv++ ;
00063       argc-- ;
00064       quiet = 1 ;
00065       if( argc < 2 ) exit(0) ;
00066    }
00067 
00068    /** Gaussian statistic **/
00069 
00070    if( argc == 2 ){
00071       pp = strtod( argv[1] , NULL ) ;
00072       if( !quiet && (pp < 1.e-20 || pp >= 0.9999999) ){
00073          fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
00074       }
00075       tt = qginv(pp) ;
00076       if( quiet ) printf("%g\n",tt) ;
00077       else        printf("p = %g  z = %g\n",pp,tt) ;
00078    }
00079 
00080    /** t statistic **/
00081 
00082    else if( argc == 3 ){
00083       pp  = strtod( argv[1] , NULL ) ;
00084       dof = strtod( argv[2] , NULL ) ;
00085       if( !quiet && (pp <= 1.e-10 || pp >= 0.9999999) ){
00086          fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
00087       }
00088       if( dof <= 1.0 ){
00089          fprintf(stderr,"command line dof value is illegal!\a\n") ; exit(-1) ;
00090       }
00091 
00092 #ifdef USE_OLD_METHOD
00093       tt = stinv( pp , dof ) ;
00094 #else
00095       bb   = lnbeta( 0.5*dof , 0.5 ) ;
00096       binv = incbeta_inverse( pp, 0.5*dof , 0.5 , bb ) ;
00097       tt   = sqrt( dof*(1.0/binv-1.0) ) ;
00098 #endif
00099 
00100       if( quiet ) printf("%g\n",tt) ;
00101       else        printf("p = %g   dof = %g   t = %g\n",pp,dof,tt) ;
00102 
00103    /** beta statistic **/
00104 
00105    } else if( argc == 5 ){
00106       pp  = strtod( argv[1] , NULL ) ;
00107       nn  = strtod( argv[2] , NULL ) ;
00108       ll  = strtod( argv[3] , NULL ) ;
00109       mm  = strtod( argv[4] , NULL ) ;
00110       if( !quiet && (pp <= 1.e-20 || pp >= 0.9999999) ){
00111          fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
00112       }
00113       if( mm < 1.0 || ll < 0.0 || nn-ll-mm < 1.0 ){
00114          fprintf(stderr,"command line N,L,M values are illegal!\a\n");exit(1) ;
00115       }
00116       bb   = lnbeta( 0.5*mm , 0.5*(nn-ll-mm) ) ;
00117       binv = incbeta_inverse( pp, 0.5*(nn-ll-mm) , 0.5*mm , bb ) ;
00118       tt   = sqrt(1.0-binv) ;
00119       if( quiet ) printf("%g\n",tt) ;
00120       else        printf("p = %g  N = %g  L = %g  M = %g  rho = %g\n",
00121                          pp,nn,ll,mm,tt) ;
00122 
00123    /** averaged t statistic **/
00124 
00125    } else if( argc == 4 ){
00126       double ww , xx , gam2,gam4 ;
00127 
00128       pp  = strtod( argv[1] , NULL ) ;
00129       dof = strtod( argv[2] , NULL ) ;
00130       nn  = strtod( argv[3] , NULL ) ;
00131 
00132       if( !quiet && (pp <= 1.e-10 || pp >= 0.9999999) ){
00133          fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
00134       }
00135       if( dof <= 6.01 || nn < 1.0 ){
00136          fprintf(stderr,"command line dof or N value is illegal!\a\n");exit(-1);
00137       }
00138 
00139       /* 4th and 6th order moments */
00140 
00141       gam2 =   6.0 / ( (dof-4.0) * nn ) ;
00142       gam4 = 240.0 / ( (dof-6.0) * (dof-4.0) * nn * nn ) ;
00143 
00144       /* Cornish-Fisher expansion */
00145 
00146       xx = qginv( 0.5 * pp ) ;  /* Gaussian approx */
00147 
00148       ww = xx + gam2 * xx * (                       xx*xx -  3.0) / 24.0
00149               + gam4 * xx * (    xx*xx*xx*xx - 10.0*xx*xx + 15.0) / 720.0
00150        - gam2 * gam2 * xx * (3.0*xx*xx*xx*xx - 24.0*xx*xx + 29.0) / 384.0 ;
00151 
00152       tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
00153 
00154       if( quiet ) printf("%g\n",tt) ;
00155       else{
00156          printf("p = %g dof = %g N = %g 4-term t = %g",pp,dof,nn,tt) ;
00157 
00158          ww = xx + gam2 * xx * ( xx*xx - 3.0) / 24.0 ;
00159          tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
00160          printf(" [2-term=%g",tt) ;
00161          ww = xx ;
00162          tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
00163          printf(" 1-term=%g]\n",tt) ;
00164       }
00165    }
00166 
00167    exit(0) ;
00168 }
00169 
00170 #ifdef USE_OLD_METHOD
00171 /*----------------------------------------------------------------------
00172    code for inverse of central t distribution
00173    Inputs: p  = double sided tail probability
00174            nu = degrees of freedom
00175    Output: T such that P( |t| > T ) = p
00176 
00177    This version is only good for nu >= 5, since it uses the
00178    approximations in Abramowitz and Stegun, Eq. 26.7.5 (p. 949).
00179 ------------------------------------------------------------------------*/
00180 
00181 double stinv( double p , double nu )
00182 {
00183    double xg , t4 ;
00184    xg = qginv(0.5*p) ;
00185    t4 = stas4( xg , nu ) ;
00186    return t4 ;
00187 }
00188 
00189 double stas4( double x , double nu)  /* this code generated by Maple */
00190 {
00191    double t1,t2,t8,t9,t14,t17,t26,t34,t37 ;
00192    t1  = x*x;
00193    t2  = t1*x;
00194    t8  = t1*t1;
00195    t9  = t8*x;
00196    t14 = nu*nu;
00197    t17 = t8*t2;
00198    t26 = t8*t8;
00199    t34 = t14*t14;
00200    t37 = x+(t2/4+x/4)/nu
00201         +(5.0/96.0*t9+t2/6+x/32)/t14
00202         +(t17/128+19.0/384.0*t9
00203         +17.0/384.0*t2-5.0/128.0*x)/t14/nu
00204         +(79.0/92160.0*t26*x+97.0/11520.0*t17+247.0/15360.0*t9
00205                                           -t2/48-21.0/2048.0*x)/t34;
00206    return t37 ;
00207 }
00208 #endif   /* USE_OLD_METHOD */
00209 
00210 
00211 #ifndef USE_OLD_METHOD
00212 
00213 /*---------------------------------------------------------------
00214   compute log of complete beta function, using the
00215   Unix math library's log gamma function.  If this is
00216   not available, tough luck.
00217 -----------------------------------------------------------------*/
00218 
00219 double lnbeta( double p , double q )
00220 {
00221    return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
00222 }
00223 
00224 /*---------------------------------------------------------------
00225      TRANSLATED FROM THE ORIGINAL FORTRAN:
00226      algorithm as 63  appl. statist. (1973), vol.22, no.3
00227 
00228      computes incomplete beta function ratio for arguments
00229      x between zero and one, p and q positive.
00230      log of complete beta function, beta, is assumed to be known
00231 -----------------------------------------------------------------*/
00232 
00233 #define ZERO 0.0
00234 #define ONE  1.0
00235 #define ACU  1.0e-15
00236 
00237 double incbeta( double x , double p , double q , double beta )
00238 {
00239    double betain , psq , cx , xx,pp,qq , term,ai , temp , rx ;
00240    int indx , ns ;
00241 
00242    if( p <= ZERO || q <= ZERO || x < ZERO || x > ONE ) return -1.0 ;
00243 
00244    if( x == ZERO ) return ZERO ;
00245    if( x == ONE  ) return ONE ;
00246 
00247    /**  change tail if necessary and determine s **/
00248 
00249    psq = p+q ;
00250    cx  = ONE-x ;
00251    if(  p < psq*x ){
00252       xx   = cx ;
00253       cx   = x ;
00254       pp   = q ;
00255       qq   = p ;
00256       indx = 1 ;
00257    } else {
00258       xx   = x ;
00259       pp   = p ;
00260       qq   = q ;
00261       indx = 0 ;
00262    }
00263 
00264    term   = ONE ;
00265    ai     = ONE ;
00266    betain = ONE ;
00267    ns     = qq + cx*psq ;
00268 
00269    /** use soper's reduction formulae **/
00270 
00271       rx = xx/cx ;
00272 
00273 lab3:
00274       temp = qq-ai ;
00275       if(ns == 0) rx = xx ;
00276 
00277 lab4:
00278       term   = term*temp*rx/(pp+ai) ;
00279       betain = betain+term ;
00280       temp   = fabs(term) ;
00281       if(temp <= ACU && temp <= ACU*betain) goto lab5 ;
00282 
00283       ai = ai+ONE ;
00284       ns = ns-1 ;
00285       if(ns >= 0) goto lab3 ;
00286       temp = psq ;
00287       psq  = psq+ONE ;
00288       goto lab4 ;
00289 
00290 lab5:
00291       betain = betain*exp(pp*log(xx)+(qq-ONE)*log(cx)-beta)/pp ;
00292       if(indx) betain=ONE-betain ;
00293 
00294    return betain ;
00295 }
00296 
00297 /*----------------------------------------------------------------
00298     algorithm as 109 appl. statist. (1977), vol.26, no.1
00299     (replacing algorithm as 64  appl. statist. (1973),
00300     vol.22, no.3)
00301 
00302     Remark AS R83 and the correction in vol40(1) p.236 have been
00303     incorporated in this version.
00304 
00305     Computes inverse of the incomplete beta function
00306     ratio for given positive values of the arguments
00307     p and q, alpha between zero and one.
00308     log of complete beta function, beta, is assumed to be known.
00309 ------------------------------------------------------------------*/
00310 
00311 
00312 #define SAE   -15.0
00313 #define TWO     2.0
00314 #define THREE   3.0
00315 #define FOUR    4.0
00316 #define FIVE    5.0
00317 #define SIX     6.0
00318 
00319 #ifndef MAX
00320 #  define MAX(a,b) (((a)<(b)) ? (b) : (a))
00321 #  define MIN(a,b) (((a)>(b)) ? (b) : (a))
00322 #endif
00323 
00324 
00325 double incbeta_inverse( double alpha , double p , double q , double beta )
00326 {
00327    int indx , iex ;
00328    double fpu , xinbta , a,pp,qq, r,y,t,s,h,w , acu ,
00329           yprev,prev,sq , g,adj,tx,xin ;
00330 
00331    fpu = pow(10.0,SAE) ;
00332 
00333    if( p <= ZERO || q <= ZERO || alpha < ZERO || alpha > ONE ) return -1.0 ;
00334 
00335    if( alpha == ZERO ) return ZERO ;
00336    if( alpha == ONE  ) return ONE ;
00337 
00338    /** change tail if necessary **/
00339 
00340    if( alpha > 0.5 ){
00341       a    = ONE-alpha ;
00342       pp   = q ;
00343       qq   = p ;
00344       indx = 1 ;
00345     } else {
00346       a    = alpha ;
00347       pp   = p ;
00348       qq   = q ;
00349       indx = 0 ;
00350    }
00351 
00352    /** calculate the initial approximation **/
00353 
00354 lab2:
00355      r = sqrt(-log(a*a)) ;
00356      y = r - (2.30753 + 0.27061*r) / (ONE+(0.99229+0.04481*r)*r) ;
00357      if(pp > ONE && qq > ONE) goto lab5 ;
00358 
00359      r = qq+qq ;
00360      t = ONE/(9.0*qq) ;
00361      t = r * pow( (ONE-t+y*sqrt(t)) , 3.0 ) ;
00362      if( t <= ZERO ) goto lab3 ;
00363 
00364      t = (FOUR*pp+r-TWO)/t ;
00365      if( t <= ONE ) goto lab4 ;
00366 
00367      xinbta = ONE-TWO/(t+ONE) ; goto lab6 ;
00368 
00369 lab3:
00370      xinbta = ONE-exp((log((ONE-a)*qq)+beta)/qq) ; goto lab6 ;
00371 
00372 lab4:
00373      xinbta = exp((log(a*pp)+beta)/pp) ; goto lab6 ;
00374 
00375 lab5:
00376      r = (y*y-THREE)/SIX ;
00377      s = ONE/(pp+pp-ONE) ;
00378      t = ONE/(qq+qq-ONE) ;
00379      h = TWO/(s+t) ;
00380      w = y*sqrt(h+r)/h-(t-s)*(r+FIVE/SIX-TWO/(THREE*h)) ;
00381      xinbta = pp/(pp+qq*exp(w+w)) ;
00382 
00383      /** solve for x by a modified newton-raphson method **/
00384 
00385 lab6:
00386     r     = ONE-pp ;
00387     t     = ONE-qq ;
00388     yprev = ZERO ;
00389     sq    = ONE ;
00390     prev  = ONE ;
00391     if(xinbta < 0.0001) xinbta = 0.0001 ;
00392     if(xinbta > 0.9999) xinbta = 0.9999 ;
00393 
00394 #if 0
00395     iex = -5.0 / (pp*pp) - 1.0/(a*a) - 13.0 ; if( iex < SAE ) iex = SAE ;
00396     acu = pow(10.0,iex) ;
00397 #else
00398     acu = fpu ;
00399 #endif
00400 
00401 lab7:
00402       y = incbeta( xinbta , pp,qq,beta ) ;
00403       if( y < ZERO ) return -1.0 ;
00404       xin = xinbta ;
00405       y = (y-a)*exp(beta+r*log(xin)+t*log(ONE-xin)) ;
00406       if(y*yprev <= ZERO) prev = MAX(sq, fpu) ;
00407       g = ONE ;
00408 
00409 lab9:
00410       adj = g*y ;
00411       sq = adj*adj ;
00412       if(sq >= prev) goto lab10 ;
00413       tx = xinbta-adj ;
00414       if(tx >= ZERO && tx <= ONE) goto lab11 ;
00415 
00416 lab10:
00417       g = g/THREE ; goto lab9 ;
00418 
00419 lab11:
00420       if(tx == ZERO  || tx == ONE ) goto lab10 ;
00421       if(prev <= acu || y*y <= acu || fabs(xinbta-tx) < fpu) goto lab12 ;
00422       xinbta = tx ;
00423       yprev = y ;
00424       goto lab7 ;
00425 
00426 lab12:
00427       xinbta = tx ;
00428       if (indx) xinbta = ONE-xinbta ;
00429 #if 0
00430 printf("alpha = %g  incbeta = %g\n",alpha, incbeta(xinbta,p,q,beta) );
00431 #endif
00432       return xinbta ;
00433 }
00434 #endif  /* USE_OLD_METHOD */
00435 
00436 /*** given p, return x such that Q(x)=p, for 0 < p < 1 ***/
00437 
00438 double qginv( double p )
00439 {
00440    double dp , dx , dt , ddq , dq ;
00441    int    newt ;
00442 
00443    dp = (p <= 0.5) ? (p) : (1.0-p) ;   /* make between 0 and 0.5 */
00444 
00445    if( dp <= 0.0 ){
00446       dx = 13.0 ;
00447       return ( (p <= 0.5) ? (dx) : (-dx) ) ;
00448    }
00449 
00450 /**  Step 1:  use 26.2.23 from Abramowitz and Stegun **/
00451 
00452       dt = sqrt( -2.0 * log(dp) ) ;
00453       dx = dt
00454            - ((.010328e+0*dt + .802853e+0)*dt + 2.525517e+0)
00455            /(((.001308e+0*dt + .189269e+0)*dt + 1.432788e+0)*dt + 1.e+0) ;
00456 
00457 /**  Step 2:  do 3 Newton steps to improve this **/
00458 
00459       for( newt=0 ; newt < 3 ; newt++ ){
00460          dq  = 0.5e+0 * erfc( dx / 1.414213562373095e+0 ) - dp ;
00461          ddq = exp( -0.5e+0 * dx * dx ) / 2.506628274631000e+0 ;
00462          dx  = dx + dq / ddq ;
00463       }
00464 
00465       return ( (p <= 0.5) ? (dx) : (-dx) ) ;  /* return with correct sign */
00466 }