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