Calculate current using inverse filter [Zp]= vb_current_calc(Bt,KW) [Zp ,post]= vb_current_calc(Bt,KW,SB_cov,Hj,sx) --- Input Bt = MEG/EEG data (Nch,Tsample) KW = VB inverse filter (Njact,Nch,1,Ntask) SB_cov = posterior sensor covariance matrix (Nch,Nch,1,Ntask) Hj = Model entropy (1,Ntask) sx = Observation noise variance (1,Ntask) --- Output Zp = current of Njact vertices (Njact,Tsample) post = posterior probability for each condition (Ntask,1) 2008-11-3 Masa-aki Sato Copyright (C) 2011, ATR All Rights Reserved. License : New BSD License(see VBMEG_LICENSE.txt)
0001 function [Zp ,post]= vb_current_calc(Bt,KW,SB_cov,Hj,sx) 0002 % Calculate current using inverse filter 0003 % [Zp]= vb_current_calc(Bt,KW) 0004 % [Zp ,post]= vb_current_calc(Bt,KW,SB_cov,Hj,sx) 0005 % --- Input 0006 % Bt = MEG/EEG data (Nch,Tsample) 0007 % KW = VB inverse filter (Njact,Nch,1,Ntask) 0008 % SB_cov = posterior sensor covariance matrix (Nch,Nch,1,Ntask) 0009 % Hj = Model entropy (1,Ntask) 0010 % sx = Observation noise variance (1,Ntask) 0011 % --- Output 0012 % Zp = current of Njact vertices (Njact,Tsample) 0013 % post = posterior probability for each condition (Ntask,1) 0014 % 0015 % 2008-11-3 Masa-aki Sato 0016 % 0017 % Copyright (C) 2011, ATR All Rights Reserved. 0018 % License : New BSD License(see VBMEG_LICENSE.txt) 0019 0020 Ntask = size(KW,4); 0021 0022 if Ntask==1, 0023 % Current 0024 Zp = (KW(:,:,1,1) * Bt); 0025 post = 1; 0026 else 0027 MinExp = -50; 0028 Njact = size(KW,1); 0029 Nt = size(Bt,2); 0030 LP = zeros(Ntask,1); 0031 0032 BB = Bt*Bt'; 0033 0034 for k=1:Ntask 0035 % Data covariance matching 0036 err = sum(sum(BB.*SB_cov(:,:,1,k))); 0037 % Data Likelihood & Free Energy 0038 LP(k) = - 0.5*(err)/sx(k) + Nt*Hj(k); 0039 end; 0040 0041 LP = max( LP - max(LP) , MinExp); 0042 post = exp(LP); 0043 post = post/sum(post); 0044 0045 Zp = zeros(Njact,Nt); 0046 0047 for k=1:Ntask 0048 % Current 0049 Zt = (KW(:,:,1,k) * Bt); 0050 % Current 0051 Zp = Zp + Zt*post(k); 0052 end 0053 end