Page 311 - Computational Statistics Handbook with MATLAB

300 Computational Statistics Handbook with MATLAB

[n,d] = size(data); % n=# pts, d=# dims
tol = 0.00001; % set up criterion for stopping EM
max_it = 100;
totprob = zeros(n,1);
We also need an initial guess at the component density parameters.
% Get the initial parameters for the model to start EM
mu(:,1) = [-1 -1]'; % each column represents a mean
mu(:,2) = [1 1]';
mix_cof = [0.3 0.7];
var_mat(:,:,1) = eye(d);
var_mat(:,:,2) = eye(d);
varup = zeros(size(var_mat));
muup = zeros(size(mu));
% Just to get started.
num_it = 1;
deltol = tol+1;% to get started
The following steps implement the EM update formulas found in
Equations 8.34 through 8.38.

while num_it <= max_it & deltol > tol
% get the posterior probabilities
totprob = zeros(n,1);
for i=1:c
posterior(:,i) = mix_cof(i)*...
csevalnorm(data,mu(:,i)',var_mat(:,:,i));
totprob = totprob+posterior(:,i);
end
den = totprob*ones(1,c);
posterior = posterior./den;
% Update the mixing coefficients.
mix_cofup = sum(posterior)/n;
% Update the means.
mut = data'*posterior;
MIX = ones(d,1)*mix_cof;
muup = mut./(MIX*n);
% Update the means and the variances.
for i=1:c
cen_data = data-ones(n,1)*mu(:,i)';
mat = cen_data'*...
diag(posterior(:,i))*cen_data;
varup(:,:,i)=mat./(mix_cof(i)*n);
end
% Get the tolerances.
delvar = max(max(max(abs(varup-var_mat))));
delmu = max(max(abs(muup-mu)));

© 2002 by Chapman & Hall/CRC