00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "mrilib.h"
00008 
00009 /*----------------------------------------------------------------------
00010    Inputs: (a,b,xc) for incomplete beta
00011    Outputs:
00012    Let Ipq = Int( x**(a-1)*(1-x)**(b-1)*ln(x)**p*ln(1-x)**q, x=0..xc ).
00013    Then
00014      bi7[0] = I00     = normalization factor
00015      bi7[1] = I10/I00 = <ln(x)>
00016      bi7[2] = I01/I00 = <ln(1-x)>
00017      bi7[3] = d(bi7[1])/da = (I20*I00-I10**2)/I00**2
00018      bi7[4] = d(bi7[1])/db = (I11*I00-I10*I01)/I00**2
00019      bi7[5] = d(bi7[2])/da = (I11*I00-I10*I01)/I00**2
00020      bi7[6] = d(bi7[2])/db = (I02*I00-I01**2)/I00**2
00021    The integrals are calculated by transforming to y=a*ln(xc/x), and
00022    then using Gauss-Laguerre quadrature:
00023 
00024    Int( x**(a-1)*(1-x)**(b-1) * f(x) , x=0..xc )
00025 
00026    transforms to
00027 
00028    xc**a
00029    ----- * Int( exp(-y)*(1-xc*exp(-y/a))**(b-1)*f(xc*exp(-y/a)), y=0..infty )
00030      a
00031 
00032    The return value of this function is -1 if an error occurred, and
00033    is 0 if all is good.
00034 -------------------------------------------------------------------------*/
00035 
00036 static int bi7func( double a , double b , double xc , double * bi7 )
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 }
00076 
00077 /*----------------------------------------------------------------------
00078    Set the range of values to search for beta distribution fit
00079 ------------------------------------------------------------------------*/
00080 
00081 #define LL   0.2
00082 #define UL   10000.0
00083 
00084 static double AL   = 0.21 ;
00085 static double AU   = 9.9 ;
00086 static double BL   = 5.9 ;
00087 static double BU   = 999.9 ;
00088 static int    NRAN = 6666 ;
00089 
00090 static void betarange( double al,double au , double bl , double bu , int nran )
00091 {
00092    if( al > 0.0 ) AL = al ;
00093    if( au > AL  ) AU = au ;
00094    if( bl > 0.0 ) BL = bl ;
00095    if( bu > BL  ) BU = bu ;
00096    if( nran > 1 ) NRAN = nran ;
00097 }
00098 
00099 static double ainit=0.0 ;
00100 static double binit=0.0 ;
00101 
00102 static void beta_init( double ai , double bi )
00103 {
00104    ainit = ai ; binit = bi ; return ;
00105 }
00106 
00107 /*--------------------------------------------------------------------
00108    Solve the two equations
00109      I10(a,b)/I00(a,b) = e0
00110      I01(a,b)/I00(a,b) = e1
00111    for (a,b), using 2D Newton's method.
00112 ----------------------------------------------------------------------*/
00113 
00114 static int betasolve( double e0, double e1, double xc, double * ap, double * bp )
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 }
00172 
00173 /*--------------------------------------------------------------------*/
00174 
00175 typedef struct {
00176   int mcount , ibot ;
00177   float * bval , * cval ;
00178 } BFIT_data ;
00179 
00180 typedef struct {
00181   int mgood , itop ;
00182   float a,b,xcut,chisq,df_chisq,q_chisq,eps ;
00183 } BFIT_result ;
00184 
00185 void BFIT_free_data( BFIT_data * bfd )
00186 {
00187    if( bfd != NULL ){
00188       if( bfd->bval != NULL ) free(bfd->bval) ;
00189       if( bfd->cval != NULL ) free(bfd->cval) ;
00190       free(bfd) ;
00191    }
00192 }
00193 
00194 void BFIT_free_result( BFIT_result * bfr ){ if( bfr != NULL ) free(bfr); }
00195 
00196 /*--------------------------------------------------------------------*/
00197 
00198 BFIT_data * BFIT_bootstrap_sample( BFIT_data * bfd )
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 }
00237 
00238 /*--------------------------------------------------------------------*/
00239 
00240 BFIT_data * BFIT_prepare_dataset(
00241                 THD_3dim_dataset * input_dset , int ival , int sqr ,
00242                 THD_3dim_dataset * mask_dset  , int miv  ,
00243                 float mask_bot  , float mask_top            )
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 }
00369 
00370 /*--------------------------------------------------------------------*/
00371 
00372 BFIT_result * BFIT_compute( BFIT_data * bfd ,
00373                             float pcut ,
00374                             float abot , float atop ,
00375                             float bbot , float btop ,
00376                             int   nran , int   nbin  )
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 }