1;


## Calculate the MLE of the class probabilities for K classes. The data
## in Y are stored one case per row.  The class densities are Gaussian
## with mean vectors in the columns of Mu, and common covariance matrix
## Cov.
function Pi = EM_known_Gaussian_densities(Y, K, Mu, Cov)

  ## Uniform starting values.
  Pi = ones(K,1)/K; 

  ## The sample size.
  n = size(Y,1);

  ## The number of observations.
  p = size(Y,2);

  CovI = inv(Cov);

  ## EM iterations
  while (1)

    ## E step
    for k=1:K
      R = Y - ones(n,1)*Mu(:,k)';
      for i=1:n
	P(i,k) = R(i,:) * CovI * R(i,:)';
      endfor
    endfor
    P = exp(-P/2) .* (ones(n,1) * Pi');
    P = P ./ (sum(P')' * ones(1,K));

    ## M step
    Pi1 = sum(P)'/n;

    ## Test for 'convergence'.
    if (max(abs(Pi1-Pi)) < 1e-4)
      break;
    endif
    Pi = Pi1;

  end

endfunction


## Calculate the MLE of the class probabilities for K classes. The data
## in Y are stored one case per row.  The class densities are Exponential
## with mean vectors in the columns of Lambda.
function Pi = EM_known_Exponential_densities(Y, K, Lambda)

  ## Uniform starting values.
  Pi = ones(K,1)/K; 

  ## The sample size.
  n = size(Y,1);

  ## The number of observations.
  p = size(Y,2);

  LR = 1./Lambda;
  PLR = prod(LR);

  ## EM iterations
  while (1)

    ## E step
    P = exp(-Y*LR) .* (ones(n,1)*Pi') .* (ones(n,1) * PLR);
    P = P ./ (sum(P')' * ones(1,K));

    ## M step
    Pi1 = sum(P)'/n;

    ## Test for 'convergence'.
    if (max(abs(Pi1-Pi)) < 1e-4)
      break;
    endif
    Pi = Pi1;

  end

endfunction


## Simulate data of smaple size n from a K-component p-dimensional
## Gaussian mixture.  Calculate the MLE of the class probabilities.
## Return the true and estimated class probabilities.
function [Pi, Pi_est] = test_EM_known_Gaussian_densities(p, K, n)

  ## Generate the population structure.
  Mu = randn(p,K);
  H = randn(p,p);
  Cov = H * H';
  Pi = rand(K,1);
  Pi = Pi / sum(Pi);
  CPi = cumsum(Pi);

  ## Generate the data.
  for i=1:n
    
    u = rand(1,1);
    kappa(i) = 1 + sum(CPi < u);
    
    Y(i,:) = Mu(:,kappa(i))' + randn(1,p) * H';
    
  endfor

  Pi_est = EM_known_Gaussian_densities(Y, K, Mu, Cov);

endfunction


## Simulate data of smaple size n from a K-component p-dimensional
## Exponential mixture.  Calculate the MLE of the class probabilities.
## Return the true and estimated class probabilities.
function [Pi, Pi_est] = test_EM_known_Exponential_densities(p, K, n)

  ## Generate the population structure.
  Lambda = rand(p,K);
  Pi = rand(K,1);
  Pi = Pi / sum(Pi);
  CPi = cumsum(Pi);

  ## Generate the data.
  for i=1:n
    
    u = rand(1,1);
    kappa(i) = 1 + sum(CPi < u);
    
    Y(i,:) = -Lambda(:,kappa(i))' .* log(rand(1,p));
    
  endfor

  Pi_est = EM_known_Exponential_densities(Y, K, Lambda);

endfunction

## Run some replicates and print the output to the screen.
for r=1:10
  [Pi, Pi_est] = test_EM_known_Exponential_densities(10, 2, 200);
  disp([Pi, Pi_est]);
  disp("\n");
endfor