Doxygen Source Code Documentation
betafit.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Data Structures | |
struct | BFIT_data |
struct | BFIT_result |
Defines | |
#define | NL 20 |
#define | LL 0.2 |
#define | UL 10000.0 |
#define | NEW_HCQ |
Functions | |
int | bi7func (double a, double b, double xc, double *bi7) |
void | betarange (double al, double au, double bl, double bu, int nran) |
void | beta_init (double ai, double bi) |
int | betasolve (double e0, double e1, double xc, double *ap, double *bp) |
void | BFIT_free_data (BFIT_data *bfd) |
void | BFIT_free_result (BFIT_result *bfr) |
BFIT_data * | BFIT_bootstrap_sample (BFIT_data *bfd) |
BFIT_data * | 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) |
BFIT_result * | BFIT_compute (BFIT_data *bfd, float pcut, float abot, float atop, float bbot, float btop, int nran, int nbin) |
Variables | |
double | AL = 0.21 |
double | AU = 9.9 |
double | BL = 5.9 |
double | BU = 999.9 |
int | NRAN = 6666 |
double | ainit = 0.0 |
double | binit = 0.0 |
Define Documentation
|
Definition at line 81 of file betafit.c. Referenced by betasolve(). |
|
|
|
|
|
Definition at line 82 of file betafit.c. Referenced by betasolve(). |
Function Documentation
|
Definition at line 102 of file betafit.c. Referenced by process_sample().
|
|
Definition at line 90 of file betafit.c. References AL, AU, BL, BU, and NRAN. Referenced by BFIT_compute().
|
|
Definition at line 114 of file betafit.c. References ainit, AL, AU, bi7func(), binit, BL, BU, LL, NRAN, UL, and xc. Referenced by BFIT_compute().
00115 { 00116 double bi7[7] , aa,bb , da,db , m11,m12,m21,m22 , r1,r2 , dd,ee ; 00117 int nite=0 , ii,jj ; 00118 00119 if( ap == NULL || bp == NULL || 00120 xc <= 0.0 || xc >= 1.0 || e0 >= 0.0 || e1 >= 0.0 ) return -1 ; 00121 00122 /* randomly search for a good starting point */ 00123 00124 dd = 1.e+20 ; aa = bb = 0.0 ; 00125 if( ainit > 0.0 && binit > 0.0 ){ 00126 ii = bi7func( ainit , binit , xc , bi7 ) ; 00127 if( ii == 0 ){ 00128 r1 = bi7[1]-e0; r2 = bi7[2]-e1; dd = fabs(r1/e0)+fabs(r2/e1); 00129 aa = ainit ; bb = binit ; 00130 } 00131 } 00132 00133 dd = 1.e+20 ; aa = bb = 0.0 ; 00134 for( jj=0 ; jj < NRAN ; jj++ ){ 00135 da = AL +(AU-AL) * drand48() ; 00136 db = BL +(BU-BL) * drand48() ; 00137 ii = bi7func( da , db , xc , bi7 ) ; if( ii ) continue ; 00138 r1 = bi7[1]-e0; r2 = bi7[2]-e1; ee = fabs(r1/e0)+fabs(r2/e1); 00139 if( ee < dd ){ aa=da ; bb=db ; dd=ee ; /*if(ee<0.05)break;*/ } 00140 } 00141 if( aa == 0.0 || bb == 0.0 ) return -1 ; 00142 #if 0 00143 fprintf(stderr,"%2d: aa=%15.10g bb=%15.10g ee=%g\n",nite,aa,bb,ee) ; 00144 #endif 00145 00146 do{ 00147 ii = bi7func( aa , bb , xc , bi7 ) ; 00148 if( ii ) return -1 ; 00149 r1 = bi7[1] - e0 ; 00150 r2 = bi7[2] - e1 ; ee = fabs(r1/e0) + fabs(r2/e1) ; 00151 m11 = bi7[3] ; m12 = bi7[4] ; m21 = bi7[5] ; m22 = bi7[6] ; 00152 dd = m11*m22 - m12*m21 ; 00153 if( dd == 0.0 ) return -1 ; 00154 da = ( m22*r1 - m12*r2 ) / dd ; 00155 db = (-m21*r1 + m11*r2 ) / dd ; 00156 nite++ ; 00157 aa -= da ; bb -=db ; 00158 #if 0 00159 if( aa < LL ) aa = LL ; else if( aa > UL ) aa = UL ; 00160 if( bb < LL ) bb = LL ; else if( bb > UL ) bb = UL ; 00161 00162 if( aa == LL || bb == LL || aa == UL || bb == UL ) return -1 ; 00163 #else 00164 if( aa < AL ) aa = AL ; else if( aa > AU ) aa = AU ; 00165 if( bb < BL ) bb = BL ; else if( bb > BU ) bb = BU ; 00166 #endif 00167 00168 } while( nite < 99 && fabs(da)+fabs(db) > 0.02 ) ; 00169 00170 *ap = aa ; *bp = bb ; return 0 ; 00171 } |
|
Definition at line 198 of file betafit.c. References BFIT_data::bval, BFIT_data::cval, BFIT_data::ibot, malloc, BFIT_data::mcount, qsort_float(), and qsort_floatfloat(). Referenced by main().
00199 { 00200 BFIT_data * nfd ; 00201 int ii , jj , mcount,ibot , nuse ; 00202 00203 if( bfd == NULL ) return NULL ; 00204 mcount = bfd->mcount ; 00205 ibot = bfd->ibot ; 00206 nuse = mcount - ibot ; 00207 00208 nfd = (BFIT_data *) malloc(sizeof(BFIT_data)) ; 00209 00210 nfd->mcount = mcount ; 00211 nfd->ibot = ibot ; 00212 nfd->bval = (float *) malloc( sizeof(float) * mcount ) ; 00213 if( bfd->cval != NULL ) 00214 nfd->cval = (float *) malloc( sizeof(float) * mcount ) ; 00215 else 00216 nfd->cval = NULL ; 00217 00218 for( ii=0 ; ii < ibot ; ii++ ){ 00219 nfd->bval[ii] = 0.0 ; 00220 if( nfd->cval != NULL ) nfd->cval[ii] = 0.0 ; 00221 } 00222 00223 for( ii=ibot ; ii < mcount ; ii++ ){ 00224 jj = ((lrand48()>>8) % nuse) + ibot ; 00225 nfd->bval[ii] = bfd->bval[jj] ; 00226 if( nfd->cval != NULL ) 00227 nfd->cval[ii] = bfd->cval[jj] ; 00228 } 00229 00230 if( nfd->cval != NULL ) 00231 qsort_floatfloat( mcount , nfd->bval , nfd->cval ) ; 00232 else 00233 qsort_float( mcount , nfd->bval ) ; 00234 00235 return nfd ; 00236 } |
|
Definition at line 372 of file betafit.c. References BFIT_result::a, BFIT_result::b, beta_p2t(), beta_t2p(), betarange(), betasolve(), BFIT_data::bval, calloc, BFIT_result::chisq, chisq_t2p(), BFIT_data::cval, BFIT_result::df_chisq, BFIT_result::eps, free, BFIT_data::ibot, BFIT_result::itop, malloc, BFIT_data::mcount, BFIT_result::mgood, mri_clear_data_pointer, mri_fix_data_pointer(), mri_free(), mri_histogram(), mri_new_vol_empty(), BFIT_result::q_chisq, SQR, xc, and BFIT_result::xcut. Referenced by BFIT_main(), main(), and process_sample().
00377 { 00378 BFIT_result * bfr ; 00379 00380 float eps,eps1 ; 00381 float *bval , *cval ; 00382 double e0,e1 , aa,bb,xc ; 00383 double chq,ccc,cdf ; 00384 int ihqbot,ihqtop ; 00385 int mcount,mgood , ii,jj , ibot,itop , sqr ; 00386 float hbot,htop,dbin ; 00387 int * hbin, * jbin ; 00388 MRI_IMAGE * flim ; 00389 00390 /* mangle inputs */ 00391 00392 if( bfd == NULL ) return NULL ; 00393 if( pcut < 20.0 || pcut > 99.0 ) return NULL ; 00394 if( abot < 0.1 || abot >= atop ) return NULL ; 00395 if( bbot < 9.9 || bbot >= btop ) return NULL ; 00396 00397 if( nran < 10 ) nran = 10 ; 00398 00399 mcount = bfd->mcount ; 00400 ibot = bfd->ibot ; 00401 bval = bfd->bval ; 00402 cval = bfd->cval ; sqr = (cval != NULL) ; 00403 00404 /* now set the cutoff value (xc) */ 00405 00406 itop = (int)( ibot + 0.01*pcut*(mcount-ibot) + 0.5 ) ; 00407 mgood = itop - ibot ; 00408 if( mgood < 999 ){ 00409 fprintf(stderr,"*** BFIT_compute: mgood=%d\n",mgood) ; 00410 return NULL ; 00411 } 00412 00413 xc = bval[itop-1] ; 00414 00415 /* compute the statistics of the values in (0,xc] */ 00416 00417 e0 = e1 = 0.0 ; 00418 for( ii=ibot ; ii < itop ; ii++ ){ 00419 e0 += log(bval[ii]) ; e1 += log(1.0-bval[ii]) ; 00420 } 00421 e0 /= mgood ; e1 /= mgood ; 00422 00423 /* and solve for the best fit parameters (aa,bb) */ 00424 00425 betarange( abot , atop , bbot , btop , nran ) ; 00426 00427 ii = betasolve( e0,e1,xc , &aa,&bb ); 00428 if( ii < 0 ) return NULL ; /* error */ 00429 00430 /*+++ At this point, could do some bootstrap to 00431 estimate how good the estimates aa and bb are +++*/ 00432 00433 /* estimate of outlier fraction */ 00434 00435 eps1 = mgood / ( (mcount-ibot)*(1.0-beta_t2p(xc,aa,bb)) ) ; 00436 eps = 1.0-eps1 ; 00437 if( eps1 > 1.0 ) eps1 = 1.0 ; 00438 eps1 = (mcount-ibot) * eps1 ; 00439 00440 /*-- compute histogram and chi-square --*/ 00441 00442 if( nbin >= 100 ){ /* don't do it if nbin is too small */ 00443 00444 #define NEW_HCQ 00445 #ifdef NEW_HCQ /* use new method */ 00446 00447 { float * xbin = (float *) malloc(sizeof(float)*nbin) ; 00448 00449 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */ 00450 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */ 00451 00452 htop = 1.0 - beta_t2p(xc,aa,bb) ; /* CDF at top */ 00453 dbin = htop / nbin ; /* d(CDF) for each bin */ 00454 ii = rint( eps1 * dbin ) ; 00455 for( jj=0 ; jj < nbin ; jj++ ){ 00456 xbin[jj] = beta_p2t( 1.0 - (jj+1)*dbin , aa , bb ) ; 00457 jbin[jj] = ii ; 00458 } 00459 xbin[nbin-1] = xc ; 00460 00461 for( ii=ibot ; ii < mcount ; ii++ ){ 00462 for( jj=0 ; jj < nbin ; jj++ ){ 00463 if( bval[ii] <= xbin[jj] ){ hbin[jj]++ ; break ; } 00464 } 00465 } 00466 00467 free(xbin) ; 00468 00469 ihqbot = 0 ; 00470 ihqtop = nbin-1 ; 00471 } 00472 00473 #else /* use old method */ 00474 00475 /* original data was already squared (e.g., R**2 values) */ 00476 00477 if( !sqr ){ 00478 hbot = 0.0 ; htop = 1.001 * xc ; 00479 dbin = (htop-hbot)/nbin ; 00480 00481 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */ 00482 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */ 00483 00484 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */ 00485 jbin[ii] = (int)( eps1 * ( beta_t2p(hbot+ii*dbin,aa,bb) 00486 -beta_t2p(hbot+ii*dbin+dbin,aa,bb) ) ) ; 00487 } 00488 00489 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ; 00490 mri_fix_data_pointer( bval+ibot , flim ) ; 00491 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ; 00492 00493 ihqbot = 0 ; 00494 ihqtop = rint( xc / dbin ) ; 00495 00496 } else { /* original data was not squared (e.g., correlations) */ 00497 00498 double hb,ht ; 00499 htop = sqrt(1.001*xc) ; 00500 hbot = -htop ; 00501 dbin = (htop-hbot)/nbin ; 00502 00503 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */ 00504 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */ 00505 00506 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */ 00507 hb = hbot+ii*dbin ; ht = hb+dbin ; 00508 hb = hb*hb ; ht = ht*ht ; 00509 if( hb > ht ){ double qq=hb ; hb=ht ; ht=qq ; } 00510 jbin[ii] = (int)( 0.5*eps1 * ( beta_t2p(hb,aa,bb) 00511 -beta_t2p(ht,aa,bb) ) ) ; 00512 } 00513 00514 ihqbot = rint( (-sqrt(xc) - hbot) / dbin ) ; 00515 ihqtop = rint( ( sqrt(xc) - hbot) / dbin ) ; 00516 00517 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ; 00518 mri_fix_data_pointer( cval+ibot , flim ) ; 00519 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ; 00520 00521 } 00522 #endif 00523 00524 /* compute upper-tail probability of chi-square */ 00525 00526 chq = cdf = 0.0 ; 00527 for( ii=ihqbot ; ii <= ihqtop ; ii++ ){ 00528 ccc = jbin[ii] ; 00529 if( ccc > 1.0 ){ 00530 chq += SQR(hbin[ii]-ccc) / ccc ; 00531 cdf++ ; 00532 } 00533 } 00534 cdf -= 3.0 ; 00535 ccc = chisq_t2p( chq , cdf ) ; 00536 00537 #ifndef NEW_HCQ 00538 mri_clear_data_pointer(flim) ; mri_free(flim) ; 00539 #endif 00540 free(hbin) ; free(jbin) ; 00541 00542 } else { 00543 chq = ccc = cdf = 0.0 ; 00544 } 00545 00546 bfr = (BFIT_result *) malloc(sizeof(BFIT_result)) ; 00547 00548 bfr->mgood = mgood ; 00549 bfr->itop = itop ; 00550 00551 bfr->a = aa ; 00552 bfr->b = bb ; 00553 bfr->xcut = xc ; 00554 bfr->chisq = chq ; 00555 bfr->q_chisq = ccc ; 00556 bfr->df_chisq = cdf ; 00557 bfr->eps = eps ; 00558 00559 return bfr ; 00560 } |
|
Definition at line 185 of file betafit.c. References BFIT_data::bval, BFIT_data::cval, and free. Referenced by BFIT_main(), and main().
|
|
Definition at line 194 of file betafit.c. References free. Referenced by BFIT_main(), main(), and process_sample().
00194 { if( bfr != NULL ) free(bfr); } |
|
Definition at line 240 of file betafit.c. References BFIT_data::bval, BFIT_data::cval, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_load, DSET_NVALS, DSET_NVOX, DSET_unload, EQUIV_DSETS, free, BFIT_data::ibot, ISVALID_DSET, malloc, BFIT_data::mcount, mmm, qsort_float(), qsort_floatfloat(), THD_countmask(), and THD_makemask(). Referenced by BFIT_main(), and main().
00244 { 00245 int mcount , ii,jj , nvox,ibot ; 00246 byte * mmm ; 00247 BFIT_data * bfd ; 00248 float * bval , * cval ; 00249 00250 /* check inputs */ 00251 00252 if( !ISVALID_DSET(input_dset) || 00253 ival < 0 || 00254 ival >= DSET_NVALS(input_dset) ) return NULL ; 00255 00256 nvox = DSET_NVOX(input_dset) ; 00257 00258 if( ISVALID_DSET(mask_dset) && 00259 (miv < 0 || 00260 miv >= DSET_NVALS(mask_dset) || 00261 DSET_NVOX(mask_dset) != nvox ) ) return NULL ; 00262 00263 DSET_load(input_dset) ; 00264 if( DSET_ARRAY(input_dset,ival) == NULL ) return NULL ; 00265 00266 /* inputs are OK */ 00267 00268 /*-- build a byte mask array --*/ 00269 00270 if( mask_dset == NULL ){ 00271 mmm = (byte *) malloc( sizeof(byte) * nvox ) ; 00272 memset( mmm , 1, nvox ) ; mcount = nvox ; 00273 } else { 00274 00275 mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ; 00276 mcount = THD_countmask( nvox , mmm ) ; 00277 00278 if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ; 00279 if( mcount < 999 ){ 00280 free(mmm) ; 00281 fprintf(stderr,"*** BFIT_prepare_dataset:\n" 00282 "*** only %d voxels survive the masking!\n", 00283 mcount ) ; 00284 return NULL ; 00285 } 00286 } 00287 00288 /*-- load values into bval --*/ 00289 00290 bval = (float *) malloc( sizeof(float) * mcount ) ; 00291 00292 switch( DSET_BRICK_TYPE(input_dset,ival) ){ 00293 00294 case MRI_short:{ 00295 short * bar = (short *) DSET_ARRAY(input_dset,ival) ; 00296 float mfac = DSET_BRICK_FACTOR(input_dset,ival) ; 00297 if( mfac == 0.0 ) mfac = 1.0 ; 00298 for( ii=jj=0 ; ii < nvox ; ii++ ) 00299 if( mmm[ii] ) bval[jj++] = mfac*bar[ii] ; 00300 } 00301 break ; 00302 00303 case MRI_byte:{ 00304 byte * bar = (byte *) DSET_ARRAY(input_dset,ival) ; 00305 float mfac = DSET_BRICK_FACTOR(input_dset,ival) ; 00306 if( mfac == 0.0 ) mfac = 1.0 ; 00307 for( ii=jj=0 ; ii < nvox ; ii++ ) 00308 if( mmm[ii] ) bval[jj++] = mfac*bar[ii] ; 00309 } 00310 break ; 00311 00312 case MRI_float:{ 00313 float * bar = (float *) DSET_ARRAY(input_dset,ival) ; 00314 float mfac = DSET_BRICK_FACTOR(input_dset,ival) ; 00315 if( mfac == 0.0 ) mfac = 1.0 ; 00316 for( ii=jj=0 ; ii < nvox ; ii++ ) 00317 if( mmm[ii] ) bval[jj++] = mfac*bar[ii] ; 00318 } 00319 break ; 00320 } 00321 00322 free(mmm) ; DSET_unload(input_dset) ; /* don't need no more */ 00323 00324 /* correlation coefficients must be squared prior to betafit, */ 00325 /* then R**2 values must be sorted. */ 00326 00327 if( sqr ){ 00328 cval = (float *) malloc( sizeof(float) * mcount ) ; 00329 for( ii=0 ; ii < mcount ; ii++ ){ 00330 cval[ii] = bval[ii] ; /* save cc values */ 00331 bval[ii] = bval[ii]*bval[ii] ; 00332 } 00333 qsort_floatfloat( mcount , bval , cval ) ; 00334 } else { /* already squared */ 00335 cval = NULL ; 00336 qsort_float( mcount , bval ) ; 00337 } 00338 00339 /* check sorted values for legality */ 00340 00341 if( bval[mcount-1] > 1.0 ){ 00342 free(bval) ; if(cval!=NULL) free(cval) ; 00343 fprintf(stderr,"*** BFIT_prepare_dataset:\n" 00344 "*** R**2 values > 1.0 exist in dataset!\n") ; 00345 return NULL ; 00346 } 00347 if( bval[0] < 0.0 ){ 00348 free(bval) ; if(cval!=NULL) free(cval) ; 00349 fprintf(stderr,"*** BFIT_prepare_dataset:\n" 00350 "*** R**2 values < 0.0 exist in dataset!\n") ; 00351 return NULL ; 00352 } 00353 00354 /* find 1st bval > 0 [we don't use 0 values] */ 00355 00356 for( ibot=0; ibot<mcount && bval[ibot]<=0.0; ibot++ ) ; /* nada */ 00357 00358 /* make output structure */ 00359 00360 bfd = (BFIT_data *) malloc( sizeof(BFIT_data) ) ; 00361 00362 bfd->mcount = mcount ; 00363 bfd->ibot = ibot ; 00364 bfd->bval = bval ; 00365 bfd->cval = cval ; 00366 00367 return bfd ; 00368 } |
|
Definition at line 36 of file betafit.c. References a, get_laguerre_table(), and xc. Referenced by betasolve().
00037 { 00038 #define NL 20 /* must be between 2 and 20 - see cs_laguerre.c */ 00039 00040 static double *yy=NULL , *ww=NULL ; 00041 double xx , s00,s10,s01,s20,s11,s02 , ff , l0,l1 ; 00042 register int ii ; 00043 00044 if( a <= 0.0 || b <= 0.0 || 00045 xc <= 0.0 || xc >= 1.0 || bi7 == NULL ) return -1 ; /* sanity check */ 00046 00047 /* initialize Laguerre integration table */ 00048 00049 if( yy == NULL ) get_laguerre_table( NL , &yy , &ww ) ; 00050 00051 s00=s10=s01=s20=s11=s02 = 0.0 ; 00052 for( ii=NL-1 ; ii >= 0 ; ii-- ){ 00053 xx = xc*exp(-yy[ii]/a) ; /* x transformed from y */ 00054 l0 = log(xx) ; l1 = log(1.0-xx) ; /* logarithms for Ipq sums */ 00055 ff = pow(1.0-xx,b-1.0) ; /* (1-x)**(b-1) */ 00056 s00 += ww[ii] * ff ; /* spq = Ipq sum */ 00057 s10 += ww[ii] * ff * l0 ; 00058 s20 += ww[ii] * ff * l0 * l0 ; 00059 s01 += ww[ii] * ff * l1 ; 00060 s02 += ww[ii] * ff * l1 * l1 ; 00061 s11 += ww[ii] * ff * l0 * l1 ; 00062 } 00063 00064 if( s00 <= 0.0 ) return -1 ; 00065 00066 bi7[0] = s00 * pow(xc,a) / a ; /* normalizer */ 00067 bi7[1] = s10/s00 ; /* R0 */ 00068 bi7[2] = s01/s00 ; /* R1 */ 00069 bi7[3] = (s20*s00-s10*s10)/(s00*s00) ; /* dR0/da */ 00070 bi7[4] = (s11*s00-s10*s01)/(s00*s00) ; /* dR0/db */ 00071 bi7[5] = bi7[4] ; /* dR1/da */ 00072 bi7[6] = (s02*s00-s01*s01)/(s00*s00) ; /* dR1/db */ 00073 00074 return 0 ; 00075 } |
Variable Documentation
|
Definition at line 99 of file betafit.c. Referenced by beta_init(), and betasolve(). |
|
Definition at line 84 of file betafit.c. Referenced by betarange(), and betasolve(). |
|
Definition at line 85 of file betafit.c. Referenced by betarange(), and betasolve(). |
|
Definition at line 100 of file betafit.c. Referenced by beta_init(), and betasolve(). |
|
Definition at line 86 of file betafit.c. Referenced by betarange(), and betasolve(). |
|
Definition at line 87 of file betafit.c. Referenced by betarange(), and betasolve(). |
|
Definition at line 88 of file betafit.c. Referenced by betarange(), and betasolve(). |