1;

% Demonstration of the EM algorithm for fitting two lines to a bunch
%  of points.  In the E step, we try to estimate the probability that
%  each point is actually on each of the lines.  Then in the M step
%  we redraw the lines to maxmimize the likelihood of the points.
%
% Rob Malouf (malouf@let.rug.nl), 25 May 2000
%
%   Alfa Informatica
%   University of Groningen
%   Postbus 716 
%   9737NT Groningen

% Set up data

x = linspace(0,1,10)';
y(1:5) = x(1:5) + 1;
y(6:10) = -x(6:10);

% Initialize random model

model = rand(2,2);
s2 = 0.1;

multiplot(1,3);

for i = 1 : 10,

% E step

r1 = model(1,1) * x + model(1,2) - y;
r2 = model(2,1) * x + model(2,2) - y;

%r = x * model(1,:)  + model(2,:) - y


w1 = exp(-r1.*r1./s2) ./ (exp(-r1.*r1./s2)+exp(-r2.*r2./s2));
w2 = exp(-r2.*r2./s2) ./ (exp(-r1.*r1./s2)+exp(-r2.*r2./s2));

% Report results 

subwindow(1,1);
mplot(x,y,"+",x,x*model(1,1)+model(1,2),x,x*model(2,1)+model(2,2));
[xb,yb] = bar(x,w1);
subwindow(1,2);
mplot(xb,yb);
[xb,yb] = bar(x,w2);
subwindow(1,3);
mplot(xb,yb);

input("Next");

% M step

A = [ sum(w1.*x.*y); sum(w1.*y)];
B = [ sum(w1.*x.*x), sum(w1.*x); sum(w1.*x), sum(w1)];
X = inv(B) * A;
model(1,:)=X';

A = [ sum(w2.*x.*y); sum(w2.*y)];
B = [ sum(w2.*x.*x), sum(w2.*x); sum(w2.*x), sum(w2)];
X = inv(B) * A;
model(2,:) = X';

end