%Generate data using GMM 
C1=[0.9 0.1;0.1 0.9]; 
[E1,D1]=eig(C1);
M1=[2 2];
C2=[0.7 0.3;0.3 0.7]; 
[E2,D2]=eig(C2);
M2=[-2 -2];
%The prior probability of the first Mixture and the second Mixtures are 
%0.2 and 0.8 respectively. 
%Generation of 10000 data with additive noise from constructed Model is as 
%given below 
figure(1) 
DATA=[]; 
INDEX=[]; 
for l=1:1:2000 
    r=rand;
    if(r<0.2) 
        temp=randn(2,1);
        temp1=E1'*D1.^(1/2)*temp+M1;
        hold on 
        plot(temp1(1),temp1(2),'ro');
        %second model 
        DATA=[DATA temp1]; 
        INDEX=[INDEX 1]; 
    else
        %first model 
        temp=randn(2,1); 
        temp2=E2'*D2.^(1/2)*temp+M2; 
        hold on 
        plot(temp2(1),temp2(2),'bo');
        DATA=[DATA temp2];
        INDEX=[INDEX 2];
    end
end
%Given data, the mean and the co-variance matrix of the %two classes are obtained as follows.
%Initialize mean
[P,Q]=kmeans(DATA',2);
M1=Q(1,:)';
M2=Q(2,:)';
[u,v]=find(P==1);
temp=DATA(:,u); 
C1=cov(temp'); 
C1=C1/2; 
p1=length(u)/length(P);
[u,v]=find(P==2);
temp=DATA(:,u);
C2=cov(temp');
p2=length(u)/length(P);
MEAN1COL=M1; 
MEAN2COL=M2; 
COV1COL{1}=C1;
COV2COL{1}=C2;
PRIOR1COL=p1;
PRIOR2COL=p2; 
GP=[];
for iteration=1:1:10 
%Expectation stage to estimate the responsibility factor. 
r1=[];
r2=[]; 
for i=1:1:length(DATA) 
    r1=[r1 p1*normalpdf(M1,DATA(:,i),C1)/(p1*normalpdf (M1,DATA(:,i),C1)+p2*normalpdf(M2,DATA(:,i),C2))];
    r2=[r2 p2*normalpdf(M2,DATA(:,i),C2)/(p1*normalpdf (M1,DATA(:,i),C1)+p2*normalpdf(M2,DATA(:,i),C2))];
end
%Maximization stage 
s1=0; s2=0; s3=0; s4=0; 
for i=1:1:length(DATA)
    s1=s1+DATA(:,i)*r1(i); 
    s2=s2+DATA(:,i)*r2(i);
end
M1=s1/sum(r1); 
M2=s2/sum(r2);
for j=1:1:length(DATA)
    s3=s3+(DATA(:,j)-M1')*(DATA(:,j)-M1')'*r1(j);
    s4=s4+(DATA(:,j)-M2')*(DATA(:,j)-M2')'*r2(j);
end
C1=s3/(2*sum(r1));
C2=s4/(2*sum(r2));
p1=sum(r1)/length(DATA);
p2=sum(r2)/length(DATA);
% T1=sum((DATA-repmat(M1,1,length(DATA))).^2); 
% T2=sum((DATA-repmat(M2,1,length(DATA))).^2); 
% [V,I]=min([T1;T2]); 
% [u,v]=find(I==1); 
% p1=length(v)/length(DATA); 
% [u,v]=find(I==2); 
% p2=length(v)/length(DATA); 
MEAN1COL=[MEAN1COL M1]; 
MEAN2COL=[MEAN2COL M2]; 
COV1COL{iteration}=C1; 
COV2COL{iteration}=C2; 
PRIOR1COL=[PRIOR1COL p1]; 
PRIOR2COL=[PRIOR2COL p2]; 
temp=0; 
for i=1:1:length(DATA) 
    temp=temp+log((p1*normalpdf(M1,DATA(:,i),C1)+p2*normalpdf(M2,DATA(:,i),C2)));
end
GP=[GP temp];
end
figure
plot(GP)
xlabel('Iteration'); 
ylabel('Generating probability'); 
figure 
%Given DATA, INDEX estimation is obtained as follows.
l=1;
INDEXEST=[];
for l=1:1:length(DATA) 
    p1=(1/2*pi)^(2/2)*(1/(det(C1)^(1/2)))*exp(-(1/2)* ((DATA(:,l)-M1')'*(inv(C1))*(DATA(:,l)-M1')));
    p2=(1/2*pi)^(2/2)*(1/(det(C2)^(1/2)))*exp(-(1/2)* ((DATA(:,l)-M2')'*(inv(C2))*(DATA(:,l)-M2')));
    if(p1>p2) INDEXEST=[INDEXEST 1]; 
        hold on 
        plot(DATA(1,l),DATA(2,l),'r*') 
    else
        INDEXEST=[INDEXEST 2];
        hold on 
        plot(DATA(1,l),DATA(2,l),'b*') 
    end
end
figure 
stem(INDEX(1:1:100)) 
hold on 
stem(INDEXEST(1:1:100),'r')
