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 /*** Not actually MRI processing routines -- just some statistics ***/
00010 
00011 /*-----------------------------------------------------------------
00012   Student t-statistic:
00013     given threshold, compute 2-sided tail probability  ("t2p"), or
00014     given statistic, compute equivalent N(0,1) deviate ("t2z"), or
00015     given 2-sided tail probability, compute threshold  ("p2t").
00016 -------------------------------------------------------------------*/
00017 
00018 double student_p2t( double pp , double dof )
00019 {
00020    double bb , binv , tt ;
00021 
00022    if( pp  <= 0.0      ) return 99.99 ;
00023    if( pp  >= 0.999999 ) return 0.0 ;
00024    if( dof < 1.0       ) return 0.0 ;
00025 
00026    bb   = lnbeta( 0.5*dof , 0.5 ) ;
00027    binv = incbeta_inverse( pp, 0.5*dof , 0.5 , bb ) ;
00028    tt   = sqrt( dof*(1.0/binv-1.0) ) ;
00029    return tt ;
00030 }
00031 
00032 double student_t2p( double tt , double dof )
00033 {
00034    double bb , xx , pp ;
00035 
00036    if( tt <= 0.0 || dof < 1.0 ) return 1.0 ;
00037 
00038    bb = lnbeta( 0.5*dof , 0.5 ) ;
00039    xx = dof/(dof + tt*tt) ;
00040    pp = incbeta( xx , 0.5*dof , 0.5 , bb ) ;
00041    return pp ;
00042 }
00043 
00044 double student_t2z( double tt , double dof )
00045 {
00046    static double bb , dof_old = -666.666 ;
00047    double xx , pp ;
00048 
00049    if( dof != dof_old ){
00050       bb      = lnbeta( 0.5*dof , 0.5 ) ;
00051       dof_old = dof ;
00052    }
00053 
00054    xx = dof/(dof + tt*tt) ;
00055    pp = incbeta( xx , 0.5*dof , 0.5 , bb ) ;
00056 
00057    if( tt > 0.0 ) pp = 1.0 - 0.5 * pp ;
00058    else           pp = 0.5 * pp ;
00059 
00060    xx = qginv(pp) ;
00061    return -xx ;
00062 }
00063 
00064 /*---------------------------------------------------------------------
00065    Correlation coefficient statistic:
00066      nsam = # of samples used to compute rho ( > nfit+nort )
00067      nfit = # of fitting parameters  ( >= 1 )
00068      nort = # of nuisance parameters ( >= 1 )
00069 -----------------------------------------------------------------------*/
00070 
00071 double correl_p2t( double pp , double nsam , double nfit , double nort )
00072 {
00073    double bb , binv , rho ;
00074 
00075    if( pp <= 0.0      ) return 0.999 ;
00076    if( pp >= 0.999999 ) return 0.0 ;
00077 
00078    if( nsam <= nfit+nort || nfit < 1.0 || nort < 1.0 ) return 0.0 ;
00079 
00080    bb   = lnbeta( 0.5*nfit , 0.5*(nsam-nfit-nort) ) ;
00081    binv = incbeta_inverse( pp, 0.5*(nsam-nfit-nort) , 0.5*nfit , bb ) ;
00082    rho  = sqrt(1.0-binv) ;
00083    return rho ;
00084 }
00085 
00086 double correl_t2p( double rho , double nsam , double nfit , double nort )
00087 {
00088    double bb , xx , pp ;
00089 
00090    if( rho <= 0.0 ||
00091        nsam <= nfit+nort || nfit < 1.0 || nort < 1.0 ) return 1.0 ;
00092 
00093    if( rho >= 0.9999999 ) return 0.0 ;
00094 
00095    bb   = lnbeta( 0.5*nfit , 0.5*(nsam-nfit-nort) ) ;
00096    xx   = 1.0 - rho*rho ;
00097    pp   = incbeta( xx , 0.5*(nsam-nfit-nort) , 0.5*nfit , bb ) ;
00098    return pp ;
00099 }
00100 
00101 /******************************/
00102 /*** added this 17 Sep 1998 ***/
00103 
00104 double correl_t2z( double rho , double nsam , double nfit , double nort )
00105 {
00106    double pp , xx ;
00107    pp = 0.5 * correl_t2p( fabs(rho) , nsam , nfit , nort ) ;
00108    xx = qginv(pp) ;
00109    return ( (rho > 0) ? xx : -xx ) ;
00110 }
00111 
00112 /*----------------------------------------------------------
00113    Averaged t-statistic
00114 ------------------------------------------------------------*/
00115 
00116 double studave_p2t( double pp , double dof , double nn )
00117 {
00118    double ww , xx , gam2,gam4 , tt ;
00119 
00120    if( pp <= 0.0      ) return 99.99 ;
00121    if( pp >= 0.999999 ) return 0.0 ;
00122 
00123    if( dof < 6.01 || nn < 1.0 ) return 0.0 ;
00124 
00125    /* 4th and 6th order moments (or scaled cumulants) */
00126 
00127    gam2 =   6.0 / ( (dof-4.0) * nn ) ;
00128    gam4 = 240.0 / ( (dof-6.0) * (dof-4.0) * nn * nn ) ;
00129 
00130    /* Cornish-Fisher expansion */
00131 
00132    xx = qginv( 0.5 * pp ) ;  /* Gaussian approx */
00133 
00134    ww = xx + gam2 * xx * (                       xx*xx -  3.0) /  24.0
00135            + gam4 * xx * (    xx*xx*xx*xx - 10.0*xx*xx + 15.0) / 720.0
00136     - gam2 * gam2 * xx * (3.0*xx*xx*xx*xx - 24.0*xx*xx + 29.0) / 384.0 ;
00137 
00138    tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
00139    return tt ;
00140 }
00141 
00142 double studave_t2p( double tt , double dof , double nn )
00143 {
00144    static int nc = 0 ;
00145    if( nc < 9 ){
00146       fprintf(stderr,"*** studave_t2p: NOT IMPLEMENTED YET!\n") ; nc++ ;
00147    }
00148    return 0.0 ;
00149 }
00150 
00151 double studave_t2z( double tt , double dof , double nn )
00152 {
00153    static int nc = 0 ;
00154    if( nc < 9 ){
00155       fprintf(stderr,"*** studave_t2z: NOT IMPLEMENTED YET!\n") ; nc++ ;
00156    }
00157    return 0.0 ;
00158 }
00159 
00160 /***********************************************************************/
00161 /*** The routines below here are wrappers for the cdflib routines    ***/
00162 /*** (cdf_*.c) from U Texas -- see file cdflib.txt for the details.  ***/
00163 /***********************************************************************/
00164 
00165 /*---------------------------------------------------------------
00166   F statistic: single sided
00167 -----------------------------------------------------------------*/
00168 
00169 double fstat_p2t( double pp , double dofnum , double dofden )
00170 {
00171    int which , status ;
00172    double p , q , f , dfn , dfd , bound ;
00173 
00174    if( pp <= 0.0      ) return 999.99 ;
00175    if( pp >= 0.999999 ) return 0.0 ;
00176 
00177    which  = 2 ;
00178    p      = 1.0 - pp ;  /* 20 Jan 1999: p and q were switched! */
00179    q      = pp ;
00180    f      = 0.0 ;
00181    dfn    = dofnum ;
00182    dfd    = dofden ;
00183 
00184    cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
00185 
00186    if( status == 0 ) return f ;
00187    else              return 0.0 ;
00188 }
00189 
00190 double fstat_t2p( double ff , double dofnum , double dofden )
00191 {
00192    int which , status ;
00193    double p , q , f , dfn , dfd , bound ;
00194 
00195    which  = 1 ;
00196    p      = 0.0 ;
00197    q      = 0.0 ;
00198    f      = ff ;
00199    dfn    = dofnum ;
00200    dfd    = dofden ;
00201 
00202    cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
00203 
00204    if( status == 0 ) return q ;
00205    else              return 1.0 ;
00206 }
00207 
00208 /******************************/
00209 /*** added this 17 Sep 1998 ***/
00210 
00211 double fstat_t2z( double ff , double dofnum , double dofden )
00212 {
00213    double pp ;
00214    pp = 0.5 * fstat_t2p( ff , dofnum , dofden ) ;
00215    return qginv(pp) ;
00216 }
00217 
00218 /*---------------------------------------------------------------
00219   compute log of complete beta function, using the
00220   Unix math library's log gamma function.  If this is
00221   not available, see the end of this file.
00222 -----------------------------------------------------------------*/
00223 
00224 double lnbeta( double p , double q )
00225 {
00226    return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
00227 }
00228 
00229 /*---------------------------------------------------------------
00230      TRANSLATED FROM THE ORIGINAL FORTRAN:
00231      algorithm as 63  appl. statist. (1973), vol.22, no.3
00232 
00233      computes incomplete beta function ratio for arguments
00234      x between zero and one, p and q positive.
00235      log of complete beta function, beta, is assumed to be known
00236 -----------------------------------------------------------------*/
00237 
00238 #define ZERO 0.0
00239 #define ONE  1.0
00240 #define ACU  1.0e-15
00241 
00242 double incbeta( double x , double p , double q , double beta )
00243 {
00244    double betain , psq , cx , xx,pp,qq , term,ai , temp , rx ;
00245    int indx , ns ;
00246 
00247    if( p <= ZERO || q <= ZERO ) return -1.0 ;  /* error! */
00248 
00249    if( x <= ZERO ) return ZERO ;
00250    if( x >= ONE  ) return ONE ;
00251 
00252    /**  change tail if necessary and determine s **/
00253 
00254    psq = p+q ;
00255    cx  = ONE-x ;
00256    if(  p < psq*x ){
00257       xx   = cx ;
00258       cx   = x ;
00259       pp   = q ;
00260       qq   = p ;
00261       indx = 1 ;
00262    } else {
00263       xx   = x ;
00264       pp   = p ;
00265       qq   = q ;
00266       indx = 0 ;
00267    }
00268 
00269    term   = ONE ;
00270    ai     = ONE ;
00271    betain = ONE ;
00272    ns     = qq + cx*psq ;
00273 
00274    /** use soper's reduction formulae **/
00275 
00276       rx = xx/cx ;
00277 
00278 lab3:
00279       temp = qq-ai ;
00280       if(ns == 0) rx = xx ;
00281 
00282 lab4:
00283       term   = term*temp*rx/(pp+ai) ;
00284       betain = betain+term ;
00285       temp   = fabs(term) ;
00286       if(temp <= ACU && temp <= ACU*betain) goto lab5 ;
00287 
00288       ai = ai+ONE ;
00289       ns = ns-1 ;
00290       if(ns >= 0) goto lab3 ;
00291       temp = psq ;
00292       psq  = psq+ONE ;
00293       goto lab4 ;
00294 
00295 lab5:
00296       betain = betain*exp(pp*log(xx)+(qq-ONE)*log(cx)-beta)/pp ;
00297       if(indx) betain=ONE-betain ;
00298 
00299    return betain ;
00300 }
00301 
00302 /*----------------------------------------------------------------
00303     algorithm as 109 appl. statist. (1977), vol.26, no.1
00304     (replacing algorithm as 64  appl. statist. (1973),
00305     vol.22, no.3)
00306 
00307     Remark AS R83 and the correction in vol40(1) p.236 have been
00308     incorporated in this version.
00309 
00310     Computes inverse of the incomplete beta function
00311     ratio for given positive values of the arguments
00312     p and q, alpha between zero and one.
00313     log of complete beta function, beta, is assumed to be known.
00314 ------------------------------------------------------------------*/
00315 
00316 #define SAE   -15.0
00317 #define TWO     2.0
00318 #define THREE   3.0
00319 #define FOUR    4.0
00320 #define FIVE    5.0
00321 #define SIX     6.0
00322 
00323 #ifndef MAX
00324 #  define MAX(a,b) (((a)<(b)) ? (b) : (a))
00325 #  define MIN(a,b) (((a)>(b)) ? (b) : (a))
00326 #endif
00327 
00328 double incbeta_inverse( double alpha , double p , double q , double beta )
00329 {
00330    int indx , iex ;
00331    double fpu , xinbta , a,pp,qq, r,y,t,s,h,w , acu ,
00332           yprev,prev,sq , g,adj,tx,xin ;
00333 
00334    fpu = pow(10.0,SAE) ;
00335 
00336    if( p <= ZERO || q <= ZERO || alpha < ZERO || alpha > ONE ) return -1.0 ;
00337 
00338    if( alpha == ZERO ) return ZERO ;
00339    if( alpha == ONE  ) return ONE ;
00340 
00341    /** change tail if necessary **/
00342 
00343    if( alpha > 0.5 ){
00344       a    = ONE-alpha ;
00345       pp   = q ;
00346       qq   = p ;
00347       indx = 1 ;
00348     } else {
00349       a    = alpha ;
00350       pp   = p ;
00351       qq   = q ;
00352       indx = 0 ;
00353    }
00354 
00355    /** calculate the initial approximation **/
00356 
00357 lab2:
00358      r = sqrt(-log(a*a)) ;
00359      y = r - (2.30753 + 0.27061*r) / (ONE+(0.99229+0.04481*r)*r) ;
00360      if(pp > ONE && qq > ONE) goto lab5 ;
00361 
00362      r = qq+qq ;
00363      t = ONE/(9.0*qq) ;
00364      t = r * pow( (ONE-t+y*sqrt(t)) , 3.0 ) ;
00365      if( t <= ZERO ) goto lab3 ;
00366 
00367      t = (FOUR*pp+r-TWO)/t ;
00368      if( t <= ONE ) goto lab4 ;
00369 
00370      xinbta = ONE-TWO/(t+ONE) ; goto lab6 ;
00371 
00372 lab3:
00373      xinbta = ONE-exp((log((ONE-a)*qq)+beta)/qq) ; goto lab6 ;
00374 
00375 lab4:
00376      xinbta = exp((log(a*pp)+beta)/pp) ; goto lab6 ;
00377 
00378 lab5:
00379      r = (y*y-THREE)/SIX ;
00380      s = ONE/(pp+pp-ONE) ;
00381      t = ONE/(qq+qq-ONE) ;
00382      h = TWO/(s+t) ;
00383      w = y*sqrt(h+r)/h-(t-s)*(r+FIVE/SIX-TWO/(THREE*h)) ;
00384      xinbta = pp/(pp+qq*exp(w+w)) ;
00385 
00386      /** solve for x by a modified newton-raphson method **/
00387 
00388 lab6:
00389     r     = ONE-pp ;
00390     t     = ONE-qq ;
00391     yprev = ZERO ;
00392     sq    = ONE ;
00393     prev  = ONE ;
00394     if(xinbta < 0.0001) xinbta = 0.0001 ;
00395     if(xinbta > 0.9999) xinbta = 0.9999 ;
00396 
00397 #if 0
00398     iex = -5.0 / (pp*pp) - 1.0/(a*a) - 13.0 ; if( iex < SAE ) iex = SAE ;
00399     acu = pow(10.0,iex) ;
00400 #else
00401     acu = fpu ;
00402 #endif
00403 
00404 lab7:
00405       y = incbeta( xinbta , pp,qq,beta ) ;
00406       if( y < ZERO ) return -1.0 ;
00407       xin = xinbta ;
00408       y = (y-a)*exp(beta+r*log(xin)+t*log(ONE-xin)) ;
00409       if(y*yprev <= ZERO) prev = MAX(sq, fpu) ;
00410       g = ONE ;
00411 
00412 lab9:
00413       adj = g*y ;
00414       sq = adj*adj ;
00415       if(sq >= prev) goto lab10 ;
00416       tx = xinbta-adj ;
00417       if(tx >= ZERO && tx <= ONE) goto lab11 ;
00418 
00419 lab10:
00420       g = g/THREE ; goto lab9 ;
00421 
00422 lab11:
00423       if(tx == ZERO  || tx == ONE ) goto lab10 ;
00424       if(prev <= acu || y*y <= acu || fabs(xinbta-tx) < fpu) goto lab12 ;
00425       xinbta = tx ;
00426       yprev = y ;
00427       goto lab7 ;
00428 
00429 lab12:
00430       xinbta = tx ;
00431       if (indx) xinbta = ONE-xinbta ;
00432 #if 0
00433    printf("alpha = %g  incbeta = %g\n",alpha, incbeta(xinbta,p,q,beta) );
00434 #endif
00435       return xinbta ;
00436 }
00437 
00438 /*******************************************************************/
00439 /****    Given p, return x such that Q(x)=p, for 0 < p < 1.     ****/
00440 /****    Q(x) = 1-P(x) = reversed cdf of N(0,1) variable.       ****/
00441 /*******************************************************************/
00442 
00443 double qg( double x ){ return 0.5*erfc(x/1.414213562373095); }
00444 
00445 double log10qg( double x )
00446 {
00447   double v = qg(x) ;
00448   if( v > 0.0 ) return log10(v) ;
00449   return -99.99 ;
00450 }
00451 
00452 double qginv( double p )
00453 {
00454    double dp , dx , dt , ddq , dq ;
00455    int    newt ;                       /* not Gingrich, but Isaac */
00456 
00457    dp = (p <= 0.5) ? (p) : (1.0-p) ;   /* make between 0 and 0.5 */
00458 
00459    if( dp <= 1.e-37 ){
00460       dx = 13.0 ;                      /* 13 sigma has p < 10**(-38) */
00461       return ( (p <= 0.5) ? (dx) : (-dx) ) ;
00462    }
00463 
00464 /**  Step 1:  use 26.2.23 from Abramowitz and Stegun **/
00465 
00466       dt = sqrt( -2.0 * log(dp) ) ;
00467       dx = dt
00468            - ((.010328*dt + .802853)*dt + 2.515517)
00469            /(((.001308*dt + .189269)*dt + 1.432788)*dt + 1.) ;
00470 
00471 /**  Step 2:  do 3 Newton steps to improve this
00472               (uses the math library erfc function) **/
00473 
00474       for( newt=0 ; newt < 3 ; newt++ ){
00475          dq  = 0.5 * erfc( dx / 1.414213562373095 ) - dp ;
00476          ddq = exp( -0.5 * dx * dx ) / 2.506628274631000 ;
00477          dx  = dx + dq / ddq ;
00478       }
00479 
00480       if( dx > 13.0 ) dx = 13.0 ;
00481       return ( (p <= 0.5) ? (dx) : (-dx) ) ;  /* return with correct sign */
00482 }
00483 
00484 /*---------------------------------------------------------------
00485   Compute double-sided tail probability for normal distribution.
00486 -----------------------------------------------------------------*/
00487 
00488 double normal_t2p( double zz )
00489 {
00490    int which , status ;
00491    double p , q , x , mean,sd,bound ;
00492 
00493    if( zz <= 0.0 ) return 1.0 ;
00494 
00495    which  = 1 ;
00496    p      = 0.0 ;
00497    q      = 0.0 ;
00498    x      = zz ;
00499    mean   = 0.0 ;
00500    sd     = 1.0 ;
00501 
00502    cdfnor( &which , &p , &q , &x , &mean , &sd , &status , &bound ) ;
00503 
00504    if( status == 0 ) return 2.0*q ;  /* double sided prob = 2 times single sided */
00505    else              return 1.0 ;
00506 }
00507 
00508 /******************************/
00509 /*** added this 17 Sep 1998 ***/
00510 
00511 double normal_p2t( double qq )
00512 {
00513    int which , status ;
00514    double p , q , x , mean,sd,bound ;
00515 
00516    if( qq <= 0.0      ) return 9.99 ;
00517    if( qq >= 0.999999 ) return 0.0 ;
00518 
00519    which  = 2 ;
00520    p      = 1.0 - 0.5 * qq ;
00521    q      = 0.5 * qq ;        /* single sided prob = 1/2 of double sided */
00522    x      = 0.0 ;
00523    mean   = 0.0 ;
00524    sd     = 1.0 ;
00525 
00526    cdfnor( &which , &p , &q , &x , &mean , &sd , &status , &bound ) ;
00527    return x ;
00528 }
00529 
00530 /*----------------------------------------------------------------
00531    Compute single-sided tail probability for Chi-square
00532 ------------------------------------------------------------------*/
00533 
00534 double chisq_t2p( double xx , double dof )
00535 {
00536    int which , status ;
00537    double p,q,x,df,bound ;
00538 
00539    which  = 1 ;
00540    p      = 0.0 ;
00541    q      = 0.0 ;
00542    x      = xx ;
00543    df     = dof ;
00544 
00545    cdfchi( &which , &p , &q , &x , &df , &status , &bound ) ;
00546 
00547    if( status == 0 ) return q ;
00548    else              return 1.0 ;
00549 }
00550 
00551 /******************************/
00552 /*** added this 17 Sep 1998 ***/
00553 
00554 double chisq_p2t( double qq , double dof )
00555 {
00556    int which , status ;
00557    double p,q,x,df,bound ;
00558 
00559    if( qq <= 0.0      ) return 999.9 ;
00560    if( qq >= 0.999999 ) return 0.0 ;
00561 
00562    which  = 2 ;
00563    p      = 1.0 - qq ;
00564    q      = qq ;
00565    x      = 0.0 ;
00566    df     = dof ;
00567 
00568    cdfchi( &which , &p , &q , &x , &df , &status , &bound ) ;
00569    return x ;
00570 }
00571 
00572 double chisq_t2z( double xx , double dof )
00573 {
00574    double pp ;
00575    pp = 0.5 * chisq_t2p( xx , dof ) ;
00576    return qginv(pp) ;
00577 }
00578 
00579 /*----------------------------------------------------------------
00580    Compute upper tail probability for incomplete beta distribution
00581 ------------------------------------------------------------------*/
00582 
00583 double beta_t2p( double xx , double aa , double bb )
00584 {
00585    int which , status ;
00586    double p,q,x,y,a,b,bound ;
00587 
00588    which  = 1 ;
00589    p      = 0.0 ;
00590    q      = 0.0 ;
00591    x      = xx ;
00592    y      = 1.0 - xx ;
00593    a      = aa ;
00594    b      = bb ;
00595 
00596    cdfbet( &which , &p , &q , &x , &y , &a , &b ,  &status , &bound ) ;
00597 
00598    if( status == 0 ) return q ;
00599    else              return 1.0 ;
00600 }
00601 
00602 /******************************/
00603 /*** added this 17 Sep 1998 ***/
00604 
00605 double beta_t2z( double xx , double aa , double bb )
00606 {
00607    double pp ;
00608    pp = 0.5 * beta_t2p( xx , aa , bb ) ;
00609    return qginv(pp) ;
00610 }
00611 
00612 double beta_p2t( double qq , double aa , double bb )
00613 {
00614    int which , status ;
00615    double p,q,x,y,a,b,bound ;
00616 
00617    if( qq <= 0.0      ) return 0.9999 ;
00618    if( qq >= 0.999999 ) return 0.0 ;
00619 
00620    which  = 2 ;
00621    p      = 1.0 - qq ;
00622    q      = qq ;
00623    x      = 0.0 ;
00624    y      = 1.0 ;
00625    a      = aa ;
00626    b      = bb ;
00627 
00628    cdfbet( &which , &p , &q , &x , &y , &a , &b ,  &status , &bound ) ;
00629 
00630    return x ;
00631 }
00632 
00633 /*----------------------------------------------------------------
00634    Compute upper tail probability for binomial distribution
00635    (that is, the probability that more than ss out of ntrial
00636     trials were successful).
00637 ------------------------------------------------------------------*/
00638 
00639 double binomial_t2p( double ss , double ntrial , double ptrial )
00640 {
00641    int which , status ;
00642    double p,q, s,xn,pr,ompr,bound ;
00643 
00644    which  = 1 ;
00645    p      = 0.0 ;
00646    q      = 0.0 ;
00647    s      = ss ;
00648    xn     = ntrial ;
00649    pr     = ptrial ;
00650    ompr   = 1.0 - ptrial ;
00651 
00652    cdfbin( &which , &p , &q , &s , &xn , &pr , &ompr , &status , &bound ) ;
00653 
00654    if( status == 0 ) return q ;
00655    else              return 1.0 ;
00656 }
00657 
00658 /******************************/
00659 /*** added this 17 Sep 1998 ***/
00660 
00661 double binomial_t2z( double ss , double ntrial , double ptrial )
00662 {
00663    double pp ;
00664    pp = 0.5 * binomial_t2p( ss , ntrial , ptrial ) ;
00665    return qginv(pp) ;
00666 }
00667 
00668 double binomial_p2t( double qq , double ntrial , double ptrial )
00669 {
00670    int which , status ;
00671    double p,q, s,xn,pr,ompr,bound ;
00672 
00673    if( qq <= 0.0      ) return 99.99 ;
00674    if( qq >= 0.999999 ) return 0.0 ;
00675 
00676    which  = 2 ;
00677    p      = 1.0 - qq ;
00678    q      = qq ;
00679    s      = 0.0 ;
00680    xn     = ntrial ;
00681    pr     = ptrial ;
00682    ompr   = 1.0 - ptrial ;
00683 
00684    cdfbin( &which , &p , &q , &s , &xn , &pr , &ompr , &status , &bound ) ;
00685 
00686    if( status == 0 ) return s ;
00687    else              return 0.0 ;
00688 }
00689 
00690 /*----------------------------------------------------------------
00691    Compute upper tail probability for the gamma distribution.
00692 ------------------------------------------------------------------*/
00693 
00694 double gamma_t2p( double xx , double sh , double sc )
00695 {
00696    int which , status ;
00697    double p,q, x,shape,scale,bound ;
00698 
00699    which  = 1 ;
00700    p      = 0.0 ;
00701    q      = 0.0 ;
00702    x      = xx ;
00703    shape  = sh ;
00704    scale  = sc ;
00705 
00706    cdfgam( &which , &p , &q , &x , &shape , &scale , &status , &bound ) ;
00707 
00708    if( status == 0 ) return q ;
00709    else              return 1.0 ;
00710 }
00711 
00712 /******************************/
00713 /*** added this 17 Sep 1998 ***/
00714 
00715 double gamma_t2z( double xx , double sh , double sc )
00716 {
00717    double pp ;
00718    pp = 0.5 * gamma_t2p( xx , sh , sc ) ;
00719    return qginv(pp) ;
00720 }
00721 
00722 double gamma_p2t( double qq , double sh , double sc )
00723 {
00724    int which , status ;
00725    double p,q, x,shape,scale,bound ;
00726 
00727    if( qq <= 0.0      ) return 999.9 ;
00728    if( qq >= 0.999999 ) return 0.0 ;
00729 
00730    which  = 2 ;
00731    p      = 1.0 - qq ;
00732    q      = qq ;
00733    x      = 0.0 ;
00734    shape  = sh ;
00735    scale  = sc ;
00736 
00737    cdfgam( &which , &p , &q , &x , &shape , &scale , &status , &bound ) ;
00738 
00739    return x ;
00740 }
00741 
00742 /*----------------------------------------------------------------
00743    Compute upper tail probability for the Poisson distribution
00744    (that is, the probability that the value is greater than xx).
00745 ------------------------------------------------------------------*/
00746 
00747 double poisson_t2p( double xx , double lambda )
00748 {
00749    int which , status ;
00750    double p,q, s,xlam,bound ;
00751 
00752    which  = 1 ;
00753    p      = 0.0 ;
00754    q      = 0.0 ;
00755    s      = xx ;
00756    xlam   = lambda ;
00757 
00758    cdfpoi( &which , &p , &q , &s , &xlam , &status , &bound ) ;
00759 
00760    if( status == 0 ) return q ;
00761    else              return 1.0 ;
00762 }
00763 
00764 /******************************/
00765 /*** added this 17 Sep 1998 ***/
00766 
00767 double poisson_t2z( double xx , double lambda )
00768 {
00769    double pp ;
00770    pp = 0.5 * poisson_t2p( xx , lambda ) ;
00771    return qginv(pp) ;
00772 }
00773 
00774 double poisson_p2t( double qq , double lambda )
00775 {
00776    int which , status ;
00777    double p,q, s,xlam,bound ;
00778 
00779    if( qq <= 0.0      ) return 999.9 ;
00780    if( qq >= 0.999999 ) return 0.0 ;
00781 
00782    which  = 2 ;
00783    p      = 1.0 - qq ;
00784    q      = qq ;
00785    s      = 0.0 ;
00786    xlam   = lambda ;
00787 
00788    cdfpoi( &which , &p , &q , &s , &xlam , &status , &bound ) ;
00789 
00790    return s ;
00791 }