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

vbmeg_multi_fmri

PURPOSE ^

Bayesian Wiener estimation for MEG.

SYNOPSIS ^

function [Model,Info] =vbmeg_multi_fmri(B, Gall, Gact, COV, vb_parm)

DESCRIPTION ^

 Bayesian Wiener estimation for MEG. 
 Estimate current variance by using time sequence of all trials. 
 Mixture of Multiple fMRI activity patterns are used for prior
 Current variance is represented by mixture of fMRI patterns
 Weight coefficient for each fMRI pattern is estimated

 -- Syntax
 function [Model,Info] = vbmeg_multi_fmri(B, Gall, Gact, COV, vb_parm),

 -- 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

 (Field of vb_parm)
 .Ntrain
 .Nskip
 .a_min
 .a_max
 .Npre_train (optional)
 .Ntrials
 .Nsession
 .Nwindow
 .Twindow
 .Ta0
 .Tv0
 .sx0
 .a0
 .v0
 .wiener
 .act
 .cont_pr
 
 -- Output
 Model.a   : Current variance (1/alpha)  
 Model.v   : Current variance for background activity
 Model.sx  : Observation noise variance

 Model.act   : Multiple fMRI activity patterns for focal area
 Model.act_v : Their estimated weights

 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. 

 --- History
 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)

 2005/12/5 M.Sato
 * New version
 * Prior is given for each time window
 2006/8/11 M.Sato
  Bayesian Wiener estimation
  Mixture of Multiple fMRI activity patterns are used for prior
 2008-11-28 Taku Yoshioka
   Use vb_disp() for displaying message

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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