1 #include "mrilib.h"
2 #include "zgaussian.c"
3
4 #define NCOR 6
5
6 static char *cor_name[NCOR] = { "Pearson" , "Spearman" , "Quadrant" , "K-Tau-b" , "TicTacToe" , "Quantile:9" } ;
7
8 /* list of functions to compute correlations [cf. thd_correlate.c] */
9
10 typedef float (*ccfun)(int,float *,float *) ;
11 static ccfun cor_func[NCOR] =
12 { THD_pearson_corr, THD_spearman_corr, THD_quadrant_corr, THD_ktaub_corr, THD_tictactoe_corr , THD_quantile_corr } ;
13
14 /* prototype for function to bootstrap the correlations */
15
16 static float_quad Corrboot( int n , float *x , float *y , float xsig , float ysig ,
17 float (*cfun)(int,float *,float *) ) ;
18
19 /* prototype for function to get normal theory confidence intervals */
20
21 static float_pair PCorrCI( int npt , float cor , float alpha ) ;
22
23 /* default number of bootstrap repetitions */
24
25 #undef NBOOT
26 #define NBOOT 4000
27 #undef NBMIN
28 #define NBMIN 400
29
30 static int nboot = NBOOT ; /* changed by the -nboot option */
31 static float alpha = 0.05f ; /* changed by the -alpha option */
32 static int doblk = 0 ; /* changed by the -block option */
33
34 /*----------------------------------------------------------------------------*/
35
usage_1dCorrelate(int detail)36 void usage_1dCorrelate( int detail ){
37 printf(
38 "Usage: 1dCorrelate [options] 1Dfile 1Dfile ...\n"
39 "------\n"
40 " * Each input 1D column is a collection of data points.\n"
41 " * The correlation coefficient between each column pair is computed, along\n"
42 " with its confidence interval (via a bias-corrected bootstrap procedure).\n"
43 " * The minimum sensible column length is 7.\n"
44 " * At least 2 columns are needed [in 1 or more .1D files].\n"
45 " * If there are N input columns, there will be N*(N-1)/2 output rows.\n"
46 " * Output appears on stdout; redirect ('>' or '>>') as needed.\n"
47 " * Only one correlation method can be used in one run of this program.\n"
48 " * This program is basically the basterd offspring of program 1ddot.\n"
49 " * Also see http://en.wikipedia.org/wiki/Confidence_interval\n"
50 "\n"
51 "-------\n"
52 "Methods [actually, only the first letter is needed to choose a method]\n"
53 "------- [and the case doesn't matter: '-P' and '-p' both = '-Pearson']\n"
54 "\n"
55 " -Pearson = Pearson correlation [the default method]\n"
56 " -Spearman = Spearman (rank) correlation [more robust vs. outliers]\n"
57 " -Quadrant = Quadrant (binarized) correlation [most robust, but weaker]\n"
58 " -Ktaub = Kendall's tau_b 'correlation' [popular somewhere, maybe]\n"
59 "\n"
60 "-------------\n"
61 "Other Options [these options cannot be abbreviated!]\n"
62 "-------------\n"
63 "\n"
64 " -nboot B = Set the number of bootstrap replicates to 'B'.\n"
65 " * The default value of B is %d.\n"
66 " * A larger number will give somewhat more accurate\n"
67 " confidence intervals, at the cost of more CPU time.\n"
68 "\n"
69 " -alpha A = Set the 2-sided confidence interval width to '100-A' percent.\n"
70 " * The default value of A is 5, giving the 2.5..97.5%% interval.\n"
71 " * The smallest allowed A is 1 (0.5%%..99.5%%) and the largest\n"
72 " allowed value of A is 20 (10%%..90%%).\n"
73 " * If you are interested assessing if the 'p-value' of a\n"
74 " correlation is smaller than 5%% (say), then you should use\n"
75 " '-alpha 10' and see if the confidence interval includes 0.\n"
76 "\n"
77 " -block = Attempt to allow for serial correlation in the data by doing\n"
78 " *OR* variable-length block resampling, rather than completely\n"
79 " -blk random resampling as in the usual bootstrap.\n"
80 " * You should NOT do this unless you believe that serial\n"
81 " correlation (along each column) is present and significant.\n"
82 " * Block resampling requires at least 20 data points in each\n"
83 " input column. Fewer than 20 will turn off this option.\n"
84 #if 0
85 "\n"
86 " -vsig ss = With this option, you specify that each individual measurement\n"
87 " variable has a nonzero standard deviation ('sig') -- that data\n"
88 " might come from some calibration studies, for example. In the\n"
89 " bootstrap analysis, the random replicates are perturbed by\n"
90 " independent zero-mean normal deviates with standard deviations\n"
91 " given in the file 'ss' -- this file should have 1 row, and that\n"
92 " row should have 1 non-negative number for every column entered\n"
93 " into the correlation analysis.\n"
94 " * In the case of the Pearson coefficient, the 'N:' normal theory\n"
95 " results do NOT include the effect of this extra randomness.\n"
96 #endif
97 "-----\n"
98 "Notes\n"
99 "-----\n"
100 "* For each pair of columns, the output include the correlation value\n"
101 " as directly calculated, plus the bias-corrected bootstrap value, and\n"
102 " the desired (100-A)%% confidence interval [also via bootstrap].\n"
103 "\n"
104 "* The primary purpose of this program is to provide an easy way to get\n"
105 " the bootstrap confidence intervals, since people almost always seem to use\n"
106 " the asymptotic normal theory to decide if a correlation is 'significant',\n"
107 " and this often seems misleading to me [especially for short columns].\n"
108 "\n"
109 "* Bootstrapping confidence intervals for the inverse correlations matrix\n"
110 " (i.e., partial correlations) would be interesting -- anyone out there\n"
111 " need this ability?\n"
112 "\n"
113 "-------------\n"
114 "Sample output [command was '1dCorrelate -alpha 10 A2.1D B2.1D']\n"
115 "-------------\n"
116 "# Pearson correlation [n=12 #col=2]\n"
117 "# Name Name Value BiasCorr 5.00%% 95.00%% N: 5.00%% N:95.00%%\n"
118 "# -------- -------- -------- -------- -------- -------- -------- --------\n"
119 " A2.1D[0] B2.1D[0] +0.57254 +0.57225 -0.03826 +0.86306 +0.10265 +0.83353\n"
120 "\n"
121 "* Bias correction of the correlation had little effect; this is very common.\n"
122 " ++ To be clear, the bootstrap bias correction is to allow for potential bias\n"
123 " in the statistical estimate of correlation when the sample size is small.\n"
124 " ++ It cannot correct for biases that result from faulty data (or faulty\n"
125 " assumptions about the data).\n"
126 "\n"
127 "* The correlation is NOT significant at this level, since the CI (confidence\n"
128 " interval) includes 0 in its range.\n"
129 "\n"
130 "* For the Pearson method ONLY, the last two columns ('N:', as above) also\n"
131 " show the widely used asymptotic normal theory confidence interval. As in\n"
132 " the example, the bootstrap interval is often (but not always) wider than\n"
133 " the theoretical interval.\n"
134 "\n"
135 "* In the example, the normal theory might indicate that the correlation is\n"
136 " significant (less than a 5%% chance that the CI includes 0), but the\n"
137 " bootstrap CI shows that is not a reasonable statistical conclusion.\n"
138 " ++ The principal reason that I wrote this program was to make it easy\n"
139 " to check if the normal (Gaussian) theory for correlation significance\n"
140 " testing is reasonable in any given case -- for small samples, it often\n"
141 " is NOT reasonable!\n"
142 "\n"
143 "* Using the same data with the '-S' option gives the table below, again\n"
144 " indicating that there is no significant correlation between the columns\n"
145 " (note also the lack of the 'N:' results for Spearman correlation):\n"
146 "\n"
147 "# Spearman correlation [n=12 #col=2]\n"
148 "# Name Name Value BiasCorr 5.00%% 95.00%%\n"
149 "# -------- -------- -------- -------- -------- --------\n"
150 " A2.1D[0] B2.1D[0] +0.46154 +0.42756 -0.23063 +0.86078\n"
151 "\n"
152 "-------------\n"
153 "SAMPLE SCRIPT\n"
154 "-------------\n"
155 "This script generates random data and correlates it until it is\n"
156 "statistically significant at some level (default=2%%). Then it\n"
157 "plots the data that looks correlated. The point is to show what\n"
158 "purely random stuff that appears correlated can look like.\n"
159 "(Like most AFNI scripts, this is written in tcsh, not bash.)\n"
160 "\n"
161 "#!/bin/tcsh\n"
162 "set npt = 20\n"
163 "set alp = 2\n"
164 "foreach fred ( `count -dig 1 1 1000` )\n"
165 " 1dcat jrandom1D:${npt},2 > qqq.1D\n"
166 " set aabb = ( `1dCorrelate -spearman -alpha $alp qqq.1D | grep qqq.1D | colrm 1 42` )\n"
167 " set ab = `ccalc -form rint \"1000 * $aabb[1] * $aabb[2]\"`\n"
168 " echo $fred $ab\n"
169 " if( $ab > 1 )then\n"
170 " 1dplot -one -noline -x qqq.1D'[0]' -xaxis -1:1:20:5 -yaxis -1:1:20:5 \\\n"
171 " -DAFNI_1DPLOT_BOXSIZE=0.012 \\\n"
172 " -plabel \"N=$npt trial#=$fred \\alpha=${alp}%% => r\\in[$aabb[1],$aabb[2]]\" \\\n"
173 " qqq.1D'[1]'\n"
174 " break\n"
175 " endif\n"
176 "end\n"
177 "\\rm qqq.1D\n"
178 "\n"
179 "----------------------------------------------------------------------\n"
180 "*** Written by RWCox (AKA Zhark the Mad Correlator) -- 19 May 2011 ***\n"
181
182 , NBOOT ) ;
183
184 PRINT_COMPILE_DATE ;
185 }
186
187 /*============================================================================*/
188
main(int argc,char * argv[])189 int main( int argc , char *argv[] )
190 {
191 int iarg , ii,jj,kk,mm , nvec , nx=0,ny , ff , vlen=4 ;
192 MRI_IMAGE *tim , *vsim=NULL ;
193 MRI_IMARR *tar ;
194 char **vecnam , *tnam ;
195 float *far , **tvec , *vsig=NULL , xsig,ysig ;
196 float_quad qcor ; float_pair pci ; float corst, cor025, cor500, cor975 ;
197 char fmt[256] ;
198 int cormeth=0 ; /* 0=Pearson, 1=Spearman, 2=Quadrant, 3=Kendall tau_b */
199 float (*corfun)(int,float *,float *) ;
200
201 /*-- start the AFNI machinery --*/
202
203 mainENTRY("1dCorrelate main") ; machdep() ;
204
205 /* check for options */
206
207 iarg = 1 ; nvec = 0 ;
208 while( iarg < argc && argv[iarg][0] == '-' ){
209 /* I get by with a little help from my friends? */
210
211 if( strcmp(argv[iarg],"-help") == 0 || strcmp(argv[iarg],"-h") == 0){
212 usage_1dCorrelate(strlen(argv[iarg])>3 ? 2:1);
213 exit(0) ;
214 }
215
216 /*--- methods ---*/
217
218 if( toupper(argv[iarg][1]) == 'P' ){ cormeth = 0 ; iarg++ ; continue ; }
219 if( toupper(argv[iarg][1]) == 'S' ){ cormeth = 1 ; iarg++ ; continue ; }
220 if( toupper(argv[iarg][1]) == 'Q' ){ cormeth = 2 ; iarg++ ; continue ; }
221 if( toupper(argv[iarg][1]) == 'K' ){ cormeth = 3 ; iarg++ ; continue ; }
222 if( toupper(argv[iarg][1]) == 'T' ){ cormeth = 4 ; iarg++ ; continue ; }
223 if( toupper(argv[iarg][1]) == 'U' ){ cormeth = 5 ; iarg++ ; continue ; }
224
225 /*--- set nboot ---*/
226
227 if( strcasecmp(argv[iarg],"-nboot") == 0 || strcasecmp(argv[iarg],"-num") == 0 ){
228 iarg++ ; if( iarg >= argc ) ERROR_exit("Need argument after '-nboot'") ;
229 nboot = (int)strtod(argv[iarg],NULL) ;
230 if( nboot < NBMIN ){
231 WARNING_message("Replacing -nboot %d with %d",nboot,NBMIN) ;
232 nboot = NBMIN ;
233 }
234 iarg++ ; continue ;
235 }
236
237 /*--- set alpha ---*/
238
239 if( strcasecmp(argv[iarg],"-alpha") == 0 ){
240 iarg++ ; if( iarg >= argc ) ERROR_exit("Need argument after '-alpha'") ;
241 alpha = (float)strtod(argv[iarg],NULL) ;
242 if( alpha < 1.0f ){
243 WARNING_message("Replacing -alpha %.1f with 1",alpha) ;
244 alpha = 0.01f ;
245 } else if( alpha > 20.0f ){
246 WARNING_message("Replacing -alpha %.1f with 20",alpha) ;
247 alpha = 0.20f ;
248 } else {
249 alpha *= 0.01f ; /* convert from percent to fraction */
250 }
251 iarg++ ; continue ;
252 }
253
254 /*--- block resampling ---*/
255
256 if( strcasecmp(argv[iarg],"-blk") == 0 || strcasecmp(argv[iarg],"-block") == 0 ){
257 doblk = 1 ; iarg++ ; continue ;
258 }
259
260 if( strcasecmp(argv[iarg],"-vsig") == 0 ){
261 if( vsim != NULL ) ERROR_exit("Can't use -vsig twice!") ;
262 if( ++iarg >= argc ) ERROR_exit("Need argument after -vsig") ;
263 vsim = mri_read_1D(argv[iarg]) ;
264 if( vsim == NULL ) ERROR_exit("Can't read -vsig file '%s'",argv[iarg]) ;
265 iarg++ ; continue ;
266 }
267
268 /*--- user should be flogged ---*/
269
270 ERROR_message("Monstrously illegal option '%s'",argv[iarg]) ;
271 suggest_best_prog_option(argv[0], argv[iarg]);
272 exit(1);
273 }
274
275 /*--- user should be flogged twice ---*/
276
277 if( argc < 2 ){
278 usage_1dCorrelate(1) ; exit(0) ;
279 }
280
281 if( iarg == argc )
282 ERROR_exit("No 1D files on command line!?\n") ;
283
284 /* the function to compute the correlation */
285
286 corfun = cor_func[cormeth] ;
287
288 /* check and assemble list of input 1D files */
289
290 ff = iarg ;
291 INIT_IMARR(tar) ;
292 for( ; iarg < argc ; iarg++ ){
293 tim = mri_read_1D( argv[iarg] ) ;
294 if( tim == NULL ) ERROR_exit("Can't read 1D file '%s'",argv[iarg]) ;
295 if( nx == 0 ){
296 nx = tim->nx ;
297 if( nx < 3 )
298 ERROR_exit("1D file '%.77s' length=%d is less than 3",argv[iarg],nx) ;
299 else if( nx < 7 )
300 WARNING_message("1D file '%.77s' length=%d is less than 7",argv[iarg],nx) ;
301 } else if( tim->nx != nx ){
302 ERROR_exit("Length of 1D file '%.77s' [%d] doesn't match first file [%d]",
303 argv[iarg] , tim->nx , nx );
304 }
305 nvec += tim->ny ;
306 ADDTO_IMARR(tar,tim) ;
307 }
308
309 /* user is really an idiot -- flogging's too good for him */
310
311 if( nvec < 2 ) ERROR_exit("Must have at least 2 input columns!") ;
312
313 if( nx < 20 && doblk ){
314 doblk = 0 ;
315 WARNING_message("Column length %d < 20 ==> cannot use block resampling",nx) ;
316 }
317
318 if( vsim != NULL ){
319 if( vsim->nvox < nvec )
320 ERROR_exit("-vsig file only has %d entries, but needs at least %d",vsim->nvox,nvec) ;
321 vsig = MRI_FLOAT_PTR(vsim) ;
322 }
323
324 /* create vectors from 1D files */
325
326 tvec = (float **)malloc( sizeof(float *)*nvec ) ;
327 vecnam = (char **)malloc( sizeof(char *)*nvec ) ;
328 for( jj=0 ; jj < nvec ; jj++ ){
329 tvec[jj] = (float *)malloc( sizeof(float)*nx ) ;
330 vecnam[jj] = (char *)malloc(sizeof(char)*THD_MAX_NAME) ;
331 }
332
333 /* copy data into new space, create output labels, check for stoopiditees */
334
335 for( kk=mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
336 tim = IMARR_SUBIM(tar,mm) ;
337 far = MRI_FLOAT_PTR(tim) ;
338 tnam = tim->name ; if( tnam == NULL ) tnam = "File" ;
339 for( jj=0 ; jj < tim->ny ; jj++,kk++ ){
340 for( ii=0 ; ii < nx ; ii++ ) tvec[kk][ii] = far[ii+jj*nx] ;
341 sprintf(vecnam[kk],"%s[%d]",THD_trailname(tnam,0),jj) ; /* vector name */
342 iarg = strlen(vecnam[kk]) ; vlen = MAX(vlen,iarg) ;
343 if( THD_is_constant(nx,tvec[kk]) )
344 ERROR_exit("Column %s is constant!",vecnam[kk]) ;
345 }
346 }
347 DESTROY_IMARR(tar) ;
348
349 /*--- Print a beeyootiful header ---*/
350
351 printf("# %s correlation [n=%d #col=%d]\n",cor_name[cormeth],nx,nvec) ;
352 sprintf(fmt,"# %%-%ds %%-%ds",vlen,vlen) ;
353 printf(fmt,"Name","Name") ;
354 printf(" Value BiasCorr %5.2f%% %5.2f%%",50.0f*alpha,100.0f-50.0f*alpha) ;
355 if( cormeth == 0 ) /* Pearson */
356 printf(" N:%5.2f%% N:%5.2f%%",50.0f*alpha,100.0f-50.0f*alpha) ;
357 printf("\n") ;
358 printf("# ") ;
359 for( ii=0 ; ii < vlen ; ii++ ) printf("-") ;
360 printf(" ") ;
361 for( ii=0 ; ii < vlen ; ii++ ) printf("-") ;
362 printf(" ") ;
363 printf(" --------") ;
364 printf(" --------") ;
365 printf(" --------") ;
366 printf(" --------") ;
367 if( cormeth == 0 ){ printf(" --------") ; printf(" --------") ; }
368 printf("\n") ;
369 if( cormeth != 0 ) /* non-Pearson */
370 sprintf(fmt," %%-%ds %%-%ds %%+8.5f %%+8.5f %%+8.5f %%+8.5f\n",vlen,vlen) ;
371 else /* Pearson */
372 sprintf(fmt," %%-%ds %%-%ds %%+8.5f %%+8.5f %%+8.5f %%+8.5f %%+8.5f %%+8.5f\n",vlen,vlen) ;
373
374 /*--- Do some actual work for a suprising change ---*/
375
376 for( jj=0 ; jj < nvec ; jj++ ){ /* loops over column pairs */
377 for( kk=jj+1 ; kk < nvec ; kk++ ){
378
379 if( vsig != NULL ){ xsig = vsig[jj]; ysig = vsig[kk]; } else { xsig = ysig = 0.0f; }
380
381 qcor = Corrboot( nx, tvec[jj], tvec[kk], xsig, ysig, corfun ) ; /* outsourced */
382
383 corst = qcor.a ; cor025 = qcor.b ; cor500 = qcor.c ; cor975 = qcor.d ;
384
385 if( cormeth == 0 ){ /* Pearson */
386 pci = PCorrCI( nx , corst , alpha ) ;
387 printf(fmt, vecnam[jj], vecnam[kk], corst, cor500, cor025, cor975, pci.a,pci.b ) ;
388 } else { /* all other methods */
389 printf(fmt, vecnam[jj], vecnam[kk], corst, cor500, cor025, cor975 ) ;
390 }
391
392 }
393 }
394
395 /* Finished -- go back to watching Star Trek reruns -- Tribbles ahoy, Cap'n! */
396
397 exit(0) ;
398 }
399
400 /*----------------------------------------------------------------------------*/
401 /* Correlate 2 vectors of length n, using the function 'cfun'.
402 Return values are in a float_quad with four elements:
403 .a = standard correlation estimate
404 .b = lower edge (alpha/2) of bias-corrected bootstrap confidence interval
405 .c = middle (50%) of bootstrap confidence interval
406 .d = upper edge (1-alpha/2) of bootstrap confidence interval
407 *//*--------------------------------------------------------------------------*/
408
Corrboot(int n,float * x,float * y,float xsig,float ysig,float (* cfun)(int,float *,float *))409 static float_quad Corrboot( int n , float *x , float *y ,
410 float xsig , float ysig ,
411 float (*cfun)(int,float *,float *) )
412 {
413 float *xar , *yar , *cbb , corst ;
414 float_quad res = {0.0f,-1.0f,0.0f,1.0f} ;
415 float_triple bci ;
416 int ii,jj,kk , nn=n ;
417
418 ENTRY("Corrboot") ;
419
420 if( nn < 3 || x == NULL || y == NULL || cfun == NULL ) RETURN(res) ;
421
422 if( xsig < 0.0f ) xsig = 0.0f ;
423 if( ysig < 0.0f ) ysig = 0.0f ;
424
425 /* workspaces */
426
427 xar = (float *)malloc(sizeof(float)*nn) ; /* resampled x vector */
428 yar = (float *)malloc(sizeof(float)*nn) ; /* resampled y vector */
429 cbb = (float *)malloc(sizeof(float)*nboot) ; /* saved correlations */
430
431 /* compute the un-resampled result */
432
433 for( ii=0 ; ii < nn ; ii++ ){ xar[ii] = x[ii] ; yar[ii] = y[ii] ; }
434
435 corst = cfun(nn,xar,yar) ; /* standard result */
436
437 /* compute bootstrap results */
438
439 for( kk=0 ; kk < nboot ; kk++ ){
440 if( !doblk ){ /* simple resampling */
441 for( ii=0 ; ii < nn ; ii++ ){
442 jj = lrand48() % nn ; xar[ii] = x[jj] ; yar[ii] = y[jj] ;
443 }
444 } else { /* block resampling */
445 int jold = lrand48() % nn ;
446 xar[0] = x[jold] ; yar[0] = y[jold] ;
447 for( ii=1 ; ii < nn ; ii++ ){
448 if( lrand48()%8 == 0 ){ /* 12.5% chance of random jump */
449 jj = lrand48() % nn ;
450 } else {
451 jj = jold+1 ; if( jj == nn ) jj = 0 ;
452 }
453 xar[ii] = x[jj] ; yar[ii] = y[jj] ; jold = jj ;
454 }
455 }
456 if( xsig > 0.0f || ysig > 0.0f ){
457 for( ii=0; ii < nn; ii++ ){
458 xar[ii] += zgaussian()*xsig; yar[ii] += zgaussian()*ysig;
459 }
460 }
461 cbb[kk] = cfun(nn,xar,yar) ; /* bootstrap result */
462 }
463
464 /* get the actual bias-corrected results [cf. thd_correlate.c] */
465
466 bci = THD_bootstrap_confinv( corst , alpha , nboot , cbb ) ;
467
468 /* empty the dustbin */
469
470 free(cbb) ; free(yar) ; free(xar) ;
471
472 /* this is bad */
473
474 if( bci.a == 0.0f && bci.b == 0.0f && bci.c == 0.0f ) RETURN(res) ;
475
476 /* this is good */
477
478 res.a = corst ; res.b = bci.a ; res.c = bci.b ; res.d = bci.c ;
479
480 RETURN(res) ;
481 }
482
483 /*----------------------------------------------------------------------------*/
484 /* Utility functions for PCorrCI (infra). */
485
486 #undef MYatanh
487 #define MYatanh(x) ( ((x)<-0.999329f) ? -4.0f \
488 :((x)>+0.999329f) ? +4.0f : atanhf(x) )
489
MYtanh(float x)490 static float MYtanh( float x )
491 {
492 register float ex , exi ;
493 if( x > 7.0f ) return 1.0f ; /* check for stupid inputs */
494 else if( x < -7.0f ) return -1.0f ;
495 ex = exp(x) ; exi = 1.0f/ex ;
496 return (ex-exi)/(ex+exi) ;
497 }
498
499 /*----------------------------------------------------------------------------*/
500 /* Get the bivariate normal theory confidence interval for Pearson */
501
PCorrCI(int npt,float cor,float alph)502 static float_pair PCorrCI( int npt , float cor , float alph )
503 {
504 float_pair ci = {-1.0f,1.0f} ;
505 float zc , zb,zt , dz ;
506
507 if( npt < 4 || cor <= -1.0f || cor >= 1.0f ) return ci ; /* bad inputs */
508
509 zc = MYatanh(cor) ; /* Piscatorial transformation */
510 dz = qginv(0.5*alph) / sqrt(npt-3.0) ; /* cf. mri_stats.c for qginv() */
511 ci.a = MYtanh(zc-dz) ; /* lower bound */
512 ci.b = MYtanh(zc+dz) ; /* upper bound */
513 return ci ;
514 }
515