function [P, F, pri, mu, sigma] = GMM(score, X, parameter, verbose)
% GMM.m         The parameters of two Gaussian distributions are estimated by the
%               EM algorithm.
%
% Synopsis:     [P, F, pri, mu, sigma] = GMM(score, X, parameter, verbose);
%
% Input:        score       = the score to be classified
%               X           = the scores of the whole class
%               parameter   = initial parameters for GMM (optional)
%               verbose     = 0 (default, no printing) (optional)
%
% Output:       P           = the probability of belonging to the leading group 
%               F           = the probability of belonging to the secondary group
%               pri         = prior of two groups
%               mu          = mean of two groups
%               sigma       = standard deviation of two groups

if nargin < 4, verbose = 0; end

% Initialization of EM algorithm
maxit       = 1000;         % maximum number of iterations
tol         = 1e-3;         % tolerance on Q - Q_old
Q_old       = -inf;         % the expectation of the complete-data log-likelihood
pri(1)      = parameter(1);
mu(1)       = parameter(2);
sigma(1)    = parameter(3);
pri(2)      = parameter(4);
mu(2)       = parameter(5);
sigma(2)    = parameter(6);

for i = 1:maxit
    if verbose
        fprintf('%3d %12.3f %12.3f %12.3f %12.3f %12.3f\n', i, Q_old, mu(1), sigma(1)^2, mu(2), sigma(2)^2);
    end

    % E-step
    % Compute Q(theta | theta_n)
    p(1,:)      = (sqrt(2*pi)*sigma(1))^(-1)*exp(-1*0.5*(X - mu(1)).^2/(sigma(1)^2));
    p(2,:)      = (sqrt(2*pi)*sigma(2))^(-1)*exp(-1*0.5*(X - mu(2)).^2/(sigma(2)^2));
    tmp         = p(1,:)*pri(1) + p(2,:)*pri(2);
    h(1,:)      = p(1,:)*pri(1)./tmp;     % the prob. that x(n) is belonging to 1st group
    h(2,:)      = p(2,:)*pri(2)./tmp;     % the prob. that x(n) is belonging to 2nd group
    Q           = sum(h(1,:).*log(p(1,:)*pri(1)) +  h(2,:).*log(p(2,:)*pri(2)));
    
    if abs(Q - Q_old) > tol
        Q_old = Q;
    else
        break;
    end
    
    % M-step
    % Find the parameters that maximize Q(theta | theta_n)
    mu(1)       = sum( h(1,:).*X )/sum( h(1,:) );
    mu(2)       = sum( h(2,:).*X )/sum( h(2,:) );
    sigma(1)    = sqrt(sum( h(1,:).*(X - mu(1)).^2)/sum( h(1,:) ));
    sigma(2)    = sqrt(sum( h(2,:).*(X - mu(2)).^2)/sum( h(2,:) ));
    pri(1)      = mean(h(1,:));
    pri(2)      = mean(h(2,:));
end

% Classify the input score
p1    = (sqrt(2*pi)*sigma(1))^(-1)*exp(-1*0.5*(score - mu(1)).^2/(sigma(1)^2));
p2    = (sqrt(2*pi)*sigma(2))^(-1)*exp(-1*0.5*(score - mu(2)).^2/(sigma(2)^2));
tmp     = p1*pri(1) + p2*pri(2);
P       = p1*pri(1)/tmp;
F       = p2*pri(2)/tmp;