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 }

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 }

00220 {
00221    return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
00222 }

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 }

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 }