%
% 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');
load forearm; x=forearm; 

% estimate a mixture density to represent this data:
% 
maxterms = 1; 
[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; ph1=plot( xinterp, fhat_1, '-og' ); hold on; grid on; 

% refine our model 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); 

% 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
figure(fh); ph2=plot( xinterp, fhat_1, '-xr' ); hold on; grid on; 
legend( [ph1,ph2], {'adaptive density estimate','refined EM model'} ); 
saveas( gcf, 'prob_8_21_pdf', 'epsc' );