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

vbmeg_ard_estimate2

PURPOSE ^

Variational Bayes estimation for MEG.

SYNOPSIS ^

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

DESCRIPTION ^

 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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