function [parameters] = EMND(x,mus,alphas,nc,steps,estimate,cvs,dim1,dim2)

%x - observations
%mu N-dimensional means for each class (need nc N-dimensional guesses)
%alpha N-dimensional weights for each class (need nc N-dimensional guesses)
%nc - number of classes. (i.e., how many Gaussians to use to fit x)

%Returns the parameters of the PDFs reached after the algorithm converges
%(i.e., after steps number of steps)


%Initialize the parameters for the probability density function
%with an initial estimate.

if(estimate==0)
  %Use input means for initial estimate
  %Start each class receiving the naive guess of Cov(x) for its covariance matrix
  K = cov(x);
  parameters.mus=mus;
  parameters.alphas=alphas;
  parameters.nc = nc;
  for class=1:nc
     parameters.covariances(:,:,class)=K;
  end
end

if(estimate==1)
  %Use k-means clustering for initial guess
  parameters.mus = estimateClassMeans(x,nc);
  parameters.covariances = estimateClassCov(mus);
  parameters.alphas=alphas;
  parameters.nc = nc;
end


if(estimate==2)
 %Use Data is two-dimensional, so plot the data and estimate manually
 handle = plot(x(:, 1), x(:, 2), 'bo');
 mus = ginput
 K = cov(x);
 parameters.nc = nc;
 parameters.mus=mus;
 parameters.alphas = alphas;
 for class=1:nc
   parameters.covariances(:,:,class)=K;
 end
end

%The dimension to display


[obs,dim] = size(x)
for step=1:steps

%E-Step
%Find the class posterior probabilities given the current parameter estimates

  %generateEllipses2(x,parameters,mus,dim1,dim2);
  postProb=posteriorProbability(x,parameters);
  %disp('Posterior Probabilities Calculated')
  %pause

  %likelihood = marginalProbability(x,parameters,2);
  %likelihood
  %big = max(likelihood);
  %small = min(likelihood);
  %mn = mean(likelihood);
  %vr = var(likelihood);
  %numBig = size(find(likelihood>mn));
  
  %big
  %small
  %mn
  %vr
  %numBig

% M-step
  
  %Calculate new estimates for the PDFs

  %Estimate new priors 
  %Sum over class responsibilities
  involvement = sum(postProb,1);
  parameters.alphas = involvement ./ obs;
  
  %Estimate new means
  parameters.mus = (postProb' * x)./(involvement'*ones(1,dim));

  %Estimate new covariance matrices
  if(cvs==2)
    for class=1:parameters.nc
      diff = x - (ones(obs, 1) * parameters.mus(class,:));
      parameters.covariances(:,:,class) = diag(sum((diff.*diff).*(postProb(:,class)*ones(1,dim)),1)./involvement(class));	
    end
  end

  if(cvs==1)
    for class=1:parameters.nc
    diff = x - (ones(obs, 1) * parameters.mus(class,:));
    diff = diff.*(sqrt(postProb(:,class))*ones(1, dim));
    parameters.covariances(:,:,class) = (diff'*diff)/involvement(class);
    end
  end
end

  %likelihood = marginalProbability(x,parameters,2);
  %likelihood
  %big = max(likelihood);
  %small = min(likelihood);
  %mn = mean(likelihood);
  %vr = var(likelihood);
  %numBig = size(find(likelihood>mn));

  %big
  %small
  %mn
  %vr
  %numBig

  %generateEllipses2(x,parameters,mus,dim1,dim2);

end

generateEllipses2(x,parameters,mus,dim1,dim2);


%Calculate the posterior probabilities
function postProb=posteriorProbability(x,parameters);

[obs,dim] = size(x);

sample = sampleNormalDistribution3(x,parameters);
%sample = sampleNormalDistribution(x,parameters);

%multiply samples from each class by the priors
responsibilities = (ones(obs, 1)*parameters.alphas).*sample;

%Sum responsibilities of a observation from each class
sumResp = sum(responsibilities, 2);

%Divide 
postProb = responsibilities./(sumResp*ones(1, parameters.nc));


function sample = sampleNormalDistribution(x,parameters)
method = 2;
[obs,dim] = size(x);
n = dim/2;
for class = 1:parameters.nc
    errors = x - (ones(obs, 1) * parameters.mus(class,:));
    K = parameters.covariances(:,:,class);
    Kinv = inv(K);
    G = det(K);
    for man = 1:obs
      sample(man,class)=(1/((2*pi)^n)*sqrt(G))*exp(-.5*errors(man,:)*Kinv*errors(man,:)');
    end
end

function sample = sampleNormalDistribution3(x,parameters)
[obs,dim] = size(x);
method=1;
n = dim/2;
normal = (2*pi)^(n);
for i = 1:parameters.nc
  diffs = x - (ones(obs, 1) * parameters.mus(i,:));
  % Use Cholesky decomposition of covariance matrix to speed computation
  c = chol(parameters.covariances(:,:,i));
  temp = diffs/c;
  sample(:,i) = exp(-0.5*sum(temp.*temp, 2))./(normal*prod(diag(c)));
end