%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating predictive posterior distributions
% from the 2 component mixture model, using different
% data sets. See KPT, exercise 15.5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

%load c:\klaus\AAEC6564\mlab\worksp\2CMMnoX_data;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lognormal sample density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load posterior results 
load c:\klaus\AAEC6564\mlab\worksp\mod10_2CMMln beta1mat beta2mat ...
    sig11mat sig22mat cpmat Zmat;
R=size(beta1mat,2); %number of draws from the GS

% obtain density evaluations for true sample density
ygrid = linspace(2,20,R)';
% R density evalautions for equally spaced points on [2,20]
ytrue = lognpdf(ygrid,log(10),sqrt(.04));

%generate PPD for predicted y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
int=randn(R,1); %so we use same set of random draws from the standard normal
% for the mixture below
% use same set for all three plots
Z=mnrnd(1,cpmat',R); %R by 2 matrix of rows with 1 zero and 1 one
yp1=beta1mat'+sqrt(sig11mat').*int;
yp2=beta2mat'+sqrt(sig22mat').*int;
yp=sum([yp1 yp2].*Z,2); %picks the y assigned by the Z matrix for each obs.
 
 %evaluate density for predicted draws
[f1,x1]=ksdensity(yp,'kernel','epanechnikov','npoints',1000);

figure(1);
subplot(3,2,1)
plot(ygrid,ytrue,'-k',x1,f1,':k');
title('PPD for 2CMM based on lognormal data');
legend('actual','predicted');
set(gca,'Xlim',[2 20]);
set(gca,'Ylim',[0 0.25]);

% Alternative approach: Predict density at each point of ygrid, averaging
% over all posterior draws (This is what KPT do generate the plots on p.
% 264)
 
 for j = 1:length(ygrid);
ydp(j,1) = mean (cpmat(1,:)'.*normpdf(ygrid(j,1),beta1mat',sqrt(sig11mat')) + ...
   cpmat(2,:)'.*normpdf(ygrid(j,1),beta2mat',sqrt(sig22mat')));
end;

subplot(3,2,2)
plot(ygrid,ytrue,'-k',ygrid,ydp,':k');
title('posterior estimates for density at different points of y, based on logn data');
legend('actual','predicted');
set(gca,'Xlim',[2 20]);
set(gca,'Ylim',[0 0.25]);

clear beta1mat beta2mat sig11mat sig22mat cpmat Zmat;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Chi2 sample density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load posterior results 
load c:\klaus\AAEC6564\mlab\worksp\mod10_2CMMchi2 beta1mat beta2mat ...
    sig11mat sig22mat cpmat Zmat;
R=size(beta1mat,2); %number of draws from the GS

% obtain density evaluations for true sample density
ygrid = linspace(0,30,R)';
% R density evalautions for equally spaced points on [0,30]
ytrue = chi2pdf(ygrid,10);

%generate PPD for predicted y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z=mnrnd(1,cpmat',R); %R by 2 matrix of rows with 1 zero and 1 one
yp1=beta1mat'+sqrt(sig11mat').*int;
yp2=beta2mat'+sqrt(sig22mat').*int;
yp=sum([yp1 yp2].*Z,2); %picks the y assigned by the Z matrix for each obs.
 
 %evaluate density for predicted draws
[f1,x1]=ksdensity(yp,'kernel','epanechnikov','npoints',1000);

subplot(3,2,3)
plot(ygrid,ytrue,'-k',x1,f1,':k');
title('PPD for 2CMM based on chi2 data');
legend('actual','predicted');
set(gca,'Xlim',[0 30]);
set(gca,'Ylim',[0 0.15]);

% Alternative approach: Predict density at each point of ygrid, averaging
% over all posterior draws (This is what KPT do generate the plots on p.
% 264)
 
 for j = 1:length(ygrid);
ydp(j,1) = mean (cpmat(1,:)'.*normpdf(ygrid(j,1),beta1mat',sqrt(sig11mat')) + ...
   cpmat(2,:)'.*normpdf(ygrid(j,1),beta2mat',sqrt(sig22mat')));
end;

subplot(3,2,4)
plot(ygrid,ytrue,'-k',ygrid,ydp,':k');
title('posterior estimates for density at different points of y, based on chi2 data');
legend('actual','predicted');
set(gca,'Xlim',[0 30]);
set(gca,'Ylim',[0 0.15]);

    
clear beta1mat beta2mat sig11mat sig22mat cpmat Zmat;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mixture sample density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load posterior results 
load c:\klaus\AAEC6564\mlab\worksp\mod10_2CMMmix beta1mat beta2mat ...
    sig11mat sig22mat cpmat Zmat;
R=size(beta1mat,2); %number of draws from the GS

% obtain density evaluations for true sample density
ygrid = linspace(-5,5,R)';
% R density evalautions for equally spaced points on [-5,5]
ytrue = 0.4*normpdf(ygrid,0,sqrt(1))+0.6*normpdf(ygrid,2,sqrt(0.5));

%generate PPD for predicted y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z=mnrnd(1,cpmat',R); %R by 2 matrix of rows with 1 zero and 1 one
yp1=beta1mat'+sqrt(sig11mat').*int;
yp2=beta2mat'+sqrt(sig22mat').*int;
yp=sum([yp1 yp2].*Z,2); %picks the y assigned by the Z matrix for each obs.
     
%evaluate density for predicted draws
[f1,x1]=ksdensity(yp,'kernel','epanechnikov','npoints',1000);

subplot(3,2,5)
plot(ygrid,ytrue,'-k',x1,f1,':k');
title('PPD for 2CMM based on mixture data');
legend('actual','predicted');
set(gca,'Xlim',[-5 5]);
set(gca,'Ylim',[0 0.5]);
    
% Alternative approach: Predict density at each point of ygrid, averaging
% over all posterior draws (This is what KPT do generate the plots on p.
% 264)
 
 for j = 1:length(ygrid);
ydp(j,1) = mean (cpmat(1,:)'.*normpdf(ygrid(j,1),beta1mat',sqrt(sig11mat')) + ...
   cpmat(2,:)'.*normpdf(ygrid(j,1),beta2mat',sqrt(sig22mat')));
end;

subplot(3,2,6)
%plot(ygrid,ytrue,'-k',x1,f1,':k');
plot(ygrid,ytrue,'-k',ygrid,ydp,':k');
title('posterior estimates for density at different points of y, based on mixture data');
legend('actual','predicted');
set(gca,'Xlim',[-5 5]);
set(gca,'Ylim',[0 0.5]);
