
clear;
close all;
load hw5em2.mat;
X = data2;
n = size(X,1);
d = size(X,2);

for m=2:5
  fprintf('m = %d\n',m);
  max_ll = -inf;
  for k=1:6
    [param,hist,ll] = em_mix(X,m);
    if (ll(end)>max_ll)
      max_ll = ll(end);
      eval(sprintf('print -depsc hw5prob2_4_fig%d.eps',m));
    end
  end
  model_complexity = m*(d + d*(d+1)/2) + (m-1);
  bic = max_ll - 0.5 * log(n) * model_complexity
end
