Variational Bayes estimation for MEG.
Estimate current variance by using time sequence of all trials.
Prune algorithm.
-- Syntax
function [Model,Info] = vbmeg_ard_estimate2(B, Gall, Gact, COV, vb_parm),
function [Model,Info] = vbmeg_ard_estimate2(B, Gall, Gact, COV, vb_parm, Modelinit)
-- Input
B{Nsession} : MEG magnetic field data matrix
size(B{i}) = [Nsensors(i), Tsample, Ntrials(i)]
Gall{Nsession}: Lead field for whole brain region as background activity
size(Gall{i}) = [Nsensors(i) Ndipole]
If Gall = [], the sparse-mode (focal dipoles only) is
executed.
Gact{Nsession} : Lead field (Current basis function)
for region with high estimation accuracy
size(Gact{i}) = [Nsensors(i) Ndipole]
If Gact = [], the iso-mode (no focal dipoles) is
executed.
COV{Nsession} : Noise covariance matrix
size(COV{i}) = [Nsensors(i) Nsensors(i)]
vb_parm : Constants parameter
-- Optional input
Modelinit : Initial Model parameters, if given.
if not given, the variables vb_parm.a0 and vb_parm.v0 are used.
The size of each variable in 'Modelinit' must be the same as that of
the output 'Model'.
-- Output
Model.a : Current variance (1/alpha)
Model.v : Current variance for background activity
Model.sx : Observation noise variance
Info : Containing convergence information such as 'FE','Ev','Err'
-- Note
Brain noise : Background brain activity is assumed
as independent Gaussian noise with same variance
"Nvact" (# of vertex for active dipoles) and
"Nt" (length of time window) are common for all sessions.
"Nsensors" (# of sensor channels in a session) and
"Ntrial" (# of trials in a session) could depend on each session.
Note that ax0_j, ax0_z, and sx are defined as variance:
the inverse of inverse-variance (i.e having the diminsion of variance).
These notations are different from those in Sato-san's note.
2004-12-5 M. Sato
2005-04-12 O.Yamashita
* The whole program are revised
so that variables have consistent names with Sato-san's note.
* The noise variance, sx, ax_z and ax_j are independently handled.
* A field 'fix_parm' is added to 'vb_parm' struct.
* The input 'COV' must be given.
2005-04-16 M. Sato
* Modified free-energy calculation
2005-05-07 M. Sato
* Modified variance update (Faster convergence)
2005-05-13 O.Yamashita
* 'Modelinit' is added as the 6th input optional argument.
* Iso-mode is implemented. (when Gact is [])
* A field 'fix_parm' is changed to two flag paramters, 'update_v' and
'update_sx'.
* lower bound of 'gain_z'
2005-05-16 O.Yamashita
* Sparse-mode is implementd. (when Gall is [])
* lower bound of 'gain_j'
2005-05-31 O.Yamashita
* constant term of free energy
2005-06-02 O.Yamashita
* constat term of free energy revised
2005/08/22 O.Yamashita --- ver.30b
* vb_parm is introduced
* Output argument is modified
* Minor bug fix
2005/09/29 O. Yamashita
* Minor bug fix (when Norinet=2 & iso-mode)
2006/04/16 M. Sato
* vb_parm.Npre_train is introduced
VB-update equation is used for iteration <= Npre_train
2006/08/20 M. Sato
* Soft normal constraint
------------ Relative variance calculation ----------------
In this implementation,
Current variance 'ax_z' and 'ax_j' are defined as ratio w.r.t sx :
ax_j * sx = <(J - Z)*(J - Z)>
ax_z * sx = < Z*Z>
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [Model,Info] = ... 0002 vbmeg_ard_estimate2(B, Gall, Gact, COV, vb_parm, Modelinit) 0003 % Variational Bayes estimation for MEG. 0004 % Estimate current variance by using time sequence of all trials. 0005 % Prune algorithm. 0006 % 0007 % -- Syntax 0008 % function [Model,Info] = vbmeg_ard_estimate2(B, Gall, Gact, COV, vb_parm), 0009 % function [Model,Info] = vbmeg_ard_estimate2(B, Gall, Gact, COV, vb_parm, Modelinit) 0010 % 0011 % -- Input 0012 % B{Nsession} : MEG magnetic field data matrix 0013 % size(B{i}) = [Nsensors(i), Tsample, Ntrials(i)] 0014 % Gall{Nsession}: Lead field for whole brain region as background activity 0015 % size(Gall{i}) = [Nsensors(i) Ndipole] 0016 % If Gall = [], the sparse-mode (focal dipoles only) is 0017 % executed. 0018 % Gact{Nsession} : Lead field (Current basis function) 0019 % for region with high estimation accuracy 0020 % size(Gact{i}) = [Nsensors(i) Ndipole] 0021 % If Gact = [], the iso-mode (no focal dipoles) is 0022 % executed. 0023 % COV{Nsession} : Noise covariance matrix 0024 % size(COV{i}) = [Nsensors(i) Nsensors(i)] 0025 % vb_parm : Constants parameter 0026 % 0027 % -- Optional input 0028 % Modelinit : Initial Model parameters, if given. 0029 % if not given, the variables vb_parm.a0 and vb_parm.v0 are used. 0030 % The size of each variable in 'Modelinit' must be the same as that of 0031 % the output 'Model'. 0032 % 0033 % -- Output 0034 % Model.a : Current variance (1/alpha) 0035 % Model.v : Current variance for background activity 0036 % Model.sx : Observation noise variance 0037 % Info : Containing convergence information such as 'FE','Ev','Err' 0038 % 0039 % -- Note 0040 % Brain noise : Background brain activity is assumed 0041 % as independent Gaussian noise with same variance 0042 % 0043 % "Nvact" (# of vertex for active dipoles) and 0044 % "Nt" (length of time window) are common for all sessions. 0045 % 0046 % "Nsensors" (# of sensor channels in a session) and 0047 % "Ntrial" (# of trials in a session) could depend on each session. 0048 % 0049 % Note that ax0_j, ax0_z, and sx are defined as variance: 0050 % the inverse of inverse-variance (i.e having the diminsion of variance). 0051 % These notations are different from those in Sato-san's note. 0052 % 0053 % 2004-12-5 M. Sato 0054 % 2005-04-12 O.Yamashita 0055 % * The whole program are revised 0056 % so that variables have consistent names with Sato-san's note. 0057 % * The noise variance, sx, ax_z and ax_j are independently handled. 0058 % * A field 'fix_parm' is added to 'vb_parm' struct. 0059 % * The input 'COV' must be given. 0060 % 2005-04-16 M. Sato 0061 % * Modified free-energy calculation 0062 % 2005-05-07 M. Sato 0063 % * Modified variance update (Faster convergence) 0064 % 2005-05-13 O.Yamashita 0065 % * 'Modelinit' is added as the 6th input optional argument. 0066 % * Iso-mode is implemented. (when Gact is []) 0067 % * A field 'fix_parm' is changed to two flag paramters, 'update_v' and 0068 % 'update_sx'. 0069 % * lower bound of 'gain_z' 0070 % 2005-05-16 O.Yamashita 0071 % * Sparse-mode is implementd. (when Gall is []) 0072 % * lower bound of 'gain_j' 0073 % 2005-05-31 O.Yamashita 0074 % * constant term of free energy 0075 % 2005-06-02 O.Yamashita 0076 % * constat term of free energy revised 0077 % 2005/08/22 O.Yamashita --- ver.30b 0078 % * vb_parm is introduced 0079 % * Output argument is modified 0080 % * Minor bug fix 0081 % 2005/09/29 O. Yamashita 0082 % * Minor bug fix (when Norinet=2 & iso-mode) 0083 % 2006/04/16 M. Sato 0084 % * vb_parm.Npre_train is introduced 0085 % VB-update equation is used for iteration <= Npre_train 0086 % 2006/08/20 M. Sato 0087 % * Soft normal constraint 0088 % 0089 % 0090 % ------------ Relative variance calculation ---------------- 0091 % In this implementation, 0092 % Current variance 'ax_z' and 'ax_j' are defined as ratio w.r.t sx : 0093 % ax_j * sx = <(J - Z)*(J - Z)> 0094 % ax_z * sx = < Z*Z> 0095 % 0096 % Copyright (C) 2011, ATR All Rights Reserved. 0097 % License : New BSD License(see VBMEG_LICENSE.txt) 0098 0099 Ntrain = vb_parm.Ntrain; 0100 Nskip = vb_parm.Nskip; 0101 a_min = vb_parm.a_min; 0102 a_max = vb_parm.a_max; 0103 0104 if isfield(vb_parm, 'Npre_train') 0105 % Number of iteration using original VB-update rule 0106 % for stable estimation 0107 Npre_train = vb_parm.Npre_train; 0108 else 0109 Npre_train = 0; 0110 end 0111 0112 fprintf('--- Initial VB-update iteration = %d\n',Npre_train) 0113 fprintf('--- Total update iteration = %d\n',Ntrain) 0114 0115 % Number parameters 0116 Ntrials = vb_parm.Ntrials; % Number of trials 0117 Nsensors = vb_parm.Nsensors; % Number of sensors 0118 Nsession = vb_parm.Nsession; % Number of sessions 0119 Nwindow = vb_parm.Nwindow; % Number of time windows 0120 Twindow = vb_parm.Twindow; % Time window index [Tstart, Tend] 0121 Ntotaltrials = sum(Ntrials); % Total number of trials in all sessions 0122 Nerror = sum(Ntrials.*Nsensors); % Total MEG channel used for error estimate 0123 0124 Nvact = vb_parm.Nvact; % Number of active current vertices 0125 Njact = vb_parm.Njact; % Number of active dipoles parameters 0126 Njall = vb_parm.Njall; % Number of whole brain vertices (dipole) 0127 % for background activity 0128 Norient = vb_parm.Norient; % Number of current orientation component 0129 Norient_var = vb_parm.Norient_var;% Number of current orientation component 0130 % for variance estimation 0131 Nvaract = Nvact * Norient_var;% Number of variance parameters 0132 % (=length(ax_z)) 0133 Ratio = Norient/Norient_var; % = (# of orientation) 0134 % /(# of orientation in variance) 0135 Normal_mode = vb_parm.soft_mode; % Soft normal constraint mode 0136 0137 % Set prior parameters 0138 Ta0 = vb_parm.Ta0 ; % confidence parameter of focal 0139 Tv0 = vb_parm.Tv0 ; % confidence parameter of global 0140 sx0 = vb_parm.sx0; % observation noise variance 0141 a0 = max(vb_parm.a0, a_min ); % initial variance of focal 0142 v0 = max(vb_parm.v0, a_min ); % initial variance of global 0143 0144 if Normal_mode == 1 & ~isempty(Gact), 0145 % Soft normal constraint 0146 % orientation variance ratio 0147 ndev = reshape(a0, [Nvact Norient_var]); 0148 ndev = ndev./repmat( sum(ndev,2) , [1 Norient_var]); 0149 ndev = ndev(:); 0150 end 0151 0152 if exist('Modelinit','var') & ~isempty(Modelinit) 0153 a0s = max(Modelinit.a, a_min ); 0154 v0s = max(Modelinit.v, a_min ); 0155 else 0156 a0s = repmat(a0, [1 Nwindow]); 0157 v0s = repmat(v0, [1 Nwindow]); 0158 end 0159 0160 % Flag cotroling updates of variance parameters 0161 update_sx = 1; 0162 update_v = 1; 0163 update_z = 1; 0164 if isfield(vb_parm, 'update_sx') 0165 update_sx = vb_parm.update_sx; 0166 end 0167 if isfield(vb_parm, 'update_v') 0168 update_v = vb_parm.update_v; 0169 end 0170 0171 %%% iso mode 0172 if isempty(Gact) 0173 Nvaract = 1; %%% avoid empty (adhoc?) 0174 for ns = 1 : Nsession 0175 Nch = Nsensors(ns); % number of channel available in each session 0176 Gact{ns} = zeros(Nch, Nvaract); 0177 end 0178 a0s = zeros(Nvaract,Nwindow); 0179 update_v = 1; 0180 update_z = 0; 0181 end 0182 0183 %%% sparse mode 0184 if isempty(Gall) 0185 for ns = 1 : Nsession 0186 Gall{ns} = 0; 0187 end 0188 v0s = zeros(1,Nwindow); 0189 update_v = 0; 0190 update_z = 1; 0191 Njall = 0; %%%% need to modify ? 0192 end 0193 0194 0195 %%% calculate BB' and GG' 0196 for ns = 1 : Nsession 0197 for nw = 1 : Nwindow 0198 B0 = 0; 0199 for i = 1 : Ntrials(ns) 0200 Bt = B{ns}(:,Twindow(nw,1):Twindow(nw,2),i); 0201 B0 = B0 + Bt * Bt'; 0202 end 0203 BBall{ns,nw} = B0; 0204 end 0205 GGall{ns} = Gall{ns}*Gall{ns}'; 0206 end 0207 0208 clear B Gall 0209 0210 % Temporal variable for variance (ratio) 0211 aall = zeros(Nvaract,Nwindow); % Modified by TY 2005-06-15 0212 vall = zeros(Nwindow,1); 0213 sxall = zeros(Nwindow,1); 0214 0215 % Free energy and Error history 0216 Info = zeros(Nwindow*Ntrain, 10); 0217 k_check = 0; 0218 0219 %%%%% Time window loop %%%%%%%% 0220 for nw=1:Nwindow 0221 0222 Nt = Twindow(nw,2) - Twindow(nw,1) + 1; 0223 Ntt = Nt*Ntotaltrials; 0224 0225 % confidence parameters 0226 gx0_z = Ta0; 0227 gx0_j = Tv0; 0228 0229 gx_z = 0.5 * Ntt * Ratio + gx0_z; 0230 gx_j = 0.5 * Ntt * Njall + gx0_j; 0231 0232 % expected variance parameters 0233 ax0_z = a0s(:,nw) ; 0234 ax0_j = v0s(nw) ; 0235 0236 % variance hyper-parameter initialization 0237 if vb_parm.cont_pr == OFF | nw == 1 0238 % use the prior variance as the initial expected value 0239 ax_j = ax0_j; 0240 ax_z = ax0_z; 0241 sx = sx0; 0242 else 0243 % use the estimated variance in the previous time window as the 0244 % initial expected value 0245 ix_zero = find(ax_z < a_min); % indices taking small variance 0246 ax_z(ix_zero) = ax0_z(ix_zero); 0247 end 0248 0249 FEconst = free_energy_constant(gx_j,gx0_j,gx_z,gx0_z,... 0250 Nt,Ntotaltrials,Nerror,Njall); 0251 0252 ix_z = [1:Nvaract]; % initial effective index 0253 0254 %%%%%% Estimation Loop %%%%%% 0255 for k=1:Ntrain 0256 0257 %%% VB E-step -- current estimation 0258 ax_z_tmp = ax_z(ix_z); % variance of effective index 0259 Nvaract_tmp = length(ix_z); % # of effective index 0260 0261 % Initialize averaging variable 0262 tr1_b = 0; % (B - G*J)^2 0263 tr2_b = 0; % <(G*J)^2 > 0264 mean_jj = 0; % (J - Z)*(J - Z) 0265 gain_j = 0; % Estimation Gain for (J - Z) 0266 mean_zz = zeros(Nvaract_tmp,1); % diag( Z*Z ) 0267 gain_z = zeros(Nvaract_tmp,1); % Estimation Gain for Z 0268 Hj = 0; % log(det(SB)) 0269 0270 for ns=1:Nsession 0271 0272 Nch = Nsensors(ns); 0273 Ntry = Ntrials(ns); 0274 0275 % Lead field for each session 0276 GG = GGall{ns}; % Gb * Gb' 0277 0278 if Ratio == 1 0279 % Each current orientation component has its own variance 0280 G = Gact{ns}(:,ix_z); 0281 else 0282 % Each current orientation component has common variance 0283 % at the same vertex 0284 G = []; 0285 for j = 1:Norient 0286 G = [G,Gact{ns}(:,ix_z+Nvact*(j-1))]; % Ga 0287 end 0288 end 0289 Gt = G'; % Ga' 0290 0291 % MEG covariance for each session 0292 BB = BBall{ns,nw}; % (B B') 0293 0294 % Noise covariance for each session 0295 Cov = COV{ns}; % Sigma_G^{-1} 0296 0297 A_z = repmat(ax_z_tmp , [Ratio Nch]); 0298 0299 % Inverse filter for each session 0300 GSG = G*(A_z.*Gt) + ax_j*GG; % Ga * Sigma_alpha *Ga' 0301 SB = GSG + Cov; % Sigma_B 0302 SB_inv = inv( SB ); % Sigam_B^{-1} 0303 GSB = Gt*SB_inv; % Ga' * Sigma_B^{-1} 0304 SBBS = SB_inv*BB*SB_inv; % Sigma_B^{-1}*BB*Sigma_B^{-1} 0305 0306 % Reconstraction error = sum( (B-G*J)'*Cov*(B-G*J) ) 0307 tr1_b = tr1_b + sum(sum(Cov.*SBBS)); 0308 % Error variance = G' * < J'*Cov*J > * G 0309 % tr2_b = tr2_b + Nt*Ntry*sum(sum(GSG.*SB_inv)); 0310 0311 % Background activity = Tr( sum( (J - Z) * (J - Z) ) ) 0312 mean_jj = mean_jj + ax_j*sum(sum(GG.*SBBS)); 0313 % Estimation gain for (J-Z) = Tr(ax_j*GG*SB_inv) 0314 gain_j = gain_j + Nt*Ntry*ax_j*sum(sum(GG.*SB_inv)); 0315 0316 % Current magnitude = diag( sum( Z * Z ) ) 0317 mean_zz = mean_zz + ax_z_tmp ... 0318 .* sum(reshape(sum((GSB*BB).*GSB,2),[Nvaract_tmp Ratio]),2); 0319 % Estimation gain for Z = diag(ax_z*Gt*SB_inv*G) 0320 gain_z = gain_z + Nt*Ntry*ax_z_tmp ... 0321 .* sum( reshape(sum(Gt.*GSB, 2),[Nvaract_tmp Ratio]), 2); 0322 0323 %if ( mod(k, Nskip)==0 || k==Ntrain ) 0324 % Model complexity for Current parameter 0325 Hj = Hj - 0.5*Nt*Ntry*vb_log_det(SB); 0326 %end 0327 end; 0328 0329 % Noise variance estimation 0330 sx_total = (tr1_b + mean_jj + sum(mean_zz)); 0331 0332 % update the observation noise variance 0333 if update_sx 0334 sx = sx_total/(Nt*Nerror); 0335 end 0336 0337 %%%%% VB M-step -- variance update 0338 0339 % Save old value for Free energy calculation 0340 ax_j_old = ax_j; 0341 ax_z_tmp_old = ax_z_tmp; 0342 ix_z_old = ix_z; 0343 % 0344 % --- Background activity variance estimation 0345 % 0346 % VB update equation 0347 % gx_j * ax_j = ... 0348 % (0.5*Ntt*Njall + gx0_j) * ax_j = ... 0349 % mean_jj_all + 0.5*ax_j*( Ntt * Njall - gain_j ) ; 0350 0351 mean_jj_all = 0.5*ax_j*mean_jj/sx + gx0_j.*ax0_j; 0352 % 0353 % --- Active current variance estimation 0354 % 0355 % VB update equation 0356 % gx_z * ax_z = ... 0357 % (0.5*Ntt*Ratio + gx0_z) * ax_z = ... 0358 % mean_zz_all + 0.5*ax_z.*( Ntt * Ratio - gain_z ); 0359 0360 mean_zz_all = 0.5*ax_z_tmp.*mean_zz/sx + gx0_z(ix_z).*ax0_z(ix_z); 0361 0362 % Update the current variance parameters 0363 if k <= Npre_train 0364 % --- VB update rule 0365 0366 % Background activity variance estimation 0367 if update_v 0368 ax_j_total = mean_jj_all + 0.5*ax_j.*(Ntt*Njall - gain_j) ; 0369 ax_j = ax_j_total ./ gx_j; 0370 end 0371 % Active current variance estimation 0372 if update_z 0373 ax_z_total = mean_zz_all + 0.5*ax_z_tmp.*(Ntt*Ratio - gain_z); 0374 ax_z_tmp = ax_z_total ./ gx_z(ix_z); 0375 end 0376 else 0377 % --- Faster convergence update rule 0378 0379 % Background activity variance estimation 0380 if update_v 0381 % elements less than eps are replaced by eps 0382 gain_j = max(gain_j, eps); 0383 ax_j = mean_jj_all ./ (0.5*gain_j + gx0_j); 0384 end 0385 % Active current variance estimation 0386 if update_z 0387 gain_z = max(gain_z, eps); 0388 ax_z_tmp = mean_zz_all ./ (0.5*gain_z + gx0_z(ix_z)); 0389 end 0390 end 0391 0392 if update_z 0393 if Normal_mode == 1, 0394 % Average orientation variance 0395 ax_z_tmp = repmat( sum(reshape(ax_z_tmp./ndev(ix_z),... 0396 [Nvaract_tmp/3 3]) ,2)/3 , [3 1]); 0397 end 0398 0399 % Prune small variance 0400 ax_z(ix_z) = ax_z_tmp; % recover ax_z 0401 ix_zero = find(ax_z < a_min); % indices taking small variance 0402 ix_z = find(ax_z >= a_min); % effective indices 0403 ax_z(ix_zero) = 0; 0404 0405 if Normal_mode == 1, 0406 ax_z(ix_z) = ax_z(ix_z) .* ndev(ix_z); 0407 end 0408 end 0409 0410 %%%%%%%% Free Energy Calculation %%%%%%%% 0411 % Model complexity for current variance 0412 if update_z 0413 rz = ax0_z(ix_z_old)./ax_z_tmp_old; 0414 Haz = sum( gx0_z(ix_z_old).*(log(rz) - rz + 1)); 0415 else 0416 Haz = 0; 0417 end 0418 0419 if gx0_j.*ax_j_old > eps 0420 Haj = gx0_j.*(log(ax0_j./ax_j_old) - ax0_j./ax_j_old + 1); 0421 else 0422 Haj = 0; 0423 end; 0424 0425 LP = - 0.5*Nt*Nerror*log(sx); 0426 0427 % Data Likelihood & Free Energy 0428 Ea = - 0.5*( sx_total/sx - Nt*Nerror ); 0429 FE = LP + Hj + Ea + Haz + Haj + FEconst; 0430 Ev = LP + Hj; % evidence (without constants) 0431 Err = tr1_b/sum(sum(BB));% Normalized Error 0432 0433 % DEBUG Info 0434 k_check = k_check + 1; 0435 Info(k_check,:) = [... 0436 FE, Ev, Hj, LP, Err,... 0437 sx, max(ax_z), ax_j, sum(gain_z), gain_j]; 0438 % end of free energy calculation 0439 0440 if mod(k, Nskip)==0 || k==Ntrain 0441 fprintf('Tn=%3d, Iter=%4d, FE=%f, Error=%e\n', nw, k, FE, Err); 0442 end; 0443 0444 end % end of learning loop 0445 0446 aall(:,nw) = ax_z; 0447 vall(nw) = ax_j; 0448 sxall(nw) = sx; 0449 end % end of time window loop 0450 0451 Model.a = aall ; % Current variance 0452 Model.v = vall ; % Current variance for background activity 0453 Model.sx = sxall; % Sensor noise variance 0454 0455 0456 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0457 %%% functions specific for this function 0458 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0459 function FEconst = free_energy_constant(gx_j, gx0_j, gx_z, gx0_z,... 0460 Nt, Ntotaltrials, Nerror, Njall) 0461 0462 gx_sx = 0.5*Nt*Nerror; 0463 L0 = 0.5*Nt*(Njall*Ntotaltrials*qG(gx_j) + Ntotaltrials*sum(qG(gx_z)) + Nerror*qG(gx_sx) - Nerror); 0464 H0 = HG0(gx_sx,0) + HG0(gx_j,gx0_j)+sum(HG0(gx_z, gx0_z)); 0465 FEconst = L0+H0; 0466 0467 %%%%% 0468 function y = qG(gx) 0469 N=length(gx); 0470 nz = find(gx ~= 0); % non-zero components 0471 y = zeros(N,1); 0472 y_tmp = psi(gx(nz)) - log(gx(nz)); 0473 y(nz) = y_tmp; 0474 0475 %%%%%% 0476 function y =HG0(gx, gx0) 0477 N=length(gx0); 0478 nz = find(gx ~= 0); % non-zero components 0479 nz0 = find(gx0 ~= 0); % non-zero components 0480 % gammma 0481 y1 = zeros(N,1); 0482 y1_tmp = gammaln(gx(nz)) - gx(nz).*psi(gx(nz)) + gx(nz); 0483 y1(nz) = y1_tmp; 0484 % gamma0 0485 y2 = zeros(N,1); 0486 y2_tmp = gammaln(gx0(nz0)) - gx0(nz0).*log(gx0(nz0)) + gx0(nz0); 0487 y2(nz0) = y2_tmp; 0488 % gamma*gamma0 0489 y3 = zeros(N,1); 0490 y3_tmp = gx0(nz0).*(psi(gx(nz0))-log(gx(nz0))); 0491 y3(nz0) = y3_tmp; 0492 0493 y = y1 - y2 + y3; 0494 0495 0496 0497