function [u,b,t,iter] = fit_mix_laplace( X,M )
%
% fit_mix_laplace - fit parameters for a mixed-laplacian distribution using EM algorithm
%
% format:   [u,b,t,iter] = fit_mix_laplacian( X,M )
%
% input:    X   - input samples, Nx1 vector
%           M   - number of gaussians which are assumed to compose the distribution
%
% output:   u   - fitted mean for each laplacian
%           b - fitted standard deviation for each laplacian
%           t   - probability of each laplacian in the complete distribution
%           iter- number of iterations done by the function
%
N           = length( X );
Z           = ones(N,M) * 1/M;                  % indicators vector
P           = zeros(N,M);                       % probabilities vector for each sample and each model
t           = ones(1,M) * 1/M;                  % distribution of the gaussian models in the samples
u           = linspace(0.2,1.4,M);        % mean vector
b           = ones(1,M) * var(X) / sqrt(M);     % variance vector
C           = 1/sqrt(2*pi);                     % just a constant
Ic          = ones(N,1);                        % - enable a row replication by the * operator
Ir          = ones(1,M);                        % - enable a column replication by the * operator
Q           = zeros(N,M);                       % user variable to determine when we have converged to a steady solution
thresh      = 1e-7;
step        = N;
last_step   = 300;         % step/last_step
iter        = 0;
min_iter    = 3000;
while ((( abs((step/last_step)-1) > thresh) & (step>(N*1e-10)) ) & (iter<min_iter) )% E step% ========Q   = Z;P   = 1./ (Ic*b) .* exp( -(1e-6+abs(X*Ir - Ic*u))./(Ic*b) );for m = 1:MZ(:,m)  = (P(:,m)*t(m))./(P*t(:));end% estimate convergence step size and update iteration numberprog_text   = sprintf(repmat( '\b',1,(iter>0)*12+ceil(log10(iter+1)) ));iter        = iter + 1;last_step   = step * (1 + eps) + eps;step        = sum(sum(abs(Q-Z)));fprintf( '%s%d iterations\n',prog_text,iter );% M step% ========Zm              = sum(Z);               % sum each columnZm(find(Zm==0)) = eps;                  % avoid devision by zerou               = sum(((X*Ir)./abs(X*Ir - Ic*u)).*Z) ./sum(1./abs(X*Ir - Ic*u).*Z) ;b               = sum((abs(X*Ir - Ic*u)).*Z) ./ Zm ;t               = Zm/N;
end
end

clc;clear all;
%generate random
xmin = -10;
xmax = 10;
Len = 10000000;
x = linspace(xmin,xmax,Len);
mu = [3,-4];
b = [0.9 0.4];
w = [0.7 0.3];
fx = w(1)/2/b(1)*exp(-abs(x-mu(1))/b(1))+ w(2)/2/b(2)*exp(-abs(x-mu(2))/b(2));
ymax = 1/b(2);
ymin = 0;
Y = (ymax-ymin)*rand(1,Len)-ymin;
data = x(Y<=fx);
%Laplace Mixture Model fitting
K = 2;
[mu_new,b_new,w_new,iter] = fit_mix_laplace( data',K);
%figure
subplot 221
hist(data,2000);
grid on;
subplot 222
numter = [xmin:.2:xmax];
plot(numter,w_new(1)/2/b_new(1)*exp(-abs(numter-mu_new(1))/b_new(1)),'r','linewidth',2);hold on;
plot(numter,w_new(2)/2/b_new(2)*exp(-abs(numter-mu_new(2))/b_new(2)),'g','linewidth',2);hold on;subplot (2,2,[3,4])
[histFreq, histXout] = hist(data, numter);
binWidth = histXout(2)-histXout(1);
%Bar
bar(histXout, histFreq/binWidth/sum(histFreq)); hold on;grid on;
plot(numter,w_new(1)/2/b_new(1)*exp(-abs(numter-mu_new(1))/b_new(1)),'r','linewidth',2);hold on;
plot(numter,w_new(2)/2/b_new(2)*exp(-abs(numter-mu_new(2))/b_new(2)),'g','linewidth',2);hold on;