00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
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
00066
00067
00068
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
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
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
00126
00127 gam2 = 6.0 / ( (dof-4.0) * nn ) ;
00128 gam4 = 240.0 / ( (dof-6.0) * (dof-4.0) * nn * nn ) ;
00129
00130
00131
00132 xx = qginv( 0.5 * pp ) ;
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
00162
00163
00164
00165
00166
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 ;
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
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
00220
00221
00222
00223
00224 double lnbeta( double p , double q )
00225 {
00226 return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
00227 }
00228
00229
00230
00231
00232
00233
00234
00235
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 ;
00248
00249 if( x <= ZERO ) return ZERO ;
00250 if( x >= ONE ) return ONE ;
00251
00252
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
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
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
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
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
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
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
00440
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 ;
00456
00457 dp = (p <= 0.5) ? (p) : (1.0-p) ;
00458
00459 if( dp <= 1.e-37 ){
00460 dx = 13.0 ;
00461 return ( (p <= 0.5) ? (dx) : (-dx) ) ;
00462 }
00463
00464
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
00472
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) ) ;
00482 }
00483
00484
00485
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 ;
00505 else return 1.0 ;
00506 }
00507
00508
00509
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 ;
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
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
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
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
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
00635
00636
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
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
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
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
00744
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
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 }