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   Outputs: ihi[0..*nhi-1] = indexes of highly correlations
00023              from rr.  ihi must be declared at least length
00024              nr to allow for the worst case.
00025 --------------------------------------------------------------*/
00026 
00027 void find_unusual_correlations( int nr    , float * rr ,
00028                                 int * nhi , int * ihi   )
00029 {
00030    int   ii , nzero,mzero ;
00031    float zmid,zsig,zmed , rmid,rcut ;
00032    float * zz , * aa ;
00033    int   * iz ;
00034 
00035    if( nhi == NULL ) return ;  /* illegal entry */
00036 
00037    *nhi = 0 ;                  /* default return */
00038 
00039    if( nr < 1000 || rr == NULL || ihi == NULL ) return ;  /* illegal entry */
00040 
00041    /*-- make workspace --*/
00042 
00043    zz = (float *) malloc(sizeof(float)*nr*2) ; aa = zz + nr ;
00044    iz = (int *)   malloc(sizeof(int)  *nr  ) ;
00045 
00046    if( zstar <= 0.0 ){               /* initialize zstar if needed */
00047       char * cp = getenv("PTAIL") ;
00048       float pp = 0.0001 ;
00049       if( cp != NULL ){
00050          float xx = strtod( cp , NULL ) ;
00051          if( xx > 0.0 && xx < 1.0 ) pp = xx ;
00052       }
00053       set_unusuality_tail( pp ) ;
00054    }
00055 
00056    /*-- copy data into workspace --*/
00057 
00058    memcpy( zz , rr , sizeof(float)*nr ) ;             /* the copy */
00059    for( ii=0 ; ii < nr ; ii++ ) iz[ii] = ii ;         /* tag location */
00060    qsort_floatint( nr , zz , iz ) ;                   /* sort now */
00061 
00062    /*- trim off 1's (perfect correlations) -*/
00063 
00064    for( ii=nr-1 ; ii > 0 && zz[ii] > 0.999 ; ii-- ) ;  /* nada */
00065    if( ii == 0 ){ free(zz); free(iz); return; }        /* shouldn't happen */
00066    nr = ii+1 ;                                         /* the trim */
00067 
00068    /*-- find median of atanh(rr) --*/
00069 
00070    if( nr%2 == 1 )              /* median */
00071       zmid = zz[nr/2] ;
00072    else
00073       zmid = 0.5 * ( zz[nr/2] + zz[nr/2-1] ) ;
00074 
00075    rmid = zmid ;        /* median of rr */
00076    zmid = atanh(zmid) ; /* median of atanh(rr) */
00077 
00078    /*-- aa = fabs( tanh( atanh(zz) - atanh(rmid) ) ) --*/
00079 
00080    for( ii=0 ; ii < nr ; ii++ )
00081       aa[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ;
00082 
00083    /*-- find MAD of atanh(rr) --*/
00084 
00085    qsort_float( nr , aa ) ;
00086    if( nr%2 == 1 )              /* MAD = median absolute deviation */
00087       zmed = aa[nr/2] ;
00088    else
00089       zmed = 0.5 * ( aa[nr/2] + aa[nr/2-1] ) ;
00090 
00091    zmed = atanh(zmed) ;         /* MAD of atanh(rr) */
00092    zsig = 1.4826 * zmed ;       /* estimate standard dev. of atanh(rr) */
00093                                 /* 1.4826 = 1/(sqrt(2)*erfinv(0.5))    */
00094 
00095    if( zsig <= 0.0 ){ free(zz); free(iz); return; } /* shouldn't happen */
00096 
00097    /*-- find values of correlation greater than threshold --*/
00098    /*-- that is, with (atanh(rr)-zmid)/zsig > zstar       --*/
00099 
00100    rcut = tanh( zsig*zstar + zmid ) ;
00101    for( ii=nr-1 ; ii > 0 ; ii-- ) if( zz[ii] < rcut ) break ;
00102 
00103    nzero = ii+1 ; mzero = nr - nzero ;
00104 
00105    *nhi = mzero ;
00106    for( ii=0 ; ii < mzero ; ii++ ) ihi[ii] = iz[nzero+ii] ;
00107 
00108    free(zz) ; free(iz) ; return ;
00109 }