1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*** Not actually MRI processing routines -- just some statistics ***/
10
11 /*-----------------------------------------------------------------
12 Student t-statistic:
13 given threshold, compute 2-sided tail probability ("t2p"), or
14 given statistic, compute equivalent N(0,1) deviate ("t2z"), or
15 given 2-sided tail probability, compute threshold ("p2t").
16 -------------------------------------------------------------------*/
17
student_p2t(double pp,double dof)18 double student_p2t( double pp , double dof )
19 {
20 double bb , binv , tt ;
21
22 if( pp <= 0.0 ) return 99.99 ;
23 if( pp >= 0.999999 ) return 0.0 ;
24 if( dof < 1.0 ) return 0.0 ;
25
26 bb = lnbeta( 0.5*dof , 0.5 ) ;
27 binv = incbeta_inverse( pp, 0.5*dof , 0.5 , bb ) ;
28 tt = sqrt( dof*(1.0/binv-1.0) ) ;
29 return tt ;
30 }
31
student_t2p(double tt,double dof)32 double student_t2p( double tt , double dof )
33 {
34 double bb , xx , pp ;
35
36 if( tt <= 0.0 || dof < 1.0 ) return 1.0 ;
37
38 bb = lnbeta( 0.5*dof , 0.5 ) ;
39 xx = dof/(dof + tt*tt) ;
40 pp = incbeta( xx , 0.5*dof , 0.5 , bb ) ;
41 return pp ;
42 }
43
student_t2z(double tt,double dof)44 double student_t2z( double tt , double dof )
45 {
46 static double bb , dof_old = -666.666 ;
47 double xx , pp ;
48
49 if( dof != dof_old ){
50 bb = lnbeta( 0.5*dof , 0.5 ) ;
51 dof_old = dof ;
52 }
53
54 xx = dof/(dof + tt*tt) ;
55 pp = incbeta( xx , 0.5*dof , 0.5 , bb ) ;
56
57 #if 1
58 if( tt > 0.0 ){
59 xx = qginv(0.5*pp) ;
60 } else {
61 xx = -qginv(0.5*pp) ;
62 }
63 return xx ;
64 #else
65 if( tt > 0.0 ) pp = 1.0 - 0.5 * pp ;
66 else pp = 0.5 * pp ;
67 xx = qginv(pp) ;
68 return -xx ;
69 #endif
70 }
71
72 /*---------------------------------------------------------------------
73 Correlation coefficient statistic:
74 nsam = # of samples used to compute rho ( > nfit+nort )
75 nfit = # of fitting parameters
76 nort = # of nuisance parameters ( nfit+nort >= 1 )
77 -----------------------------------------------------------------------*/
78
correl_p2t(double pp,double nsam,double nfit,double nort)79 double correl_p2t( double pp , double nsam , double nfit , double nort )
80 {
81 double bb , binv , rho ;
82
83 if( pp <= 0.0 ) return 0.999 ;
84 if( pp >= 0.999999 ) return 0.0 ;
85
86 if( nsam <= nfit+nort || nfit < 1.0 || nort < 1.0 ) return 0.0 ;
87
88 bb = lnbeta( 0.5*nfit , 0.5*(nsam-nfit-nort) ) ;
89 binv = incbeta_inverse( pp, 0.5*(nsam-nfit-nort) , 0.5*nfit , bb ) ;
90 rho = sqrt(1.0-binv) ;
91 return rho ;
92 }
93
correl_t2p(double rho,double nsam,double nfit,double nort)94 double correl_t2p( double rho , double nsam , double nfit , double nort )
95 {
96 double bb , xx , pp ;
97
98 if( rho <= 0.0 || nsam <= nfit+nort || nfit+nort < 1.0 ) return 1.0 ;
99
100 if( rho >= 0.9999999 ) return 0.0 ;
101
102 bb = lnbeta( 0.5*nfit , 0.5*(nsam-nfit-nort) ) ;
103 xx = 1.0 - rho*rho ;
104 pp = incbeta( xx , 0.5*(nsam-nfit-nort) , 0.5*nfit , bb ) ;
105 return pp ;
106 }
107
108 /******************************/
109 /*** added this 17 Sep 1998 ***/
110
correl_t2z(double rho,double nsam,double nfit,double nort)111 double correl_t2z( double rho , double nsam , double nfit , double nort )
112 {
113 double pp , xx ;
114 pp = 0.5 * correl_t2p( fabs(rho) , nsam , nfit , nort ) ;
115 xx = qginv(pp) ;
116 return ( (rho > 0) ? xx : -xx ) ;
117 }
118
119 /*----------------------------------------------------------
120 Averaged t-statistic
121 ------------------------------------------------------------*/
122
studave_p2t(double pp,double dof,double nn)123 double studave_p2t( double pp , double dof , double nn )
124 {
125 double ww , xx , gam2,gam4 , tt ;
126
127 if( pp <= 0.0 ) return 99.99 ;
128 if( pp >= 0.999999 ) return 0.0 ;
129
130 if( dof < 6.01 || nn < 1.0 ) return 0.0 ;
131
132 /* 4th and 6th order moments (or scaled cumulants) */
133
134 gam2 = 6.0 / ( (dof-4.0) * nn ) ;
135 gam4 = 240.0 / ( (dof-6.0) * (dof-4.0) * nn * nn ) ;
136
137 /* Cornish-Fisher expansion */
138
139 xx = qginv( 0.5 * pp ) ; /* Gaussian approx */
140
141 ww = xx + gam2 * xx * ( xx*xx - 3.0) / 24.0
142 + gam4 * xx * ( xx*xx*xx*xx - 10.0*xx*xx + 15.0) / 720.0
143 - gam2 * gam2 * xx * (3.0*xx*xx*xx*xx - 24.0*xx*xx + 29.0) / 384.0 ;
144
145 tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
146 return tt ;
147 }
148
studave_t2p(double tt,double dof,double nn)149 double studave_t2p( double tt , double dof , double nn )
150 {
151 static int nc = 0 ;
152 if( nc < 9 ){
153 fprintf(stderr,"*** studave_t2p: NOT IMPLEMENTED YET!\n") ; nc++ ;
154 }
155 return 0.0 ;
156 }
157
studave_t2z(double tt,double dof,double nn)158 double studave_t2z( double tt , double dof , double nn )
159 {
160 static int nc = 0 ;
161 if( nc < 3 ){
162 fprintf(stderr,"*** studave_t2z: NOT IMPLEMENTED YET!\n") ; nc++ ;
163 }
164 return 0.0 ;
165 }
166
167 /***********************************************************************/
168 /*** The routines below here are wrappers for the cdflib routines ***/
169 /*** (cdf_*.c) from U Texas -- see file cdflib.txt for the details. ***/
170 /***********************************************************************/
171
172 /*---------------------------------------------------------------
173 F statistic: single sided
174 -----------------------------------------------------------------*/
175
fstat_p2t(double pp,double dofnum,double dofden)176 double fstat_p2t( double pp , double dofnum , double dofden )
177 {
178 int which , status ;
179 double p , q , f , dfn , dfd , bound ;
180
181 if( pp <= 0.0 ) return 999.99 ;
182 if( pp >= 0.999999 ) return 0.0 ;
183
184 which = 2 ;
185 p = 1.0 - pp ; /* 20 Jan 1999: p and q were switched! */
186 q = pp ;
187 f = 0.0 ;
188 dfn = dofnum ;
189 dfd = dofden ;
190
191 cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
192
193 if( status == 0 ) return f ;
194 else return 0.0 ;
195 }
196
fstat_t2p(double ff,double dofnum,double dofden)197 double fstat_t2p( double ff , double dofnum , double dofden )
198 {
199 int which , status ;
200 double p , q , f , dfn , dfd , bound ;
201
202 if( ff <= 0.0 ) return 1.0 ;
203
204 which = 1 ;
205 p = 0.0 ;
206 q = 0.0 ;
207 f = ff ;
208 dfn = dofnum ;
209 dfd = dofden ;
210
211 cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
212
213 if( status == 0 ) return q ;
214 else return 1.0 ;
215 }
216
217 /******************************/
218 /*** added this 17 Sep 1998 ***/
219
fstat_t2z(double ff,double dofnum,double dofden)220 double fstat_t2z( double ff , double dofnum , double dofden )
221 {
222 double pp ;
223 pp = 0.5 * fstat_t2p( ff , dofnum , dofden ) ;
224 return qginv(pp) ;
225 }
226
227 /*---------------------------------------------------------------
228 compute log of complete beta function, using the
229 Unix math library's log gamma function. If this is
230 not available, see the end of this file.
231 -----------------------------------------------------------------*/
232
lnbeta(double p,double q)233 double lnbeta( double p , double q )
234 {
235 return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
236 }
237
238 /*---------------------------------------------------------------
239 TRANSLATED FROM THE ORIGINAL FORTRAN:
240 algorithm as 63 appl. statist. (1973), vol.22, no.3
241
242 computes incomplete beta function ratio for arguments
243 x between zero and one, p and q positive.
244 log of complete beta function, beta, is assumed to be known
245 -----------------------------------------------------------------*/
246
247 #define ZERO 0.0
248 #define ONE 1.0
249 #define ACU 1.0e-15
250
incbeta(double x,double p,double q,double beta)251 double incbeta( double x , double p , double q , double beta )
252 {
253 double betain , psq , cx , xx,pp,qq , term,ai , temp , rx ;
254 int indx , ns ;
255
256 if( p <= ZERO || q <= ZERO ) return -1.0 ; /* error! */
257
258 if( x <= ZERO ) return ZERO ;
259 if( x >= ONE ) return ONE ;
260
261 /** change tail if necessary and determine s **/
262
263 psq = p+q ;
264 cx = ONE-x ;
265 if( p < psq*x ){
266 xx = cx ;
267 cx = x ;
268 pp = q ;
269 qq = p ;
270 indx = 1 ;
271 } else {
272 xx = x ;
273 pp = p ;
274 qq = q ;
275 indx = 0 ;
276 }
277
278 term = ONE ;
279 ai = ONE ;
280 betain = ONE ;
281 ns = qq + cx*psq ;
282
283 /** use soper's reduction formulae **/
284
285 rx = xx/cx ;
286
287 lab3:
288 temp = qq-ai ;
289 if(ns == 0) rx = xx ;
290
291 lab4:
292 term = term*temp*rx/(pp+ai) ;
293 betain = betain+term ;
294 temp = fabs(term) ;
295 if(temp <= ACU && temp <= ACU*betain) goto lab5 ;
296
297 ai = ai+ONE ;
298 ns = ns-1 ;
299 if(ns >= 0) goto lab3 ;
300 temp = psq ;
301 psq = psq+ONE ;
302 goto lab4 ;
303
304 lab5:
305 betain = betain*exp(pp*log(xx)+(qq-ONE)*log(cx)-beta)/pp ;
306 if(indx) betain=ONE-betain ;
307
308 return betain ;
309 }
310
311 /*----------------------------------------------------------------
312 algorithm as 109 appl. statist. (1977), vol.26, no.1
313 (replacing algorithm as 64 appl. statist. (1973),
314 vol.22, no.3)
315
316 Remark AS R83 and the correction in vol40(1) p.236 have been
317 incorporated in this version.
318
319 Computes inverse of the incomplete beta function
320 ratio for given positive values of the arguments
321 p and q, alpha between zero and one.
322 log of complete beta function, beta, is assumed to be known.
323 ------------------------------------------------------------------*/
324
325 #define SAE -15.0
326 #define TWO 2.0
327 #define THREE 3.0
328 #define FOUR 4.0
329 #define FIVE 5.0
330 #define SIX 6.0
331
332 #ifndef MAX
333 # define MAX(a,b) (((a)<(b)) ? (b) : (a))
334 # define MIN(a,b) (((a)>(b)) ? (b) : (a))
335 #endif
336
incbeta_inverse(double alpha,double p,double q,double beta)337 double incbeta_inverse( double alpha , double p , double q , double beta )
338 {
339 int indx , iex ;
340 double fpu , xinbta , a,pp,qq, r,y,t,s,h,w , acu ,
341 yprev,prev,sq , g,adj,tx,xin ;
342
343 fpu = pow(10.0,SAE) ;
344
345 if( p <= ZERO || q <= ZERO || alpha < ZERO || alpha > ONE ) return -1.0 ;
346
347 if( alpha == ZERO ) return ZERO ;
348 if( alpha == ONE ) return ONE ;
349
350 /** change tail if necessary **/
351
352 if( alpha > 0.5 ){
353 a = ONE-alpha ;
354 pp = q ;
355 qq = p ;
356 indx = 1 ;
357 } else {
358 a = alpha ;
359 pp = p ;
360 qq = q ;
361 indx = 0 ;
362 }
363
364 /** calculate the initial approximation **/
365
366 lab2:
367 r = sqrt(-log(a*a)) ;
368 y = r - (2.30753 + 0.27061*r) / (ONE+(0.99229+0.04481*r)*r) ;
369 if(pp > ONE && qq > ONE) goto lab5 ;
370
371 r = qq+qq ;
372 t = ONE/(9.0*qq) ;
373 t = r * pow( (ONE-t+y*sqrt(t)) , 3.0 ) ;
374 if( t <= ZERO ) goto lab3 ;
375
376 t = (FOUR*pp+r-TWO)/t ;
377 if( t <= ONE ) goto lab4 ;
378
379 xinbta = ONE-TWO/(t+ONE) ; goto lab6 ;
380
381 lab3:
382 xinbta = ONE-exp((log((ONE-a)*qq)+beta)/qq) ; goto lab6 ;
383
384 lab4:
385 xinbta = exp((log(a*pp)+beta)/pp) ; goto lab6 ;
386
387 lab5:
388 r = (y*y-THREE)/SIX ;
389 s = ONE/(pp+pp-ONE) ;
390 t = ONE/(qq+qq-ONE) ;
391 h = TWO/(s+t) ;
392 w = y*sqrt(h+r)/h-(t-s)*(r+FIVE/SIX-TWO/(THREE*h)) ;
393 xinbta = pp/(pp+qq*exp(w+w)) ;
394
395 /** solve for x by a modified newton-raphson method **/
396
397 lab6:
398 r = ONE-pp ;
399 t = ONE-qq ;
400 yprev = ZERO ;
401 sq = ONE ;
402 prev = ONE ;
403 if(xinbta < 0.0001) xinbta = 0.0001 ;
404 if(xinbta > 0.9999) xinbta = 0.9999 ;
405
406 #if 0
407 iex = -5.0 / (pp*pp) - 1.0/(a*a) - 13.0 ; if( iex < SAE ) iex = SAE ;
408 acu = pow(10.0,iex) ;
409 #else
410 acu = fpu ;
411 #endif
412
413 lab7:
414 y = incbeta( xinbta , pp,qq,beta ) ;
415 if( y < ZERO ) return -1.0 ;
416 xin = xinbta ;
417 y = (y-a)*exp(beta+r*log(xin)+t*log(ONE-xin)) ;
418 if(y*yprev <= ZERO) prev = MAX(sq, fpu) ;
419 g = ONE ;
420
421 lab9:
422 adj = g*y ;
423 sq = adj*adj ;
424 if(sq >= prev) goto lab10 ;
425 tx = xinbta-adj ;
426 if(tx >= ZERO && tx <= ONE) goto lab11 ;
427
428 lab10:
429 g = g/THREE ; goto lab9 ;
430
431 lab11:
432 if(tx == ZERO || tx == ONE ) goto lab10 ;
433 if(prev <= acu || y*y <= acu || fabs(xinbta-tx) < fpu) goto lab12 ;
434 xinbta = tx ;
435 yprev = y ;
436 goto lab7 ;
437
438 lab12:
439 xinbta = tx ;
440 if (indx) xinbta = ONE-xinbta ;
441 #if 0
442 printf("alpha = %g incbeta = %g\n",alpha, incbeta(xinbta,p,q,beta) );
443 #endif
444 return xinbta ;
445 }
446
447 /*******************************************************************/
448 /**** Given p, return x such that Q(x)=p, for 0 < p < 1. ****/
449 /**** Q(x) = 1-P(x) = reversed cdf of N(0,1) variable. ****/
450 /*******************************************************************/
451
qg(double x)452 double qg( double x ){ return 0.5*erfc(x/1.414213562373095); }
453
log10qg(double x)454 double log10qg( double x )
455 {
456 double v = qg(x) ;
457 if( v > 0.0 ) return log10(v) ;
458 return -99.99 ;
459 }
460
qginv(double p)461 double qginv( double p )
462 {
463 double dp , dx , dt , ddq , dq ;
464 int newt ; /* not Gingrich, but Isaac */
465
466 dp = (p <= 0.5) ? (p) : (1.0-p) ; /* make between 0 and 0.5 */
467
468 if( dp <= 6.1172e-39 ){ /* cut off at 13 sigma */
469 dx = 13.0 ;
470 return ( (p <= 0.5) ? (dx) : (-dx) ) ;
471 }
472
473 /** Step 1: use 26.2.23 from Abramowitz and Stegun **/
474
475 dt = sqrt( -2.0 * log(dp) ) ;
476 dx = dt
477 - ((.010328*dt + .802853)*dt + 2.515517)
478 /(((.001308*dt + .189269)*dt + 1.432788)*dt + 1.) ;
479
480 /** Step 2: do 3 Newton steps to improve this
481 (uses the math library erfc function) **/
482
483 for( newt=0 ; newt < 3 ; newt++ ){
484 dq = 0.5 * erfc( dx / 1.414213562373095 ) - dp ;
485 ddq = exp( -0.5 * dx * dx ) / 2.506628274631000 ;
486 dx = dx + dq / ddq ;
487 }
488
489 if( dx > 13.0 ) dx = 13.0 ;
490 return ( (p <= 0.5) ? (dx) : (-dx) ) ; /* return with correct sign */
491 }
492
493 /*---------------------------------------------------------------
494 Compute double-sided tail probability for normal distribution.
495 -----------------------------------------------------------------*/
496
normal_t2p(double zz)497 double normal_t2p( double zz )
498 {
499 int which , status ;
500 double p , q , x , mean,sd,bound ;
501
502 if( zz == 0.0 ) return 1.0 ;
503 if( zz < 0.0 ) zz = -zz ; /* 19 Oct 2010 */
504
505 which = 1 ;
506 p = 0.0 ;
507 q = 0.0 ;
508 x = zz ;
509 mean = 0.0 ;
510 sd = 1.0 ;
511
512 cdfnor( &which , &p , &q , &x , &mean , &sd , &status , &bound ) ;
513
514 if( status == 0 ) return 2.0*q ; /* double sided prob = 2 times single sided */
515 else return 1.0 ;
516 }
517
518 /******************************/
519 /*** added this 17 Sep 1998 ***/
520
normal_p2t(double qq)521 double normal_p2t( double qq )
522 {
523 int which , status ;
524 double p , q , x , mean,sd,bound ;
525
526 if( qq <= 0.0 ) return 9.99 ;
527 if( qq >= 0.999999 ) return 0.0 ;
528
529 which = 2 ;
530 p = 1.0 - 0.5 * qq ;
531 q = 0.5 * qq ; /* single sided prob = 1/2 of double sided */
532 x = 0.0 ;
533 mean = 0.0 ;
534 sd = 1.0 ;
535
536 cdfnor( &which , &p , &q , &x , &mean , &sd , &status , &bound ) ;
537 return x ;
538 }
539
540 /*----------------------------------------------------------------
541 Compute single-sided tail probability for Chi-square
542 ------------------------------------------------------------------*/
543
chisq_t2p(double xx,double dof)544 double chisq_t2p( double xx , double dof )
545 {
546 int which , status ;
547 double p,q,x,df,bound ;
548
549 if( xx <= 0.0 ) return 1.0 ;
550
551 which = 1 ;
552 p = 0.0 ;
553 q = 0.0 ;
554 x = xx ;
555 df = dof ;
556
557 cdfchi( &which , &p , &q , &x , &df , &status , &bound ) ;
558
559 if( status == 0 ) return q ;
560 else return 1.0 ;
561 }
562
563 /******************************/
564 /*** added this 17 Sep 1998 ***/
565
chisq_p2t(double qq,double dof)566 double chisq_p2t( double qq , double dof )
567 {
568 int which , status ;
569 double p,q,x,df,bound ;
570
571 if( qq <= 0.0 ) return 999.9 ;
572 if( qq >= 0.999999 ) return 0.0 ;
573
574 which = 2 ;
575 p = 1.0 - qq ;
576 q = qq ;
577 x = 0.0 ;
578 df = dof ;
579
580 cdfchi( &which , &p , &q , &x , &df , &status , &bound ) ;
581 return x ;
582 }
583
chisq_t2z(double xx,double dof)584 double chisq_t2z( double xx , double dof )
585 {
586 double pp ;
587 pp = 0.5 * chisq_t2p( xx , dof ) ;
588 return qginv(pp) ;
589 }
590
591 /*----------------------------------------------------------------
592 Compute upper tail probability for incomplete beta distribution
593 ------------------------------------------------------------------*/
594
beta_t2p(double xx,double aa,double bb)595 double beta_t2p( double xx , double aa , double bb )
596 {
597 int which , status ;
598 double p,q,x,y,a,b,bound ;
599
600 if( xx <= 0.0 ) return 1.0 ;
601
602 which = 1 ;
603 p = 0.0 ;
604 q = 0.0 ;
605 x = xx ;
606 y = 1.0 - xx ;
607 a = aa ;
608 b = bb ;
609
610 cdfbet( &which , &p , &q , &x , &y , &a , &b , &status , &bound ) ;
611
612 if( status == 0 ) return q ;
613 else return 1.0 ;
614 }
615
616 /******************************/
617 /*** added this 17 Sep 1998 ***/
618
beta_t2z(double xx,double aa,double bb)619 double beta_t2z( double xx , double aa , double bb )
620 {
621 double pp ;
622 pp = 0.5 * beta_t2p( xx , aa , bb ) ;
623 return qginv(pp) ;
624 }
625
beta_p2t(double qq,double aa,double bb)626 double beta_p2t( double qq , double aa , double bb )
627 {
628 int which , status ;
629 double p,q,x,y,a,b,bound ;
630
631 if( qq <= 0.0 ) return 0.9999 ;
632 if( qq >= 0.999999 ) return 0.0 ;
633
634 which = 2 ;
635 p = 1.0 - qq ;
636 q = qq ;
637 x = 0.0 ;
638 y = 1.0 ;
639 a = aa ;
640 b = bb ;
641
642 cdfbet( &which , &p , &q , &x , &y , &a , &b , &status , &bound ) ;
643
644 return x ;
645 }
646
647 /*----------------------------------------------------------------
648 Compute upper tail probability for binomial distribution
649 (that is, the probability that more than ss out of ntrial
650 trials were successful).
651 ------------------------------------------------------------------*/
652
binomial_t2p(double ss,double ntrial,double ptrial)653 double binomial_t2p( double ss , double ntrial , double ptrial )
654 {
655 int which , status ;
656 double p,q, s,xn,pr,ompr,bound ;
657
658 which = 1 ;
659 p = 0.0 ;
660 q = 0.0 ;
661 s = ss ;
662 xn = ntrial ;
663 pr = ptrial ;
664 ompr = 1.0 - ptrial ;
665
666 cdfbin( &which , &p , &q , &s , &xn , &pr , &ompr , &status , &bound ) ;
667
668 if( status == 0 ) return q ;
669 else return 1.0 ;
670 }
671
672 /******************************/
673 /*** added this 17 Sep 1998 ***/
674
binomial_t2z(double ss,double ntrial,double ptrial)675 double binomial_t2z( double ss , double ntrial , double ptrial )
676 {
677 double pp ;
678 pp = 0.5 * binomial_t2p( ss , ntrial , ptrial ) ;
679 return qginv(pp) ;
680 }
681
binomial_p2t(double qq,double ntrial,double ptrial)682 double binomial_p2t( double qq , double ntrial , double ptrial )
683 {
684 int which , status ;
685 double p,q, s,xn,pr,ompr,bound ;
686
687 if( qq <= 0.0 ) return 99.99 ;
688 if( qq >= 0.999999 ) return 0.0 ;
689
690 which = 2 ;
691 p = 1.0 - qq ;
692 q = qq ;
693 s = 0.0 ;
694 xn = ntrial ;
695 pr = ptrial ;
696 ompr = 1.0 - ptrial ;
697
698 cdfbin( &which , &p , &q , &s , &xn , &pr , &ompr , &status , &bound ) ;
699
700 if( status == 0 ) return s ;
701 else return 0.0 ;
702 }
703
704 /*----------------------------------------------------------------
705 Compute upper tail probability for the gamma distribution.
706 ------------------------------------------------------------------*/
707
gamma_t2p(double xx,double sh,double sc)708 double gamma_t2p( double xx , double sh , double sc )
709 {
710 int which , status ;
711 double p,q, x,shape,scale,bound ;
712
713 if( xx <= 0.0 ) return 1.0 ;
714
715 which = 1 ;
716 p = 0.0 ;
717 q = 0.0 ;
718 x = xx ;
719 shape = sh ;
720 scale = sc ;
721
722 cdfgam( &which , &p , &q , &x , &shape , &scale , &status , &bound ) ;
723
724 if( status == 0 ) return q ;
725 else return 1.0 ;
726 }
727
728 /******************************/
729 /*** added this 17 Sep 1998 ***/
730
gamma_t2z(double xx,double sh,double sc)731 double gamma_t2z( double xx , double sh , double sc )
732 {
733 double pp ;
734 pp = 0.5 * gamma_t2p( xx , sh , sc ) ;
735 return qginv(pp) ;
736 }
737
gamma_p2t(double qq,double sh,double sc)738 double gamma_p2t( double qq , double sh , double sc )
739 {
740 int which , status ;
741 double p,q, x,shape,scale,bound ;
742
743 if( qq <= 0.0 ) return 999.9 ;
744 if( qq >= 0.999999 ) return 0.0 ;
745
746 which = 2 ;
747 p = 1.0 - qq ;
748 q = qq ;
749 x = 0.0 ;
750 shape = sh ;
751 scale = sc ;
752
753 cdfgam( &which , &p , &q , &x , &shape , &scale , &status , &bound ) ;
754
755 return x ;
756 }
757
758 /*----------------------------------------------------------------
759 Compute upper tail probability for the Poisson distribution
760 (that is, the probability that the value is greater than xx).
761 ------------------------------------------------------------------*/
762
poisson_t2p(double xx,double lambda)763 double poisson_t2p( double xx , double lambda )
764 {
765 int which , status ;
766 double p,q, s,xlam,bound ;
767
768 which = 1 ;
769 p = 0.0 ;
770 q = 0.0 ;
771 s = xx ;
772 xlam = lambda ;
773
774 cdfpoi( &which , &p , &q , &s , &xlam , &status , &bound ) ;
775
776 if( status == 0 ) return q ;
777 else return 1.0 ;
778 }
779
780 /******************************/
781 /*** added this 17 Sep 1998 ***/
782
poisson_t2z(double xx,double lambda)783 double poisson_t2z( double xx , double lambda )
784 {
785 double pp ;
786 pp = 0.5 * poisson_t2p( xx , lambda ) ;
787 return qginv(pp) ;
788 }
789
poisson_p2t(double qq,double lambda)790 double poisson_p2t( double qq , double lambda )
791 {
792 int which , status ;
793 double p,q, s,xlam,bound ;
794
795 if( qq <= 0.0 ) return 999.9 ;
796 if( qq >= 0.999999 ) return 0.0 ;
797
798 which = 2 ;
799 p = 1.0 - qq ;
800 q = qq ;
801 s = 0.0 ;
802 xlam = lambda ;
803
804 cdfpoi( &which , &p , &q , &s , &xlam , &status , &bound ) ;
805
806 return s ;
807 }
808