% This code generates the data points by Gaussian mixture model.

% 20 Gaussians generate 5 pieces of data, 

% 8 for the ring; 1 for the center of ring; 5 for the half ring; 5 for the other half ring; 1 extra.  



clear;

N=1000; % total number of data points



% ring data

centers(1:9,:) = [1,0; 1/sqrt(2),1/sqrt(2); 0,1; -1/sqrt(2),1/sqrt(2); ...

                 -1,0; -1/sqrt(2),-1/sqrt(2); 0,-1; 1/sqrt(2),-1/sqrt(2); 0,0];

    

covars(:,:,1) = [0.01,0;0,0.05];  covars(:,:,2) = [0.03,-0.02;-0.02,0.03]; 

covars(:,:,3) = [0.05,0;0,0.01];  covars(:,:,4) = [0.03,0.02;0.02,0.03]; 

covars(:,:,5) = [0.01,0;0,0.05];  covars(:,:,6) = [0.03,-0.02;-0.02,0.03]; 

covars(:,:,7) = [0.05,0;0,0.01];  covars(:,:,8) = [0.03,0.02;0.02,0.03]; 

covars(:,:,9) = [0.02,0;0,0.02];  

      

% two half ring data

centers(10:19,:) = [1,0; 1/sqrt(2),1/sqrt(2); 0,1; -1/sqrt(2),1/sqrt(2); -1,0; ...

                   -1,0; -1/sqrt(2),-1/sqrt(2); 0,-1; 1/sqrt(2),-1/sqrt(2); 1,0];

centers(10:14,:)=centers(10:14,:)+ones(5,1)*[5,5];

centers(15:19,:)=centers(15:19,:)+ones(5,1)*[4,4.8];



covars(:,:,10) = [0.01,0;0,0.05];  covars(:,:,11) = [0.03,-0.02;-0.02,0.03]; 

covars(:,:,12) = [0.05,0;0,0.01];  covars(:,:,13) = [0.03,0.02;0.02,0.03]; 

covars(:,:,14) = [0.01,0;0,0.05];  covars(:,:,15) = [0.01,0;0,0.05];

covars(:,:,16) = [0.03,-0.02;-0.02,0.03]; covars(:,:,17) = [0.05,0;0,0.01];

covars(:,:,18) = [0.03,0.02;0.02,0.03]; covars(:,:,19) = [0.01,0;0,0.05];    



%Gaussian data  

centers(20,:) = [6,3];

covars(:,:,20) = [0.02,0;0,0.06];



% prior

p=[0.05*ones(1,8), 0.1, 0.04*ones(1,10), 0.1];



% cumulation of priors

cump=cumsum(p);



% generate the data by given mixture models

for i=1:N

   s=rand;

   supk=find(cump-s>0);

   k=supk(1);

   x(i,:) = mvnrnd(centers(k,:),covars(:,:,k),1);

end



% likelihood

for k=1:20

   cx = x-ones(N,1)*centers(k,:); 

   likeli(:,k) = 1/(2*pi*sqrt(det(covars(:,:,k))))*exp(-0.5*sum(cx*inv(covars(:,:,k)).*cx,2));

end



% posterior   

pp=likeli.*(ones(N,1)*p)./(sum(likeli.*(ones(N,1)*p),2)*ones(1,20));


[dummy,idx0]=max(pp,[],1);   groups=[1,9,10,15,20]; idx=idx0(groups);  % Added, AL: 2/28/04
plot(x(:,1),x(:,2),'b.',x(idx,1),x(idx,2),'rx')

save mix_data pp x idx

break,

