clear;
close all;

load hw5em1.mat;
X = data1;
n = size(X,1);
d = size(X,2);

max_loglik = -inf;
P_y = [];
P_xj_given_y = [];
niter=30;
ntrial=6;
L = zeros(niter,ntrial);
for trial = 1:ntrial % run trials, take best model
  fprintf('trial %d\n',trial);
  [P_y_trial,P_xj_given_y_trial,loglik_trial] = discrete_em(niter,X,2);
  L(:,trial)=loglik_trial;
  if (loglik_trial(end)>max_loglik)
    P_y = P_y_trial;
    P_xj_given_y = P_xj_given_y_trial;
    max_loglik = loglik_trial(end);
  end
end

display(max_loglik);
display(P_y);
display(P_xj_given_y{1});
display(P_xj_given_y{2});
display(P_xj_given_y{3});
display(P_xj_given_y{4});
plot(L);
xlabel('EM iterations');
ylabel('avg. log-likelihood');
refresh;
print -deps hw5prob2_2.eps;
