1 #include "mrilib.h"
2
3 /*------------------------------------------------------------
4 Set the one-sided tail probability at which we will cutoff
5 the unusuality test.
6 --------------------------------------------------------------*/
7
8 static float zstar = 0.0 ; /* the actual cutoff */
9 static float pstar = 0.0 ; /* tail probability */
10
set_unusuality_tail(float p)11 void set_unusuality_tail( float p )
12 {
13 if( p > 0.0 && p < 1.0 ){
14 zstar = qginv(p) ;
15 pstar = p ;
16 }
17 return ;
18 }
19
20 #define NEAR1 0.999
21
22 /*------------------------------------------------------------
23 Inputs: rr[0..nr-1] = array of correlation coefficients.
24 Outputs: ihi[0..*nhi-1] = indexes of highly correlations
25 from rr. ihi must be declared at least length
26 nr to allow for the worst case.
27 --------------------------------------------------------------*/
28
find_unusual_correlations(int nrin,float * rr,int * nhi,int * ihi)29 void find_unusual_correlations( int nrin , float * rr ,
30 int * nhi , int * ihi )
31 {
32 register int ii,jj,nr ;
33 register float *zz ;
34 float zmid,zsig,zmed , rmid,rcut ;
35
36 if( nhi == NULL ) return ; /* illegal entry */
37
38 *nhi = 0 ; /* default return */
39
40 if( nrin < 1000 || rr == NULL || ihi == NULL ) return ; /* illegal entry */
41
42 /*-- make workspace --*/
43
44 zz = (float *) malloc(sizeof(float)*nrin) ;
45
46 if( zstar <= 0.0 ){ /* initialize zstar if needed */
47 char * cp = getenv("PTAIL") ;
48 float pp = 0.0001 ;
49 if( cp != NULL ){
50 float xx = strtod( cp , NULL ) ;
51 if( xx > 0.0 && xx < 1.0 ) pp = xx ;
52 }
53 set_unusuality_tail( pp ) ;
54 }
55
56 /*-- copy data into workspace (eliding 1's) --*/
57
58 for( ii=jj=0 ; ii < nrin ; ii++ )
59 if( rr[ii] <= NEAR1 ) zz[jj++] = rr[ii] ;
60 nr = jj ;
61 if( nr < 2 ){ free(zz); return; } /* shouldn't happen */
62
63 rmid = qmed_float( nr , zz ) ; /* median of correlations */
64 zmid = atanh(rmid) ; /* median of atanh(rr) */
65
66 /*-- zz <- fabs( tanh( atanh(zz) - atanh(rmid) ) ) --*/
67
68 for( ii=0 ; ii < nr ; ii++ )
69 zz[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ;
70
71 /*-- find MAD of atanh(rr) --*/
72
73 zmed = qmed_float( nr , zz ) ; /* MAD = median absolute deviation of rr */
74 zmed = atanh(zmed) ; /* MAD of atanh(rr) */
75 zsig = 1.4826 * zmed ; /* estimate standard dev. of atanh(rr) */
76 /* 1.4826 = 1/(sqrt(2)*erfinv(0.5)) */
77 free(zz) ;
78
79 if( zsig <= 0.0 ) return ; /* shouldn't happen */
80
81 /*-- find values of correlation greater than threshold --*/
82 /*-- that is, with (atanh(rr)-zmid)/zsig > zstar --*/
83
84 rcut = tanh( zsig*zstar + zmid ) ;
85
86 for( ii=jj=0 ; ii < nrin ; ii++ )
87 if( rr[ii] >= rcut && rr[ii] <= NEAR1 ) ihi[jj++] = ii ;
88
89 *nhi = jj ; return ;
90 }
91