00001 #include "mrilib.h"
00002 
00003 /*------------------------------------------------------------
00004   Set the one-sided tail probability at which we will cutoff
00005   the unusuality test.
00006 --------------------------------------------------------------*/
00007 
00008 static float zstar = 0.0 ;            /* the actual cutoff */
00009 static float pstar = 0.0 ;            /* tail probability  */
00010 
00011 void set_unusuality_tail( float p )
00012 {
00013    if( p > 0.0 && p < 1.0 ){
00014       zstar = qginv(p) ;
00015       pstar = p ;
00016    }
00017    return ;
00018 }
00019 
00020 /*------------------------------------------------------------
00021   Inputs: rr[0..nr-1] = array of correlation coefficients.
00022 --------------------------------------------------------------*/
00023 
00024 #undef TANHALL
00025 
00026 float unusuality( int nr , float * rr )
00027 {
00028    int ii , nzero , mzero ;
00029    float * zz , * aa ;
00030    float zmid,zsig,zmed, uval, fac, zrat, ff,fp, ss,ds,pp,ee , sigma ;
00031 #ifndef TANHALL
00032    float rmid , rcut ;
00033 #endif
00034 
00035    if( nr < 1000 || rr == NULL ) return 0.0 ;
00036 
00037    /*-- make workspace --*/
00038 
00039    zz = (float *) malloc(sizeof(float)*nr*2) ; aa = zz + nr ;
00040 
00041    if( zstar <= 0.0 ){
00042       char * cp = getenv("PTAIL") ;
00043       float pp = 0.01 ;
00044       if( cp != NULL ){
00045          float xx = strtod( cp , NULL ) ;
00046          if( xx > 0.0 && xx < 1.0 ) pp = xx ;
00047       }
00048       set_unusuality_tail( pp ) ;
00049    }
00050 
00051    /*-- copy data into workspace, converting to atanh --*/
00052 
00053    memcpy( zz , rr , sizeof(float)*nr ) ;
00054    qsort_float( nr , zz ) ;                           /* sort now */
00055 
00056    /*- trim off 1's (perfect correlations) -*/
00057 
00058    for( ii=nr-1 ; ii > 0 && zz[ii] > 0.999 ; ii-- ) ; /* nada */
00059    if( ii == 0 ){ free(zz) ; return 0.0 ; }           /* shouldn't happen */
00060    nr = ii+1 ;                                        /* the trim */
00061 
00062 #ifdef TANHALL
00063    for( ii=0 ; ii < nr ; ii++ ) zz[ii] = atanh(rr[ii]) ;
00064 #endif
00065 
00066    /*-- find median of zz [brute force sort] --*/
00067 
00068    if( nr%2 == 1 )              /* median */
00069       zmid = zz[nr/2] ;
00070    else
00071       zmid = 0.5 * ( zz[nr/2] + zz[nr/2-1] ) ;
00072 
00073 #ifdef TANHALL
00074    for( ii=0 ; ii < nr ; ii++ ) aa[ii] = fabs(zz[ii]-zmid) ;
00075 #else
00076    rmid = zmid ; zmid = atanh(zmid) ;
00077    for( ii=0 ; ii < nr ; ii++ )
00078       aa[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ;
00079 #endif
00080 
00081    /*-- find MAD of zz --*/
00082 
00083    qsort_float( nr , aa ) ;
00084    if( nr%2 == 1 )              /* MAD = median absolute deviation */
00085       zmed = aa[nr/2] ;
00086    else
00087       zmed = 0.5 * ( aa[nr/2] + aa[nr/2-1] ) ;
00088 
00089 #ifndef TANHALL
00090    zmed = atanh(zmed) ;
00091 #endif
00092 
00093    zsig = 1.4826 * zmed ;           /* estimate standard deviation of zz */
00094                                     /* 1/1.4826 = sqrt(2)*erfinv(0.5)    */
00095 
00096    if( zsig <= 0.0 ){               /* should not happen */
00097       free(zz) ; return 0.0 ;
00098    }
00099 
00100    /*-- normalize zz (is already sorted) --*/
00101    /*-- then, find values >= zstar       --*/
00102 
00103    fac = 1.0 / zsig ;
00104 #ifdef TANHALL
00105    for( ii=0 ; ii < nr ; ii++ ) zz[ii] = fac * ( zz[ii] - zmid ) ;
00106    for( ii=nr-1 ; ii > 0 ; ii-- ) if( zz[ii] < zstar ) break ;
00107    nzero = ii+1 ; mzero = nr - nzero ;
00108 #else
00109    rcut = tanh( zsig * zstar + zmid ) ;
00110    for( ii=nr-1 ; ii > 0 ; ii-- ){
00111       if( zz[ii] < rcut ) break ;
00112       else                zz[ii] = fac * ( atanh(zz[ii]) - zmid ) ;
00113    }
00114    nzero = ii+1 ; mzero = nr - nzero ;
00115 
00116 #if 0
00117    fprintf(stderr,"uuu: nr=%d rcut=%g mzero=%d\n",nr,rcut,mzero) ;
00118 #endif
00119 #endif
00120 
00121    if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){          /* too weird */
00122       free(zz) ; return 0.0 ;
00123    }
00124 
00125    /*-- compute sigma-tilde squared --*/
00126 
00127    zsig = 0.0 ;
00128    for( ii=nzero ; ii < nr ; ii++ ) zsig += zz[ii]*zz[ii] ;
00129    zsig = zsig / mzero ;
00130 
00131    /*-- set up to compute f(s) --*/
00132 
00133 #define SQRT_2PI 2.5066283
00134 
00135    zrat = zstar*zstar / zsig ;
00136    fac  = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
00137 
00138    ss   = zstar ;          /* initial guess for s = zstar/sigma */
00139 
00140    /*-- Newton's method [almost] --*/
00141 
00142 #undef  PHI
00143 #define PHI(s) (1.0-0.5*normal_t2p(ss))           /* N(0,1) cdf */
00144 
00145    for( ii=0 ; ii < 5 ; ii++ ){
00146       pp = PHI(ss) ;                              /* Phi(ss) \approx 1 */
00147       ee = exp(-0.5*ss*ss) ;
00148 
00149       ff = ss*ss - (fac/pp) * ss * ee - zrat ;    /* f(s) */
00150 
00151       fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */
00152 
00153       ds = ff / fp ;                              /* Newton step */
00154 
00155 #if 0
00156       fprintf(stderr,"Newton: ss=%g ds=%g ff=%g fp=%g pp=%g\n",ss,ds,ff,fp,pp) ;
00157 #endif
00158 
00159       ss -= ds ;                                  /* update */
00160    }
00161 
00162    sigma = zstar / ss ;                           /* actual estimate of sigma */
00163                                                   /* from the upper tail data */
00164 
00165    if( sigma <= 1.0 ){                            /* the boring case */
00166       free(zz) ; return 0.0 ;
00167    }
00168 
00169    /*-- compute the log-likelihood difference next --*/
00170 
00171    uval =  nzero * log( PHI(ss)/(1.0-pstar) )
00172          - mzero * ( log(sigma) + 0.5 * zsig * (1.0/(sigma*sigma)-1.0) ) ;
00173 
00174    /*-- done! --*/
00175 
00176    free(zz) ; return uval ;
00177 }