Group2EMActivity

From Computational Statistics Course Wiki
Jump to navigation Jump to search
x = [5;9;8;4;7];
thetaInit = [0.6;0.4];

% EM method
theta = thetaInit;
while(1)
    thetaPrime = theta;
    % E-step
    p = zeros(5,2);
    for(i = 1:5)
        for(j = 1:2)
            p(i,j) = thetaPrime(j)^x(i)*(1 - thetaPrime(j))^(10 - x(i));
        end
        p(i,:) = p(i,:)/sum(p(i,:));
    end
    
    % M-step
    for(i = 1:2)
        theta(i) = sum(p(:,i).*x)/(10*sum(p(:,i)));
    end
    
    % termination
    if(norm(theta - thetaPrime) < 1e-6)
        break;
    end
end
theta

% naive method
theta = thetaInit;
while(1)
    thetaPrime = theta;
    % E-step
    p = zeros(5,2);
    for(i = 1:5)
        for(j = 1:2)
            p(i,j) = thetaPrime(j)^x(i)*(1 - thetaPrime(j))^(10 - x(i));
        end
        p(i,:) = p(i,:)/sum(p(i,:));
    end
    p = round(p);
    
    % M-step
    for(i = 1:2)
        theta(i) = sum(p(:,i).*x)/(10*sum(p(:,i)));
    end
    
    % termination
    if(norm(theta - thetaPrime) < 1e-6)
        break;
    end
end
theta

Running this results in the following

theta =

    0.7968
    0.5196


theta =

    0.8000
    0.4500