function o = tutorial2()

    pts = [0,20,32];

% Initialization 

sigma_cur = 8;
mean1_cur = 10;
mean2_cur = 32;

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 = [mean1_cur, mean2_cur];