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 }

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 }

00606 {
00607    double pp ;
00608    pp = 0.5 * beta_t2p( xx , aa , bb ) ;
00609    return qginv(pp) ;
00610 }

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 }

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 }

00662 {
00663    double pp ;
00664    pp = 0.5 * binomial_t2p( ss , ntrial , ptrial ) ;
00665    return qginv(pp) ;
00666 }

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 }

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 }

00573 {
00574    double pp ;
00575    pp = 0.5 * chisq_t2p( xx , dof ) ;
00576    return qginv(pp) ;
00577 }

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 }

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 }

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 }

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 }

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 }

00212 {
00213    double pp ;
00214    pp = 0.5 * fstat_t2p( ff , dofnum , dofden ) ;
00215    return qginv(pp) ;
00216 }

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 }

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 }

00716 {
00717    double pp ;
00718    pp = 0.5 * gamma_t2p( xx , sh , sc ) ;
00719    return qginv(pp) ;
00720 }

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 }

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 }

00225 {
00226    return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
00227 }

00446 {
00447   double v = qg(x) ;
00448   if( v > 0.0 ) return log10(v) ;
00449   return -99.99 ;
00450 }

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 }

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 }

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 }

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 }

00768 {
00769    double pp ;
00770    pp = 0.5 * poisson_t2p( xx , lambda ) ;
00771    return qginv(pp) ;
00772 }

00443 { return 0.5*erfc(x/1.414213562373095); }

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 }

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 }

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 }

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 }

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 }

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 }

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 }