Doxygen Source Code Documentation
uuu.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | SQRT_2PI 2.5066283 |
#define | PHI(s) (1.0-0.5*normal_t2p(ss)) |
Functions | |
void | set_unusuality_tail (float p) |
float | unusuality (int nr, float *rr) |
Variables | |
float | zstar = 0.0 |
float | pstar = 0.0 |
Define Documentation
|
|
|
|
Function Documentation
|
Definition at line 11 of file uuu.c. References p, pstar, qginv(), and zstar.
|
|
Definition at line 26 of file uuu.c. References free, getenv(), malloc, MAX, nr, pstar, qsort_float(), set_unusuality_tail(), strtod(), and zstar. Referenced by CORREL_main(), evaluate_bitvec(), and UC_unusuality().
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 } |
Variable Documentation
|
Definition at line 9 of file uuu.c. Referenced by set_unusuality_tail(), and unusuality(). |
|
Definition at line 8 of file uuu.c. Referenced by set_unusuality_tail(), and unusuality(). |