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 #undef USE_OLD_METHOD
10
11 #ifdef USE_OLD_METHOD
12 double stas4( double , double ) ;
13 double stinv( double , double ) ;
14 #else
15 double lnbeta( double , double ) ;
16 double incbeta( double,double,double,double ) ;
17 double incbeta_inverse( double,double,double,double ) ;
18 #endif
19
20 double qginv( double ) ; /* prototype */
21
main(int argc,char * argv[])22 int main( int argc , char * argv[] )
23 {
24 double pp , dof , tt , bb , binv ;
25 double nn,ll,mm ;
26 int quiet = 0 ;
27
28 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ){
29 printf("\n*** NOTE: This program has been superseded by program 'cdf' ***\n\n") ;
30
31 printf("Usage #1: p2t p dof\n"
32 " where p = double sided tail probability for t-distribution\n"
33 " dof = number of degrees of freedom to use\n"
34 " OUTPUT = t value that matches the input p\n"
35 "\n"
36 "Usage #2: p2t p N L M\n"
37 " where p = double sided tail probability of beta distribution\n"
38 " N = number of measured data points\n"
39 " L = number of nuisance parameters (orts)\n"
40 " M = number of fit parameters\n"
41 " OUTPUT = threshold for correlation coefficient\n"
42 "\n"
43 "Usage #3: p2t p\n"
44 " where p = one sided tail probability of Gaussian distribution\n"
45 " OUTPUT = z value for which P(x>z) = p\n"
46 "\n"
47 "Usage #4: p2t p dof N\n"
48 " where p = double sided tail probability for distribution of\n"
49 " the mean of N iid zero-mean t-variables\n"
50 " dof = number of degrees of freedom of each t-variable\n"
51 " N = number of t variables averaged\n"
52 " OUTPUT = threshold for the t average statistic\n"
53 " N.B.: The method used for this calculation is the Cornish-\n"
54 " Fisher expansion in N, and is only an approximation.\n"
55 " This also requires dof > 6, and the results will be\n"
56 " less accurate as dof approaches 6 from above!\n"
57 ) ;
58 exit(0) ;
59 }
60
61 if( strcmp(argv[1],"-q") == 0 ){
62 argv++ ;
63 argc-- ;
64 quiet = 1 ;
65 if( argc < 2 ) exit(0) ;
66 }
67
68 /** Gaussian statistic **/
69
70 if( argc == 2 ){
71 pp = strtod( argv[1] , NULL ) ;
72 if( !quiet && (pp < 1.e-20 || pp >= 0.9999999) ){
73 fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
74 }
75 tt = qginv(pp) ;
76 if( quiet ) printf("%g\n",tt) ;
77 else printf("p = %g z = %g\n",pp,tt) ;
78 }
79
80 /** t statistic **/
81
82 else if( argc == 3 ){
83 pp = strtod( argv[1] , NULL ) ;
84 dof = strtod( argv[2] , NULL ) ;
85 if( !quiet && (pp <= 1.e-10 || pp >= 0.9999999) ){
86 fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
87 }
88 if( dof <= 1.0 ){
89 fprintf(stderr,"command line dof value is illegal!\a\n") ; exit(-1) ;
90 }
91
92 #ifdef USE_OLD_METHOD
93 tt = stinv( pp , dof ) ;
94 #else
95 bb = lnbeta( 0.5*dof , 0.5 ) ;
96 binv = incbeta_inverse( pp, 0.5*dof , 0.5 , bb ) ;
97 tt = sqrt( dof*(1.0/binv-1.0) ) ;
98 #endif
99
100 if( quiet ) printf("%g\n",tt) ;
101 else printf("p = %g dof = %g t = %g\n",pp,dof,tt) ;
102
103 /** beta statistic **/
104
105 } else if( argc == 5 ){
106 pp = strtod( argv[1] , NULL ) ;
107 nn = strtod( argv[2] , NULL ) ;
108 ll = strtod( argv[3] , NULL ) ;
109 mm = strtod( argv[4] , NULL ) ;
110 if( !quiet && (pp <= 1.e-20 || pp >= 0.9999999) ){
111 fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
112 }
113 if( mm < 1.0 || ll < 0.0 || nn-ll-mm < 1.0 ){
114 fprintf(stderr,"command line N,L,M values are illegal!\a\n");exit(1) ;
115 }
116 bb = lnbeta( 0.5*mm , 0.5*(nn-ll-mm) ) ;
117 binv = incbeta_inverse( pp, 0.5*(nn-ll-mm) , 0.5*mm , bb ) ;
118 tt = sqrt(1.0-binv) ;
119 if( quiet ) printf("%g\n",tt) ;
120 else printf("p = %g N = %g L = %g M = %g rho = %g\n",
121 pp,nn,ll,mm,tt) ;
122
123 /** averaged t statistic **/
124
125 } else if( argc == 4 ){
126 double ww , xx , gam2,gam4 ;
127
128 pp = strtod( argv[1] , NULL ) ;
129 dof = strtod( argv[2] , NULL ) ;
130 nn = strtod( argv[3] , NULL ) ;
131
132 if( !quiet && (pp <= 1.e-10 || pp >= 0.9999999) ){
133 fprintf(stderr,"command line p value is illegal!\a\n") ; exit(-1) ;
134 }
135 if( dof <= 6.01 || nn < 1.0 ){
136 fprintf(stderr,"command line dof or N value is illegal!\a\n");exit(-1);
137 }
138
139 /* 4th and 6th order moments */
140
141 gam2 = 6.0 / ( (dof-4.0) * nn ) ;
142 gam4 = 240.0 / ( (dof-6.0) * (dof-4.0) * nn * nn ) ;
143
144 /* Cornish-Fisher expansion */
145
146 xx = qginv( 0.5 * pp ) ; /* Gaussian approx */
147
148 ww = xx + gam2 * xx * ( xx*xx - 3.0) / 24.0
149 + gam4 * xx * ( xx*xx*xx*xx - 10.0*xx*xx + 15.0) / 720.0
150 - gam2 * gam2 * xx * (3.0*xx*xx*xx*xx - 24.0*xx*xx + 29.0) / 384.0 ;
151
152 tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
153
154 if( quiet ) printf("%g\n",tt) ;
155 else{
156 printf("p = %g dof = %g N = %g 4-term t = %g",pp,dof,nn,tt) ;
157
158 ww = xx + gam2 * xx * ( xx*xx - 3.0) / 24.0 ;
159 tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
160 printf(" [2-term=%g",tt) ;
161 ww = xx ;
162 tt = sqrt( dof/(dof-2.0)/nn ) * ww ;
163 printf(" 1-term=%g]\n",tt) ;
164 }
165 }
166
167 exit(0) ;
168 }
169
170 #ifdef USE_OLD_METHOD
171 /*----------------------------------------------------------------------
172 code for inverse of central t distribution
173 Inputs: p = double sided tail probability
174 nu = degrees of freedom
175 Output: T such that P( |t| > T ) = p
176
177 This version is only good for nu >= 5, since it uses the
178 approximations in Abramowitz and Stegun, Eq. 26.7.5 (p. 949).
179 ------------------------------------------------------------------------*/
180
stinv(double p,double nu)181 double stinv( double p , double nu )
182 {
183 double xg , t4 ;
184 xg = qginv(0.5*p) ;
185 t4 = stas4( xg , nu ) ;
186 return t4 ;
187 }
188
stas4(double x,double nu)189 double stas4( double x , double nu) /* this code generated by Maple */
190 {
191 double t1,t2,t8,t9,t14,t17,t26,t34,t37 ;
192 t1 = x*x;
193 t2 = t1*x;
194 t8 = t1*t1;
195 t9 = t8*x;
196 t14 = nu*nu;
197 t17 = t8*t2;
198 t26 = t8*t8;
199 t34 = t14*t14;
200 t37 = x+(t2/4+x/4)/nu
201 +(5.0/96.0*t9+t2/6+x/32)/t14
202 +(t17/128+19.0/384.0*t9
203 +17.0/384.0*t2-5.0/128.0*x)/t14/nu
204 +(79.0/92160.0*t26*x+97.0/11520.0*t17+247.0/15360.0*t9
205 -t2/48-21.0/2048.0*x)/t34;
206 return t37 ;
207 }
208 #endif /* USE_OLD_METHOD */
209
210
211 #ifndef USE_OLD_METHOD
212
213 /*---------------------------------------------------------------
214 compute log of complete beta function, using the
215 Unix math library's log gamma function. If this is
216 not available, tough luck.
217 -----------------------------------------------------------------*/
218
lnbeta(double p,double q)219 double lnbeta( double p , double q )
220 {
221 return (lgamma(p) + lgamma(q) - lgamma(p+q)) ;
222 }
223
224 /*---------------------------------------------------------------
225 TRANSLATED FROM THE ORIGINAL FORTRAN:
226 algorithm as 63 appl. statist. (1973), vol.22, no.3
227
228 computes incomplete beta function ratio for arguments
229 x between zero and one, p and q positive.
230 log of complete beta function, beta, is assumed to be known
231 -----------------------------------------------------------------*/
232
233 #define ZERO 0.0
234 #define ONE 1.0
235 #define ACU 1.0e-15
236
incbeta(double x,double p,double q,double beta)237 double incbeta( double x , double p , double q , double beta )
238 {
239 double betain , psq , cx , xx,pp,qq , term,ai , temp , rx ;
240 int indx , ns ;
241
242 if( p <= ZERO || q <= ZERO || x < ZERO || x > ONE ) return -1.0 ;
243
244 if( x == ZERO ) return ZERO ;
245 if( x == ONE ) return ONE ;
246
247 /** change tail if necessary and determine s **/
248
249 psq = p+q ;
250 cx = ONE-x ;
251 if( p < psq*x ){
252 xx = cx ;
253 cx = x ;
254 pp = q ;
255 qq = p ;
256 indx = 1 ;
257 } else {
258 xx = x ;
259 pp = p ;
260 qq = q ;
261 indx = 0 ;
262 }
263
264 term = ONE ;
265 ai = ONE ;
266 betain = ONE ;
267 ns = qq + cx*psq ;
268
269 /** use soper's reduction formulae **/
270
271 rx = xx/cx ;
272
273 lab3:
274 temp = qq-ai ;
275 if(ns == 0) rx = xx ;
276
277 lab4:
278 term = term*temp*rx/(pp+ai) ;
279 betain = betain+term ;
280 temp = fabs(term) ;
281 if(temp <= ACU && temp <= ACU*betain) goto lab5 ;
282
283 ai = ai+ONE ;
284 ns = ns-1 ;
285 if(ns >= 0) goto lab3 ;
286 temp = psq ;
287 psq = psq+ONE ;
288 goto lab4 ;
289
290 lab5:
291 betain = betain*exp(pp*log(xx)+(qq-ONE)*log(cx)-beta)/pp ;
292 if(indx) betain=ONE-betain ;
293
294 return betain ;
295 }
296
297 /*----------------------------------------------------------------
298 algorithm as 109 appl. statist. (1977), vol.26, no.1
299 (replacing algorithm as 64 appl. statist. (1973),
300 vol.22, no.3)
301
302 Remark AS R83 and the correction in vol40(1) p.236 have been
303 incorporated in this version.
304
305 Computes inverse of the incomplete beta function
306 ratio for given positive values of the arguments
307 p and q, alpha between zero and one.
308 log of complete beta function, beta, is assumed to be known.
309 ------------------------------------------------------------------*/
310
311
312 #define SAE -15.0
313 #define TWO 2.0
314 #define THREE 3.0
315 #define FOUR 4.0
316 #define FIVE 5.0
317 #define SIX 6.0
318
319 #ifndef MAX
320 # define MAX(a,b) (((a)<(b)) ? (b) : (a))
321 # define MIN(a,b) (((a)>(b)) ? (b) : (a))
322 #endif
323
324
incbeta_inverse(double alpha,double p,double q,double beta)325 double incbeta_inverse( double alpha , double p , double q , double beta )
326 {
327 int indx , iex ;
328 double fpu , xinbta , a,pp,qq, r,y,t,s,h,w , acu ,
329 yprev,prev,sq , g,adj,tx,xin ;
330
331 fpu = pow(10.0,SAE) ;
332
333 if( p <= ZERO || q <= ZERO || alpha < ZERO || alpha > ONE ) return -1.0 ;
334
335 if( alpha == ZERO ) return ZERO ;
336 if( alpha == ONE ) return ONE ;
337
338 /** change tail if necessary **/
339
340 if( alpha > 0.5 ){
341 a = ONE-alpha ;
342 pp = q ;
343 qq = p ;
344 indx = 1 ;
345 } else {
346 a = alpha ;
347 pp = p ;
348 qq = q ;
349 indx = 0 ;
350 }
351
352 /** calculate the initial approximation **/
353
354 lab2:
355 r = sqrt(-log(a*a)) ;
356 y = r - (2.30753 + 0.27061*r) / (ONE+(0.99229+0.04481*r)*r) ;
357 if(pp > ONE && qq > ONE) goto lab5 ;
358
359 r = qq+qq ;
360 t = ONE/(9.0*qq) ;
361 t = r * pow( (ONE-t+y*sqrt(t)) , 3.0 ) ;
362 if( t <= ZERO ) goto lab3 ;
363
364 t = (FOUR*pp+r-TWO)/t ;
365 if( t <= ONE ) goto lab4 ;
366
367 xinbta = ONE-TWO/(t+ONE) ; goto lab6 ;
368
369 lab3:
370 xinbta = ONE-exp((log((ONE-a)*qq)+beta)/qq) ; goto lab6 ;
371
372 lab4:
373 xinbta = exp((log(a*pp)+beta)/pp) ; goto lab6 ;
374
375 lab5:
376 r = (y*y-THREE)/SIX ;
377 s = ONE/(pp+pp-ONE) ;
378 t = ONE/(qq+qq-ONE) ;
379 h = TWO/(s+t) ;
380 w = y*sqrt(h+r)/h-(t-s)*(r+FIVE/SIX-TWO/(THREE*h)) ;
381 xinbta = pp/(pp+qq*exp(w+w)) ;
382
383 /** solve for x by a modified newton-raphson method **/
384
385 lab6:
386 r = ONE-pp ;
387 t = ONE-qq ;
388 yprev = ZERO ;
389 sq = ONE ;
390 prev = ONE ;
391 if(xinbta < 0.0001) xinbta = 0.0001 ;
392 if(xinbta > 0.9999) xinbta = 0.9999 ;
393
394 #if 0
395 iex = -5.0 / (pp*pp) - 1.0/(a*a) - 13.0 ; if( iex < SAE ) iex = SAE ;
396 acu = pow(10.0,iex) ;
397 #else
398 acu = fpu ;
399 #endif
400
401 lab7:
402 y = incbeta( xinbta , pp,qq,beta ) ;
403 if( y < ZERO ) return -1.0 ;
404 xin = xinbta ;
405 y = (y-a)*exp(beta+r*log(xin)+t*log(ONE-xin)) ;
406 if(y*yprev <= ZERO) prev = MAX(sq, fpu) ;
407 g = ONE ;
408
409 lab9:
410 adj = g*y ;
411 sq = adj*adj ;
412 if(sq >= prev) goto lab10 ;
413 tx = xinbta-adj ;
414 if(tx >= ZERO && tx <= ONE) goto lab11 ;
415
416 lab10:
417 g = g/THREE ; goto lab9 ;
418
419 lab11:
420 if(tx == ZERO || tx == ONE ) goto lab10 ;
421 if(prev <= acu || y*y <= acu || fabs(xinbta-tx) < fpu) goto lab12 ;
422 xinbta = tx ;
423 yprev = y ;
424 goto lab7 ;
425
426 lab12:
427 xinbta = tx ;
428 if (indx) xinbta = ONE-xinbta ;
429 #if 0
430 printf("alpha = %g incbeta = %g\n",alpha, incbeta(xinbta,p,q,beta) );
431 #endif
432 return xinbta ;
433 }
434 #endif /* USE_OLD_METHOD */
435
436 /*** given p, return x such that Q(x)=p, for 0 < p < 1 ***/
437
qginv(double p)438 double qginv( double p )
439 {
440 double dp , dx , dt , ddq , dq ;
441 int newt ;
442
443 dp = (p <= 0.5) ? (p) : (1.0-p) ; /* make between 0 and 0.5 */
444
445 if( dp <= 0.0 ){
446 dx = 13.0 ;
447 return ( (p <= 0.5) ? (dx) : (-dx) ) ;
448 }
449
450 /** Step 1: use 26.2.23 from Abramowitz and Stegun **/
451
452 dt = sqrt( -2.0 * log(dp) ) ;
453 dx = dt
454 - ((.010328e+0*dt + .802853e+0)*dt + 2.525517e+0)
455 /(((.001308e+0*dt + .189269e+0)*dt + 1.432788e+0)*dt + 1.e+0) ;
456
457 /** Step 2: do 3 Newton steps to improve this **/
458
459 for( newt=0 ; newt < 3 ; newt++ ){
460 dq = 0.5e+0 * erfc( dx / 1.414213562373095e+0 ) - dp ;
461 ddq = exp( -0.5e+0 * dx * dx ) / 2.506628274631000e+0 ;
462 dx = dx + dq / ddq ;
463 }
464
465 return ( (p <= 0.5) ? (dx) : (-dx) ) ; /* return with correct sign */
466 }
467