% demo of GMM clustering using MATLAB statistic toolbox 

function demoGMM_statisticTB


close all;
clc;
clear;
%% Generate 2D data using function mvnrnd refering to Multivariate normal random numbers
mu1 = [1 2];
sigma1 = [3 .2; .2 2];
mu2 = [-1 -2];
sigma2 = [2 0; 0 1];

X = [mvnrnd(mu1,sigma1,200);mvnrnd(mu2,sigma2,100)];

scatter(X(:,1),X(:,2),10,'ko');

%% Training GMM 
options = statset('Display','final');%,'MaxIter',10);
% options = statset('Display','iter');
gm = gmdistribution.fit(X,2,'Options',options,'Replicates',1, 'CovType','full');

%AIC http://en.wikipedia.org/wiki/Akaike_information_criterion
%BIC http://en.wikipedia.org/wiki/Bayesian_information_criterion

%%
hold on;
colourmp = hsv(gm.NComponents);  
for i = 1 : gm.NComponents  
plot(gm.mu(i,1),gm.mu(i,2),'x','lineWidth', 10,'color', colourmp(i,:)')
end

%% Plot GMM components
hold on
ezcontour(@(x,y)pdf(gm,[x y]),[-8 6],[-8 6]);
for i = 1 : gm.NComponents  
plot(gm.mu(i,1),gm.mu(i,2),'x','lineWidth', 10,'color', colourmp(i,:)')
end
hold off

%% Partition the Data
idx = cluster(gm,X);
cluster1 = (idx == 1);
cluster2 = (idx == 2);

scatter(X(cluster1,1),X(cluster1,2),10,'r+');
hold on
scatter(X(cluster2,1),X(cluster2,2),10,'bo');
hold off
legend('Cluster 1','Cluster 2','Location','NW')

%% Compute posterior for data points

P = posterior(gm,X);

scatter(X(cluster1,1),X(cluster1,2),10,P(cluster1,1),'+')
hold on
scatter(X(cluster2,1),X(cluster2,2),10,P(cluster2,1),'o')
hold off
legend('Cluster 1','Cluster 2','Location','NW')
clrmap = jet(80); colormap(clrmap(9:72,:))
ylabel(colorbar,'Component 1 Posterior Probability')

%% Test on new data
Y = [mvnrnd(mu1,sigma1,50);mvnrnd(mu2,sigma2,25)];

scatter(Y(:,1),Y(:,2),10,'ko');

idx = cluster(gm,Y);

cluster1 = (idx == 1);
cluster2 = (idx == 2);

scatter(Y(cluster1,1),Y(cluster1,2),10,'ro');
hold on;
scatter(Y(cluster2,1),Y(cluster2,2),10,'bo');

ezcontour(@(x,y)pdf(gm,[x y]),[-8 6],[-8 6]);
hold on
for i = 1 : gm.NComponents  
plot(gm.mu(i,1),gm.mu(i,2),'x','lineWidth', 10,'color', colourmp(i,:)')
end

hold off
legend('Cluster 1','Cluster 2','Location','NW')






