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 /*----------------------------------------------------------------------
10 Inputs: (a,b,xc) for incomplete beta
11 Outputs:
12 Let Ipq = Int( x**(a-1)*(1-x)**(b-1)*ln(x)**p*ln(1-x)**q, x=0..xc ).
13 Then
14 bi7[0] = I00 = normalization factor
15 bi7[1] = I10/I00 = <ln(x)>
16 bi7[2] = I01/I00 = <ln(1-x)>
17 bi7[3] = d(bi7[1])/da = (I20*I00-I10**2)/I00**2
18 bi7[4] = d(bi7[1])/db = (I11*I00-I10*I01)/I00**2
19 bi7[5] = d(bi7[2])/da = (I11*I00-I10*I01)/I00**2
20 bi7[6] = d(bi7[2])/db = (I02*I00-I01**2)/I00**2
21 The integrals are calculated by transforming to y=a*ln(xc/x), and
22 then using Gauss-Laguerre quadrature:
23
24 Int( x**(a-1)*(1-x)**(b-1) * f(x) , x=0..xc )
25
26 transforms to
27
28 xc**a
29 ----- * Int( exp(-y)*(1-xc*exp(-y/a))**(b-1)*f(xc*exp(-y/a)), y=0..infty )
30 a
31
32 The return value of this function is -1 if an error occurred, and
33 is 0 if all is good.
34 -------------------------------------------------------------------------*/
35
bi7func(double a,double b,double xc,double * bi7)36 static int bi7func( double a , double b , double xc , double * bi7 )
37 {
38 #define NL 20 /* must be between 2 and 20 - see cs_laguerre.c */
39
40 static double *yy=NULL , *ww=NULL ;
41 double xx , s00,s10,s01,s20,s11,s02 , ff , l0,l1 ;
42 register int ii ;
43
44 if( a <= 0.0 || b <= 0.0 ||
45 xc <= 0.0 || xc >= 1.0 || bi7 == NULL ) return -1 ; /* sanity check */
46
47 /* initialize Laguerre integration table */
48
49 if( yy == NULL ) get_laguerre_table( NL , &yy , &ww ) ;
50
51 s00=s10=s01=s20=s11=s02 = 0.0 ;
52 for( ii=NL-1 ; ii >= 0 ; ii-- ){
53 xx = xc*exp(-yy[ii]/a) ; /* x transformed from y */
54 l0 = log(xx) ; l1 = log(1.0-xx) ; /* logarithms for Ipq sums */
55 ff = pow(1.0-xx,b-1.0) ; /* (1-x)**(b-1) */
56 s00 += ww[ii] * ff ; /* spq = Ipq sum */
57 s10 += ww[ii] * ff * l0 ;
58 s20 += ww[ii] * ff * l0 * l0 ;
59 s01 += ww[ii] * ff * l1 ;
60 s02 += ww[ii] * ff * l1 * l1 ;
61 s11 += ww[ii] * ff * l0 * l1 ;
62 }
63
64 if( s00 <= 0.0 ) return -1 ;
65
66 bi7[0] = s00 * pow(xc,a) / a ; /* normalizer */
67 bi7[1] = s10/s00 ; /* R0 */
68 bi7[2] = s01/s00 ; /* R1 */
69 bi7[3] = (s20*s00-s10*s10)/(s00*s00) ; /* dR0/da */
70 bi7[4] = (s11*s00-s10*s01)/(s00*s00) ; /* dR0/db */
71 bi7[5] = bi7[4] ; /* dR1/da */
72 bi7[6] = (s02*s00-s01*s01)/(s00*s00) ; /* dR1/db */
73
74 return 0 ;
75 }
76
77 /*----------------------------------------------------------------------
78 Set the range of values to search for beta distribution fit
79 ------------------------------------------------------------------------*/
80
81 #define LL 0.2
82 #define UL 10000.0
83
84 static double AL = 0.21 ;
85 static double AU = 9.9 ;
86 static double BL = 5.9 ;
87 static double BU = 999.9 ;
88 static int NRAN = 6666 ;
89
betarange(double al,double au,double bl,double bu,int nran)90 static void betarange( double al,double au , double bl , double bu , int nran )
91 {
92 if( al > 0.0 ) AL = al ;
93 if( au > AL ) AU = au ;
94 if( bl > 0.0 ) BL = bl ;
95 if( bu > BL ) BU = bu ;
96 if( nran > 1 ) NRAN = nran ;
97 }
98
99 static double ainit=0.0 ;
100 static double binit=0.0 ;
101
beta_init(double ai,double bi)102 static void beta_init( double ai , double bi )
103 {
104 ainit = ai ; binit = bi ; return ;
105 }
106
107 /*--------------------------------------------------------------------
108 Solve the two equations
109 I10(a,b)/I00(a,b) = e0
110 I01(a,b)/I00(a,b) = e1
111 for (a,b), using 2D Newton's method.
112 ----------------------------------------------------------------------*/
113
betasolve(double e0,double e1,double xc,double * ap,double * bp)114 static int betasolve( double e0, double e1, double xc, double * ap, double * bp )
115 {
116 double bi7[7] , aa,bb , da,db , m11,m12,m21,m22 , r1,r2 , dd,ee ;
117 int nite=0 , ii,jj ;
118
119 if( ap == NULL || bp == NULL ||
120 xc <= 0.0 || xc >= 1.0 || e0 >= 0.0 || e1 >= 0.0 ) return -1 ;
121
122 /* randomly search for a good starting point */
123
124 dd = 1.e+20 ; aa = bb = 0.0 ;
125 if( ainit > 0.0 && binit > 0.0 ){
126 ii = bi7func( ainit , binit , xc , bi7 ) ;
127 if( ii == 0 ){
128 r1 = bi7[1]-e0; r2 = bi7[2]-e1; dd = fabs(r1/e0)+fabs(r2/e1);
129 aa = ainit ; bb = binit ;
130 }
131 }
132
133 dd = 1.e+20 ; aa = bb = 0.0 ;
134 for( jj=0 ; jj < NRAN ; jj++ ){
135 da = AL +(AU-AL) * drand48() ;
136 db = BL +(BU-BL) * drand48() ;
137 ii = bi7func( da , db , xc , bi7 ) ; if( ii ) continue ;
138 r1 = bi7[1]-e0; r2 = bi7[2]-e1; ee = fabs(r1/e0)+fabs(r2/e1);
139 if( ee < dd ){ aa=da ; bb=db ; dd=ee ; /*if(ee<0.05)break;*/ }
140 }
141 if( aa == 0.0 || bb == 0.0 ) return -1 ;
142 #if 0
143 fprintf(stderr,"%2d: aa=%15.10g bb=%15.10g ee=%g\n",nite,aa,bb,ee) ;
144 #endif
145
146 do{
147 ii = bi7func( aa , bb , xc , bi7 ) ;
148 if( ii ) return -1 ;
149 r1 = bi7[1] - e0 ;
150 r2 = bi7[2] - e1 ; ee = fabs(r1/e0) + fabs(r2/e1) ;
151 m11 = bi7[3] ; m12 = bi7[4] ; m21 = bi7[5] ; m22 = bi7[6] ;
152 dd = m11*m22 - m12*m21 ;
153 if( dd == 0.0 ) return -1 ;
154 da = ( m22*r1 - m12*r2 ) / dd ;
155 db = (-m21*r1 + m11*r2 ) / dd ;
156 nite++ ;
157 aa -= da ; bb -=db ;
158 #if 0
159 if( aa < LL ) aa = LL ; else if( aa > UL ) aa = UL ;
160 if( bb < LL ) bb = LL ; else if( bb > UL ) bb = UL ;
161
162 if( aa == LL || bb == LL || aa == UL || bb == UL ) return -1 ;
163 #else
164 if( aa < AL ) aa = AL ; else if( aa > AU ) aa = AU ;
165 if( bb < BL ) bb = BL ; else if( bb > BU ) bb = BU ;
166 #endif
167
168 } while( nite < 99 && fabs(da)+fabs(db) > 0.02 ) ;
169
170 *ap = aa ; *bp = bb ; return 0 ;
171 }
172
173 /*--------------------------------------------------------------------*/
174
175 typedef struct {
176 int mcount , ibot ;
177 float * bval , * cval ;
178 } BFIT_data ;
179
180 typedef struct {
181 int mgood , itop ;
182 float a,b,xcut,chisq,df_chisq,q_chisq,eps ;
183 } BFIT_result ;
184
BFIT_free_data(BFIT_data * bfd)185 void BFIT_free_data( BFIT_data * bfd )
186 {
187 if( bfd != NULL ){
188 if( bfd->bval != NULL ) free(bfd->bval) ;
189 if( bfd->cval != NULL ) free(bfd->cval) ;
190 free(bfd) ;
191 }
192 }
193
BFIT_free_result(BFIT_result * bfr)194 void BFIT_free_result( BFIT_result * bfr ){ if( bfr != NULL ) free(bfr); }
195
196 /*--------------------------------------------------------------------*/
197
BFIT_bootstrap_sample(BFIT_data * bfd)198 BFIT_data * BFIT_bootstrap_sample( BFIT_data * bfd )
199 {
200 BFIT_data * nfd ;
201 int ii , jj , mcount,ibot , nuse ;
202
203 if( bfd == NULL ) return NULL ;
204 mcount = bfd->mcount ;
205 ibot = bfd->ibot ;
206 nuse = mcount - ibot ;
207
208 nfd = (BFIT_data *) malloc(sizeof(BFIT_data)) ;
209
210 nfd->mcount = mcount ;
211 nfd->ibot = ibot ;
212 nfd->bval = (float *) malloc( sizeof(float) * mcount ) ;
213 if( bfd->cval != NULL )
214 nfd->cval = (float *) malloc( sizeof(float) * mcount ) ;
215 else
216 nfd->cval = NULL ;
217
218 for( ii=0 ; ii < ibot ; ii++ ){
219 nfd->bval[ii] = 0.0 ;
220 if( nfd->cval != NULL ) nfd->cval[ii] = 0.0 ;
221 }
222
223 for( ii=ibot ; ii < mcount ; ii++ ){
224 jj = ((lrand48()>>8) % nuse) + ibot ;
225 nfd->bval[ii] = bfd->bval[jj] ;
226 if( nfd->cval != NULL )
227 nfd->cval[ii] = bfd->cval[jj] ;
228 }
229
230 if( nfd->cval != NULL )
231 qsort_floatfloat( mcount , nfd->bval , nfd->cval ) ;
232 else
233 qsort_float( mcount , nfd->bval ) ;
234
235 return nfd ;
236 }
237
238 /*--------------------------------------------------------------------*/
239
BFIT_prepare_dataset(THD_3dim_dataset * input_dset,int ival,int sqr,THD_3dim_dataset * mask_dset,int miv,float mask_bot,float mask_top)240 BFIT_data * BFIT_prepare_dataset(
241 THD_3dim_dataset * input_dset , int ival , int sqr ,
242 THD_3dim_dataset * mask_dset , int miv ,
243 float mask_bot , float mask_top )
244 {
245 int mcount , ii,jj , nvox,ibot ;
246 byte * mmm ;
247 BFIT_data * bfd ;
248 float * bval , * cval ;
249
250 /* check inputs */
251
252 if( !ISVALID_DSET(input_dset) ||
253 ival < 0 ||
254 ival >= DSET_NVALS(input_dset) ) return NULL ;
255
256 nvox = DSET_NVOX(input_dset) ;
257
258 if( ISVALID_DSET(mask_dset) &&
259 (miv < 0 ||
260 miv >= DSET_NVALS(mask_dset) ||
261 DSET_NVOX(mask_dset) != nvox ) ) return NULL ;
262
263 DSET_load(input_dset) ;
264 if( DSET_ARRAY(input_dset,ival) == NULL ) return NULL ;
265
266 /* inputs are OK */
267
268 /*-- build a byte mask array --*/
269
270 if( mask_dset == NULL ){
271 mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
272 memset( mmm , 1, nvox ) ; mcount = nvox ;
273 } else {
274
275 mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
276 mcount = THD_countmask( nvox , mmm ) ;
277
278 if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
279 if( mcount < 999 ){
280 free(mmm) ;
281 fprintf(stderr,"*** BFIT_prepare_dataset:\n"
282 "*** only %d voxels survive the masking!\n",
283 mcount ) ;
284 return NULL ;
285 }
286 }
287
288 /*-- load values into bval --*/
289
290 bval = (float *) malloc( sizeof(float) * mcount ) ;
291
292 switch( DSET_BRICK_TYPE(input_dset,ival) ){
293
294 case MRI_short:{
295 short * bar = (short *) DSET_ARRAY(input_dset,ival) ;
296 float mfac = DSET_BRICK_FACTOR(input_dset,ival) ;
297 if( mfac == 0.0 ) mfac = 1.0 ;
298 for( ii=jj=0 ; ii < nvox ; ii++ )
299 if( mmm[ii] ) bval[jj++] = mfac*bar[ii] ;
300 }
301 break ;
302
303 case MRI_byte:{
304 byte * bar = (byte *) DSET_ARRAY(input_dset,ival) ;
305 float mfac = DSET_BRICK_FACTOR(input_dset,ival) ;
306 if( mfac == 0.0 ) mfac = 1.0 ;
307 for( ii=jj=0 ; ii < nvox ; ii++ )
308 if( mmm[ii] ) bval[jj++] = mfac*bar[ii] ;
309 }
310 break ;
311
312 case MRI_float:{
313 float * bar = (float *) DSET_ARRAY(input_dset,ival) ;
314 float mfac = DSET_BRICK_FACTOR(input_dset,ival) ;
315 if( mfac == 0.0 ) mfac = 1.0 ;
316 for( ii=jj=0 ; ii < nvox ; ii++ )
317 if( mmm[ii] ) bval[jj++] = mfac*bar[ii] ;
318 }
319 break ;
320 }
321
322 free(mmm) ; DSET_unload(input_dset) ; /* don't need no more */
323
324 /* correlation coefficients must be squared prior to betafit, */
325 /* then R**2 values must be sorted. */
326
327 if( sqr ){
328 cval = (float *) malloc( sizeof(float) * mcount ) ;
329 for( ii=0 ; ii < mcount ; ii++ ){
330 cval[ii] = bval[ii] ; /* save cc values */
331 bval[ii] = bval[ii]*bval[ii] ;
332 }
333 qsort_floatfloat( mcount , bval , cval ) ;
334 } else { /* already squared */
335 cval = NULL ;
336 qsort_float( mcount , bval ) ;
337 }
338
339 /* check sorted values for legality */
340
341 if( bval[mcount-1] > 1.0 ){
342 free(bval) ; if(cval!=NULL) free(cval) ;
343 fprintf(stderr,"*** BFIT_prepare_dataset:\n"
344 "*** R**2 values > 1.0 exist in dataset!\n") ;
345 return NULL ;
346 }
347 if( bval[0] < 0.0 ){
348 free(bval) ; if(cval!=NULL) free(cval) ;
349 fprintf(stderr,"*** BFIT_prepare_dataset:\n"
350 "*** R**2 values < 0.0 exist in dataset!\n") ;
351 return NULL ;
352 }
353
354 /* find 1st bval > 0 [we don't use 0 values] */
355
356 for( ibot=0; ibot<mcount && bval[ibot]<=0.0; ibot++ ) ; /* nada */
357
358 /* make output structure */
359
360 bfd = (BFIT_data *) malloc( sizeof(BFIT_data) ) ;
361
362 bfd->mcount = mcount ;
363 bfd->ibot = ibot ;
364 bfd->bval = bval ;
365 bfd->cval = cval ;
366
367 return bfd ;
368 }
369
370 /*--------------------------------------------------------------------*/
371
BFIT_compute(BFIT_data * bfd,float pcut,float abot,float atop,float bbot,float btop,int nran,int nbin)372 BFIT_result * BFIT_compute( BFIT_data * bfd ,
373 float pcut ,
374 float abot , float atop ,
375 float bbot , float btop ,
376 int nran , int nbin )
377 {
378 BFIT_result * bfr ;
379
380 float eps,eps1 ;
381 float *bval , *cval ;
382 double e0,e1 , aa,bb,xc ;
383 double chq,ccc,cdf ;
384 int ihqbot,ihqtop ;
385 int mcount,mgood , ii,jj , ibot,itop , sqr ;
386 float hbot,htop,dbin ;
387 int * hbin, * jbin ;
388 MRI_IMAGE * flim ;
389
390 /* mangle inputs */
391
392 if( bfd == NULL ) return NULL ;
393 if( pcut < 20.0 || pcut > 99.0 ) return NULL ;
394 if( abot < 0.1 || abot >= atop ) return NULL ;
395 if( bbot < 9.9 || bbot >= btop ) return NULL ;
396
397 if( nran < 10 ) nran = 10 ;
398
399 mcount = bfd->mcount ;
400 ibot = bfd->ibot ;
401 bval = bfd->bval ;
402 cval = bfd->cval ; sqr = (cval != NULL) ;
403
404 /* now set the cutoff value (xc) */
405
406 itop = (int)( ibot + 0.01*pcut*(mcount-ibot) + 0.5 ) ;
407 mgood = itop - ibot ;
408 if( mgood < 999 ){
409 fprintf(stderr,"*** BFIT_compute: mgood=%d\n",mgood) ;
410 return NULL ;
411 }
412
413 xc = bval[itop-1] ;
414
415 /* compute the statistics of the values in (0,xc] */
416
417 e0 = e1 = 0.0 ;
418 for( ii=ibot ; ii < itop ; ii++ ){
419 e0 += log(bval[ii]) ; e1 += log(1.0-bval[ii]) ;
420 }
421 e0 /= mgood ; e1 /= mgood ;
422
423 /* and solve for the best fit parameters (aa,bb) */
424
425 betarange( abot , atop , bbot , btop , nran ) ;
426
427 ii = betasolve( e0,e1,xc , &aa,&bb );
428 if( ii < 0 ) return NULL ; /* error */
429
430 /*+++ At this point, could do some bootstrap to
431 estimate how good the estimates aa and bb are +++*/
432
433 /* estimate of outlier fraction */
434
435 eps1 = mgood / ( (mcount-ibot)*(1.0-beta_t2p(xc,aa,bb)) ) ;
436 eps = 1.0-eps1 ;
437 if( eps1 > 1.0 ) eps1 = 1.0 ;
438 eps1 = (mcount-ibot) * eps1 ;
439
440 /*-- compute histogram and chi-square --*/
441
442 if( nbin >= 100 ){ /* don't do it if nbin is too small */
443
444 #define NEW_HCQ
445 #ifdef NEW_HCQ /* use new method */
446
447 { float * xbin = (float *) malloc(sizeof(float)*nbin) ;
448
449 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */
450 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */
451
452 htop = 1.0 - beta_t2p(xc,aa,bb) ; /* CDF at top */
453 dbin = htop / nbin ; /* d(CDF) for each bin */
454 ii = rint( eps1 * dbin ) ;
455 for( jj=0 ; jj < nbin ; jj++ ){
456 xbin[jj] = beta_p2t( 1.0 - (jj+1)*dbin , aa , bb ) ;
457 jbin[jj] = ii ;
458 }
459 xbin[nbin-1] = xc ;
460
461 for( ii=ibot ; ii < mcount ; ii++ ){
462 for( jj=0 ; jj < nbin ; jj++ ){
463 if( bval[ii] <= xbin[jj] ){ hbin[jj]++ ; break ; }
464 }
465 }
466
467 free(xbin) ;
468
469 ihqbot = 0 ;
470 ihqtop = nbin-1 ;
471 }
472
473 #else /* use old method */
474
475 /* original data was already squared (e.g., R**2 values) */
476
477 if( !sqr ){
478 hbot = 0.0 ; htop = 1.001 * xc ;
479 dbin = (htop-hbot)/nbin ;
480
481 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */
482 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */
483
484 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */
485 jbin[ii] = (int)( eps1 * ( beta_t2p(hbot+ii*dbin,aa,bb)
486 -beta_t2p(hbot+ii*dbin+dbin,aa,bb) ) ) ;
487 }
488
489 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ;
490 mri_fix_data_pointer( bval+ibot , flim ) ;
491 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
492
493 ihqbot = 0 ;
494 ihqtop = rint( xc / dbin ) ;
495
496 } else { /* original data was not squared (e.g., correlations) */
497
498 double hb,ht ;
499 htop = sqrt(1.001*xc) ;
500 hbot = -htop ;
501 dbin = (htop-hbot)/nbin ;
502
503 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */
504 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */
505
506 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */
507 hb = hbot+ii*dbin ; ht = hb+dbin ;
508 hb = hb*hb ; ht = ht*ht ;
509 if( hb > ht ){ double qq=hb ; hb=ht ; ht=qq ; }
510 jbin[ii] = (int)( 0.5*eps1 * ( beta_t2p(hb,aa,bb)
511 -beta_t2p(ht,aa,bb) ) ) ;
512 }
513
514 ihqbot = rint( (-sqrt(xc) - hbot) / dbin ) ;
515 ihqtop = rint( ( sqrt(xc) - hbot) / dbin ) ;
516
517 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ;
518 mri_fix_data_pointer( cval+ibot , flim ) ;
519 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
520
521 }
522 #endif
523
524 /* compute upper-tail probability of chi-square */
525
526 chq = cdf = 0.0 ;
527 for( ii=ihqbot ; ii <= ihqtop ; ii++ ){
528 ccc = jbin[ii] ;
529 if( ccc > 1.0 ){
530 chq += SQR(hbin[ii]-ccc) / ccc ;
531 cdf++ ;
532 }
533 }
534 cdf -= 3.0 ;
535 ccc = chisq_t2p( chq , cdf ) ;
536
537 #ifndef NEW_HCQ
538 mri_clear_data_pointer(flim) ; mri_free(flim) ;
539 #endif
540 free(hbin) ; free(jbin) ;
541
542 } else {
543 chq = ccc = cdf = 0.0 ;
544 }
545
546 bfr = (BFIT_result *) malloc(sizeof(BFIT_result)) ;
547
548 bfr->mgood = mgood ;
549 bfr->itop = itop ;
550
551 bfr->a = aa ;
552 bfr->b = bb ;
553 bfr->xcut = xc ;
554 bfr->chisq = chq ;
555 bfr->q_chisq = ccc ;
556 bfr->df_chisq = cdf ;
557 bfr->eps = eps ;
558
559 return bfr ;
560 }
561