Doxygen Source Code Documentation
Deconvolve.c File Reference
#include "RegAna.c"Go to the source code of this file.
Defines | |
| #define | SPOL ((legendre_polort) ? "P_" : "t^") |
| #define | IBOT(ss) min_lag[ss] |
| #define | ITOP(ss) max_lag[ss] |
Functions | |
| double | legendre (double x, int m) |
| int | init_indep_var_matrix (int p, int qp, int polort, int nt, int N, int *good_list, int *block_list, int num_blocks, int num_stimts, float **stimulus, int *stim_length, int *min_lag, int *max_lag, int *nptr, int *stim_base, matrix *xgood) |
| int | init_regression_analysis (int p, int qp, int num_stimts, int *baseline, int *min_lag, int *max_lag, matrix xdata, matrix *x_full, matrix *xtxinv_full, matrix *xtxinvxt_full, matrix *x_base, matrix *xtxinvxt_base, matrix *x_rdcd, matrix *xtxinvxt_rdcd) |
| int | init_glt_analysis (matrix xtxinv, int glt_num, matrix *glt_cmat, matrix *glt_amat, matrix *cxtxinvct) |
| void | regression_analysis (int N, int p, int q, int num_stimts, int *min_lag, int *max_lag, matrix x_full, matrix xtxinv_full, matrix xtxinvxt_full, matrix x_base, matrix xtxinvxt_base, matrix *x_rdcd, matrix *xtxinvxt_rdcd, vector y, float rms_min, float *mse, vector *coef_full, vector *scoef_full, vector *tcoef_full, float *fpart, float *rpart, float *ffull, float *rfull, int *novar, float *fitts, float *errts) |
| void | glt_analysis (int N, int p, matrix x, vector y, float ssef, vector coef, int novar, matrix *cxtxinvct, int glt_num, int *glt_rows, matrix *glt_cmat, matrix *glt_amat, vector *glt_coef, vector *glt_tcoef, float *fglt, float *rglt) |
| double | student_t2pp (double tt, double dof) |
| double | fstat_t2pp (double ff, double dofnum, double dofden) |
| void | report_results (int N, int qp, int q, int p, int polort, int *block_list, int num_blocks, int num_stimts, char **stim_label, int *baseline, int *min_lag, int *max_lag, vector coef, vector tcoef, float *fpart, float *rpart, float ffull, float rfull, float mse, int glt_num, char **glt_label, int *glt_rows, vector *glt_coef, vector *glt_tcoef, float *fglt, float *rglt, char **label) |
Variables | |
| int | legendre_polort = 1 |
| int | demean_base = 1 |
| char | lbuf [65536] |
| char | sbuf [512] |
Define Documentation
|
|
if here, m > 20 ==> use recurrence relation * Definition at line 229 of file Deconvolve.c. Referenced by init_indep_var_matrix(), init_regression_analysis(), regression_analysis(), and report_results(). |
|
|
Definition at line 230 of file Deconvolve.c. Referenced by init_indep_var_matrix(), init_regression_analysis(), regression_analysis(), and report_results(). |
|
|
Definition at line 103 of file Deconvolve.c. Referenced by report_results(). |
Function Documentation
|
||||||||||||||||
|
Definition at line 801 of file Deconvolve.c. Referenced by report_results().
00802 {
00803 int which , status ;
00804 double p , q , f , dfn , dfd , bound ;
00805
00806 if (ff >= 1000.0) return 0.0;
00807
00808 which = 1 ;
00809 p = 0.0 ;
00810 q = 0.0 ;
00811 f = ff ;
00812 dfn = dofnum ;
00813 dfd = dofden ;
00814
00815 cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
00816
00817 if( status == 0 ) return q ;
00818 else return 1.0 ;
00819 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 695 of file Deconvolve.c. References calc_freg(), calc_lcoef(), calc_rcoef(), calc_rsqr(), calc_sse(), calc_tcoef(), ENTRY, p, q, vector_create(), vector_destroy(), and vector_initialize(). Referenced by calculate_results(), and do_xrestore_stuff().
00703 : C(1/(X'X))C' for GLT */ 00704 int glt_num, /* number of general linear tests */ 00705 int * glt_rows, /* number of linear constraints in glt */ 00706 matrix * glt_cmat, /* general linear test matrices */ 00707 matrix * glt_amat, /* constant matrices */ 00708 vector * glt_coef, /* linear combinations from GLT matrices */ 00709 vector * glt_tcoef, /* t-statistics for GLT linear combinations */ 00710 float * fglt, /* F-statistics for the general linear tests */ 00711 float * rglt /* R^2 statistics for the general linear tests */ 00712 ) 00713 00714 { 00715 int iglt; /* index for general linear test */ 00716 int q; /* number of parameters in the rdcd model */ 00717 float sser; /* error sum of squares, reduced model */ 00718 vector rcoef; /* regression parameters for reduced model */ 00719 vector scoef; /* std. devs. for regression parameters */ 00720 00721 ENTRY("glt_analysis") ; 00722 00723 /*----- Initialization -----*/ 00724 vector_initialize (&rcoef); 00725 vector_initialize (&scoef); 00726 00727 00728 /*----- Loop over multiple general linear tests -----*/ 00729 for (iglt = 0; iglt < glt_num; iglt++) 00730 { 00731 /*----- Test for insufficient variation in data -----*/ 00732 if (novar) 00733 { 00734 vector_create (glt_rows[iglt], &glt_coef[iglt]); 00735 vector_create (glt_rows[iglt], &glt_tcoef[iglt]); 00736 fglt[iglt] = 0.0; 00737 rglt[iglt] = 0.0; 00738 } 00739 else 00740 { 00741 /*----- Calculate the GLT linear combinations -----*/ 00742 calc_lcoef (glt_cmat[iglt], coef, &glt_coef[iglt]); 00743 00744 /*----- Calculate t-statistics for GLT linear combinations -----*/ 00745 calc_tcoef (N, p, ssef, cxtxinvct[iglt], 00746 glt_coef[iglt], &scoef, &glt_tcoef[iglt]); 00747 00748 /*----- Calculate regression parameters for the reduced model -----*/ 00749 /*----- (that is, the model in the column space of X but ) -----*/ 00750 /*----- (orthogonal to the restricted column space of XC') -----*/ 00751 calc_rcoef (glt_amat[iglt], coef, &rcoef); 00752 00753 /*----- Calculate error sum of squares for the reduced model -----*/ 00754 sser = calc_sse (x, rcoef, y); 00755 00756 /*----- Calculate the F-statistic for this GLT -----*/ 00757 q = p - glt_rows[iglt]; 00758 fglt[iglt] = calc_freg (N, p, q, ssef, sser); 00759 00760 /*----- Calculate the R^2 statistic for this GLT -----*/ 00761 rglt[iglt] = calc_rsqr (ssef, sser); 00762 00763 } 00764 } 00765 00766 00767 /*----- Dispose of vectors -----*/ 00768 vector_destroy (&rcoef); 00769 vector_destroy (&scoef); 00770 00771 EXRETURN ; 00772 } |
|
||||||||||||||||||||||||
|
Definition at line 523 of file Deconvolve.c. References calc_glt_matrix(), ENTRY, and RETURN. Referenced by calculate_results(), and do_xrestore_stuff().
00524 : 1/(X'X) for full model */ 00525 int glt_num, /* number of general linear tests */ 00526 matrix * glt_cmat, /* general linear test matrices */ 00527 matrix * glt_amat, /* constant GLT matrices for later use */ 00528 matrix * cxtxinvct /* matrices: C(1/(X'X))C' for GLT */ 00529 00530 ) 00531 00532 { 00533 int iglt; /* index for general linear test */ 00534 int ok; /* flag for successful matrix inversion */ 00535 00536 00537 ENTRY("init_glt_analysis") ; 00538 00539 for (iglt = 0; iglt < glt_num; iglt++) 00540 { 00541 ok = calc_glt_matrix (xtxinv, glt_cmat[iglt], &(glt_amat[iglt]), 00542 &(cxtxinvct[iglt])); 00543 if (! ok) RETURN (0); 00544 } 00545 00546 RETURN (1); 00547 } |
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 239 of file Deconvolve.c. References DC_error(), demean_base, matrix::elts, ENTRY, IBOT, ITOP, legendre(), legendre_polort, matrix_create(), matrix_destroy(), matrix_extract_rows(), matrix_initialize(), MRI_FLOAT_PTR, p, RETURN, and STATUS.
00258 {
00259 int m; /* X matrix column index */
00260 int n; /* X matrix row index */
00261 int is; /* input stimulus index */
00262 int ilag; /* time lag index */
00263 int ib; /* block (run) index */
00264 int nfirst, nlast; /* time boundaries of a block (run) */
00265 int mfirst, mlast; /* column boundaries of baseline parameters
00266 for a block (run) */
00267
00268 float * stim_array; /* stimulus function time series */
00269 matrix x; /* X matrix */
00270
00271 int mold ; /* 12 Aug 2004 */
00272 int ibot,itop ;
00273
00274 ENTRY("init_indep_var_matrix") ;
00275
00276
00277 /*----- Initialize X matrix -----*/
00278
00279 STATUS("create x matrix" ) ;
00280 matrix_initialize (&x);
00281 matrix_create (nt, p, &x);
00282
00283
00284 /*----- Set up columns of X matrix corresponding to
00285 the baseline (null hypothesis) signal model -----*/
00286
00287 STATUS("loop over blocks") ;
00288
00289 for (ib = 0; ib < num_blocks; ib++)
00290 {
00291 nfirst = block_list[ib]; /* start time index for this run */
00292 if (ib+1 < num_blocks)
00293 nlast = block_list[ib+1]; /* last+1 time index for this run */
00294 else
00295 nlast = nt;
00296
00297 if(PRINT_TRACING){
00298 char str[256] ;
00299 sprintf(str," block #%d = %d .. %d",ib,nfirst,nlast-1) ; STATUS(str) ;
00300 }
00301
00302 for (n = nfirst; n < nlast; n++)
00303 {
00304 mfirst = ib * (polort+1); /* first column index */
00305 mlast = (ib+1) * (polort+1); /* last+1 column index */
00306
00307 if( !legendre_polort ){ /* the old way: powers */
00308 for (m = mfirst; m < mlast; m++)
00309 x.elts[n][m] = pow ((double)(n-nfirst), (double)(m-mfirst));
00310
00311 } else { /* 15 Jul 2004: the new way: Legendre - RWCox */
00312
00313 double xx , aa=2.0/(nlast-nfirst-1.0) ; /* map nfirst..nlast-1 */
00314 for( m=mfirst ; m < mlast ; m++ ){ /* to interval [-1,1] */
00315 xx = aa*(n-nfirst) - 1.0 ;
00316 x.elts[n][m] = legendre( xx , m-mfirst ) ;
00317 }
00318 }
00319 }
00320
00321 if( mfirst+1 < mlast && demean_base ){ /* 12 Aug 2004: remove means? */
00322 float sum ;
00323 for( m=mfirst+1 ; m < mlast ; m++ ){
00324 sum = 0.0f ;
00325 for( n=nfirst ; n < nlast ; n++ ) sum += x.elts[n][m] ;
00326 sum /= (nlast-nfirst) ;
00327 for( n=nfirst ; n < nlast ; n++ ) x.elts[n][m] -= sum ;
00328 }
00329 }
00330 }
00331
00332
00333 /*----- Set up columns of X matrix corresponding to
00334 time delayed versions of the input stimulus -----*/
00335
00336 STATUS("loop over stimulus time series") ;
00337
00338 m = qp;
00339 for (is = 0; is < num_stimts; is++){
00340 #ifdef USE_BASIS
00341 if( basis_vect[is] != NULL ){ /* 16 Aug 2004 */
00342 float *bv=MRI_FLOAT_PTR(basis_vect[is]) ;
00343 int nf=basis_vect[is]->ny , jj ;
00344 if( PRINT_TRACING ){
00345 char str[256] ;
00346 sprintf(str," stim #%d: expanding into %d basis vectors",is,nf) ;
00347 STATUS(str) ;
00348 }
00349 for( jj=0 ; jj < nf ; jj++ ){
00350 for( n=0 ; n < nt ; n++ ) x.elts[n][m] = bv[n+jj*nt] ;
00351 m++ ;
00352 }
00353 }
00354 else {
00355 #endif
00356 if (stim_length[is] < nt*nptr[is])
00357 {
00358 DC_error ("Input stimulus time series is too short");
00359 RETURN (0);
00360 }
00361 stim_array = stimulus[is]; mold = m ; /* mold = col index we start at */
00362 ibot = IBOT(is) ; itop = ITOP(is) ;
00363 if( PRINT_TRACING ){
00364 char str[256] ;
00365 sprintf(str," stim #%d: ibot=%d itop=%d",is,ibot,itop) ;
00366 STATUS(str) ;
00367 }
00368 for( ilag=ibot ; ilag <= itop ; ilag++ )
00369 {
00370 for (n = 0; n < nt; n++)
00371 {
00372 if (n*nptr[is] < ilag)
00373 x.elts[n][m] = 0.0;
00374 else
00375 x.elts[n][m] = stim_array[n*nptr[is]-ilag];
00376 }
00377 m++;
00378 }
00379
00380 /* 12 Aug 2004: remove mean of baseline columns? */
00381 /* 07 Feb 2005: Oops -- used [m] instead of [mm] in the for(n) loops! */
00382
00383 if( stim_base != NULL && stim_base[is] && demean_base ){
00384 int mm ; float sum ;
00385 STATUS(" remove baseline mean") ;
00386 for( mm=mold ; mm < m ; mm++ ){
00387 sum = 0.0f ;
00388 for( n=0 ; n < nt ; n++ ) sum += x.elts[n][mm] ;
00389 sum /= nt ;
00390 for( n=0 ; n < nt ; n++ ) x.elts[n][mm] -= sum ;
00391 }
00392 }
00393 #ifdef USE_BASIS
00394 }
00395 #endif
00396 }
00397
00398
00399 /*----- Keep only those rows of the X matrix which correspond to
00400 usable time points -----*/
00401
00402 STATUS("extract xgood matrix") ;
00403
00404 matrix_extract_rows (x, N, good_list, xgood);
00405 matrix_destroy (&x);
00406
00407
00408 RETURN (1);
00409
00410 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 419 of file Deconvolve.c. References calc_matrices(), ENTRY, free, IBOT, ITOP, malloc, matrix_destroy(), matrix_initialize(), MTEST, p, and RETURN.
00428 : 1/(X'X) for full model */ 00429 matrix * xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */ 00430 matrix * x_base, /* extracted X matrix for baseline model */ 00431 matrix * xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */ 00432 matrix * x_rdcd, /* extracted X matrices for reduced models */ 00433 matrix * xtxinvxt_rdcd /* matrix: (1/(X'X))X' for reduced models */ 00434 ) 00435 00436 { 00437 int * plist = NULL; /* list of model parameters */ 00438 int ip, it; /* parameter indices */ 00439 int is, js; /* stimulus indices */ 00440 int im, jm; /* lag index */ 00441 int ok; /* flag for successful matrix calculation */ 00442 matrix xtxinv_temp; /* intermediate results */ 00443 int ibot,itop ; 00444 00445 00446 ENTRY("init_regression_analysis") ; 00447 00448 /*----- Initialize matrix -----*/ 00449 matrix_initialize (&xtxinv_temp); 00450 00451 00452 /*----- Initialize matrices for the baseline model -----*/ 00453 plist = (int *) malloc (sizeof(int) * p); MTEST(plist); 00454 for (ip = 0; ip < qp; ip++) 00455 plist[ip] = ip; 00456 it = ip = qp; 00457 for (is = 0; is < num_stimts; is++) 00458 { 00459 ibot = IBOT(is) ; itop = ITOP(is) ; 00460 for (im = ibot; im <= itop; im++) 00461 { 00462 if (baseline[is]) 00463 { 00464 plist[ip] = it; 00465 ip++; 00466 } 00467 it++; 00468 } 00469 } 00470 ok = calc_matrices (xdata, ip, plist, x_base, &xtxinv_temp, xtxinvxt_base); 00471 if (!ok) { matrix_destroy (&xtxinv_temp); RETURN (0); }; 00472 00473 00474 /*----- Initialize matrices for stimulus functions -----*/ 00475 for (is = 0; is < num_stimts; is++) 00476 { 00477 for (ip = 0; ip < qp; ip++) 00478 { 00479 plist[ip] = ip; 00480 } 00481 00482 it = ip = qp; 00483 00484 for (js = 0; js < num_stimts; js++) 00485 { 00486 ibot = IBOT(js) ; itop = ITOP(js) ; 00487 for (jm = ibot; jm <= itop; jm++) 00488 { 00489 if (is != js){ plist[ip] = it; ip++; } 00490 it++; 00491 } 00492 } 00493 00494 ibot = IBOT(is) ; itop = ITOP(is) ; 00495 ok = calc_matrices (xdata, p-(itop-ibot+1), 00496 plist, &(x_rdcd[is]), &xtxinv_temp, &(xtxinvxt_rdcd[is])); 00497 if (!ok) { matrix_destroy (&xtxinv_temp); RETURN (0); }; 00498 } 00499 00500 00501 /*----- Initialize matrices for full model -----*/ 00502 for (ip = 0; ip < p; ip++) 00503 plist[ip] = ip; 00504 ok = calc_matrices (xdata, p, plist, x_full, xtxinv_full, xtxinvxt_full); 00505 if (!ok) { matrix_destroy (&xtxinv_temp); RETURN (0); }; 00506 00507 00508 /*----- Destroy matrix -----*/ 00509 matrix_destroy (&xtxinv_temp); 00510 00511 if (plist != NULL) { free(plist); plist = NULL; } 00512 00513 RETURN (1); 00514 } |
|
||||||||||||
|
Definition at line 105 of file Deconvolve.c. References legendre().
00106 {
00107 if( m < 0 ) return 1.0 ; /* bad input */
00108
00109 switch( m ){ /*** P_m(x) for m=0..20 ***/
00110 case 0: return 1.0 ;
00111 case 1: return x ;
00112 case 2: return (3.0*x*x-1.0)/2.0 ;
00113 case 3: return (5.0*x*x-3.0)*x/2.0 ;
00114 case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
00115 case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
00116 case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
00117 case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
00118 case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
00119
00120 /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
00121
00122 case 9:
00123 return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
00124 + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
00125
00126 case 10:
00127 return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
00128 (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
00129 * x * x) * x * x) * x * x) * x * x;
00130
00131 case 11:
00132 return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
00133 (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
00134 * x * x) * x * x) * x * x) * x * x) * x;
00135
00136 case 12:
00137 return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
00138 (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
00139 + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
00140 * x * x;
00141
00142 case 13:
00143 return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
00144 (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
00145 + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00146 * x * x) * x;
00147
00148 case 14:
00149 return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
00150 (0.236808837890625e4 + (-0.710426513671875e4 +
00151 (0.1089320654296875e5 + (-0.825242919921875e4 +
00152 0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00153 * x * x) * x * x;
00154
00155 case 15:
00156 return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
00157 + (0.710426513671875e4 + (-0.1815534423828125e5 +
00158 (0.2475728759765625e5 + (-0.1713966064453125e5 +
00159 0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
00160 * x * x) * x * x) * x * x) * x;
00161
00162 case 16:
00163 return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
00164 + (-0.4972985595703125e4 + (0.2042476226806641e5 +
00165 (-0.4538836059570312e5 + (0.5570389709472656e5 +
00166 (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
00167 * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00168
00169 case 17:
00170 return (0.3338470458984375e1 + (-0.169149169921875e3 +
00171 (0.2486492797851562e4 + (-0.1633980981445312e5 +
00172 (0.5673545074462891e5 + (-0.1114077941894531e6 +
00173 (0.1242625396728516e6 + (-0.7337407104492188e5 +
00174 0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
00175 * x * x) * x * x) * x * x) * x * x) * x;
00176
00177 case 18:
00178 return -0.1854705810546875e0 + (0.3171546936035156e2 +
00179 (-0.8880331420898438e3 + (0.9531555725097656e4 +
00180 (-0.5106190567016602e5 + (0.153185717010498e6 +
00181 (-0.2692355026245117e6 + (0.275152766418457e6 +
00182 (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
00183 * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00184
00185 case 19:
00186 return (-0.3523941040039062e1 + (0.2220082855224609e3 +
00187 (-0.4084952453613281e4 + (0.3404127044677734e5 +
00188 (-0.153185717010498e6 + (0.4038532539367676e6 +
00189 (-0.6420231216430664e6 + (0.6053360861206055e6 +
00190 (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
00191 * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
00192
00193 case 20:
00194 return 0.1761970520019531e0 + (-0.3700138092041016e2 +
00195 (0.127654764175415e4 + (-0.1702063522338867e5 +
00196 (0.1148892877578735e6 + (-0.4442385793304443e6 +
00197 (0.1043287572669983e7 + (-0.1513340215301514e7 +
00198 (0.1324172688388824e7 + (-0.6404495355606079e6 +
00199 0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
00200 * x * x) * x * x) * x * x) * x * x) * x * x;
00201 }
00202
00203 #if 0
00204 /* order out of range: return Chebyshev instead (it's easy) */
00205
00206 if( x >= 1.0 ) x = 0.0 ;
00207 else if ( x <= -1.0 ) x = 3.14159265358979323846 ;
00208 else x = acos(x) ;
00209 return cos(m*x) ;
00210 #else
00211 /** if here, m > 20 ==> use recurrence relation **/
00212
00213 { double pk, pkm1, pkm2 ; int k ;
00214 pkm2 = legendre( x , 19 ) ;
00215 pkm1 = legendre( x , 20 ) ;
00216 for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
00217 pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
00218 return pk ;
00219 }
00220 #endif
00221 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 556 of file Deconvolve.c. References calc_coef(), calc_freg(), calc_rsqr(), calc_sse(), calc_sse_fit(), calc_tcoef(), ENTRY, IBOT, ITOP, p, q, vector_create(), vector_destroy(), and vector_initialize().
00564 : 1/(X'X) for full model */ 00565 matrix xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */ 00566 matrix x_base, /* extracted X matrix for baseline model */ 00567 matrix xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */ 00568 matrix * x_rdcd, /* extracted X matrices for reduced models */ 00569 matrix * xtxinvxt_rdcd, /* matrix: (1/(X'X))X' for reduced models */ 00570 vector y, /* vector of measured data */ 00571 float rms_min, /* minimum variation in data to fit full model */ 00572 float * mse, /* mean square error from full model */ 00573 vector * coef_full, /* regression parameters */ 00574 vector * scoef_full, /* std. devs. for regression parameters */ 00575 vector * tcoef_full, /* t-statistics for regression parameters */ 00576 float * fpart, /* partial F-statistics for the stimuli */ 00577 float * rpart, /* partial R^2 stats. for the stimuli */ 00578 float * ffull, /* full model F-statistics */ 00579 float * rfull, /* full model R^2 stats. */ 00580 int * novar, /* flag for insufficient variation in data */ 00581 float * fitts, /* full model fitted time series */ 00582 float * errts /* full model residual error time series */ 00583 ) 00584 00585 { 00586 int is; /* input stimulus index */ 00587 float sse_base; /* error sum of squares, baseline model */ 00588 float sse_rdcd; /* error sum of squares, reduced model */ 00589 float sse_full; /* error sum of squares, full model */ 00590 vector coef_temp; /* intermediate results */ 00591 int ibot,itop ; 00592 00593 00594 ENTRY("regression_analysis") ; 00595 00596 /*----- Initialization -----*/ 00597 vector_initialize (&coef_temp); 00598 00599 00600 /*----- Calculate regression coefficients for baseline model -----*/ 00601 calc_coef (xtxinvxt_base, y, &coef_temp); 00602 00603 00604 /*----- Calculate the error sum of squares for the baseline model -----*/ 00605 sse_base = calc_sse (x_base, coef_temp, y); 00606 00607 00608 /*----- Stop here if variation about baseline is sufficiently low -----*/ 00609 if (sqrt(sse_base/N) < rms_min) 00610 { 00611 int it; 00612 00613 *novar = 1; 00614 vector_create (p, coef_full); 00615 vector_create (p, scoef_full); 00616 vector_create (p, tcoef_full); 00617 for (is = 0; is < num_stimts; is++) 00618 { 00619 fpart[is] = 0.0; 00620 rpart[is] = 0.0; 00621 } 00622 for (it = 0; it < N; it++) 00623 { 00624 fitts[it] = 0.0; 00625 errts[it] = 0.0; 00626 } 00627 *mse = 0.0; 00628 *rfull = 0.0; 00629 *ffull = 0.0; 00630 vector_destroy (&coef_temp); 00631 EXRETURN; 00632 } 00633 else 00634 *novar = 0; 00635 00636 00637 /*----- Calculate regression coefficients for the full model -----*/ 00638 calc_coef (xtxinvxt_full, y, coef_full); 00639 00640 00641 /*----- Calculate the error sum of squares for the full model -----*/ 00642 sse_full = calc_sse_fit (x_full, *coef_full, y, fitts, errts); 00643 *mse = sse_full / (N-p); 00644 00645 00646 /*----- Calculate t-statistics for the regression coefficients -----*/ 00647 calc_tcoef (N, p, sse_full, xtxinv_full, 00648 *coef_full, scoef_full, tcoef_full); 00649 00650 00651 /*----- Determine significance of the individual stimuli -----*/ 00652 for (is = 0; is < num_stimts; is++) 00653 { 00654 00655 /*----- Calculate regression coefficients for reduced model -----*/ 00656 calc_coef (xtxinvxt_rdcd[is], y, &coef_temp); 00657 00658 00659 /*----- Calculate the error sum of squares for the reduced model -----*/ 00660 sse_rdcd = calc_sse (x_rdcd[is], coef_temp, y); 00661 00662 00663 /*----- Calculate partial F-stat for significance of the stimulus -----*/ 00664 ibot = IBOT(is) ; itop = ITOP(is) ; 00665 fpart[is] = calc_freg (N, p, p-(itop-ibot+1), sse_full, sse_rdcd); 00666 00667 00668 /*----- Calculate partial R^2 for this stimulus -----*/ 00669 rpart[is] = calc_rsqr (sse_full, sse_rdcd); 00670 00671 } 00672 00673 00674 /*----- Calculate coefficient of multiple determination R^2 -----*/ 00675 *rfull = calc_rsqr (sse_full, sse_base); 00676 00677 00678 /*----- Calculate the total regression F-statistic -----*/ 00679 *ffull = calc_freg (N, p, q, sse_full, sse_base); 00680 00681 00682 /*----- Dispose of vector -----*/ 00683 vector_destroy (&coef_temp); 00684 00685 EXRETURN ; 00686 } |
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 828 of file Deconvolve.c. References fstat_t2pp(), IBOT, ITOP, lbuf, p, q, r, sbuf, SPOL, and student_t2pp().
00858 {
00859 const int MAXBUF = 65000; /* maximum buffer string length */
00860 int m; /* coefficient index */
00861 int is; /* stimulus index */
00862 int ilag; /* time lag index */
00863
00864 int iglt; /* general linear test index */
00865 int ilc; /* linear combination index */
00866
00867 double pvalue; /* p-value corresponding to F-value */
00868 int r; /* number of parameters in the reduced model */
00869
00870 int ib; /* block (run) index */
00871 int mfirst, mlast; /* column boundaries of baseline parameters
00872 for a block (run) */
00873 int ibot,itop ;
00874
00875
00876 lbuf[0] = '\0' ; /* make this a 0 length string to start */
00877
00878 /** for each reference, make a string into sbuf **/
00879
00880
00881 /*----- Statistical results for baseline fit -----*/
00882 if (num_blocks == 1)
00883 {
00884 sprintf (sbuf, "\nBaseline: \n");
00885 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf); else goto finisher ;
00886 for (m=0; m < qp; m++)
00887 {
00888 sprintf (sbuf, "%s%d coef = %10.4f ", SPOL,m, coef.elts[m]);
00889 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
00890 sprintf (sbuf, "%s%d t-st = %10.4f ", SPOL,m, tcoef.elts[m]);
00891 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
00892 pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
00893 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
00894 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00895 }
00896 }
00897 else
00898 {
00899 for (ib = 0; ib < num_blocks; ib++)
00900 {
00901 sprintf (sbuf, "\nBaseline for Run #%d: \n", ib+1);
00902 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf); else goto finisher ;
00903
00904 mfirst = ib * (polort+1);
00905 mlast = (ib+1) * (polort+1);
00906 for (m = mfirst; m < mlast; m++)
00907 {
00908 sprintf (sbuf, "%s%d coef = %10.4f ",
00909 SPOL,m - mfirst, coef.elts[m]);
00910 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
00911 sprintf (sbuf, "%s%d t-st = %10.4f ",
00912 SPOL,m - mfirst, tcoef.elts[m]);
00913 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
00914 pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
00915 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
00916 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00917 }
00918 }
00919 }
00920
00921
00922 /*----- Statistical results for stimulus response -----*/
00923 m = qp;
00924 for (is = 0; is < num_stimts; is++)
00925 {
00926 if (baseline[is])
00927 sprintf (sbuf, "\nBaseline: %s \n", stim_label[is]);
00928 else
00929 sprintf (sbuf, "\nStimulus: %s \n", stim_label[is]);
00930 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf); else goto finisher ;
00931 ibot = IBOT(is) ; itop = ITOP(is) ;
00932 for (ilag = ibot; ilag <= itop; ilag++)
00933 {
00934 sprintf (sbuf,"h[%2d] coef = %10.4f ", ilag, coef.elts[m]);
00935 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
00936 sprintf (sbuf,"h[%2d] t-st = %10.4f ", ilag, tcoef.elts[m]);
00937 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00938 pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
00939 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
00940 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00941 m++;
00942 }
00943
00944 sprintf (sbuf, " R^2 = %10.4f ", rpart[is]);
00945 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00946 r = p - (itop-ibot+1);
00947 sprintf (sbuf, "F[%2d,%3d] = %10.4f ", p-r, N-p, fpart[is]);
00948 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00949 pvalue = fstat_t2pp ((double)fpart[is], (double)(p-r), (double)(N-p));
00950 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
00951 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00952 }
00953
00954
00955 /*----- Statistical results for full model -----*/
00956 sprintf (sbuf, "\nFull Model: \n");
00957 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00958
00959 sprintf (sbuf, " MSE = %10.4f \n", mse);
00960 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00961
00962 sprintf (sbuf, " R^2 = %10.4f ", rfull);
00963 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00964
00965 sprintf (sbuf, "F[%2d,%3d] = %10.4f ", p-q, N-p, ffull);
00966 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00967 pvalue = fstat_t2pp ((double)ffull, (double)(p-q), (double)(N-p));
00968 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
00969 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00970
00971
00972 /*----- Statistical results for general linear test -----*/
00973 if (glt_num > 0)
00974 {
00975 for (iglt = 0; iglt < glt_num; iglt++)
00976 {
00977 sprintf (sbuf, "\nGeneral Linear Test: %s \n", glt_label[iglt]);
00978 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
00979 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
00980 {
00981 sprintf (sbuf, "LC[%d] coef = %10.4f ",
00982 ilc, glt_coef[iglt].elts[ilc]);
00983 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
00984 sprintf (sbuf, "LC[%d] t-st = %10.4f ",
00985 ilc, glt_tcoef[iglt].elts[ilc]);
00986 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
00987 pvalue = student_t2pp ((double)glt_tcoef[iglt].elts[ilc],
00988 (double)(N-p));
00989 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
00990 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00991 }
00992
00993 sprintf (sbuf, " R^2 = %10.4f ", rglt[iglt]);
00994 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
00995
00996 r = p - glt_rows[iglt];
00997 sprintf (sbuf, "F[%2d,%3d] = %10.4f ", p-r, N-p, fglt[iglt]);
00998 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
00999 pvalue = fstat_t2pp ((double)fglt[iglt],
01000 (double)(p-r), (double)(N-p));
01001 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
01002 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
01003 }
01004 }
01005
01006 finisher:
01007 if (strlen(lbuf) >= MAXBUF)
01008 strcat (lbuf, "\n\nWarning: Screen output buffer is full. \n");
01009
01010 *label = lbuf ; /* send address of lbuf back in what label points to */
01011
01012 }
|
|
||||||||||||
|
Definition at line 784 of file Deconvolve.c. References incbeta(), lnbeta(), and tt. Referenced by report_results().
|
Variable Documentation
|
|
Definition at line 101 of file Deconvolve.c. Referenced by init_indep_var_matrix(). |
|
|
Definition at line 823 of file Deconvolve.c. Referenced by report_results(). |
|
|
Definition at line 100 of file Deconvolve.c. Referenced by init_indep_var_matrix(). |
|
|
Definition at line 824 of file Deconvolve.c. Referenced by report_results(). |
