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