function o = tutorial;
sigma = 5;
mean1 = 100;
mean2 = 130;

pts1 = sigma*randn(1,20) + mean1;
pts2 = sigma*randn(1,20) + mean2;
pts = [pts1,pts2];

% Initialization 

sigma_cur = 50;
mean1_cur = 255*rand;
mean2_cur = 255*rand;

for i = 1:10
    
    % Expectation step
    
    prob1 = exp(- ((pts - mean1_cur).^2 / sigma_cur^2));
    prob2 = exp(- ((pts - mean2_cur).^2 / sigma_cur^2));
    weights1 = prob1./(prob1+prob2);
    weights2 = prob2./(prob1+prob2);
    
    plot_em(pts, weights1, mean1_cur, mean2_cur, sigma_cur);
    if i == 10
        'last iteration'
    end
    keyboard;
    
    % Maximization step
    
    mean1_cur = sum(pts.*weights1)/sum(weights1);
    mean2_cur = sum(pts.*weights2)/sum(weights2);
    sigma_cur = sqrt( (sum(((pts - mean1_cur).^2).*weights1) +  ...
        sum(((pts - mean2_cur).^2).*weights2)) /length(pts));
    
end
o = 1;