%% Compute the Multivariate Normal pdf
%

% Copyright 2015 The MathWorks, Inc.


%%
% Compute and plot the pdf of a multivariate normal distribution.
mu = [0 0];
Sigma = [1 0; 0 1 ]
%Sigma = [2 0; 0 0.5 ]
% Sigma = [1 1; 1 5 ]
%   Sigma = [.25 0; 0 1];
%    Sigma = [4 .3; .3 1];
x1 = -4:.2:5; x2 = -4:.2:5;
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x2),length(x1));
figure(11); clf;
surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
axis([min(x1) max(x1) min(x2) max(x2) 0 1])
xlabel('x1'); ylabel('x2'); zlabel('Probability Density');

pause

Sigma = [0.5 0; 0 4 ]
%Sigma = [2 0; 0 0.5 ]
% Sigma = [1 1; 1 5 ]
%   Sigma = [.25 0; 0 1];
%    Sigma = [4 .3; .3 1];
x1 = -4:.2:5; x2 = -4:.2:5;
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x2),length(x1));
figure(12); clf;
surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
axis([min(x1) max(x1) min(x2) max(x2) 0 1])
xlabel('x1'); ylabel('x2'); zlabel('Probability Density');


pause

Sigma = [0.5 0.7; 0.7 4 ]
%Sigma = [2 0; 0 0.5 ]
% Sigma = [1 1; 1 5 ]
%   Sigma = [.25 0; 0 1];
%    Sigma = [4 .3; .3 1];
x1 = -4:.2:5; x2 = -4:.2:5;
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x2),length(x1));
figure(13); clf;
surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
axis([min(x1) max(x1) min(x2) max(x2) 0 1])
xlabel('x1'); ylabel('x2'); zlabel('Probability Density');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pause;

% first distribution is centered at (0,0), second at (-1,3)
mu = [0 0; 3 3];
% covariance of first is identity matrix, second diagonal
%  sigma = cat(3, eye(2), [5 0; 0 1]);
% sigma = cat(3, eye(2), [4 0; 0 0.9]);
  sigma = cat(3, eye(2), [3 0.8; 0.8 1.5])
% again I'm using equal priors
pim = [0.7; 0.3];


% build GMM
gmm = gmdistribution(mu, sigma, pim);
figure(21);
% 2D projection
ezcontourf(@(x,y) pdf(gmm,[x y]));

% view PDF surface
figure(22);
ezsurfc(@(x,y) pdf(gmm,[x y]));
pause();



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 mu = [0 0; 5 -1; -1 4];
 sigma = cat(3, eye(2), [5 0.25; 0.25 1], [1 0; 0 1]);
 pim = [0.5; 0.1; 0.4];


 %% HW problem
 mu = [3 3; 7 4; 5 7];
 sigma = cat(3, [1 0; 0 3], [1 0.5; 0.5 1], [1 0.2; 0.2 1]);

 pim = [0.2; 0.3; 0.5];

figure(31);clf;

% build GMM
gmm2 = gmdistribution(mu, sigma, pim);

% 2D projection
ezcontourf(@(x,y) pdf(gmm2,[x y]), [0 10 1 10]);

% view PDF surface
figure(32);clf;
ezsurfc(@(x,y) pdf(gmm2,[x y]), [0 10 1 10]);
pause();




%pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%N=10000;
%X = random(gmm2,N);
%X = [mvnrnd(mu(1,:), sigma(:,:,1), 2000) ; mvnrnd(mu(2,:), sigma(:,:,2), 3000) ; mvnrnd(mu(3,:), sigma(:,:,3), 5000) ]; 
X=load ('HW3/3gaussian.txt');
[N,D] = size(X)

whos
mu_1 = [0 1];
mu_2 = [ 2 -1]; 
mu_3 = [-2 3];
w = [0.33; 0.33; 0.33];
K=3
D=2
sigma_1 = eye(2); 
sigma_2 = eye(2);
sigma_3 = eye(2);
Pi = rand(N,K);
for r=1:1000
	r
	%% E-step
	F1 = mvnpdf(X, mu_1, sigma_1);	
	F2 = mvnpdf(X, mu_2, sigma_2);	
	F3 = mvnpdf(X, mu_3, sigma_3);
	Pi(:,1) = F1*w(1);
	Pi(:,2) = F2*w(2);
	Pi(:,3) = F3*w(3);
	Pi = Pi ./ sum(Pi,2);
	
  	
	
	
	%% M-step	
	n_ = sum(Pi,1);
	w = n_ / N	
	
	mu_1 = sum(X .* Pi(:,1),1) / n_(1)
	mu_2 = sum(X .* Pi(:,2),1)  / n_(2)
	mu_3 = sum(X .* Pi(:,3),1)  / n_(3)
	
	
	sigma_1 = zeros(D);
	M = bsxfun(@minus, X, mu_1);
	for i=1:N
		sigma_1  = sigma_1 + M(i,:)' * M(i,:) * Pi(i,1);
	end;


	sigma_2 = zeros(D);
	M = bsxfun(@minus, X, mu_2);
	for i=1:N
		sigma_2  = sigma_2 + M(i,:)' * M(i,:) * Pi(i,2);
	end;

	sigma_3 = zeros(D);
	M = bsxfun(@minus, X, mu_3);
	for i=1:N
		sigma_3  = sigma_3 + M(i,:)' * M(i,:) * Pi(i,3);
	end;
	
	sigma_1 = sigma_1 / n_(1)
	sigma_2 = sigma_2 / n_(2)
	sigma_3 = sigma_3 / n_(3)
	%whos;	
	pause	
end;
		
		

