%%% Expectation Maximization %%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K=9; %%% Components of learning clusters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=360; %%% Number of samples
KURIKAESHI=1000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% make samples %%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
for j=1:1:9
X0(1,j)=floor((j-1)/3);
X0(2,j)=mod(j-1,3);
end
for i=1:1:n
XX(:,i)=0.3*randn(2,1)+X0(:,mod(i-1,9)+1);
end
%}
XX=2*rand(2,n);
%%%%%%%%%%%%%%%%%%%%%%% plot samples %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,4,1);
plot(XX(1,:),XX(2,:),'b.'); hold on
title('Samples'); hold on
xlim([-0.5,2.5]);
ylim([-0.5,2.5]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gauss=@(x1,x2,b1,b2,s)(exp(-((x1-b1).^2+(x2-b2).^2)./(2*s))./(2*pi*s));
%%%%%%%%%%%%%%%%%%%%%%% plot samples %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ttt=1:1:11
%%%%%%%%%% Initialize
aaa=1/K*ones(1,K);
bbb1=mean(XX(1,:))+0.2*randn(1,K);
bbb2=mean(XX(2,:))+0.2*randn(1,K);
sigma=(var(XX(:,1))+var(XX(:,2)))/2*ones(1,K);
%%%%%%%%%% Recursive VB Start
for kuri=1:1:KURIKAESHI
for i=1:1:n
each=aaa.*gauss(XX(1,i),XX(2,i),bbb1,bbb2,sigma);
wa=sum(each);
YYY(:,i)=each/wa;
end
for k=1:1:K
sumy=sum(YYY(k,:));
aaanew(k)=sumy/n;
bbb1new(k)=sum(YYY(k,:).*XX(1,:))/sumy;
bbb2new(k)=sum(YYY(k,:).*XX(2,:))/sumy;
sigmanew(k)=sum(YYY(k,:).*((XX(1,:)-bbb1(k)).^2+(XX(2,:)-bbb2(k)).^2))/(2*sumy);
end
aaa=aaanew;
bbb1=bbb1new;
bbb2=bbb2new;
sigma=sigmanew;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[x1,y1]=meshgrid(-0.5:0.1:2.5,-0.5:0.1:2.5);
zzz=0*x1;
for k=1:1:K
zzz = zzz + aaa(k)*gauss(x1,y1,bbb1(k),bbb2(k),sigma(k));
end
loglike=0;
for i=1:1:n
loglike=loglike+sum(aaa.*gauss(XX(1,i),XX(2,i),bbb1,bbb2,sigma));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('%f\n',loglike);
%%%%%%%%%%%%%%%%%%%% plot true and estimated %%%%%%%%%%%%%
subplot(3,4,ttt+1);
contour(x1,y1,zzz);
xlim([-0.5,2.5]);
ylim([-0.5,2.5]);
end