Home > vbmeg > functions > estimation > bayes > vb_current_calc.m

vb_current_calc

PURPOSE ^

Calculate current using inverse filter

SYNOPSIS ^

function [Zp ,post]= vb_current_calc(Bt,KW,SB_cov,Hj,sx)

DESCRIPTION ^

 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005