%
% Repeat Example 8.12. epage 316
% 
% compute the ISE between the models ... before an application of the EM algorithm 
% 
% use EM on both models to refine the estimates ... 
% 
% recompute the ISE between the models 
% 
% The problem is on epage 327
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');

% generate data from a KNOWN mixture model: 
% 
pi_true = [0.3, 0.3, 0.4];
n = 100; 
x = zeros(n,1); 
r = rand(1,n); 
ind1 = length(find(r<=0.3)); 
ind2 = length(find(r>0.3 & r<=0.6));
ind3 = length(find(r>0.6));
x(1:ind1)           = randn(ind1,1) - 3; 
x((ind1+1):(ind2+ind1)) = randn(ind2,1); 
x((ind1+ind2+1):n)    = randn(ind3,1)*sqrt(0.5)+2; 

% estimate a mixture density to represent this data:
% 
maxterms = 25; 
[pihat_m1,muhat_m1,varhat_m1] = csadpmix(x,maxterms);

% plot this density estimate (model #1): 
%
xinterp = linspace( min(x)-1.0, max(x)+1.0 ); 
fhat_1 = zeros(1,length(xinterp)); 
for ii=1:length(pihat_m1),
  fhat_1 = fhat_1 + pihat_m1(ii)*normpdf(xinterp,muhat_m1(ii),sqrt(varhat_m1(ii))); 
end
fh=figure; plot( xinterp, fhat_1, '-og' ); hold on; grid on; 

% reorder the data and estimates a second density:
%
ind = randperm(n); 
x_2 = x(ind);
[pihat_m2,muhat_m2,varhat_m2] = csadpmix(x_2,maxterms);

% plot this density estimate (model #2): 
%
xinterp = linspace( min(x_2)-1.0, max(x_2)+1.0 ); 
fhat_2 = zeros(1,length(xinterp)); 
for ii=1:length(pihat_m1),
  fhat_2 = fhat_2 + pihat_m1(ii)*normpdf(xinterp,muhat_m1(ii),sqrt(varhat_m1(ii))); 
end
figure(fh); plot( xinterp, fhat_2, '-xr' ); hold on; grid on; 
saveas( gcf, 'prob_8_17_initial_den_estimate', 'epsc' ); 

% compute the error between these two models
%
err=trapz(xinterp,(fhat_1-fhat_2).^2);
fprintf('the ISE between the two density estimates is %10.5g\n',err); 
err=trapz(xinterp,abs(fhat_1-fhat_2));
fprintf('the IAE between the two density estimates is %10.5g\n',err); 

% refine both models using the EM algorithm: 
%
max_its = 100000; tol = 1e-9; 
[pi_m1_em,mu_m1_em,var_m1_em] = csfinmix(x,muhat_m1,varhat_m1,pihat_m1,max_its,tol); 
[pi_m2_em,mu_m2_em,var_m2_em] = csfinmix(x_2,muhat_m2,varhat_m2,pihat_m2,max_its,tol); 
%[pi_m2_em,mu_m2_em,var_m2_em] = csfinmix(x,muhat_m2,varhat_m2,pihat_m2,max_its,tol); 

clear fhat_1 fhat_2; 

% plot this NEW density estimate (model #1): 
%
xinterp = linspace( min(x)-1.0, max(x)+1.0 ); 
fhat_1 = zeros(1,length(xinterp)); 
for ii=1:length(pi_m1_em),
  fhat_1 = fhat_1 + pi_m1_em(ii)*normpdf(xinterp,mu_m1_em(ii),sqrt(var_m1_em(ii))); 
end
fh_em=figure; plot( xinterp, fhat_1, '-og' ); hold on; grid on; 

% plot this NEW density estimate (model #2): 
%
xinterp = linspace( min(x_2)-1.0, max(x_2)+1.0 ); 
fhat_2 = zeros(1,length(xinterp)); 
for ii=1:length(pi_m2_em),
  fhat_2 = fhat_2 + pi_m2_em(ii)*normpdf(xinterp,mu_m2_em(ii),sqrt(var_m2_em(ii))); 
end
figure(fh_em); plot( xinterp, fhat_2, '-xr' ); hold on; grid on; 
saveas( gcf, 'prob_8_17_em_refined_den_estimate', 'epsc' ); 

% compute the error between these two models
%
err=trapz(xinterp,(fhat_1-fhat_2).^2);
fprintf('the ISE between the two density estimates is %10.5g\n',err); 
err=trapz(xinterp,abs(fhat_1-fhat_2));
fprintf('the IAE between the two density estimates is %10.5g\n',err);