1 #include "mrilib.h"
2
3 static int acorr_demean = 1 ;
4
5 #undef PHI
6 #define PHI(x) (1.0-qg(x)) /* CDF of N(0,1) */
7
8 #undef PHINV
9 #define PHINV(x) qginv(1.0-(x)) /* inverse CDF of N(0,1) */
10
11 /*----------------------------------------------------------------------------*/
12 /*! atanh() of Pearson correlation of x[] and y[] */
13
acorrfun(int n,float * x,float * y)14 static INLINE float acorrfun( int n, float *x, float *y )
15 {
16 float xv=0.0f , yv=0.0f , xy=0.0f , vv,ww ;
17 float xm=0.0f , ym=0.0f ;
18 register int jj ;
19
20 if( acorr_demean ){
21 for( jj=0 ; jj < n ; jj++ ){ xm += x[jj]; ym += y[jj]; }
22 xm /= n ; ym /= n ;
23 }
24 for( jj=0 ; jj < n ; jj++ ){
25 vv = x[jj]-xm; ww = y[jj]-ym; xv += vv*vv; yv += ww*ww; xy += vv*ww;
26 }
27
28 if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
29 return atanhf(xy/sqrtf(xv*yv)) ;
30 }
31
32 /*----------------------------------------------------------------------------*/
33 /*! atanh() of Correlation of x[] and y[], omitting index qq.
34 For doing jackknife statistics on the correlation. */
35
acorrfun_jack(int n,float * x,float * y,int qq)36 static INLINE float acorrfun_jack( int n, float *x, float *y, int qq )
37 {
38 float xv=0.0f , yv=0.0f , xy=0.0f , vv,ww ;
39 float xm=0.0f , ym=0.0f ;
40 register int jj ;
41
42 if( acorr_demean ){
43 for( jj=0 ; jj < n ; jj++ ){ if( jj != qq ){ xm += x[jj]; ym += y[jj]; } }
44 xm /= (n-1.0f) ; ym /= (n-1.0f) ;
45 }
46 for( jj=0 ; jj < n ; jj++ ){
47 if( jj != qq ){
48 vv = x[jj]-xm; ww = y[jj]-ym; xv += vv*vv; yv += ww*ww; xy += vv*ww;
49 }
50 }
51
52 if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
53 return atanhf(xy/sqrtf(xv*yv)) ;
54 }
55
56 /*----------------------------------------------------------------------------*/
57
58 static float bca_alpha_1 = 0.05f ;
59 static float bca_phial_1 = 1.64485 ;
60
61 static float bca_alpha_2 = 0.1f ;
62 static float bca_phial_2 = 1.28155f ;
63
THD_bstrap_set_bca_alphas(float al1,float al2)64 void THD_bstrap_set_bca_alphas( float al1 , float al2 )
65 {
66 if( al1 < 0.01f ) al1 = 0.01f ; else if( al1 > 0.1f ) al1 = 0.1f ;
67 if( al2 <= 1.2f*al1 ) al2 = 2.0f*al1 ; else if( al2 > 0.2f ) al2 = 0.2f ;
68 bca_alpha_1 = al1 ; bca_phial_1 = (float)qginv(al1) ;
69 bca_alpha_2 = al2 ; bca_phial_2 = (float)qginv(al2) ;
70 }
71
72 /*----------------------------------------------------------------------------*/
73 /* Bias-Corrected accelerated (BCa) boostrap estimation of correlation
74 and confidence intervals: M Mudelsee, Mathematical Geology 35:651-665.
75 *//*--------------------------------------------------------------------------*/
76
THD_bstrap_corr(int npt,float * xx,float * yy,float tau,int nboot)77 float_triple THD_bstrap_corr( int npt, float *xx, float *yy, float tau, int nboot )
78 {
79 float_triple retval = {0.0f,0.0f,0.0f} ;
80 int kb , ii , qq ;
81 float rxy , pboot , *rboot , *xb , *yb ; int qboot ;
82 float rjbar , ahat , anum,aden , val , z0hat , alpha_b ;
83 float rb_al1 , rb_al2 , rb_au1 , rb_au2 , sig1,sig2,sig , med ;
84
85 ENTRY("THD_bstrap_corr") ;
86
87 if( npt < 19 || xx == NULL || yy == NULL ) RETURN(retval) ;
88
89 if( nboot < 666 ) nboot = 666 ;
90
91 #pragma omp critical (MALLOC)
92 { rboot = (float *)malloc(sizeof(float)*(nboot+2*npt)) ; }
93 xb = rboot + nboot ;
94 yb = xb + npt ;
95
96 rxy = acorrfun( npt , xx , yy ) ; /* "raw" correlation */
97
98 pboot = (tau <= 0.5f) ? 0.5f : 0.25f/tau ; /* prob of random jump */
99 qboot = (int)(pboot * ((1u << 31)-1u)) ; /* scaled to int range */
100
101 /** bootstrap work **/
102
103 for( kb=0 ; kb < nboot ; kb++ ){
104 qq = (lrand48() >> 3) % npt ; /* random starting point */
105 xb[0] = xx[qq] ; yb[0] = yy[qq] ;
106 for( ii=1 ; ii < npt ; ii++ ){
107 if( lrand48() < qboot ) qq = (lrand48() >> 3) % npt ; /* random jump */
108 else qq = (qq+1) % npt ; /* next point */
109 xb[ii] = xx[qq] ; yb[ii] = yy[qq] ;
110 }
111 rboot[kb] = acorrfun( npt , xb , yb ) ;
112 }
113
114 /* find location of rxy in the sorted rboot array,
115 then use that to estimate the bias correction parameter z0hat */
116
117 qsort_float( nboot , rboot ) ;
118 for( kb=0 ; kb < nboot && rboot[kb] < rxy ; kb++ ) ; /*nada*/
119 z0hat = PHINV(kb/(double)nboot) ; /* kb = number of rboot < rxy */
120 if( z0hat < -1.0f ) z0hat = -1.0f ; else if( z0hat > 1.0f ) z0hat = 1.0f ;
121
122 /** jackknife work to get "acceleration" parameter ahat **/
123
124 rjbar = 0.0f ;
125 for( kb=0 ; kb < npt ; kb++ ){
126 xb[kb] = acorrfun_jack( npt, xx,yy , kb ) ; rjbar += xb[kb] ;
127 }
128 rjbar /= npt ;
129 anum = aden = 0.0f ;
130 for( kb=0 ; kb < npt ; kb++ ){
131 val = rjbar - xb[kb] ; anum += val*val*val ; aden += val*val ;
132 }
133 ahat = (aden > 0.0f ) ? anum / (6.0f*aden*sqrtf(aden)) : 0.0f ;
134 if( ahat < -0.2f ) ahat = -0.2f ; else if( ahat > 0.2f ) ahat = 0.2f ;
135
136 /* compute bootstrap value at alpha_1 level, adjusted */
137
138 val = z0hat + (z0hat - bca_phial_1) / ( 1.0f - ahat*(z0hat - bca_phial_1) ) ;
139 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
140 rb_al1 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
141
142 /* compute bootstrap value at alpha_2 level, adjusted */
143
144 val = z0hat + (z0hat - bca_phial_2) / ( 1.0f - ahat*(z0hat - bca_phial_2) ) ;
145 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
146 rb_al2 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
147
148 /* compute bootstrap value at 1-alpha_2 level, adjusted */
149
150 val = z0hat + (z0hat + bca_phial_2) / ( 1.0f - ahat*(z0hat + bca_phial_2) ) ;
151 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
152 rb_au2 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
153
154 /* compute bootstrap value at 1-alpha_1 level, adjusted */
155
156 val = z0hat + (z0hat + bca_phial_1) / ( 1.0f - ahat*(z0hat + bca_phial_1) ) ;
157 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
158 rb_au1 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
159
160 /* compute bootstrap value at alpha_1 level, adjusted */
161
162 val = z0hat + (z0hat + bca_alpha_1) / ( 1.0f - ahat*(z0hat + bca_alpha_1) ) ;
163 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
164 rb_al1 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
165
166 /* compute bootstrap value at alpha_2 level, adjusted */
167
168 val = z0hat + (z0hat + bca_alpha_2) / ( 1.0f - ahat*(z0hat + bca_alpha_2) ) ;
169 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
170 rb_al2 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
171
172 /* compute bootstrap value at 1-alpha_1 level, adjusted */
173
174 val = z0hat + (z0hat + 1.0f-bca_alpha_1) / ( 1.0f - ahat*(z0hat + 1.0f-bca_alpha_1) ) ;
175 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
176 rb_au1 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
177
178 /* compute bootstrap value at 1-alpha_2 level, adjusted */
179
180 val = z0hat + (z0hat + 1.0f-bca_alpha_2) / ( 1.0f - ahat*(z0hat + 1.0f-bca_alpha_2) ) ;
181 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
182 rb_au2 = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
183
184 /* now estimate sigma both ways, then combine them */
185
186 sig1 = (rb_au1 - rb_al1) / (2.0f*bca_phial_1) ;
187 sig2 = (rb_au2 - rb_al2) / (2.0f*bca_phial_2) ;
188 sig = 0.5f * (sig1+sig2) ;
189
190 /* and estimate middle of distribution */
191
192 val = z0hat + (z0hat + 0.5f) / ( 1.0f - ahat*(z0hat + 0.5f) ) ;
193 val = nboot * PHI(val) ; kb = (int)val ; val = val-kb ;
194 med = (1.0f-val)*rboot[kb] + val*rboot[kb+1] ;
195
196 retval.a = med ; retval.b = sig ; retval.c = rxy ;
197
198 #if 0
199 /* convert rboot into alpha values */
200
201 for( kb=0 ; kb < nboot ; kb++ ){
202 phab = PHINV( (kb+0.5f)/nboot ) - z0hat ;
203 val = phab / ( 1.0f + aa*phab ) - z0hat ;
204 if( val < -5.0f ) val = -5.0f ; else if ( val > 5.0f ) val = 5.0f ;
205 phab = PHI(val) ;
206 printf("%f %f\n",rboot[kb],val) ;
207 }
208 #endif
209
210 /***/
211
212 #pragma omp critical (MALLOC)
213 { free(rboot) ; }
214 RETURN(retval) ;
215 }
216
217 /*----------------------------------------------------------------------------*/
218
219 #undef A
220 #define A(i,j) asym[(i)+(j)*nsym]
221 #undef XPT
222 #define XPT(q) ( (xtyp<=0) ? xx+(q)*nn : xar[q] )
223
224 /*----------------------------------------------------------------------------*/
225
mean_vector(int n,int m,int xtyp,void * xp,float * uvec)226 static void mean_vector( int n , int m , int xtyp , void *xp , float *uvec )
227 {
228 int nn=n , mm=m , jj ; register int ii ;
229 register float *xj , fac ; float *xx=NULL , **xar=NULL ;
230
231 if( nn < 1 || mm < 1 || xp == NULL || uvec == NULL ) return ;
232
233 if( xtyp <= 0 ) xx = (float * )xp ;
234 else xar = (float **)xp ;
235
236 for( ii=0 ; ii < nn ; ii++ ) uvec[ii] = 0.0f ;
237
238 for( jj=0 ; jj < mm ; jj++ ){
239 xj = XPT(jj) ;
240 for( ii=0 ; ii < nn ; ii++ ) uvec[ii] += xj[ii] ;
241 }
242
243 fac = 1.0f / nn ;
244 for( ii=0 ; ii < nn ; ii++ ) uvec[ii] *= fac ;
245 return ;
246 }
247
248 /*----------------------------------------------------------------------------*/
249 /*! Compute the principal singular vector of a set of m columns, each
250 of length n, stored in array xx[i+j*n] for i=0..n-1, j=0..m-1.
251 * The singular value is returned, and the vector is stored into uvec[].
252 * tvec is a vector so that the sign of uvec dot tvec will be non-negative.
253 * If the return value is not positive, something bad happened.
254 *//*--------------------------------------------------------------------------*/
255
principal_vector(int n,int m,int xtyp,void * xp,float * uvec,float * tvec)256 static float principal_vector( int n , int m , int xtyp , void *xp ,
257 float *uvec , float *tvec )
258 {
259 int nn=n , mm=m , nsym , ii,jj,kk,qq ;
260 double *asym , *deval ;
261 register double sum , qsum ; register float *xj , *xk ;
262 float sval , fsum , *xx=NULL , **xar=NULL ;
263
264 nsym = MIN(nn,mm) ; /* size of the symmetric matrix to create */
265
266 if( nsym < 1 || xp == NULL || uvec == NULL ) return -666.0f ;
267
268 #pragma omp critical (MALLOC)
269 { asym = (double *)malloc(sizeof(double)*nsym*nsym) ; /* symmetric matrix */
270 deval = (double *)malloc(sizeof(double)*nsym) ; /* its eigenvalues */
271 }
272
273 if( xtyp <= 0 ) xx = (float * )xp ;
274 else xar = (float **)xp ;
275
276 /** setup matrix to eigensolve: choose smaller of [X]'[X] and [X][X]' **/
277 /** since [X] is n x m, [X]'[X] is m x m and [X][X]' is n x n **/
278
279 if( nn > mm ){ /* more rows than columns: */
280 /* so [A] = [X]'[X] = m x m */
281 int n1 = nn-1 ;
282 for( jj=0 ; jj < mm ; jj++ ){
283 xj = XPT(jj) ;
284 for( kk=0 ; kk <= jj ; kk++ ){
285 sum = 0.0 ; xk = XPT(kk) ;
286 for( ii=0 ; ii < n1 ; ii+=2 ) sum += xj[ii]*xk[ii] + xj[ii+1]*xk[ii+1];
287 if( ii == n1 ) sum += xj[ii]*xk[ii] ;
288 A(jj,kk) = sum ; if( kk < jj ) A(kk,jj) = sum ;
289 }
290 }
291
292 } else { /* more columns than rows: */
293 /* so [A] = [X][X]' = n x n */
294 float *xt ; int m1=mm-1 ;
295 #pragma omp critical (MALLOC)
296 xt = (float *)malloc(sizeof(float)*nn*mm) ;
297 for( jj=0 ; jj < mm ; jj++ ){ /* form [X]' into array xt */
298 if( xtyp <= 0 )
299 for( ii=0 ; ii < nn ; ii++ ) xt[jj+ii*mm] = xx[ii+jj*nn] ;
300 else
301 for( ii=0 ; ii < nn ; ii++ ) xt[jj+ii*mm] = xar[jj][ii] ;
302 }
303
304 for( jj=0 ; jj < nn ; jj++ ){
305 xj = xt + jj*mm ;
306 for( kk=0 ; kk <= jj ; kk++ ){
307 sum = 0.0 ; xk = xt + kk*mm ;
308 for( ii=0 ; ii < m1 ; ii+=2 ) sum += xj[ii]*xk[ii] + xj[ii+1]*xk[ii+1];
309 if( ii == m1 ) sum += xj[ii]*xk[ii] ;
310 A(jj,kk) = sum ; if( kk < jj ) A(kk,jj) = sum ;
311 }
312 }
313
314 #pragma omp critical (MALLOC)
315 free(xt) ; /* don't need this no more */
316 }
317
318 /** compute the first eigenvector corresponding to largest eigenvalue **/
319 /** this eigenvector is stored on top of first column of asym **/
320
321 ii = symeig_irange( nsym, asym, deval, nsym-1, nsym-1, 0 ) ;
322
323 if( ii != 0 || deval[0] < 0.0 ){
324 #pragma omp critical (MALLOC)
325 { free(deval) ; free(asym) ; }
326 return -999.0f ; /* eigensolver failed!? */
327 }
328
329 /** Store singular value (sqrt of eigenvalue) **/
330
331 sval = (float)sqrt(deval[0]) ;
332
333 /** SVD is [X] = [U] [S] [V]', where [U] = desired output vectors
334
335 case n <= m: [A] = [X][X]' = [U] [S][S]' [U]'
336 so [A][U] = [U] [S][S]'
337 so eigenvectors of [A] are just [U]
338
339 case n > m: [A] = [X]'[X] = [V] [S]'[S] [V]'
340 so [A][V] = [V] [S'][S]
341 so eigenvectors of [A] are [V], but we want [U]
342 note that [X][V] = [U] [S]
343 so pre-multiplying each column vector in [V] by matrix [X]
344 will give the corresponding column in [U], but scaled;
345 below, just L2-normalize the column to get output vector **/
346
347 if( nn <= mm ){ /* copy eigenvector into output directly */
348 /* (e.g., more vectors than time points) */
349
350 for( ii=0 ; ii < nn ; ii++ ) uvec[ii] = (float)asym[ii] ;
351
352 } else { /* n > m: transform eigenvector to get left singular vector */
353 /* (e.g., more time points than vectors) */
354
355 qsum = 0.0 ;
356 for( ii=0 ; ii < nn ; ii++ ){
357 sum = 0.0 ;
358 if( xtyp <= 0 )
359 for( kk=0 ; kk < mm ; kk++ ) sum += xx[ii+kk*nn] * asym[kk] ;
360 else
361 for( kk=0 ; kk < mm ; kk++ ) sum += xar[kk][ii] * asym[kk] ;
362 uvec[ii] = sum ; qsum += sum*sum ;
363 }
364 if( qsum > 0.0 ){ /* L2 normalize */
365 register float fac ;
366 fac = (float)(1.0/sqrt(qsum)) ;
367 for( ii=0 ; ii < nn ; ii++ ) uvec[ii] *= fac ;
368 }
369 }
370
371 /** make it so that uvec dotted into the mean vector is positive **/
372
373 fsum = 0.0f ;
374 for( ii=0 ; ii < nn ; ii++ ) fsum += tvec[ii]*uvec[ii] ;
375 if( fsum < 0.0f ){
376 for( ii=0 ; ii < nn ; ii++ ) uvec[ii] = -uvec[ii] ;
377 }
378
379 /** free at last!!! **/
380
381 #pragma omp critical (MALLOC)
382 { free(deval) ; free(asym) ; }
383
384 return sval ;
385 }
386
387 #undef XPT
388 #undef A
389
390 /*----------------------------------------------------------------------------*/
391 /* Return bootstrap BCa estimates for atanh(correlation) and stdev of that,
392 plus the un-corrected 'raw' correlation, in a triple of floats.
393 *//*--------------------------------------------------------------------------*/
394
395 static int use_pca = 0 ;
396
THD_bstrap_vectcorr(int nlen,int nboot,int xnum,float * xx,int ynum,float * yy)397 float_triple THD_bstrap_vectcorr( int nlen , int nboot ,
398 int xnum , float *xx , int ynum , float *yy )
399 {
400 float_triple retval = {0.0f,0.0f,0.0f} ;
401 float **xar, **yar, *rboot, *rjack, *xbar, *ybar, *uvec, *vvec, *zvec ;
402 float rxy , rjbar , anum,aden,ahat , val , z0hat , zval ;
403 float bmed , bsig , bcor ;
404 int kb,ii,jj , njack ;
405
406 ENTRY("THD_bstrap_vectcorr") ;
407
408 if( nlen < 3 || xnum < 1 || ynum < 1 || xx == NULL || yy == NULL )
409 RETURN(retval) ;
410
411 /* trivial case */
412
413 if( xnum == 1 && ynum == 1 ){
414 retval = THD_bstrap_corr( nlen , xx , yy , 0.0f , nboot ) ;
415 RETURN(retval) ;
416 }
417
418 /* compute mean vectors */
419
420 #pragma omp critical (MALLOC)
421 { xbar = (float *)malloc(sizeof(float)*nlen) ;
422 ybar = (float *)malloc(sizeof(float)*nlen) ; }
423
424 mean_vector( nlen , xnum , 0 , xx , xbar ) ;
425 mean_vector( nlen , ynum , 0 , yy , ybar ) ;
426
427 /* almost trivial case */
428
429 if( xnum < 7 || ynum < 7 ){
430 retval = THD_bstrap_corr( nlen , xbar , ybar , 0.0f , nboot ) ;
431 #pragma omp critical (MALLOC)
432 { free(ybar) ; free(xbar) ; }
433 RETURN(retval) ;
434 }
435
436 /** allocate memory for bootstrapping **/
437
438 if( nboot < 66 ) nboot = 66 ;
439 njack = xnum + ynum ;
440
441 #pragma omp critical (MALLOC)
442 { xar = (float **)malloc(sizeof(float *)*xnum ) ;
443 yar = (float **)malloc(sizeof(float *)*ynum ) ;
444 rboot = (float * )malloc(sizeof(float) *nboot) ;
445 rjack = (float * )malloc(sizeof(float) *njack) ;
446 uvec = (float * )malloc(sizeof(float) *nlen ) ;
447 vvec = (float * )malloc(sizeof(float) *nlen ) ;
448 zvec = (float * )malloc(sizeof(float) *nlen ) ;
449 }
450
451 /* correlation with no resampling at all */
452
453 if( use_pca ){
454 (void)principal_vector( nlen , xnum , 0 , xx , uvec , xbar ) ;
455 (void)principal_vector( nlen , ynum , 0 , yy , vvec , ybar ) ;
456 } else {
457 #pragma omp critical (MEMCPY)
458 { memcpy(uvec,xbar,sizeof(float)*nlen) ;
459 memcpy(vvec,ybar,sizeof(float)*nlen) ; }
460 }
461 rxy = acorrfun( nlen , uvec , vvec ) ;
462 INFO_message("rxy = %g",rxy) ;
463
464 /* bootstrap correlations */
465
466 INFO_message("bootstrap %d times",nboot) ;
467 for( kb=0 ; kb < nboot ; kb++ ){
468 for( ii=0 ; ii < xnum ; ii++ ){ /* resample xx vectors */
469 jj = (lrand48() >> 3) % xnum ; xar[ii] = xx + jj*nlen ;
470 }
471 for( ii=0 ; ii < ynum ; ii++ ){ /* resample yy vectors */
472 jj = (lrand48() >> 3) % ynum ; yar[ii] = yy + jj*nlen ;
473 }
474 if( use_pca ){
475 (void)principal_vector( nlen , xnum , 1 , xar , uvec , xbar ) ;
476 (void)principal_vector( nlen , ynum , 1 , yar , vvec , ybar ) ;
477 } else {
478 mean_vector( nlen , xnum , 1 , xar , uvec ) ;
479 mean_vector( nlen , ynum , 1 , yar , vvec ) ;
480 }
481 rboot[kb] = acorrfun( nlen , uvec , vvec ) ;
482 }
483
484 /* statistics */
485
486 qmedmadbmv_float( nboot , rboot , &bmed , NULL , &bsig ) ;
487
488 bcor = 2*rxy - bmed ;
489
490 val = fabsf(rxy) ; if( val < 0.1f ) val = 0.1f ;
491 val = fabsf(bcor-rxy) / val ; if( val > 1.0f ) val = 1.0f ;
492 bsig = bsig * sqrtf( 1.0f + 4.0f*val*val ) ;
493
494 /* cleanup the trash */
495
496 #pragma omp critical (MALLOC)
497 { free(zvec); free(vvec); free(uvec); free(rjack); free(rboot);
498 free(yar); free(xar);
499 }
500
501 retval.a = bcor ; retval.b = bsig ; retval.c = rxy ;
502
503 RETURN(retval) ;
504 }
505
506 /*----------------------------------------------------------------------------*/
507
main(int argc,char * argv[])508 int main( int argc , char *argv[] )
509 {
510 int iarg=1 , nboot=666 ;
511 MRI_IMAGE *aim , *bim ;
512 float *aar , *bar ;
513 float_triple rval ;
514
515 if( argc < 3 ) ERROR_exit("need 2 file args") ;
516
517 if( AFNI_yesenv("NO_DEMEAN") ) acorr_demean = 0 ;
518 SET_RAND_SEED ;
519 THD_bstrap_set_bca_alphas( 0.1f , 0.2f ) ;
520
521 aim = mri_read_1D(argv[1]) ; if( aim == NULL ) ERROR_exit("Can't read aim") ;
522 bim = mri_read_1D(argv[2]) ; if( bim == NULL ) ERROR_exit("Can't read bim") ;
523 if( aim->nx != aim->nx ) ERROR_exit("nx not same") ;
524
525 if( argc > 3 ){
526 nboot = (int)strtod(argv[3],NULL) ;
527 if( nboot < 66 || nboot > 999999 ) nboot = 666 ;
528 }
529
530 use_pca = 0 ;
531 fprintf(stderr,"----------- MEAN ---------\n") ;
532 rval = THD_bstrap_vectcorr( aim->nx , nboot,
533 aim->ny , MRI_FLOAT_PTR(aim) ,
534 bim->ny , MRI_FLOAT_PTR(bim) ) ;
535 fprintf(stderr,"final = %.5f %.5f %.5f\n",rval.a,rval.b,rval.c) ;
536
537 #if 1
538 fprintf(stderr,"----------- SVD ---------\n") ;
539 use_pca = 1 ;
540 rval = THD_bstrap_vectcorr( aim->nx , nboot,
541 aim->ny , MRI_FLOAT_PTR(aim) ,
542 bim->ny , MRI_FLOAT_PTR(bim) ) ;
543 fprintf(stderr,"final = %.5f %.5f %.5f\n",rval.a,rval.b,rval.c) ;
544 #endif
545 exit(0) ;
546 }
547