function [Pi,Mu1,Sigma1,Mu2, L] = EMGauExp(X,ltol,maxiter,pflag,Init, Verbose)
% 
% EM algorithm for the mixture of Gaussian and Exponential estimation
%
% Inputs:
%   X - input data
%   ltol - percentage of the log likelihood difference between 2 iterations ([] for none)
%   maxiter - maximum number of iteration allowed ([] for none)
%   pflag - 1 for plotting GM for 1D or 2D cases only, 0 otherwise ([] for none)
%   Init - structure of initial Pi, Mu1, Sigma1, Mu2: Init.Pi, Init.Mu1, Init.Sigma1, Init.Mu2 ([] for none)
%
% Ouputs:
%   Pi - estimated weights of mixture components
%   Mu1 - estimated the mean of Gaussians
%   Sigma1 - estimated the standard deviation the Gaussians
%   Mu2 - estimated the mean for Exponential distribution
%   L - log likelihood of estimates
%
% April 2009, Keshi Dai
% 
% Credits
% EM_GM - Patrick P. C. Tsui
% http://www.mathworks.com/matlabcentral/fileexchange/8636-emgm
%

if nargin < 6
    Verbose = false;
end

if isempty(ltol)
    ltol = 0.01;
end

if isempty(maxiter)
    maxiter = 500;
end

%%%% Initialize Pi, Mu1, Sigma1, Mu2 %%%%
t = cputime;
if isempty(Init),  
    [Pi, Mu1, Sigma1, Mu2] = Init_EM(X);     
else
    Pi = Init.Pi;
    Mu1 = Init.Mu1;
    Sigma1 = Init.Sigma1;
    Mu2 = Init.Mu2;
end
%[newL, E] = Expectation(X,Pi,Mu1,Sigma1,Mu2); % Initialize log likelihood
newL = -Inf;
oldL = 2*newL;

%%%% EM algorithm %%%%
niter = 0;
while niter<=maxiter
    [newL, E] = Expectation(X,Pi,Mu1,Sigma1,Mu2); % E-step    
    [Pi,Mu1,Sigma1,Mu2] = Maximization(X,E);  % M-step
    if abs(newL-oldL)<ltol
        break;
    end
    oldL = newL;
    %newL = Likelihood(X,Pi,Mu1,Sigma1,Mu2);
    niter = niter + 1;
    if Verbose
        disp(sprintf('iterations: %d, log-likelihood: %f',niter-1, newL));
    end
end 
L = newL;

%%%% Plot %%%%
if pflag==1,
    if size(X,2)>2
        disp('Can only plot 1 or 2 dimensional applications!/n');
    else
        Plot_GM(X,Pi,Mu1,Sigma1,Mu2);
    end
    elapsed_time = sprintf('CPU time used for EM: %5.2fs',cputime-t);
    if Verbose
        disp(elapsed_time);
        disp(sprintf('Number of iterations: %d',niter-1));
    end
end
%%%%%%%%%%%%%%%%%%%%%%
%%%% End of EM    %%%%
%%%%%%%%%%%%%%%%%%%%%%

function [L, E] = Expectation(X,Pi,Mu1,Sigma1,Mu2)
E = zeros(size(X,1), 2);
% for n=1:size(X,1)
%     E(n,1) = ComputeGamma(X(n), 1, Pi, Mu1, Sigma1, Mu2);
%     E(n,2) = ComputeGamma(X(n), 2, Pi, Mu1, Sigma1, Mu2);
% end
p1 = Pi(1) * normpdf(X,Mu1,Sigma1);
p2 = Pi(2) * exppdf(X,Mu2);
E(:,1) = p1./(p1+p2);
E(:,2) = p2./(p1+p2);

% Compute Likelihood
L = sum(E(:,1) .* ( repmat(log(Pi(1)), size(X)) - repmat(log(2*pi)/2, size(X)) - ...
    repmat(log(Sigma1), size(X)) - (X-Mu1).^2/2*Sigma1^2 ) + ...
    E(:,2) .* ( repmat(log(Pi(2)), size(X)) - X/Mu2 - repmat(log(Mu2), size(X)) ));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End of Expectation %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Pi,Mu1,Sigma1,Mu2] = Maximization(X,E)
Pi(1) = sum(E(:,1))/size(X,1);
Pi(2) = sum(E(:,2))/size(X,1);

Mu1 = sum(E(:,1).*X)/sum(E(:,1));
Sigma1 = sqrt(sum(E(:,1).*((X-Mu1).^2))/sum(E(:,1)));
Mu2 = sum(E(:,2).*X)/sum(E(:,2));
if Sigma1 < 0.0001
    Sigma1 = rand();
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End of Maximization %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function L = Likelihood(X,Pi,Mu1,Sigma1,Mu2)
% % Compute L based on K. V. Mardia, "Multivariate Analysis", Academic Press, 1979, PP. 96-97
% % to enchance computational speed
% L = 0;
% for n=1:size(X,1)
%     gamma1 = ComputeGamma(X(n),1,Pi, Mu1, Sigma1, Mu2);
%     gamma2 = ComputeGamma(X(n),2,Pi, Mu1, Sigma1, Mu2);
%     L1 = gamma1 * (log(Pi(1)) - log(2*pi)/2 - log(Sigma1) - (X(n)-Mu1)^2/2*Sigma1^2);
%     L2 = gamma2 * (log(Pi(2)) - X(n)/Mu2 - log(Mu2));
%     L = L + L1 + L2;
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End of Likelihood %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function gamma = ComputeGamma(xn, i, Pi, Mu1, Sigma1, Mu2)
% gamma = 0;
% p1 = Pi(1)*normpdf(xn,Mu1,Sigma1);
% p2 = Pi(2)*exppdf(xn,Mu2);
% if i==1
%     gamma = p1/(p1+p2);
% end
% if i==2
%     gamma = p2/(p1+p2);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End of ComputeGamma %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Pi, Mu1, Sigma1, Mu2] = Init_EM(X)
% [Ci,C] = kmeans(X,2,'Start','cluster', ...
%     'Maxiter',100, ...
%     'EmptyAction','drop', ...
%     'Display','off'); % Ci(nx1) - cluster indeices; C(k,d) - cluster centroid (i.e. mean)
% while sum(isnan(C))>0,
%     [Ci,C] = kmeans(X,2,'Start','cluster', ...
%         'Maxiter',100, ...
%         'EmptyAction','drop', ...
%         'Display','off');
% end
% Pi(1) = length(find(Ci==1))/size(X,1);
% Pi(2) = length(find(Ci==2))/size(X,1);
% Mu1 = C(1);
% Mu2 = C(2);
% g1 = X(Ci==1);
[p,m1,s1,m2,s2] = EMGauGau(X, [], 500, 0, [],false);
if m1 < m2
    Mu1 = m1;
    Sigma1 = s1;
    Mu2 = m2;
else
    Mu1 = m2;
    Sigma1 = s2;
    Mu2 = m1;
end
if Sigma1 < 0.0001
    Sigma1 = rand();
end
Pi = p;
%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End of Init_EM %%%%
%%%%%%%%%%%%%%%%%%%%%%%%

function Plot_GM(X, Pi, Mu1, Sigma1, Mu2) %W,M,V
[f,r] = ecdf(X);
ecdfhist(f,r)
hold on

Q = zeros(size(r));

P = Pi(1)*normpdf(r,Mu1,Sigma1);
Q = Q + P;
plot(r,P,'r--', 'linewidth',1.5); 
hold on

P = Pi(2)*exppdf(r,Mu2);
Q = Q + P;
plot(r,P,'r--', 'linewidth',1.5); 
hold on

plot(r,Q,'g-','linewidth',2, 'EraseMode', 'none');
grid on
xlabel('X');
ylabel('Probability density');

title('Mixture Models estimated by EM');
%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End of Plot_GM %%%%
%%%%%%%%%%%%%%%%%%%%%%%%
