% EM_EXP1 Example on using EM algorithm.
%

% Modifications:
% 25-aug-2004, VF


% ground truth model
% -------------------------------------------------------
true_model.Mean = [[-3;3],[3;3],[3;-3],[-3;-3]];
true_model.Cov(:,:,1) = [1.5 0.9; 0.9 1.3];
true_model.Cov(:,:,2) = [1.4 -0.3; -0.3 0.7];
true_model.Cov(:,:,3) = [0.9 0; 0 1.4];
true_model.Cov(:,:,4) = [0.9 -0.3; -0.3 1.7];
true_model.Prior = [0.2 0.3 0.2 0.3];

% sample training and testing data
% -------------------------------------------------------
trn_data = gmmsamp(true_model, 400);

% estimate GMM with 4 components
options.ncomp = 4;
options.init = 'kmeans';
options.verb = 1;

gmm_model = emgmm(trn_data.X, options);

% evaluate log-likelihood of the training sdata
% -------------------------------------------------------
logL = sum(log( pdfgmm( trn_data.X, gmm_model)))

% classify training data to classes (clusters).
% -------------------------------------------------------
[dummy,trn_data.y] = max(gmm_model.Alpha);


% displayes solution (works only for 2D data)
% -------------------------------------------------------
for i=1:length(gmm_model.Prior),
  bayes_model.Pclass{i}.Mean = gmm_model.Mean(:,i);
  bayes_model.Pclass{i}.Cov = gmm_model.Cov(:,:,i);
  bayes_model.Pclass{i}.fun = 'pdfgauss';
end
bayes_model.Prior = gmm_model.Prior;
bayes_model.fun = 'bayescls';

figure; title('Input unlabeled data'); ppatterns(trn_data.X); 

figure; title('Clustered data and estimated components of GMM'); 
ppatterns(trn_data); pboundary(bayes_model); pgauss(gmm_model)

figure; title('Plot of optimized Log-likelihood');
plot(gmm_model.logL); xlabel('iterations'); ylabel('logL');

figure; title('Estimate distribution function');
pgmm(gmm_model,{'visual','surf'}); grid on;

% EOF
