Variational Bayes for GLM-AR modelling in a slice of fMRI
FORMAT [slice] = spm_vb_glmar (Y,slice)
Y [T x N] time series with T time points, N voxels
slice data structure containing the following fields:
.X [T x k] the design matrix
.p order of AR model
.D [N x N] spatial precision matrix
(see spm_vb_set_priors.m)
The above fields are mandatory. The fields below are
optional or are filled in by this function.
OPTIMISIATION PARAMETERS:
.tol termination tolerance (default = 0.01% increase in F)
.maxits maximum number of iterations (default=4)
.verbose '1' for description of actions (default=1)
.update_??? set to 1 to update parameter ??? (set to 0 to fix)
eg. update_alpha=1; % update prior precision on W
ESTIMATED REGRESSION COEFFICIENTS:
.wk_mean [k x N] VB regression coefficients
.wk_ols [k x N] OLS " "
.w_cov N-element cell array with entries [k x k]
.w_dev [k x N] standard deviation of regression coeffs
ESTIMATED AR COEFFICIENTS:
.ap_mean [p x N] VB AR coefficients
.ap_ols [p x N] OLS AR coefficients
.a_cov N-element cell array with entries [p x p]
ESTIMATED NOISE PRECISION:
.b_lambda [N x 1] temporal noise precisions
.c_lambda
.mean_lambda
MODEL COMPARISON AND COEFFICIENT RESELS:
.gamma_tot [k x 1] Coefficient RESELS
.F Negative free energy (used for model selection)
.F_record [its x 1] record of F at each iteration
.elapsed_seconds estimation time
PRIORS:
.b_alpha [k x 1] spatial prior precisions for W
.c_alpha
.mean_alpha
.b_beta [p x 1] spatial prior precisions for AR
.c_beta
.mean_beta
HYPERPRIORS:
.b_alpha_prior priors on alpha
.c_alpha_prior
.b_beta_prior priors on beta
.c_beta_prior
.b_lambda_prior priors on temporal noise precisions
.c_lambda_prior
There are other fields that are used internally
%W% Will Penny and Nelson Trujillo-Barreto %E%0001 function [slice] = spm_vb_glmar (Y,slice) 0002 % Variational Bayes for GLM-AR modelling in a slice of fMRI 0003 % FORMAT [slice] = spm_vb_glmar (Y,slice) 0004 % 0005 % Y [T x N] time series with T time points, N voxels 0006 % 0007 % slice data structure containing the following fields: 0008 % 0009 % .X [T x k] the design matrix 0010 % .p order of AR model 0011 % .D [N x N] spatial precision matrix 0012 % (see spm_vb_set_priors.m) 0013 % 0014 % The above fields are mandatory. The fields below are 0015 % optional or are filled in by this function. 0016 % 0017 % OPTIMISIATION PARAMETERS: 0018 % 0019 % .tol termination tolerance (default = 0.01% increase in F) 0020 % .maxits maximum number of iterations (default=4) 0021 % .verbose '1' for description of actions (default=1) 0022 % .update_??? set to 1 to update parameter ??? (set to 0 to fix) 0023 % eg. update_alpha=1; % update prior precision on W 0024 % 0025 % ESTIMATED REGRESSION COEFFICIENTS: 0026 % 0027 % .wk_mean [k x N] VB regression coefficients 0028 % .wk_ols [k x N] OLS " " 0029 % .w_cov N-element cell array with entries [k x k] 0030 % .w_dev [k x N] standard deviation of regression coeffs 0031 % 0032 % ESTIMATED AR COEFFICIENTS: 0033 % 0034 % .ap_mean [p x N] VB AR coefficients 0035 % .ap_ols [p x N] OLS AR coefficients 0036 % .a_cov N-element cell array with entries [p x p] 0037 % 0038 % ESTIMATED NOISE PRECISION: 0039 % 0040 % .b_lambda [N x 1] temporal noise precisions 0041 % .c_lambda 0042 % .mean_lambda 0043 % 0044 % MODEL COMPARISON AND COEFFICIENT RESELS: 0045 % 0046 % .gamma_tot [k x 1] Coefficient RESELS 0047 % .F Negative free energy (used for model selection) 0048 % .F_record [its x 1] record of F at each iteration 0049 % .elapsed_seconds estimation time 0050 % PRIORS: 0051 % 0052 % .b_alpha [k x 1] spatial prior precisions for W 0053 % .c_alpha 0054 % .mean_alpha 0055 % 0056 % .b_beta [p x 1] spatial prior precisions for AR 0057 % .c_beta 0058 % .mean_beta 0059 % 0060 % HYPERPRIORS: 0061 % 0062 % .b_alpha_prior priors on alpha 0063 % .c_alpha_prior 0064 % 0065 % .b_beta_prior priors on beta 0066 % .c_beta_prior 0067 % 0068 % .b_lambda_prior priors on temporal noise precisions 0069 % .c_lambda_prior 0070 % 0071 % There are other fields that are used internally 0072 % 0073 % %W% Will Penny and Nelson Trujillo-Barreto %E% 0074 0075 tic; 0076 0077 % Check arguments and set defaults 0078 if nargin < 2, 0079 disp('Error in in spm_vbglmar_slice: function needs two arguments'); 0080 return 0081 end 0082 [T,N]=size(Y); 0083 slice.T=T; 0084 slice.N=N; 0085 if ~isfield(slice,'X') 0086 disp('Error in spm_vbglmar_slice: mandatory field X missing'); 0087 return 0088 else 0089 X=slice.X; 0090 end 0091 [tmp,k]=size(X); 0092 if ~(tmp==T) 0093 disp('Error in spm_vbglmar_slice: X is not of compatible size to Y'); 0094 return 0095 end 0096 slice.k=k; 0097 p=slice.p; 0098 0099 if ~isfield(slice,'D') 0100 disp('Error in spm_vbglmar_slice: mandatory field D missing'); 0101 return 0102 end 0103 0104 F = 0; 0105 last_F=0; 0106 slice=spm_vb_init_slice(Y,slice); 0107 0108 if slice.verbose 0109 disp(' '); 0110 disp('Starting VB-GLM-AR-SLICE'); 0111 end 0112 0113 for it = 1:slice.maxits, % Loop over iterations 0114 0115 if slice.update_w 0116 slice=spm_vb_w (Y,slice); 0117 end 0118 if (slice.p>0) & (slice.update_a) 0119 slice=spm_vb_a (Y,slice); 0120 end 0121 if slice.update_lambda 0122 slice=spm_vb_lambda (Y,slice); 0123 end 0124 if slice.update_alpha 0125 slice=spm_vb_alpha (Y,slice); 0126 end 0127 if (slice.p>0) & (slice.update_beta) 0128 slice=spm_vb_beta (Y,slice); 0129 end 0130 if slice.update_F 0131 F=spm_vb_F (Y,slice); 0132 end 0133 if slice.verbose 0134 disp(sprintf('Iteration %d, F=%1.2f',it,F)); 0135 end 0136 0137 if slice.update_F 0138 slice.F_record(it)=F; 0139 delta_F=F-last_F; 0140 if it > 2 0141 if delta_F < 0 0142 disp(sprintf('********** Warning: decrease in F of %1.4f per cent *************',100*(delta_F/F))); 0143 break; 0144 elseif abs(delta_F/F) < slice.tol, 0145 break; 0146 end; 0147 end 0148 last_F=F; 0149 end 0150 end; 0151 0152 if slice.update_F 0153 slice.F=F; 0154 end 0155 0156 slice=spm_vb_gamma(Y,slice); 0157 0158 slice.elapsed_seconds=toc; 0159