Bayesian Wiener estimation for MEG.
Estimate current variance by using time sequence of all trials.
Mixture of Multiple fMRI activity patterns are used for prior
Current variance is represented by mixture of fMRI patterns
Weight coefficient for each fMRI pattern is estimated
-- Syntax
function [Model,Info] = vbmeg_multi_fmri(B, Gall, Gact, COV, vb_parm),
-- Input
B{Nsession} : MEG magnetic field data matrix
size(B{i}) = [Nsensors(i), Tsample, Ntrials(i)]
Gall{Nsession}: Lead field for whole brain region as background activity
size(Gall{i}) = [Nsensors(i) Ndipole]
If Gall = [], the sparse-mode (focal dipoles only) is
executed.
Gact{Nsession} : Lead field (Current basis function)
for region with high estimation accuracy
size(Gact{i}) = [Nsensors(i) Ndipole]
If Gact = [], the iso-mode (no focal dipoles) is
executed.
COV{Nsession} : Noise covariance matrix
size(COV{i}) = [Nsensors(i) Nsensors(i)]
vb_parm : Constants parameter
(Field of vb_parm)
.Ntrain
.Nskip
.a_min
.a_max
.Npre_train (optional)
.Ntrials
.Nsession
.Nwindow
.Twindow
.Ta0
.Tv0
.sx0
.a0
.v0
.wiener
.act
.cont_pr
-- Output
Model.a : Current variance (1/alpha)
Model.v : Current variance for background activity
Model.sx : Observation noise variance
Model.act : Multiple fMRI activity patterns for focal area
Model.act_v : Their estimated weights
Info : Containing convergence information such as 'FE','Ev','Err'
-- Note
Brain noise : Background brain activity is assumed
as independent Gaussian noise with same variance
"Nvact" (# of vertex for active dipoles) and
"Nt" (length of time window) are common for all sessions.
"Nsensors" (# of sensor channels in a session) and
"Ntrial" (# of trials in a session) could depend on each session.
Note that ax0_j, ax0_z, and sx are defined as variance:
the inverse of inverse-variance (i.e having the diminsion of variance).
These notations are different from those in Sato-san's note.
--- History
2004-12-5 M. Sato
2005-04-12 O.Yamashita
* The whole program are revised
so that variables have consistent names with Sato-san's note.
* The noise variance, sx, ax_z and ax_j are independently handled.
* A field 'fix_parm' is added to 'vb_parm' struct.
* The input 'COV' must be given.
2005-04-16 M. Sato
* Modified free-energy calculation
2005-05-07 M. Sato
* Modified variance update (Faster convergence)
2005-05-13 O.Yamashita
* 'Modelinit' is added as the 6th input optional argument.
* Iso-mode is implemented. (when Gact is [])
* A field 'fix_parm' is changed to two flag paramters, 'update_v' and
'update_sx'.
* lower bound of 'gain_z'
2005-05-16 O.Yamashita
* Sparse-mode is implementd. (when Gall is [])
* lower bound of 'gain_j'
2005-05-31 O.Yamashita
* constant term of free energy
2005-06-02 O.Yamashita
* constat term of free energy revised
2005/08/22 O.Yamashita --- ver.30b
* vb_parm is introduced
* Output argument is modified
* Minor bug fix
2005/09/29 O. Yamashita
* Minor bug fix (when Norinet=2 & iso-mode)
2005/12/5 M.Sato
* New version
* Prior is given for each time window
2006/8/11 M.Sato
Bayesian Wiener estimation
Mixture of Multiple fMRI activity patterns are used for prior
2008-11-28 Taku Yoshioka
Use vb_disp() for displaying message
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [Model,Info] = ... 0002 vbmeg_multi_fmri(B, Gall, Gact, COV, vb_parm) 0003 % Bayesian Wiener estimation for MEG. 0004 % Estimate current variance by using time sequence of all trials. 0005 % Mixture of Multiple fMRI activity patterns are used for prior 0006 % Current variance is represented by mixture of fMRI patterns 0007 % Weight coefficient for each fMRI pattern is estimated 0008 % 0009 % -- Syntax 0010 % function [Model,Info] = vbmeg_multi_fmri(B, Gall, Gact, COV, vb_parm), 0011 % 0012 % -- Input 0013 % B{Nsession} : MEG magnetic field data matrix 0014 % size(B{i}) = [Nsensors(i), Tsample, Ntrials(i)] 0015 % Gall{Nsession}: Lead field for whole brain region as background activity 0016 % size(Gall{i}) = [Nsensors(i) Ndipole] 0017 % If Gall = [], the sparse-mode (focal dipoles only) is 0018 % executed. 0019 % Gact{Nsession} : Lead field (Current basis function) 0020 % for region with high estimation accuracy 0021 % size(Gact{i}) = [Nsensors(i) Ndipole] 0022 % If Gact = [], the iso-mode (no focal dipoles) is 0023 % executed. 0024 % COV{Nsession} : Noise covariance matrix 0025 % size(COV{i}) = [Nsensors(i) Nsensors(i)] 0026 % vb_parm : Constants parameter 0027 % 0028 % (Field of vb_parm) 0029 % .Ntrain 0030 % .Nskip 0031 % .a_min 0032 % .a_max 0033 % .Npre_train (optional) 0034 % .Ntrials 0035 % .Nsession 0036 % .Nwindow 0037 % .Twindow 0038 % .Ta0 0039 % .Tv0 0040 % .sx0 0041 % .a0 0042 % .v0 0043 % .wiener 0044 % .act 0045 % .cont_pr 0046 % 0047 % -- Output 0048 % Model.a : Current variance (1/alpha) 0049 % Model.v : Current variance for background activity 0050 % Model.sx : Observation noise variance 0051 % 0052 % Model.act : Multiple fMRI activity patterns for focal area 0053 % Model.act_v : Their estimated weights 0054 % 0055 % Info : Containing convergence information such as 'FE','Ev','Err' 0056 % 0057 % -- Note 0058 % Brain noise : Background brain activity is assumed 0059 % as independent Gaussian noise with same variance 0060 % 0061 % "Nvact" (# of vertex for active dipoles) and 0062 % "Nt" (length of time window) are common for all sessions. 0063 % 0064 % "Nsensors" (# of sensor channels in a session) and 0065 % "Ntrial" (# of trials in a session) could depend on each session. 0066 % 0067 % Note that ax0_j, ax0_z, and sx are defined as variance: 0068 % the inverse of inverse-variance (i.e having the diminsion of variance). 0069 % These notations are different from those in Sato-san's note. 0070 % 0071 % --- History 0072 % 2004-12-5 M. Sato 0073 % 2005-04-12 O.Yamashita 0074 % * The whole program are revised 0075 % so that variables have consistent names with Sato-san's note. 0076 % * The noise variance, sx, ax_z and ax_j are independently handled. 0077 % * A field 'fix_parm' is added to 'vb_parm' struct. 0078 % * The input 'COV' must be given. 0079 % 2005-04-16 M. Sato 0080 % * Modified free-energy calculation 0081 % 2005-05-07 M. Sato 0082 % * Modified variance update (Faster convergence) 0083 % 2005-05-13 O.Yamashita 0084 % * 'Modelinit' is added as the 6th input optional argument. 0085 % * Iso-mode is implemented. (when Gact is []) 0086 % * A field 'fix_parm' is changed to two flag paramters, 'update_v' and 0087 % 'update_sx'. 0088 % * lower bound of 'gain_z' 0089 % 2005-05-16 O.Yamashita 0090 % * Sparse-mode is implementd. (when Gall is []) 0091 % * lower bound of 'gain_j' 0092 % 2005-05-31 O.Yamashita 0093 % * constant term of free energy 0094 % 2005-06-02 O.Yamashita 0095 % * constat term of free energy revised 0096 % 2005/08/22 O.Yamashita --- ver.30b 0097 % * vb_parm is introduced 0098 % * Output argument is modified 0099 % * Minor bug fix 0100 % 2005/09/29 O. Yamashita 0101 % * Minor bug fix (when Norinet=2 & iso-mode) 0102 % 0103 % 2005/12/5 M.Sato 0104 % * New version 0105 % * Prior is given for each time window 0106 % 2006/8/11 M.Sato 0107 % Bayesian Wiener estimation 0108 % Mixture of Multiple fMRI activity patterns are used for prior 0109 % 2008-11-28 Taku Yoshioka 0110 % Use vb_disp() for displaying message 0111 % 0112 % Copyright (C) 2011, ATR All Rights Reserved. 0113 % License : New BSD License(see VBMEG_LICENSE.txt) 0114 0115 const = vb_define_verbose; 0116 VERBOSE_LEVEL_NOTICE = const.VERBOSE_LEVEL_NOTICE; 0117 VERBOSE_LEVEL_INFO = const.VERBOSE_LEVEL_INFO; 0118 0119 Ntrain = vb_parm.Ntrain; 0120 Nskip = vb_parm.Nskip; 0121 a_min = vb_parm.a_min; 0122 a_max = vb_parm.a_max; 0123 0124 if isfield(vb_parm, 'Npre_train') 0125 % Number of iteration using original VB-update rule 0126 % for stable estimation 0127 Npre_train = vb_parm.Npre_train; 0128 else 0129 Npre_train = 0; 0130 end 0131 0132 vb_disp(['--- Initial VB-update iteration = ' num2str(Npre_train) ], ... 0133 VERBOSE_LEVEL_NOTICE); 0134 vb_disp(['--- Total update iteration = ' num2str(Ntrain) ], ... 0135 VERBOSE_LEVEL_NOTICE); 0136 0137 % Number parameters 0138 Ntrials = vb_parm.Ntrials; % Number of trials 0139 Nsensors = vb_parm.Nsensors; % Number of sensors 0140 Nsession = vb_parm.Nsession; % Number of sessions 0141 Nwindow = vb_parm.Nwindow; % Number of time windows 0142 Twindow = vb_parm.Twindow; % Time window index (Tstart, Tend) 0143 0144 Ntotaltrials = sum(Ntrials); % Total number of trials in all sessions 0145 Nerror = sum(Ntrials.*Nsensors); % Total MEG channel used for error estimate 0146 0147 % Set prior parameters 0148 Ta0 = vb_parm.Ta0 ; % confidence parameter of focal 0149 Tv0 = vb_parm.Tv0 ; % confidence parameter of global 0150 sx0 = vb_parm.sx0; % observation noise variance 0151 a0 = max(vb_parm.a0, a_min ); % initial variance of focal 0152 v0 = max(vb_parm.v0, a_min ); % initial weight of global & multiple fMRI 0153 0154 % Sensor noise update flag 0155 if isfield(vb_parm, 'wiener') 0156 update_sx = vb_parm.wiener.update_sx; 0157 else 0158 update_sx = 1; 0159 end 0160 0161 % Set fMRI activity pattern 0162 if isfield(vb_parm, 'act') 0163 act = vb_parm.act; 0164 Nact_ptn = size(act,2); 0165 else 0166 act = []; 0167 Nact_ptn = 0; 0168 end 0169 0170 % Nvarall = 0 (No Global variance estimation) 0171 % = 1 ( Global variance estimation) 0172 0173 if isempty(Gall) 0174 % No global area case 0175 Njall = []; 0176 Tv0 = []; 0177 v0 = []; 0178 Nvarall = 0; 0179 update_v = 0; 0180 else 0181 % Global area variance estimation case 0182 Njall = size(Gall{1},2); % Number of global area dipoles 0183 Nvarall = 1; 0184 update_v = 1; 0185 end 0186 0187 % Nvarea = # of variance parameters 0188 Narea = (Nvarall + Nact_ptn); 0189 0190 if isempty(Gact) 0191 % No focal area case 0192 Njact = 1; % avoid empt 0193 act_id = []; 0194 act = 0; 0195 update_z = 0; 0196 elseif Nact_ptn > 0 0197 % Focal area variance estimation case 0198 Njact = size(Gact{1},2); % Number of active dipoles 0199 Nvaract = Nact_ptn; 0200 act_id = (Nvarall+1):Narea; % variance index for focal areas 0201 Njall = [Njall , Njact * ones(1,Nvaract)]; 0202 Tv0 = [Tv0 , Ta0 * ones(1,Nvaract)]; 0203 v0 = [ v0 , a0 * ones(1,Nvaract)]; 0204 0205 update_z = 1; 0206 if isfield(vb_parm, 'wiener') 0207 update_v = vb_parm.wiener.update_v; 0208 else 0209 update_v = 0; 0210 end 0211 else 0212 error('No active patterns are given') 0213 end 0214 0215 mesg = {'is fixed', 'is estimated'}; 0216 vb_disp(['Sensor noise variance ' mesg{update_sx + 1}], ... 0217 VERBOSE_LEVEL_NOTICE); 0218 vb_disp(['Background current variance ' mesg{update_v + 1}], ... 0219 VERBOSE_LEVEL_NOTICE); 0220 0221 % Calculate (G * Az(j) * G') for multiple patterns 0222 GGall = cell(Nsession,1); 0223 0224 for ns = 1 : Nsession 0225 Nch = Nsensors(ns); % number of channel available in each session 0226 0227 if isempty(Gall) 0228 GGall{ns} = []; 0229 else 0230 GG = Gall{ns}*Gall{ns}'; 0231 GGall{ns} = GG(:); 0232 end 0233 0234 if ~isempty(Gact) 0235 if Njact ~= size(act,1), 0236 error('Activity pattern does not match with focal area\n') 0237 end 0238 0239 G = Gact{ns}; 0240 Gt = G'; 0241 for j = 1:Nact_ptn, 0242 Az = repmat(act(:,j) , [1 Nch]); 0243 GG = G*(Az.*Gt); 0244 0245 GGall{ns} = [GGall{ns} GG(:)]; 0246 end 0247 end 0248 end 0249 0250 %%% calculate BB' 0251 for ns = 1 : Nsession 0252 for nw = 1 : Nwindow 0253 B0 = 0; 0254 for i = 1 : Ntrials(ns) 0255 Bt = B{ns}(:,Twindow(nw,1):Twindow(nw,2),i); 0256 B0 = B0 + Bt * Bt'; 0257 end 0258 BBall{ns,nw} = B0; 0259 end 0260 end 0261 0262 clear B Gall 0263 0264 [Nv1, Nv2 ]= size(v0); 0265 0266 if Nv2 ~= Narea, error('Dimension mismatch for V0'); end; 0267 0268 % Temporal variable for variance (ratio) 0269 aall = zeros(Njact, Nwindow); 0270 vall = zeros(Narea,Nwindow); 0271 sxall = zeros(1,Nwindow); 0272 0273 % Free energy and Error history 0274 Info = zeros(Nwindow*Ntrain, 8); 0275 k_check = 0; 0276 0277 %%%%% Time window loop %%%%%%%% 0278 for nw=1:Nwindow 0279 0280 Nt = Twindow(nw,2) - Twindow(nw,1) + 1; 0281 Ntt = Nt*Ntotaltrials; 0282 0283 % confidence parameters 0284 gx0_j = Tv0; 0285 0286 gx_j = 0.5 * Ntt * Njall + gx0_j; 0287 0288 % expected variance parameters 0289 ax0_j = v0; 0290 sx = sx0; 0291 0292 % variance hyper-parameter initialization 0293 if vb_parm.cont_pr == OFF | nw == 1 0294 % use the prior variance as the initial expected value 0295 ax_j = ax0_j; 0296 sx = sx0; 0297 else 0298 % use the estimated variance in the previous time window as the 0299 % initial expected value 0300 end 0301 0302 %%%%%% Estimation Loop %%%%%% 0303 for k=1:Ntrain 0304 0305 %%% VB E-step -- current estimation 0306 0307 % Initialize averaging variable 0308 tr1_b = 0; % Reconstruction error 0309 tr2_b = 0; % Error variance 0310 mean_jj = zeros(1,Narea); % Z*Z 0311 gain_j = zeros(1,Narea); % Estimation Gain for Z 0312 Hj = 0; % log(det(SB)) 0313 0314 for ns=1:Nsession 0315 0316 Nch = Nsensors(ns); 0317 Nch2 = Nch^2; 0318 Ntry = Ntrials(ns); 0319 0320 % Lead field for each session 0321 GG = GGall{ns}; % ( G * Az(j) * G' ) 0322 0323 % MEG covariance for each session 0324 BB = BBall{ns,nw}; % (B B') 0325 0326 % Noise covariance for each session 0327 Cov = COV{ns}; % Sigma_G^{-1} 0328 0329 % Inverse filter for each session 0330 GSG = reshape(GG*ax_j(:), [Nch Nch]); % Ga * Sigma_alpha *Ga' 0331 SB = GSG + Cov; % Sigma_B 0332 SB_inv = inv( SB ); % Sigam_B^{-1} 0333 SBBS = SB_inv*BB*SB_inv; % Sigma_B^{-1}*BB*Sigma_B^{-1} 0334 0335 % Reconstraction error = sum( (B-G*J)'*Cov*(B-G*J) ) 0336 tr1_b = tr1_b + sum(sum(Cov.*SBBS)); 0337 % Error variance 0338 tr2_b = tr2_b + sum(sum(BB.*SB_inv)); 0339 0340 % Current magnitude = diag( sum( Z * Z ) ) 0341 mean_jj = mean_jj + ( reshape(SBBS,[1 Nch2]) * GG ).*ax_j; 0342 % Estimation gain for Z = diag(ax_z*Gt*SB_inv*G) 0343 gain_j = gain_j + Nt*Ntry*( reshape(SB_inv,[1 Nch2]) * GG ).*ax_j; 0344 0345 if ( mod(k, Nskip)==0 || k==Ntrain ) 0346 % Model complexity for Current parameter 0347 Hj = Hj - 0.5*Nt*Ntry*vb_log_det(SB); 0348 end 0349 end; 0350 0351 % Noise variance estimation 0352 % sx_total = (tr1_b + mean_jj + sum(mean_zz)); 0353 sx_total = tr2_b; 0354 0355 % update the observation noise variance 0356 if update_sx 0357 sx = sx_total/(Nt*Nerror); 0358 end 0359 0360 %%%%% VB M-step -- variance update 0361 0362 % Save old value for Free energy calculation 0363 ax_j_old = ax_j; 0364 0365 % variance estimation 0366 mean_jj_all = 0.5*ax_j.*mean_jj/sx + gx0_j.*ax0_j; 0367 0368 % Update the current variance parameters 0369 % VB update equation 0370 % gx_j * ax_j = ... 0371 % (0.5*Ntt*Njall + gx0_j) * ax_j = ... 0372 % mean_jj_all + 0.5*ax_j*( Ntt * Njall - gain_j ) ; 0373 0374 % Update the current variance parameters 0375 if k <= Npre_train 0376 % --- VB update rule 0377 ax_j_total = mean_jj_all + 0.5*ax_j.*(Ntt*Njall - gain_j) ; 0378 0379 % Background activity variance estimation 0380 if update_v 0381 ax_j(1) = ax_j_total(1) ./ gx_j(1); 0382 end 0383 % Active current variance estimation 0384 if update_z 0385 ax_j(act_id) = ax_j_total(act_id) ./ gx_j(act_id); 0386 end 0387 else 0388 % --- Faster convergence update rule 0389 gain_j = max(gain_j, eps); 0390 0391 % Background activity variance estimation 0392 if update_v 0393 ax_j(1) = mean_jj_all(1) ./ (0.5*gain_j(1) + gx0_j(1)); 0394 end 0395 % Active current variance estimation 0396 if update_z 0397 ax_j(act_id) = mean_jj_all(act_id) ... 0398 ./ (0.5*gain_j(act_id) + gx0_j(act_id)); 0399 end 0400 end 0401 0402 %%%%%%%% Free Energy Calculation %%%%%%%% 0403 0404 % Model complexity for current variance 0405 ix = find( gx0_j > eps & ax0_j > eps ); 0406 0407 if ~isempty(ix) 0408 Haj = sum(gx0_j(ix).*( log(ax0_j(ix)./ax_j_old(ix)) ... 0409 - ax0_j(ix)./ax_j_old(ix) + 1)); 0410 else 0411 Haj = 0; 0412 end; 0413 0414 LP = - 0.5*Nt*Nerror*log(sx); 0415 0416 % Data Likelihood & Free Energy 0417 Ea = - 0.5*( sx_total/sx - Nt*Nerror ); 0418 FE = LP + Hj + Ea + Haj; 0419 Ev = LP + Hj; % evidence (without constants) 0420 Err = tr1_b/sum(sum(BB));% Normalized Error 0421 0422 % DEBUG Info 0423 k_check = k_check + 1; 0424 Info(k_check,:) = [... 0425 FE, Ev, Hj, LP, Err,... 0426 sx, max(ax_j), sum(gain_j)]; 0427 % end of free energy calculation 0428 0429 if mod(k, Nskip)==0 || k==Ntrain 0430 vb_disp(['Tn=' num2str(nw,'%3d') ', Iter=' num2str(k,'%4d') ... 0431 ', FE=' num2str(FE,'%6.2f') ', Error=' ... 0432 num2str(Err,'%6.2e')],VERBOSE_LEVEL_INFO); 0433 end; 0434 0435 end % end of learning loop 0436 0437 sxall(nw) = sx; 0438 vall(:,nw) = ax_j(:); 0439 0440 if ~isempty(act_id) 0441 aall(:,nw) = act*ax_j(act_id)'; 0442 end 0443 end % end of time window loop 0444 0445 0446 Model.sx = sxall; % Sensor noise variance 0447 Model.a = aall ; % Current variance 0448 0449 % Current variance for background activity 0450 if Njall == 0, 0451 Model.v = zeros(1,Nwindow) ; 0452 else 0453 Model.v = vall(1,:) ; 0454 end 0455 0456 % Model.act : Multiple fMRI activity patterns for focal area 0457 % Model.act_v : Their estimated weights 0458 if ~isempty(act_id) 0459 Model.act = act; 0460 Model.act_v = vall(act_id,:); 0461 end 0462 0463 %fprintf('Wiener estimation for (multi) fmri patterns\n') 0464 0465 return 0466 0467 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0468 %%% functions specific for this function 0469 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0470 function FEconst = free_energy_constant(gx_j, gx0_j, gx_z, gx0_z,... 0471 Nt, Ntotaltrials, Nerror, Njall) 0472 0473 gx_sx = 0.5*Nt*Nerror; 0474 L0 = 0.5*Nt*(Njall*Ntotaltrials*qG(gx_j) + Ntotaltrials*sum(qG(gx_z)) + Nerror*qG(gx_sx) - Nerror); 0475 H0 = HG0(gx_sx,0) + HG0(gx_j,gx0_j)+sum(HG0(gx_z, gx0_z)); 0476 FEconst = L0+H0; 0477 0478 %%%%% 0479 function y = qG(gx) 0480 N=length(gx); 0481 nz = find(gx ~= 0); % non-zero components 0482 y = zeros(N,1); 0483 y_tmp = psi(gx(nz)) - log(gx(nz)); 0484 y(nz) = y_tmp; 0485 0486 %%%%%% 0487 function y =HG0(gx, gx0) 0488 N=length(gx0); 0489 nz = find(gx ~= 0); % non-zero components 0490 nz0 = find(gx0 ~= 0); % non-zero components 0491 % gammma 0492 y1 = zeros(N,1); 0493 y1_tmp = gammaln(gx(nz)) - gx(nz).*psi(gx(nz)) + gx(nz); 0494 y1(nz) = y1_tmp; 0495 % gamma0 0496 y2 = zeros(N,1); 0497 y2_tmp = gammaln(gx0(nz0)) - gx0(nz0).*log(gx0(nz0)) + gx0(nz0); 0498 y2(nz0) = y2_tmp; 0499 % gamma*gamma0 0500 y3 = zeros(N,1); 0501 y3_tmp = gx0(nz0).*(psi(gx(nz0))-log(gx(nz0))); 0502 y3(nz0) = y3_tmp; 0503 0504 y = y1 - y2 + y3; 0505 0506 0507 0508