Estimate parameters of hierarchical Bayse model using variational Bayes
from both MEG and EEG (VBMEG public function)
(modified of vb_job_vb.m)
[syntax]
vb_job_vb_meeg(bayes_parm)
vb_job_vb_meeg(proj_root,bayes_parm) [old style]
[input]
proj_root : <<string>> VBMEG project root directory.
bayes_parm: <<struct>> Parameters for variational Bayes method.
--- fields of bayes_parm
(filename parameters)
brainfile : <<string>> Cortical surface model file (.brain.mat).
actfile : <<string>> Cortical activity map file (.act.mat).
areafile : <<string>> Cortical area file (.area.mat).
megfile : <<cell>> Cell of MEG-MAT file (.meg.mat).
megfile_baseline: <<cell>> Cell of MEG-MAT file (.meg.mat).
basisfile_meg : <<cell>> Cell of leadfield file (.basis.mat).
eegfile : <<cell>> Cell of MEG-MAT file (.meg.mat).
eegfile_baseline: <<cell>> Cell of MEG-MAT file (.meg.mat).
basisfile_eeg : <<cell>> Cell of leadfield file (.basis.mat).
bayesfile : <<string>> Current variance file (.bayes.mat).
(noise model parameters)
noise_model: <<int>> Specify observation noise covariance matrix. It
should be noted that the covariance matrix is normalized
so that sum(diag(Cov)) is equal to the number of
sensors. In VB algorithm, a scale parameter of
observation noise variance (beta^{-1}) is estimated with
the covariance matrix (constant).
noise_model=1: Identity matrix. MEG data is not used.
noise_model=2: Only diagonal elements are estimated from MEG data.
noise_model=3: Full covariance matrix is estimated from MEG data.
noise_reg : <<float>> Regularization parameter for observarion noise
covariance matrix. This value is used to suppress
singularity of the covariance matrix. It is effective
when noise_mode=2 or 3.
(prior parameters: see Yoshioka et al., 2008, Appendix A)
act_key: <<string>> ID of activity map, which is used to specify prior
information for current variance, corresponding to
\hat{t}_{i}. The activity map is automatically normalized
(maximum=1) in estimation procedure.
a0 : <<double>>
a0_act : <<double>> Maximum value of current variance parameter. If
a0=1, this value is the same with a variance magnification
parameter m_{0}.
prior_weight : Relative influence of prior information on data
(from 0 to 1 real).
<<The following two parameters are dependent variables of prior_weight>>
Ta0 : <<double>> Confidence parameter of prior current variance
parameters (gamma distribution) at which the activity map
('act_key') takes the minimum value.
Ta0_act: <<double>> Confidence parameter of prior current variance
parameters (gamma distribution) at which the activity map
('act_key') takes the maximum value.
v0 : <<double>>
Tv0 : <<double>>
(time window parameters)
twin_noise : <<vector>> Observation noise covariance matrix is
calculated using MEG data "bayes_parm.megfile_baseline"
within this time window, specified by start and end.
twin_baseline: <<vector>> Baseline observation noise variance
(rho_{base}) and baseline current variance parameter
(nu_{0,base}) are estimated from MEG data
"megfile_baseline" within the time window specified by
this parameter. See Yoshioka et al., 2008, Appendix A.
twin_meg : <<vector>> Current variance parameters (alpha^{-1}) are
estimated from MEG data "bayes_parm.megfile" within the
time window specified by this parameter.
Tperiod : <<int/vector>> Partial time window length of MEG data.
Tnext : <<int/vector>> Moving step length of partial time
window of MEG data.
(forward model parameters: see Yoshioka et al., 2008, Appendix B)
forward_model : <<string>> Set to 'focal'.
area_key : <<string>> ID of cortical area at which current
dipoles are assumed.
reduce : <<float>> Current dipole reduction ratio (0-1).
Rfilt : <<float>> Smoothing filter radius in [m]. Elements of
matrix W is determined from this parameter and
inter-dipole distance.
.area_key_global: <<string>> ID of cortical area at whith current
dipoles are assumed in estimation of baseline
variance parameter (nu_{0,base}).
.reduce_global : <<float>> Current dipole reduction ratio in
estimation of baseline variance parameter
(nu_{0,base}) (0-1).
.Rfilt_global : <<float>> Smoothing filter radius in [m] in
estimation of baseline variance parameter
(nu_{0,base}).
.patch_norm : <<bool>> If true, patch area normalization is
applied. That is, dipole current moment (Am) is
replaced by its (area) density (Am/m^2).
.expand_spatial_filter: <<bool>> Set to 'true'.
.spatial_smoother : <<string>> smoothing leadfield
= 'std' : using Standard brain(nextIX, nextDD, xx)
= 'subj': using Individual brain(subj.nextIX,
subj.nextDD, subj.xx)
.remove_crossed_area : <<bool>> If true, cortical area specified by
'area_key' is removed in estimation of baseline
variance parameter (nu_{0,base}).
(extra-dipole parameters)
extra : <optional> <<struct>> Parameter set of
extra-dipole.
extra.basisfile : <<string>> Filename of leadfield file (.basis.mat)
for extra-dipole.
extra.a0 : <<vector>> Alpha_{0}^{-1} for all extra-dipoles.
extra.Ta0 : <<vector>> Gamma_{0} for all extra-dipoles.
extra.a0_physical: <<bool>> If true, a0 is interpreted as physical
unit (Am^2, current variance). Specifically,
vb_parm.a0 is set to extra.a0/sx0, where sx0 is the
baseline of observation noise variance.
(algorithm parameters)
Ntrain : <<int>> Iteration number of VB algorithm.
Npre_train: <<int>> Iteration number of fast VB algorithm.
skip : <<int>> Calculate free energy at each 'skip' steps.
update_sx : <<bool>> If true, observarion noise variance is
simultaneously updated with prior current variance
parameter (alpha_{0}).
update_v : (bool)
Fdmin : <<float>> Constant for convergence test of VB algorithm.
a_min : <<float>> Minimum value of alpha_{0}^{-1} for stability of
the VB alrogithm.
a_max : <<float>> Maximum value of alpha_{0}^{-1} for stability of
the VB alrogithm.
cont_pr : <<bool>> If true, estimated current variance parameters
(alpha^{-1}) are used as the initial value of current
variance parameters in the next time window.
(Parameters for soft normal constraints (not publicly supported in
current version of VBMEG))
soft_mode : <<bool>>
variance_orientation: <<bool>>
var_max : <<float>>
tan_var : <<float>>
cosval : <<float>>
(preprocessing parameters (not recommended to use))
.trial_average : <<bool>> If true, trial average is applied to the
original MEG data.
.temporal_filter: <<bool>> If true, temporal smoothing filter is
applied to the original MEG data.
[output]
Estimation result will be saved to the file specified by
bayes_parm.bayesfile with the following variables:
Model.a : <<vector>> (alpha_{i}^{-1})
Model.sx : <<float>> (beta^{-1})
Model.v : <<float>> (nu_{0,base})
Model.a0 : <<vector>> (alpha_{0i}^{-1}).
Model.sx0 : <<float>> (rho_{base}=beta_{base}^{-1})
Model.v0 : <<float?>> (nu_{0,base})
Model.ix : <<vector>> Dipole indices of Model.a.
Model.Cov : <<cell?>> Noise covariance matrix.
bayes_parm: <<struct>> All estimation parameter.
vb_parm : <<struct>> Parameters for VB algorithm.
Info : <<struct>> Free energy, etc.
Pointlist : <<struct>> Extra dipole information.
[note]
1. 'bayes_parm.reduce', J-current and Z-current
Cortical current J is represented by an auxillary parameter Z using
equation J=WZ, where W is Gaussian smoothing matrix (see Yoshioka et
al., 2008, Appendix A). The number of current dipoles I (=size(J,1))
is reduced to ('bayes_parm.reduce' times I), which is the same with
the number of current variance parameters (size(Model.a,1)). J and Z
are referred to as J-current and Z-current, respectively.
2. Time window parameters
Firstly, it is noted that the unit of time window parameters (e.g.,
'twin_meg') is the sampling point depending on the sampling frequency
of M/EEG measurement. In addition, these values are NOT relative to
the pre-trigger, but relative to the start of each trial. For
simplicity, we assume 1 ms sampling data in the following
explanation.
MEG data is separated into a number of time windows and model
parameters are estimated for each time window. The following setting
is a simple example:
>> bayes_parm.twin_meg = [501 1000];
>> bayes_parm.Tperiod = 500;
>> bayes_parm.Tnext = 500;
In this case, model parameters are estimated from (500 x N) samples
501 ms after each trial, where N denotes the number of trials. As a
result, size(Model.a) becomes [I 1], where I denotes the number of
current dipoles of Z-current. In the next example, whole time window
is separated into 5 time windows:
>> bayes_parm.twin_meg = [501 1000];
>> bayes_parm.Tperiod = [100 100 100 100 100];
>> bayes_parm.Tnext = [1 101 201 301 401];
The start of time windows are 1, 101, ..., 401 ms and the length is
100 ms. The following parameter set is the same with the above
example:
>> bayes_parm.twin_meg = [501 1000];
>> bayes_parm.Tperiod = 100;
>> bayes_parm.Tnext = 100;
'Tperiod' and 'Tnext' are interpreted as the length and moving steps
of the time window, respectively. Time windows can be overlapped:
>> bayes_parm.twin_meg = [501 1000];
>> bayes_parm.Tperiod = [200 200 200 200];
>> bayes_parm.Tnext = [1 101 201 301];
Within the overlapped periods, dipole currents are calculated by the
average of dipole currents using estimated parameters corresponding to
overlapped time windows. This is implemented not in job_vb.m, but in
current estimation functions.
Different length of time windows can be specified:
>> bayes_parm.twin_meg = [401 1000];
>> bayes_parm.Tperiod = [100 500];
>> bayes_parm.Tnext = [1 101];
In this case, model parameters are estimated separately for
pre-trigger (100 ms) and post-trigger (500 ms) of MEG data.
[reference]
Sato et al. (2004) NeuroImage 23, 806-826.
Yoshioka et al. (2008) NeuroImage 42, 1397-1413.
[history]
2017-06-15 Yusuke Fujiwara refined
2017-03-16 rhayashi
spatial_smoother is added.
2017-03-10 Yusuke Fujiwara
MEG+EEG
2008-06-27 Taku Yoshioka
2008-07-07 Masa-aki Sato
Model.Cov is added
2008-08-08 Taku Yoshioka
Support for physical unit for extra dipole.
2008-10-06 Taku Yoshioka
Bug fix (physical unit for extra dipole)
2008-12-03 Taku Yoshioka
Use vb_disp() for displaying messages
2017-11-14 oyamashi
prior_weight implementation.
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function vb_job_vb_meeg(varargin) 0002 % Estimate parameters of hierarchical Bayse model using variational Bayes 0003 % from both MEG and EEG (VBMEG public function) 0004 % (modified of vb_job_vb.m) 0005 % 0006 % [syntax] 0007 % vb_job_vb_meeg(bayes_parm) 0008 % vb_job_vb_meeg(proj_root,bayes_parm) [old style] 0009 % 0010 % [input] 0011 % proj_root : <<string>> VBMEG project root directory. 0012 % bayes_parm: <<struct>> Parameters for variational Bayes method. 0013 % --- fields of bayes_parm 0014 % (filename parameters) 0015 % brainfile : <<string>> Cortical surface model file (.brain.mat). 0016 % actfile : <<string>> Cortical activity map file (.act.mat). 0017 % areafile : <<string>> Cortical area file (.area.mat). 0018 % megfile : <<cell>> Cell of MEG-MAT file (.meg.mat). 0019 % megfile_baseline: <<cell>> Cell of MEG-MAT file (.meg.mat). 0020 % basisfile_meg : <<cell>> Cell of leadfield file (.basis.mat). 0021 % eegfile : <<cell>> Cell of MEG-MAT file (.meg.mat). 0022 % eegfile_baseline: <<cell>> Cell of MEG-MAT file (.meg.mat). 0023 % basisfile_eeg : <<cell>> Cell of leadfield file (.basis.mat). 0024 % 0025 % bayesfile : <<string>> Current variance file (.bayes.mat). 0026 % 0027 % (noise model parameters) 0028 % noise_model: <<int>> Specify observation noise covariance matrix. It 0029 % should be noted that the covariance matrix is normalized 0030 % so that sum(diag(Cov)) is equal to the number of 0031 % sensors. In VB algorithm, a scale parameter of 0032 % observation noise variance (beta^{-1}) is estimated with 0033 % the covariance matrix (constant). 0034 % noise_model=1: Identity matrix. MEG data is not used. 0035 % noise_model=2: Only diagonal elements are estimated from MEG data. 0036 % noise_model=3: Full covariance matrix is estimated from MEG data. 0037 % noise_reg : <<float>> Regularization parameter for observarion noise 0038 % covariance matrix. This value is used to suppress 0039 % singularity of the covariance matrix. It is effective 0040 % when noise_mode=2 or 3. 0041 % 0042 % (prior parameters: see Yoshioka et al., 2008, Appendix A) 0043 % act_key: <<string>> ID of activity map, which is used to specify prior 0044 % information for current variance, corresponding to 0045 % \hat{t}_{i}. The activity map is automatically normalized 0046 % (maximum=1) in estimation procedure. 0047 % a0 : <<double>> 0048 % a0_act : <<double>> Maximum value of current variance parameter. If 0049 % a0=1, this value is the same with a variance magnification 0050 % parameter m_{0}. 0051 % 0052 % prior_weight : Relative influence of prior information on data 0053 % (from 0 to 1 real). 0054 % <<The following two parameters are dependent variables of prior_weight>> 0055 % Ta0 : <<double>> Confidence parameter of prior current variance 0056 % parameters (gamma distribution) at which the activity map 0057 % ('act_key') takes the minimum value. 0058 % Ta0_act: <<double>> Confidence parameter of prior current variance 0059 % parameters (gamma distribution) at which the activity map 0060 % ('act_key') takes the maximum value. 0061 % v0 : <<double>> 0062 % Tv0 : <<double>> 0063 % 0064 % (time window parameters) 0065 % twin_noise : <<vector>> Observation noise covariance matrix is 0066 % calculated using MEG data "bayes_parm.megfile_baseline" 0067 % within this time window, specified by start and end. 0068 % twin_baseline: <<vector>> Baseline observation noise variance 0069 % (rho_{base}) and baseline current variance parameter 0070 % (nu_{0,base}) are estimated from MEG data 0071 % "megfile_baseline" within the time window specified by 0072 % this parameter. See Yoshioka et al., 2008, Appendix A. 0073 % twin_meg : <<vector>> Current variance parameters (alpha^{-1}) are 0074 % estimated from MEG data "bayes_parm.megfile" within the 0075 % time window specified by this parameter. 0076 % Tperiod : <<int/vector>> Partial time window length of MEG data. 0077 % Tnext : <<int/vector>> Moving step length of partial time 0078 % window of MEG data. 0079 % 0080 % (forward model parameters: see Yoshioka et al., 2008, Appendix B) 0081 % forward_model : <<string>> Set to 'focal'. 0082 % area_key : <<string>> ID of cortical area at which current 0083 % dipoles are assumed. 0084 % reduce : <<float>> Current dipole reduction ratio (0-1). 0085 % Rfilt : <<float>> Smoothing filter radius in [m]. Elements of 0086 % matrix W is determined from this parameter and 0087 % inter-dipole distance. 0088 % .area_key_global: <<string>> ID of cortical area at whith current 0089 % dipoles are assumed in estimation of baseline 0090 % variance parameter (nu_{0,base}). 0091 % .reduce_global : <<float>> Current dipole reduction ratio in 0092 % estimation of baseline variance parameter 0093 % (nu_{0,base}) (0-1). 0094 % .Rfilt_global : <<float>> Smoothing filter radius in [m] in 0095 % estimation of baseline variance parameter 0096 % (nu_{0,base}). 0097 % .patch_norm : <<bool>> If true, patch area normalization is 0098 % applied. That is, dipole current moment (Am) is 0099 % replaced by its (area) density (Am/m^2). 0100 % .expand_spatial_filter: <<bool>> Set to 'true'. 0101 % .spatial_smoother : <<string>> smoothing leadfield 0102 % = 'std' : using Standard brain(nextIX, nextDD, xx) 0103 % = 'subj': using Individual brain(subj.nextIX, 0104 % subj.nextDD, subj.xx) 0105 % .remove_crossed_area : <<bool>> If true, cortical area specified by 0106 % 'area_key' is removed in estimation of baseline 0107 % variance parameter (nu_{0,base}). 0108 % 0109 % (extra-dipole parameters) 0110 % extra : <optional> <<struct>> Parameter set of 0111 % extra-dipole. 0112 % extra.basisfile : <<string>> Filename of leadfield file (.basis.mat) 0113 % for extra-dipole. 0114 % extra.a0 : <<vector>> Alpha_{0}^{-1} for all extra-dipoles. 0115 % extra.Ta0 : <<vector>> Gamma_{0} for all extra-dipoles. 0116 % extra.a0_physical: <<bool>> If true, a0 is interpreted as physical 0117 % unit (Am^2, current variance). Specifically, 0118 % vb_parm.a0 is set to extra.a0/sx0, where sx0 is the 0119 % baseline of observation noise variance. 0120 % 0121 % (algorithm parameters) 0122 % Ntrain : <<int>> Iteration number of VB algorithm. 0123 % Npre_train: <<int>> Iteration number of fast VB algorithm. 0124 % skip : <<int>> Calculate free energy at each 'skip' steps. 0125 % update_sx : <<bool>> If true, observarion noise variance is 0126 % simultaneously updated with prior current variance 0127 % parameter (alpha_{0}). 0128 % update_v : (bool) 0129 % Fdmin : <<float>> Constant for convergence test of VB algorithm. 0130 % a_min : <<float>> Minimum value of alpha_{0}^{-1} for stability of 0131 % the VB alrogithm. 0132 % a_max : <<float>> Maximum value of alpha_{0}^{-1} for stability of 0133 % the VB alrogithm. 0134 % cont_pr : <<bool>> If true, estimated current variance parameters 0135 % (alpha^{-1}) are used as the initial value of current 0136 % variance parameters in the next time window. 0137 % 0138 % (Parameters for soft normal constraints (not publicly supported in 0139 % current version of VBMEG)) 0140 % soft_mode : <<bool>> 0141 % variance_orientation: <<bool>> 0142 % var_max : <<float>> 0143 % tan_var : <<float>> 0144 % cosval : <<float>> 0145 % 0146 % (preprocessing parameters (not recommended to use)) 0147 % .trial_average : <<bool>> If true, trial average is applied to the 0148 % original MEG data. 0149 % .temporal_filter: <<bool>> If true, temporal smoothing filter is 0150 % applied to the original MEG data. 0151 % 0152 % [output] 0153 % Estimation result will be saved to the file specified by 0154 % bayes_parm.bayesfile with the following variables: 0155 % Model.a : <<vector>> (alpha_{i}^{-1}) 0156 % Model.sx : <<float>> (beta^{-1}) 0157 % Model.v : <<float>> (nu_{0,base}) 0158 % Model.a0 : <<vector>> (alpha_{0i}^{-1}). 0159 % Model.sx0 : <<float>> (rho_{base}=beta_{base}^{-1}) 0160 % Model.v0 : <<float?>> (nu_{0,base}) 0161 % Model.ix : <<vector>> Dipole indices of Model.a. 0162 % Model.Cov : <<cell?>> Noise covariance matrix. 0163 % bayes_parm: <<struct>> All estimation parameter. 0164 % vb_parm : <<struct>> Parameters for VB algorithm. 0165 % Info : <<struct>> Free energy, etc. 0166 % Pointlist : <<struct>> Extra dipole information. 0167 % 0168 % [note] 0169 % 1. 'bayes_parm.reduce', J-current and Z-current 0170 % Cortical current J is represented by an auxillary parameter Z using 0171 % equation J=WZ, where W is Gaussian smoothing matrix (see Yoshioka et 0172 % al., 2008, Appendix A). The number of current dipoles I (=size(J,1)) 0173 % is reduced to ('bayes_parm.reduce' times I), which is the same with 0174 % the number of current variance parameters (size(Model.a,1)). J and Z 0175 % are referred to as J-current and Z-current, respectively. 0176 % 0177 % 2. Time window parameters 0178 % Firstly, it is noted that the unit of time window parameters (e.g., 0179 % 'twin_meg') is the sampling point depending on the sampling frequency 0180 % of M/EEG measurement. In addition, these values are NOT relative to 0181 % the pre-trigger, but relative to the start of each trial. For 0182 % simplicity, we assume 1 ms sampling data in the following 0183 % explanation. 0184 % MEG data is separated into a number of time windows and model 0185 % parameters are estimated for each time window. The following setting 0186 % is a simple example: 0187 % >> bayes_parm.twin_meg = [501 1000]; 0188 % >> bayes_parm.Tperiod = 500; 0189 % >> bayes_parm.Tnext = 500; 0190 % In this case, model parameters are estimated from (500 x N) samples 0191 % 501 ms after each trial, where N denotes the number of trials. As a 0192 % result, size(Model.a) becomes [I 1], where I denotes the number of 0193 % current dipoles of Z-current. In the next example, whole time window 0194 % is separated into 5 time windows: 0195 % >> bayes_parm.twin_meg = [501 1000]; 0196 % >> bayes_parm.Tperiod = [100 100 100 100 100]; 0197 % >> bayes_parm.Tnext = [1 101 201 301 401]; 0198 % The start of time windows are 1, 101, ..., 401 ms and the length is 0199 % 100 ms. The following parameter set is the same with the above 0200 % example: 0201 % >> bayes_parm.twin_meg = [501 1000]; 0202 % >> bayes_parm.Tperiod = 100; 0203 % >> bayes_parm.Tnext = 100; 0204 % 'Tperiod' and 'Tnext' are interpreted as the length and moving steps 0205 % of the time window, respectively. Time windows can be overlapped: 0206 % >> bayes_parm.twin_meg = [501 1000]; 0207 % >> bayes_parm.Tperiod = [200 200 200 200]; 0208 % >> bayes_parm.Tnext = [1 101 201 301]; 0209 % Within the overlapped periods, dipole currents are calculated by the 0210 % average of dipole currents using estimated parameters corresponding to 0211 % overlapped time windows. This is implemented not in job_vb.m, but in 0212 % current estimation functions. 0213 % Different length of time windows can be specified: 0214 % >> bayes_parm.twin_meg = [401 1000]; 0215 % >> bayes_parm.Tperiod = [100 500]; 0216 % >> bayes_parm.Tnext = [1 101]; 0217 % In this case, model parameters are estimated separately for 0218 % pre-trigger (100 ms) and post-trigger (500 ms) of MEG data. 0219 % 0220 % [reference] 0221 % Sato et al. (2004) NeuroImage 23, 806-826. 0222 % Yoshioka et al. (2008) NeuroImage 42, 1397-1413. 0223 % 0224 % [history] 0225 % 2017-06-15 Yusuke Fujiwara refined 0226 % 2017-03-16 rhayashi 0227 % spatial_smoother is added. 0228 % 2017-03-10 Yusuke Fujiwara 0229 % MEG+EEG 0230 % 2008-06-27 Taku Yoshioka 0231 % 2008-07-07 Masa-aki Sato 0232 % Model.Cov is added 0233 % 2008-08-08 Taku Yoshioka 0234 % Support for physical unit for extra dipole. 0235 % 2008-10-06 Taku Yoshioka 0236 % Bug fix (physical unit for extra dipole) 0237 % 2008-12-03 Taku Yoshioka 0238 % Use vb_disp() for displaying messages 0239 % 2017-11-14 oyamashi 0240 % prior_weight implementation. 0241 % Copyright (C) 2011, ATR All Rights Reserved. 0242 % License : New BSD License(see VBMEG_LICENSE.txt) 0243 0244 if length(varargin) == 1 0245 proj_root = []; 0246 bayes_parm = varargin{1}; 0247 elseif length(varargin) == 2 0248 proj_root = varargin{1}; 0249 bayes_parm = varargin{2}; 0250 end 0251 0252 proj_root = vb_rm_trailing_slash(proj_root); 0253 0254 const = vb_define_verbose; 0255 VERBOSE_LEVEL_INFO = const.VERBOSE_LEVEL_INFO; 0256 0257 % 0258 % Change path to absolute path 0259 % 0260 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0261 vb_disp('----- New VBMEG -----',VERBOSE_LEVEL_INFO) 0262 if ~isfield(bayes_parm, 'spatial_smoother') 0263 dp = vb_set_bayes_default_parameters; 0264 bayes_parm.spatial_smoother = dp.spatial_smoother; 0265 end 0266 bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm); 0267 0268 Nsession = length(bayes_parm_abs.eegfile_baseline); 0269 for nf = 1:Nsession 0270 bayes_parm_abs.eegfile_baseline{nf} = fullfile(proj_root, bayes_parm_abs.eegfile_baseline{nf}); 0271 bayes_parm_abs.eegfile{nf} = fullfile(proj_root, bayes_parm_abs.eegfile{nf}); 0272 end 0273 bayes_parm_abs.basisfile_meg = fullfile(proj_root, bayes_parm_abs.basisfile_meg); 0274 bayes_parm_abs.basisfile_eeg = fullfile(proj_root, bayes_parm_abs.basisfile_eeg); 0275 0276 0277 % 0278 % Preparation of lead fields 0279 % 0280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0281 vb_disp('----- Leadfield matrix preparation',VERBOSE_LEVEL_INFO); 0282 0283 lf_parm.brainfile = bayes_parm_abs.brainfile; 0284 lf_parm.areafile = bayes_parm_abs.areafile; 0285 lf_parm.patch_norm = bayes_parm_abs.patch_norm; 0286 lf_parm.expand_spatial_filter = bayes_parm_abs.expand_spatial_filter; 0287 lf_parm.spatial_smoother = bayes_parm_abs.spatial_smoother; 0288 0289 % % Global window 0290 % vb_disp('--- Lead field matrix of global window ',VERBOSE_LEVEL_INFO); 0291 % lf_parm.basisfile = bayes_parm_abs.basisfile_global; 0292 % lf_parm.area_key = bayes_parm_abs.area_key_global; 0293 % lf_parm.reduce = bayes_parm_abs.reduce_global; 0294 % lf_parm.Rfilt = bayes_parm_abs.Rfilt_global; 0295 % if strcmp(bayes_parm_abs.forward_model,'focal+global')==1 ... 0296 % && bayes_parm_abs.remove_crossed_area==ON, 0297 % % Remove focal area from global area 0298 % lf_parm.remove_area_key = bayes_parm_abs.area_key; 0299 % else 0300 % lf_parm.remove_area_key = []; 0301 % end 0302 % [Gall_meg, ix_all] = vb_leadfield_preparation(lf_parm); 0303 0304 0305 %%%%%%%%%%%%%%%%%%% MEG %%%%%%%%%%%%%%%%%%%%% 0306 % Focal window 0307 vb_disp('--- Lead field matrix of focal window ',VERBOSE_LEVEL_INFO); 0308 lf_parm.basisfile = bayes_parm_abs.basisfile_meg; 0309 lf_parm.area_key = bayes_parm_abs.area_key; 0310 lf_parm.reduce = bayes_parm_abs.reduce; 0311 lf_parm.Rfilt = bayes_parm_abs.Rfilt; 0312 lf_parm.remove_area_key = []; 0313 [Gact_meg, ix_act_meg] = vb_leadfield_preparation(lf_parm); 0314 0315 0316 %%%%%%%%%%%%%%%%%%% EEG %%%%%%%%%%%%%%%%%%%%% 0317 lf_parm.basisfile = bayes_parm_abs.basisfile_eeg; 0318 [Gact_eeg, ix_act_eeg] = vb_leadfield_preparation(lf_parm); 0319 0320 ix_act = ix_act_meg; 0321 0322 0323 % 0324 % Normalization of lead fields 0325 % 0326 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0327 0328 % Normalization constant of leadfield 0329 Nmeg = size(Gact_meg{1},1); 0330 Neeg = size(Gact_eeg{1},1); 0331 0332 for n = 1: Nsession 0333 GGall_meg = sum(Gact_meg{n}.^2, 2); 0334 Bnorm_meg(n) = sum(GGall_meg)/length(GGall_meg); 0335 GGall_eeg = sum(Gact_eeg{n}.^2, 2); 0336 Bnorm_eeg(n) = sum(GGall_eeg)/length(GGall_eeg); 0337 end 0338 bsnorm_meg = sqrt(mean(Bnorm_meg)); 0339 bsnorm_eeg = sqrt(mean(Bnorm_eeg)); 0340 0341 % Normalization of leadfield 0342 for n=1:Nsession 0343 nGact_meg{n} = Gact_meg{n}/bsnorm_meg; 0344 nGact_eeg{n} = Gact_eeg{n}/bsnorm_eeg; 0345 nGact{n} = cat(1,nGact_meg{n},nGact_eeg{n}); 0346 end 0347 nGall_meg = nGact_meg; 0348 nGall_eeg = nGact_eeg; 0349 nGall = nGact; 0350 0351 0352 % 0353 % Observation noise estimation 0354 % 0355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0356 vb_disp('----- Observation noise model specification', ... 0357 VERBOSE_LEVEL_INFO); 0358 0359 %Cov Normalized noise covariance matrix 0360 0361 %MEG covariance 0362 tmp_parm = bayes_parm_abs; 0363 %tmp_parm.basisfile = tmp_parm.basisfile_meg; 0364 [Cov_meg,sx0_meg] = vb_observation_noise_specification(tmp_parm); 0365 0366 %EEG covariance 0367 %tmp_parm.basisfile = tmp_parm.basisfile_eeg; 0368 tmp_parm.megfile_baseline = tmp_parm.eegfile_baseline; 0369 [Cov_eeg,sx0_eeg] = vb_observation_noise_specification(tmp_parm); 0370 0371 sx0_meg2 = sx0_meg/(bsnorm_meg.^2); 0372 sx0_eeg2 = sx0_eeg/(bsnorm_eeg.^2); 0373 0374 %ix_meg = 1:size(Cov_meg{1},1); 0375 %ix_eeg = ix_meg(end)+1:ix_meg(end)+size(Cov_eeg{1},1); 0376 0377 for nf = 1:Nsession 0378 Cov{nf} = blkdiag(Cov_meg{nf},(sx0_eeg2/sx0_meg2).*Cov_eeg{nf}); 0379 end 0380 0381 sx0 = sx0_meg2; 0382 0383 % 0384 % MEG data preparation 0385 % 0386 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0387 [B0, Ntrial, tmp, Tsample] = vb_megdata_preparation(bayes_parm_abs); 0388 0389 % EEG data 0390 megfile = bayes_parm_abs.megfile; 0391 bayes_parm_abs.megfile = bayes_parm_abs.eegfile; 0392 [E0] = vb_megdata_preparation(bayes_parm_abs); 0393 0394 for n=1:Nsession 0395 B0{n} = B0{n}/bsnorm_meg; 0396 E0{n} = E0{n}/bsnorm_eeg; 0397 B{n} = cat(1,B0{n},E0{n}); 0398 end 0399 bayes_parm_abs.megfile = megfile; 0400 0401 % 0402 % Check variable consistency 0403 % 0404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0405 vb_check_variable(B, Cov, nGact, nGall); 0406 0407 % 0408 % Estimation of baseline current variance 0409 % 0410 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0411 vb_disp('----- Baseline current variance estimation',VERBOSE_LEVEL_INFO); 0412 0413 % current baseline estimation at pretrigger period 0414 % together with spherical sensor noise estimation 0415 %MEG 0416 tmp_parm = bayes_parm_abs; 0417 tmp_parm.megfile_baseline = tmp_parm.megfile_baseline; 0418 tmp_parm.basisfile{1} = tmp_parm.basisfile_meg; 0419 tmp_parm.basisfile_global{1} = tmp_parm.basisfile_meg; 0420 [v0m, sx0b] = vb_baseline_variance_acquisition(tmp_parm, nGall_meg); 0421 0422 %EEG 0423 tmp_parm.megfile_baseline = tmp_parm.eegfile_baseline; 0424 tmp_parm.basisfile{1} = tmp_parm.basisfile_eeg; 0425 tmp_parm.basisfile_global{1} = tmp_parm.basisfile_eeg; 0426 [v0e, sx0e] = vb_baseline_variance_acquisition(tmp_parm, nGall_eeg); 0427 0428 v0 = (v0m + v0e)/2; 0429 0430 % 0431 % Set VB parameters 0432 % 0433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0434 0435 Njact = size(nGact{1},2); 0436 Njall = size(nGall{1},2); 0437 %vb_parm = vb_set_vb_parm(bayes_parm_abs, B, Njact, Njall, v0, sx0); 0438 vb_parm = vb_set_vb_parm(tmp_parm, B, Njact, Njall, v0, sx0); 0439 0440 % 0441 % Set prior parameters 0442 % 0443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0444 prior_parm.ix = ix_act; 0445 prior_parm.v0 = v0; 0446 prior_parm.a0 = bayes_parm.a0; 0447 prior_parm.a0_act = bayes_parm.a0_act; 0448 0449 w = bayes_parm.prior_weight; 0450 prior_parm.Ta0 = Ntrial * Tsample * w / (2*(1-w)); 0451 prior_parm.Ta0_act = prior_parm.Ta0; 0452 0453 prior_parm.act = vb_get_act(bayes_parm_abs.actfile, ... 0454 bayes_parm_abs.act_key); 0455 [a0,Ta0] = vb_set_vb_prior(prior_parm); 0456 vb_parm.a0 = a0; 0457 vb_parm.Ta0 = Ta0; 0458 0459 % 0460 % Modify prior parameters for soft normal constraint 0461 % 0462 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0463 0464 if bayes_parm.variance_orientation == ON, 0465 soft_parm.ix = ix_act; 0466 soft_parm.a0 = a0; 0467 soft_parm.Ta0 = Ta0; 0468 soft_parm.Norient = vb_parm.Norient; 0469 soft_parm.soft_mode = vb_parm.soft_mode; 0470 soft_parm.brainfile = bayes_parm_abs.brainfile; 0471 soft_parm.var_max = vb_parm.var_max; 0472 soft_parm.tan_var = vb_parm.tan_var; 0473 [a0,Ta0] = set_soft_prior(soft_parm); 0474 vb_parm.a0 = a0; 0475 vb_parm.Ta0 = Ta0; 0476 end 0477 0478 % 0479 % Add extra dipole 0480 % 0481 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0482 0483 if isfield(bayes_parm_abs,'extra') & ~isempty(bayes_parm_abs.extra), 0484 vb_struct2vars(bayes_parm_abs,{'extra'}); 0485 0486 % prior for extra dipole 0487 if extra.a0_physical == true, e0 = extra.a0/sx0.*bsnorm^2; 0488 else e0 = extra.a0*v0; end 0489 vb_parm.a0 = [vb_parm.a0; e0]; 0490 vb_parm.Ta0 = [vb_parm.Ta0; extra.Ta0]; 0491 0492 for i=1:length(nGact) 0493 extra_basis = vb_load_basis(extra.basisfile{i}); 0494 nGact{i} = [nGact{i} (extra_basis./bsnorm)']; 0495 end 0496 load(extra.basisfile{1},'Pointlist'); 0497 Nextra = length(Pointlist); 0498 vb_parm.Njact = vb_parm.Njact+Nextra*3; 0499 vb_parm.Nvact = vb_parm.Nvact+Nextra*3; 0500 else 0501 Pointlist = []; 0502 Nextra = 0; 0503 end 0504 0505 % 0506 % Wiener filter (removed from the previous version) 0507 % 0508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0509 0510 Model = []; 0511 0512 % 0513 % Estimation with focal (and global) window 0514 % 0515 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0516 0517 switch bayes_parm_abs.forward_model 0518 case 'focal+global', 0519 [Model,Info] = vbmeg_ard_estimate2(B, nGall, nGact, ... 0520 Cov, vb_parm, Model); 0521 case 'focal', 0522 vb_parm.v0 = 0; 0523 vb_parm.Tv0 = 0; 0524 [Model,Info] = vbmeg_ard_estimate4(B, [], nGact, ... 0525 Cov, vb_parm, Model); 0526 0527 otherwise, 0528 error(['---- no such forward model : ' ... 0529 bayes_parm_abs.forward_model '----']); 0530 end 0531 0532 Model.sx0 = sx0; 0533 0534 % 0535 % Denormalization 0536 % 0537 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0538 0539 if isfield(bayes_parm,'DEBUG') && bayes_parm.DEBUG==1 0540 % For DEBUG 0541 fprintf('Old saving mode\n') 0542 else 0543 % Alpha is scaled back by bsnorm 0544 %fprintf('Alpha is scaled back by bsnorm\n') 0545 %Model.a = (Model.a)./bsnorm^2; 0546 %Model.v = (Model.v)./bsnorm^2; 0547 %Model.a0 = (vb_parm.a0)./bsnorm^2; 0548 %Model.v0 = (vb_parm.v0)./bsnorm^2; 0549 Model.Cov = Cov; 0550 Model.ix = ix_act; 0551 Model.bsnorm_meg = bsnorm_meg; 0552 Model.bsnorm_eeg = bsnorm_eeg; 0553 end 0554 0555 % 0556 % Split variance parameters into cortial and extra dipoles 0557 % 0558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0559 0560 if Nextra>0, % Nextra is the number of extra points 0561 Model_ext.sx = Model.sx; 0562 Model_ext.sx0 = Model.sx0; 0563 Model_ext.a = Model.a(end-Nextra*3+1:end,:); 0564 Model_ext.a0 = Model.a0(end-Nextra*3+1:end,:); 0565 Model.a = Model.a(1:end-Nextra*3,:); 0566 Model.a0 = Model.a0(1:end-Nextra*3,:); 0567 else 0568 Model_ext = []; 0569 end 0570 0571 0572 % 0573 % Save result 0574 % 0575 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0576 fprintf('----- Save estimation result in %s\n',bayes_parm_abs.bayesfile) 0577 0578 %bayes_parm = bayes_parm_abs; 0579 vb_fsave([bayes_parm_abs.bayesfile], ... 0580 'Model','Model_ext','bayes_parm','vb_parm',... 0581 'Info','Pointlist'); 0582 0583 % project_file save 0584 proj_file = get_project_filename; 0585 if isempty(proj_file) 0586 return; 0587 end 0588 0589 project_file_mgr('load', proj_file); 0590 project_file_mgr('add', 'bayes_parm', bayes_parm);