%
% Script to generate example figures demonstrating the B-spline series
% probability density estimation package.
%

% 
% Copyright 2014-2018 Biometrics Northwest LLC
% 
%

%
% Modified: 2018-09-12
% Author : Kevin R. Gehringer
%
%    Added examples comparing the reflection boundary correction with the
%    histogram boundary correction.
%
% Modified: 2016-11-04
% Author : Kevin R. Gehringer
%
%    Added 1-d and 2-d examples for bounded domains with active boundaries.
%
% Modified: 2016-06-25
% Author : Kevin R. Gehringer
%
%    Added probability functions and other examples.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

%
% Clear the workspace
%
clear

%
% reset the randome number generator
%
rng('default');

%
% Define the output file ID 
%
fid_out = stdout;

%
%-------------------------------------------------------------------------
%

%
% Simple example in the user's guide.
%

% Generate data from a 2-D standard normal distribution

data = randn(1000,2);

% Estimate a pdf from data

bss1 = bsspdfest(data);

% Evaluate the estimated PDF at the data points to get
% their likelihood values

lh = bsseval(data,bss1,'pdf');

% Plot the likelihood values in 3-D to show the
% approximation to the PDF.

figure;
plot3(data(:,1),data(:,2),lh,'+');
axis tight;
grid
xlabel('X1')
ylabel('X2')
zlabel('Likelihood')

%
%-------------------------------------------------------------------------
%

%
% 1-D example figures.
%

%
% Generate data from a standard normal distribution and estimate
% the probability density function
%

nsample   = 1000;
ndata     = nsample;
evalbnds  = [-5.5 5.5];
xx1       = linspace(evalbnds(1),evalbnds(2),501)';
xtest1    = randn(nsample,1);
bss1test  = bsspdfest(xtest1);
y1        = bsseval(xx1,bss1test,'pdf');
yn        = normpdf(xx1,0,1);

figure;
subplot(2,1,1)
plot(xx1,[y1 yn])
axis tight;
legend('Estimate','N(0,1)');
title('Probability density functions');
subplot(2,1,2)
plot(xx1,yn-y1)
axis tight;
title('Error')

fprintf(fid_out,'Test bsspdfest and bsseval using 1D data:\n');
fprintf(fid_out,'   Distribution                 = N(0,1)\n');
fprintf(fid_out,'   Evaluation interval          = (%f,%f)\n',evalbnds);
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx1));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y1)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(yn)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(yn-y1))*(xx1(2)-xx1(1)));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((yn-y1).^2)/length(yn)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((yn-y1).^2)*(xx1(2)-xx1(1))));
fprintf(fid_out,'\n');

y1cdf = bsseval(xx1,bss1test,'cdf');
yncdf    = normcdf(xx1,0,1);
figure;
subplot(2,1,1)
plot(xx1,[y1cdf yncdf])
axis tight;
legend('Estimate','N(0,1)');
title('Cumulative distribution functions');
subplot(2,1,2)
plot(xx1,yncdf-y1cdf)
axis tight;
title('Error')

y1surv = bsseval(xx1,bss1test,'survivor');
ynsurv = 1-normcdf(xx1,0,1);
figure;
subplot(2,1,1)
plot(xx1,[y1surv ynsurv])
axis tight;
legend('Estimate','N(0,1)');
title('Survivor functions');
subplot(2,1,2)
plot(xx1,ynsurv-y1surv)
axis tight;
title('Error')

y1ch = bsseval(xx1,bss1test,'cumhazard');
ynch = 1-normcdf(xx1,0,1);
gt0 = ynch>0;
ynch(~gt0) = nan;
ynch(gt0)  = -log(ynch(gt0));
figure;
subplot(2,1,1)
plot(xx1,[y1ch ynch])
axis tight;
legend('Estimate','N(0,1)','location','northwest');
title('Cumulative hazard functions');
subplot(2,1,2)
plot(xx1,ynch-y1ch)
axis tight;
title('Error')

xx1p = linspace(0,1,101)';
y1icdf = bsseval(xx1p,bss1test,'icdf');
ynicdf = norminv(xx1p,0,1);
figure;
subplot(2,1,1)
plot(xx1p,[y1icdf ynicdf])
axis tight;
legend('Estimate','N(0,1)','location','northwest');
title('Inverse cumulative distribution functions');
subplot(2,1,2)
plot(xx1p,ynicdf-y1icdf)
axis tight;
title('Error')

%
% Generate data from a mixture of 2 1D normal distributions:
%   1/3 * N(0,1) + 2/3 * N(3,0.5)
% and estimate the density function.
%
nsample = 500;
ndata     = 3*nsample;
evalbnds  = [-2 6];
xx1       = linspace(evalbnds(1),evalbnds(2),501)';
xtest1    = [0.25*randn(nsample,1); 3+0.5*randn(2*nsample,1)];
bss1test  = bsspdfest(xtest1);
y1        = bsseval(xx1,bss1test,'pdf');
yn        = (1/3)*normpdf(xx1,0,0.25) + (2/3)*normpdf(xx1,3,0.5);

figure;
subplot(2,1,1)
plot(xx1,[y1 yn])
axis tight;
legend('Estimate','Mixture');
title('Probability density functions');
subplot(2,1,2)
plot(xx1,yn-y1)
axis tight;
title('Error')

fprintf(fid_out,'Test bsspdfest and bsseval using 1D data:\n');
fprintf(fid_out,'   Distribution                 = 1/3 * N(0,0.25) + 2/3 * N(3,0.5)\n');
fprintf(fid_out,'   Evaluation interval          = (%f,%f)\n',evalbnds);
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx1));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y1)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(yn)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(yn-y1))*(xx1(2)-xx1(1)));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((yn-y1).^2)/length(yn)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((yn-y1).^2)*(xx1(2)-xx1(1))));
fprintf(fid_out,'\n');

y1cdf = bsseval(xx1,bss1test,'cdf');
yncdf = (1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5);
figure;
subplot(2,1,1)
plot(xx1,[y1cdf yncdf])
axis tight;
legend('Estimate','Mixture');
title('Cumulative distribution functions');
subplot(2,1,2)
plot(xx1,yncdf-y1cdf)
axis tight;
title('Error')

y1surv = bsseval(xx1,bss1test,'survivor');
ynsurv = 1-((1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5));
figure;
subplot(2,1,1)
plot(xx1,[y1surv ynsurv])
axis tight;
legend('Estimate','Mixture');
title('Survivor functions');
subplot(2,1,2)
plot(xx1,ynsurv-y1surv)
axis tight;
title('Error')

y1ch = bsseval(xx1,bss1test,'cumhazard');
ynch = 1-((1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5));
gt0 = ynch>0;
ynch(~gt0) = nan;
ynch(gt0)  = -log(ynch(gt0));
figure;
subplot(2,1,1)
plot(xx1,[y1ch ynch])
axis tight;
legend('Estimate','Mixture)','location','northwest');
title('Cumulative hazard functions');
subplot(2,1,2)
plot(xx1,ynch-y1ch)
axis tight;
title('Error')

xx1p = linspace(0,1,501)';
y1icdf = bsseval(xx1p,bss1test,'icdf');
xx = linspace(-2,6,1001)';
ynpdf = (1/3)*normpdf(xx,0,0.25) + (2/3)*normpdf(xx,3,0.5);
yncdf = cumtrapz(xx,ynpdf);
ztol = 2*max([min(yncdf) min(abs(1-yncdf))]);
idx0 = find(yncdf<ztol,1,'last');
idx1 = find(abs(1-yncdf)<ztol,1,'first');
idxbad  = find(diff(yncdf)==0);
idxbad  = idxbad((idxbad>=idx0 & idxbad <=idx1));
idxgood = (idx0:idx1)';
idxgood = setdiff(idxgood,idxbad);
ynicdf  = interp1(yncdf(idxgood),xx(idxgood),xx1p,'linear','extrap');
figure;
subplot(2,1,1)
plot(xx1p,[y1icdf ynicdf])
axis tight;
legend('Estimate','Mixture','location','northwest');
title('Inverse cumulative distribution functions');
subplot(2,1,2)
plot(xx1p,ynicdf-y1icdf)
axis tight;
title('Error')

%
% Repeat the mixture example using more data.
%
nsample   = 2500;
ndata     = 3*nsample;
evalbnds  = [-2 6];
xx1       = linspace(evalbnds(1),evalbnds(2),501)';
xtest1    = [0.25*randn(nsample,1); 3+0.5*randn(2*nsample,1)];
bss1test  = bsspdfest(xtest1);
y1        = bsseval(xx1,bss1test,'pdf');
yn        = (1/3)*normpdf(xx1,0,0.25) + (2/3)*normpdf(xx1,3,0.5);

figure;
subplot(2,1,1)
plot(xx1,[y1 yn])
axis tight;
legend('Estimate','Mixture');
title('Probability density functions');
subplot(2,1,2)
plot(xx1,yn-y1)
axis tight;
title('Error')

fprintf(fid_out,'Test bsspdfest and bsseval using 1D data:\n');
fprintf(fid_out,'   Distribution                 = 1/3 * N(0,0.25) + 2/3 * N(3,0.5)\n');
fprintf(fid_out,'   Evaluation interval          = (%f,%f)\n',evalbnds);
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx1));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y1)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(yn)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(yn-y1))*(xx1(2)-xx1(1)));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((yn-y1).^2)/length(yn)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((yn-y1).^2)*(xx1(2)-xx1(1))));
fprintf(fid_out,'\n');

y1cdf = bsseval(xx1,bss1test,'cdf');
yncdf = (1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5);
figure;
subplot(2,1,1)
plot(xx1,[y1cdf yncdf])
axis tight;
legend('Estimate','Mixture');
title('Cumulative distribution functions');
subplot(2,1,2)
plot(xx1,yncdf-y1cdf)
axis tight;
title('Error')

y1surv = bsseval(xx1,bss1test,'survivor');
ynsurv = 1-((1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5));
figure;
subplot(2,1,1)
plot(xx1,[y1surv ynsurv])
axis tight;
legend('Estimate','Mixture');
title('Survivor functions');
subplot(2,1,2)
plot(xx1,ynsurv-y1surv)
axis tight;
title('Error')

y1ch = bsseval(xx1,bss1test,'cumhazard');
ynch = 1-((1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5));
gt0 = ynch>0;
ynch(~gt0) = nan;
ynch(gt0)  = -log(ynch(gt0));
figure;
subplot(2,1,1)
plot(xx1,[y1ch ynch])
axis tight;
legend('Estimate','Mixture)','location','northwest');
title('Cumulative hazard functions');
subplot(2,1,2)
plot(xx1,ynch-y1ch)
axis tight;
title('Error')

xx1p = linspace(0,1,501)';
y1icdf = bsseval(xx1p,bss1test,'icdf');
xx = linspace(-2,6,1001)';
ynpdf = (1/3)*normpdf(xx,0,0.25) + (2/3)*normpdf(xx,3,0.5);
yncdf = cumtrapz(xx,ynpdf);
ztol = 2*max([min(yncdf) min(abs(1-yncdf))]);
idx0 = find(yncdf<ztol,1,'last');
idx1 = find(abs(1-yncdf)<ztol,1,'first');
idxbad  = find(diff(yncdf)==0);
idxbad  = idxbad((idxbad>=idx0 & idxbad <=idx1));
idxgood = (idx0:idx1)';
idxgood = setdiff(idxgood,idxbad);
ynicdf  = interp1(yncdf(idxgood),xx(idxgood),xx1p,'linear','extrap');
figure;
subplot(2,1,1)
plot(xx1p,[y1icdf ynicdf])
axis tight;
legend('Estimate','Mixture','location','northwest');
title('Inverse cumulative distribution functions');
subplot(2,1,2)
plot(xx1p,ynicdf-y1icdf)
axis tight;
title('Error')

%
% Generate data from a mixture of 2 1D normal distributions:
%   1/3 * N(0,1) + 2/3 * N(3,0.5)
% and estimate the density function as a histogram.
%
nsample = 500;
ndata     = 3*nsample;
evalbnds  = [-2 6];
xx1       = linspace(evalbnds(1),evalbnds(2),501)';
xtest1    = [0.25*randn(nsample,1); 3+0.5*randn(2*nsample,1)];
bss1test  = bsspdfest(xtest1,[],[],1);
y1        = bsseval(xx1,bss1test,'pdf');
yn        = (1/3)*normpdf(xx1,0,0.25) + (2/3)*normpdf(xx1,3,0.5);

figure;
subplot(2,1,1)
plot(xx1,[y1 yn])
axis tight;
legend('Estimate','Mixture');
title('Probability density functions');
subplot(2,1,2)
plot(xx1,yn-y1)
axis tight;
title('Error')

fprintf(fid_out,'Test bsspdfest and bsseval using 1D data:\n');
fprintf(fid_out,'   Distribution                 = 1/3 * N(0,.25) + 2/3 * N(3,0.5)\n');
fprintf(fid_out,'   Evaluation interval          = (%f,%f)\n',evalbnds);
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx1));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y1)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(yn)*(xx1(2)-xx1(1)));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(yn-y1))*(xx1(2)-xx1(1)));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((yn-y1).^2)/length(yn)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((yn-y1).^2)*(xx1(2)-xx1(1))));
fprintf(fid_out,'\n');

y1cdf = bsseval(xx1,bss1test,'cdf');
yncdf = (1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5);
figure;
subplot(2,1,1)
plot(xx1,[y1cdf yncdf])
axis tight;
legend('Estimate','Mixture');
title('Cumulative distribution functions');
subplot(2,1,2)
plot(xx1,yncdf-y1cdf)
axis tight;
title('Error')

y1surv = bsseval(xx1,bss1test,'survivor');
ynsurv = 1-((1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5));
figure;
subplot(2,1,1)
plot(xx1,[y1surv ynsurv])
axis tight;
legend('Estimate','Mixture');
title('Survivor functions');
subplot(2,1,2)
axis tight;
plot(xx1,ynsurv-y1surv)
title('Error')

y1ch = bsseval(xx1,bss1test,'cumhazard');
ynch = 1-((1/3)*normcdf(xx1,0,0.25) + (2/3)*normcdf(xx1,3,0.5));
gt0 = ynch>0;
ynch(~gt0) = nan;
ynch(gt0)  = -log(ynch(gt0));
figure;
subplot(2,1,1)
plot(xx1,[y1ch ynch])
axis tight;
legend('Estimate','Mixture)','location','northwest');
title('Cumulative hazard functions');
subplot(2,1,2)
plot(xx1,ynch-y1ch)
axis tight;
title('Error')

%xx1p = linspace(0,1,501)';
%y1icdf = bsseval(xx1p,bss1test,'icdf');
%xx = linspace(-2,6,1001)';
%ynpdf = (1/3)*normpdf(xx,0,0.25) + (2/3)*normpdf(xx,3,0.5);
%yncdf = cumtrapz(xx,ynpdf);
%ztol = 2*max([min(yncdf) min(abs(1-yncdf))]);
%idx0 = find(yncdf<ztol,1,'last');
%idx1 = find(abs(1-yncdf)<ztol,1,'first');
%idxbad  = find(diff(yncdf)==0);
%idxbad  = idxbad((idxbad>=idx0 & idxbad <=idx1));
%idxgood = (idx0:idx1)';
%idxgood = setdiff(idxgood,idxbad);
%ynicdf  = interp1(yncdf(idxgood),xx(idxgood),xx1p,'linear','extrap');
%figure;
%subplot(2,1,1)
%plot(xx1p,[y1icdf ynicdf])
%legend('Estimate','Mixture','location','northwest');
%title('Inverse cumulative distribution functions');
%subplot(2,1,2)
%plot(xx1p,ynicdf-y1icdf)
%title('Error')

%
%-------------------------------------------------------------------------
%

%
% 2-D example figures
%

%
% Generate data from a standard 2D normal distribution and estimate
% the probability density function using bsspdfest and bsseval.
%
% N(0,I)
%
neval = 81;
nsample  = 5000;
ndata    = nsample;
evalbnds = repmat([-5.5 5.5],2,1);
xtest2   = randn(nsample,2);
bss2test = bsspdfest(xtest2);
xx2      = [
    linspace(evalbnds(1,1),evalbnds(1,2),neval)
    linspace(evalbnds(2,1),evalbnds(2,2),neval)
    ]';
[xx2mg1, xx2mg2] = ndgrid(xx2(:,1),xx2(:,2));
xx2mg            = [xx2mg1(:) xx2mg2(:)];
y2mg     = bsseval(xx2mg,bss2test,'pdf');
y2       = reshape(y2mg,size(xx2mg1));
y2nmg    = mvnpdf(xx2mg,[0 0],ones(1,2));
y2n      = reshape(y2nmg,size(xx2mg1));

figure;
plot(xtest2(:,1),xtest2(:,2),'.');
axis tight;
title('Data');
xlabel('X1');
ylabel('X2');
axis([-5 5 -5 5]);

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Error');
xlabel('X1');
ylabel('X2');

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

tmpx = [xx2(1,1) xx2(2,1); xx2(1,2) xx2(2,2)]';
h = diff(tmpx);

fprintf(fid_out,'Test bsspdfest and bsseval using 2D data:\n');
fprintf(fid_out,'   Distribution                 = N(0,I)\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx2mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y2(:))*prod(h));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(y2n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y2n(:)-y2(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y2n(:)-y2(:)).^2)/length(y2n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y2n(:)-y2(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

y2mg     = bsseval(xx2mg,bss2test,'cdf');
y2       = reshape(y2mg,size(xx2mg1));
y2nmg    = mvncdf(xx2mg,[0 0],ones(1,2));
y2n      = reshape(y2nmg,size(xx2mg1));

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated CDF');
xlabel('X1');
ylabel('X2');
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('CDF Error');
xlabel('X1');
ylabel('X2');

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated CDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('CDF Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

y2mg     = bsseval(xx2mg,bss2test,'survivor');
y2       = reshape(y2mg,size(xx2mg1));
y2nmg    = 1-mvncdf(xx2mg,[0 0],ones(1,2));
y2n      = reshape(y2nmg,size(xx2mg1));

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated Survivor');
xlabel('X1');
ylabel('X2');
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True Survivor');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Survivor Error');
xlabel('X1');
ylabel('X2');

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated Survivor');
xlabel('X1');
ylabel('X2');
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True Survivor');
xlabel('X1');
ylabel('X2');
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Survivor Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

%
% Repeat the previous example with 2 times the data.
%
% N(0,I)
%
neval = 81;
nsample  = 10000;
ndata    = nsample;
evalbnds = repmat([-5.5 5.5],2,1);
xtest2   = randn(nsample,2);
bss2test = bsspdfest(xtest2);
xx2      = [
    linspace(evalbnds(1,1),evalbnds(1,2),neval)
    linspace(evalbnds(2,1),evalbnds(2,2),neval)
    ]';
[xx2mg1, xx2mg2] = ndgrid(xx2(:,1),xx2(:,2));
xx2mg            = [xx2mg1(:) xx2mg2(:)];
y2mg     = bsseval(xx2mg,bss2test,'pdf');
y2       = reshape(y2mg,size(xx2mg1));
y2nmg    = mvnpdf(xx2mg,[0 0],ones(1,2));
y2n      = reshape(y2nmg,size(xx2mg1));

figure;
plot(xtest2(:,1),xtest2(:,2),'.');
axis tight;
title('Data');
xlabel('X1');
ylabel('X2');

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Error');
xlabel('X1');
ylabel('X2');

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

tmpx = [xx2(1,1) xx2(2,1); xx2(1,2) xx2(2,2)]';
h = diff(tmpx);

fprintf(fid_out,'Test bsspdfest and bsseval using 2D data:\n');
fprintf(fid_out,'   Distribution                 = N(0,I)\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx2mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y2(:))*prod(h));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(y2n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y2n(:)-y2(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y2n(:)-y2(:)).^2)/length(y2n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y2n(:)-y2(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

y2mg     = bsseval(xx2mg,bss2test,'cdf');
y2       = reshape(y2mg,size(xx2mg1));
y2nmg    = mvncdf(xx2mg,[0 0],ones(1,2));
y2n      = reshape(y2nmg,size(xx2mg1));

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated CDF');
xlabel('X1');
ylabel('X2');
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('CDF Error');
xlabel('X1');
ylabel('X2');

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated CDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('CDF Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

y2mg     = bsseval(xx2mg,bss2test,'survivor');
y2       = reshape(y2mg,size(xx2mg1));
y2nmg    = 1-mvncdf(xx2mg,[0 0],ones(1,2));
y2n      = reshape(y2nmg,size(xx2mg1));

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated Survivor');
xlabel('X1');
ylabel('X2');
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True Survivor');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Survivor Error');
xlabel('X1');
ylabel('X2');

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated Survivor');
xlabel('X1');
ylabel('X2');
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True Survivor');
xlabel('X1');
ylabel('X2');
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('Survivor Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

%
%  Now generate data for a mixture of 5 2d normal distributions with
%  elliptical shapes and different orientations. The PDF is estimated and
%  then compared with the actual distribution.
%
neval    = 81;
nsample  = 10000;
ndata    = 5*nsample;
xmu      = [0 0; -3 -3; 3 3; 3 -3; -3 3];
xcov0    = [0.9 -0.4; -0.4 0.3];
xcov3    = [0.9  0.4;  0.4 0.3];
chxcov0  = cholcov(xcov0);
chxcov3  = cholcov(xcov3);
evalbnds = [-7 7;-6 6];

%
% Generate the data and estimate pdf
%
xtest2  = [
    repmat(xmu(1,:),nsample,1) + randn(nsample,2)*chxcov0;
    repmat(xmu(2,:),nsample,1) + randn(nsample,2)*chxcov3;
    repmat(xmu(3,:),nsample,1) + randn(nsample,2)*chxcov3;
    repmat(xmu(4,:),nsample,1) + randn(nsample,2)*chxcov3;
    repmat(xmu(5,:),nsample,1) + randn(nsample,2)*chxcov3;
    ];
bss2test = bsspdfest(xtest2);
xx2      = [
    linspace(evalbnds(1,1),evalbnds(1,2),neval)
    linspace(evalbnds(2,1),evalbnds(2,2),neval)
    ]';
[xx2mg1, xx2mg2] = ndgrid(xx2(:,1),xx2(:,2));
xx2mg            = [xx2mg1(:) xx2mg2(:)];
y2mg     = bsseval(xx2mg,bss2test,'pdf');
y2       = reshape(y2mg,size(xx2mg1));

%
% Compute the true density function
%
y2nmg        = zeros(length(xx2mg),1);
for i = 1:5
   if (i == 1)
      y2nmg = y2nmg + (1/5)*mvnpdf(xx2mg,xmu(i,:),xcov0);
   else
      y2nmg = y2nmg + (1/5)*mvnpdf(xx2mg,xmu(i,:),xcov3);
   end
end
y2n = reshape(y2nmg,size(xx2mg1));

tmpx = [xx2(1,1) xx2(2,1); xx2(1,2) xx2(2,2)]';
h = diff(tmpx);

figure;
plot(xtest2(:,1),xtest2(:,2),'.');
axis tight;
title('Data');
xlabel('X1');
ylabel('X2');

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('PDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2)
axis tight;
title('PDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

fprintf(fid_out,'Test bsspdfest and bsseval using 2D data:\n');
fprintf(fid_out,'   Distribution                 = mixture of 5 2d normal distributions\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx2mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y2(:))*prod(h));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(y2n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y2n(:)-y2(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y2n(:)-y2(:)).^2)/length(y2n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y2n(:)-y2(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

y2mg     = bsseval(xx2mg,bss2test,'cdf');
y2       = reshape(y2mg,size(xx2mg1));

y2nmg        = zeros(length(xx2mg),1);
for i = 1:5
   if (i == 1)
      y2nmg = y2nmg + (1/5)*mvncdf(xx2mg,xmu(i,:),xcov0);
   else
      y2nmg = y2nmg + (1/5)*mvncdf(xx2mg,xmu(i,:),xcov3);
   end
end
y2ncdf = reshape(y2nmg,size(xx2mg1));

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated CDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2ncdf)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
surf(xx2mg1,xx2mg2,y2ncdf-y2)
axis tight;
title('CDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated CDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2ncdf)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2ncdf-y2)
axis tight;
title('CDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

y2mg     = bsseval(xx2mg,bss2test,'survivor');
y2       = reshape(y2mg,size(xx2mg1));
y2nsurv   = 1-y2ncdf;

figure;
subplot(2,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated Survivor');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(2,1,2)
surf(xx2mg1,xx2mg2,y2nsurv)
axis tight;
title('True Survivor');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
surf(xx2mg1,xx2mg2,y2nsurv-y2)
axis tight;
title('Survivor Error');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
subplot(3,1,1)
axis tight;
surf(xx2mg1,xx2mg2,y2)
title('Estimated Survivor');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2nsurv)
axis tight;
title('True Survivor');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2nsurv-y2)
axis tight;
title('Survivor Error');
xlabel('X1');
ylabel('X2');
view([-14 52])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

%
%-------------------------------------------------------------------------
%

%
% 3-D example figures
%

%
% Generate data from a standard 3D normal distribution and estimate
% the probability density function using the direct 3D functions:
% bsspdfest3d and bsseval3d.
%
% N(0,I)
%
neval    = 121;
sliceidx = 1+floor(neval/2);
nsample  = 50000;
ndata    = nsample;
evalbnds = repmat([-5.5 5.5],3,1);
xx3      = [
    linspace(evalbnds(1,1),evalbnds(1,2),neval);
    linspace(evalbnds(2,1),evalbnds(2,2),neval);
    linspace(evalbnds(3,1),evalbnds(3,2),neval);
    ]';
[xx3mg1, xx3mg2, xx3mg3] = ndgrid(xx3(:,1),xx3(:,2),xx3(:,3));
xx3mg            = [xx3mg1(:) xx3mg2(:) xx3mg3(:)];
xtest3           = randn(nsample,3);
bss3test         = bsspdfest(xtest3);
y3mg             = bsseval(xx3mg,bss3test,'pdf');
y3               = reshape(y3mg,size(xx3mg1));
y3nmg            = mvnpdf(xx3mg,[0 0 0],ones(1,3));
y3n              = reshape(y3nmg,size(xx3mg1));

figure;
plot3(xtest3(:,1),xtest3(:,2),xtest3(:,3),'.');
title('Data');
xlabel('X1');
ylabel('X2');
zlabel('X3');
grid;

figure;
plotmatrix(xtest3);
 
figure;
subplot(2,1,1)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3(:,:,sliceidx)))
title('Estimated PDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx)))
title('True PDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3(:,:,sliceidx)))
title('Estimated PDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx)))
title('True PDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx))-squeeze(y3(:,:,sliceidx)));
axis tight
title('PDF Error, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3(:,sliceidx,:)))
title('Estimated PDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:)))
title('True PDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3(:,sliceidx,:)))
title('Estimated PDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:)))
title('True PDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:))-squeeze(y3(:,sliceidx,:)));
axis tight;
title('PDF Error, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3(sliceidx,:,:)))
title('Estimated PDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:)))
title('True PDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3(sliceidx,:,:)))
title('Estimated PDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:)))
title('True PDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:))-squeeze(y3(sliceidx,:,:)));
axis tight
title('PDF Error, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

tmpx = xx3(1:2,:);
h = diff(tmpx)';

fprintf(fid_out,'Test bsspdfest and bsseval using 3D data:\n');
fprintf(fid_out,'   Distribution                 = N(0,I)\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Evaluation interval, x3      = (%f,%f)\n',evalbnds(3,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx3mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y3(:))*prod(h));
fprintf(fid_out,'   Area under the N(0,1) pdf    = %f\n',sum(y3n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y3n(:)-y3(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y3n(:)-y3(:)).^2)/length(y3n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y3n(:)-y3(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

y3mg             = bsseval(xx3mg,bss3test,'cdf');
y3               = reshape(y3mg,size(xx3mg1));
y3nmg            = mvncdf(xx3mg,[0 0 0],ones(1,3));
y3n              = reshape(y3nmg,size(xx3mg1));
 
figure;
subplot(2,1,1)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3(:,:,sliceidx)))
title('Estimated CDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx)))
title('True CDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3(:,:,sliceidx)))
title('Estimated CDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx)))
title('True CDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx))-squeeze(y3(:,:,sliceidx)));
axis tight
title('CDF Error, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3(:,sliceidx,:)))
title('Estimated CDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:)))
title('True CDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3(:,sliceidx,:)))
title('Estimated CDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:)))
title('True CDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:))-squeeze(y3(:,sliceidx,:)));
axis tight
title('CDF Error, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3(sliceidx,:,:)))
title('Estimated CDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:)))
title('True CDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3(sliceidx,:,:)))
title('Estimated CDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:)))
title('True CDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:))-squeeze(y3(sliceidx,:,:)));
axis tight
title('CDF Error, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

y3mg             = bsseval(xx3mg,bss3test,'survivor');
y3               = reshape(y3mg,size(xx3mg1));
y3nmg            = 1-mvncdf(xx3mg,[0 0 0],ones(1,3));
y3n              = reshape(y3nmg,size(xx3mg1));
 
figure;
subplot(2,1,1)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3(:,:,sliceidx)))
title('Estimated Survivor, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx)))
title('True Survivor, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3(:,:,sliceidx)))
title('Estimated Survivor, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx)))
title('True Survivor, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg1(:,:,sliceidx)),squeeze(xx3mg2(:,:,sliceidx)),squeeze(y3n(:,:,sliceidx))-squeeze(y3(:,:,sliceidx)));
axis tight
title('Survivor Error, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3(:,sliceidx,:)))
title('Estimated Survivor, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:)))
title('True Survivor, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3(:,sliceidx,:)))
title('Estimated Survivor, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:)))
title('True Survivor, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg1(:,sliceidx,:)),squeeze(xx3mg3(:,sliceidx,:)),squeeze(y3n(:,sliceidx,:))-squeeze(y3(:,sliceidx,:)));
axis tight
title('Survivor Error, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3(sliceidx,:,:)))
title('Estimated Survivor, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(2,1,2)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:)))
title('True Survivor, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight

figure;
subplot(3,1,1)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3(sliceidx,:,:)))
title('Estimated Survivor, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:)))
title('True Survivor, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,3)
surf(squeeze(xx3mg2(sliceidx,:,:)),squeeze(xx3mg3(sliceidx,:,:)),squeeze(y3n(sliceidx,:,:))-squeeze(y3(sliceidx,:,:)));
axis tight
title('Survivor Error, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

%
%  Now generate data for a mixture of 5 3d normal distributions with
%  N(0,I) at (0,0,0) and elliptical shapes for the others located at
%  (+-3,+-3,0). The PDF is estimated and  then compared with the actual
%  distribution.
%
neval    = 151;
%sliceidx = 1+floor(neval/2);
nsample  = 5000;
ndata    = 5*nsample;
xmu      = [0 0 0; -3 -3 0; 3 3 0; 3 -3 0; -3 3 0];
xcov0    = [1 0 0; 0 1 0; 0 0 1];
xcov3    = [0.9  0 0;  0 0.3 0; 0 0 1];
chxcov0  = cholcov(xcov0);
chxcov3  = cholcov(xcov3);
evalbnds = [-10 10;-10 10;-10 10];
xx3      = [
    linspace(evalbnds(1,1),evalbnds(1,2),neval);
    linspace(evalbnds(2,1),evalbnds(2,2),neval);
    linspace(evalbnds(3,1),evalbnds(3,2),neval);
    ]';
[xx3mg1, xx3mg2, xx3mg3] = ndgrid(xx3(:,1),xx3(:,2),xx3(:,3));
xx3mg            = [xx3mg1(:) xx3mg2(:) xx3mg3(:)];

%
% Compute the true density function
%
y3nmg        = zeros(length(xx3mg),1);
for i = 1:5
   if (i == 1)
      y3nmg = y3nmg + (1/5)*mvnpdf(xx3mg,xmu(i,:),xcov0);
   else
      y3nmg = y3nmg + (1/5)*mvnpdf(xx3mg,xmu(i,:),xcov3);
   end
end
y3n = reshape(y3nmg,size(xx3mg1));

%
% Generate the data and estimate pdf
%
xtest3  = [
    repmat(xmu(1,:),nsample,1) + randn(nsample,3)*chxcov0;
    repmat(xmu(2,:),nsample,1) + randn(nsample,3)*chxcov3;
    repmat(xmu(3,:),nsample,1) + randn(nsample,3)*chxcov3;
    repmat(xmu(4,:),nsample,1) + randn(nsample,3)*chxcov3;
    repmat(xmu(5,:),nsample,1) + randn(nsample,3)*chxcov3;
    ];
bss3test = bsspdfest(xtest3);
y3mg     = bsseval(xx3mg,bss3test,'pdf');
y3       = reshape(y3mg,size(xx3mg1));

tmpx = xx3(1:2,:);
h = diff(tmpx)';

figure;
plot3(xtest3(:,1),xtest3(:,2),xtest3(:,3),'.');
axis tight;
view([-14 51]);
title('Data');
xlabel('X1');
ylabel('X2');
zlabel('X3');
grid;

figure;
[~,AX,~,~,PAx]=plotmatrix(xtest3);
subplot(AX(1,1));
axis([-7 7 -7 7]);
subplot(AX(1,2));
axis([-7 7 -7 7]);
subplot(AX(1,3));
axis([-7 7 -7 7]);
subplot(AX(2,1));
axis([-7 7 -7 7]);
subplot(AX(2,2));
axis([-7 7 -7 7]);
subplot(AX(2,3));
axis([-7 7 -7 7]);
subplot(AX(3,1));
axis([-7 7 -7 7]);
subplot(AX(3,2));
axis([-7 7 -7 7]);
subplot(AX(3,3));
axis([-7 7 -7 7]);
subplot(PAx(1));
a = axis;
axis([-7 7 a(3:4)]);
subplot(PAx(2));
a = axis;
axis([-7 7 a(3:4)]);
subplot(PAx(3));
a = axis;
axis([-7 7 a(3:4)]);
 
fprintf(fid_out,'Test bsspdfest and bsseval using 3D data:\n');
fprintf(fid_out,'   Distribution                 = 3d mixture\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Evaluation interval, x3      = (%f,%f)\n',evalbnds(3,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx3mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y3(:))*prod(h));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(y3n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y3n(:)-y3(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y3n(:)-y3(:)).^2)/length(y3n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y3n(:)-y3(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

%
% Repeat the previous example with twice the data.
%
neval    = 151;
%sliceidx = 1+floor(neval/2);
nsample  = 10000;
ndata    = 5*nsample;
xmu      = [0 0 0; -3 -3 0; 3 3 0; 3 -3 0; -3 3 0];
xcov0    = [0.9 0 0; 0 0.3 0; 0 0 1];
xcov3    = [0.9  0 0;  0 0.3 0; 0 0 1];
chxcov0  = cholcov(xcov0);
chxcov3  = cholcov(xcov3);
evalbnds = [-10 10;-10 10;-10 10];
xx3      = [
    linspace(evalbnds(1,1),evalbnds(1,2),neval);
    linspace(evalbnds(2,1),evalbnds(2,2),neval);
    linspace(evalbnds(3,1),evalbnds(3,2),neval);
    ]';
[xx3mg1, xx3mg2, xx3mg3] = ndgrid(xx3(:,1),xx3(:,2),xx3(:,3));
xx3mg            = [xx3mg1(:) xx3mg2(:) xx3mg3(:)];

%
% Compute the true density function
%
y3nmg        = zeros(length(xx3mg),1);
for i = 1:5
   if (i == 1)
      y3nmg = y3nmg + (1/5)*mvnpdf(xx3mg,xmu(i,:),xcov0);
   else
      y3nmg = y3nmg + (1/5)*mvnpdf(xx3mg,xmu(i,:),xcov3);
   end
end
y3n = reshape(y3nmg,size(xx3mg1));

%
% Generate the data and estimate pdf
%
xtest3  = [
    repmat(xmu(1,:),nsample,1) + randn(nsample,3)*chxcov0;
    repmat(xmu(2,:),nsample,1) + randn(nsample,3)*chxcov3;
    repmat(xmu(3,:),nsample,1) + randn(nsample,3)*chxcov3;
    repmat(xmu(4,:),nsample,1) + randn(nsample,3)*chxcov3;
    repmat(xmu(5,:),nsample,1) + randn(nsample,3)*chxcov3;
    ];
bss3test = bsspdfest(xtest3);
y3mg     = bsseval(xx3mg,bss3test,'pdf');
y3       = reshape(y3mg,size(xx3mg1));

tmpx = xx3(1:2,:);
h = diff(tmpx)';

figure;
plot3(xtest3(:,1),xtest3(:,2),xtest3(:,3),'.');
axis tight;
view([-14 51]);
title('Data');
xlabel('X1');
ylabel('X2');
zlabel('X3');
grid;

figure;
[~,AX,~,~,PAx]=plotmatrix(xtest3);
subplot(AX(1,1));
axis([-7 7 -7 7]);
subplot(AX(1,2));
axis([-7 7 -7 7]);
subplot(AX(1,3));
axis([-7 7 -7 7]);
subplot(AX(2,1));
axis([-7 7 -7 7]);
subplot(AX(2,2));
axis([-7 7 -7 7]);
subplot(AX(2,3));
axis([-7 7 -7 7]);
subplot(AX(3,1));
axis([-7 7 -7 7]);
subplot(AX(3,2));
axis([-7 7 -7 7]);
subplot(AX(3,3));
axis([-7 7 -7 7]);
subplot(PAx(1));
a = axis;
axis([-7 7 a(3:4)]);
subplot(PAx(2));
a = axis;
axis([-7 7 a(3:4)]);
subplot(PAx(3));
a = axis;
axis([-7 7 a(3:4)]);

fprintf(fid_out,'Test bsspdfest and bsseval using 3D data:\n');
fprintf(fid_out,'   Distribution                 = 3d mixture\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Evaluation interval, x3      = (%f,%f)\n',evalbnds(3,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx3mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y3(:))*prod(h));
fprintf(fid_out,'   Area under the true pdf      = %f\n',sum(y3n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y3n(:)-y3(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y3n(:)-y3(:)).^2)/length(y3n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y3n(:)-y3(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

%
%-------------------------------------------------------------------------
%

%
% N-D example figures
%

%
% Generate data from a standard 4D normal distribution and estimate
% the probability density function using the functions:
% bsspdfest and bsseval.
%
% N(0,I)
%
neval            = 51;
sliceidx = 1+floor(neval/2);
nsample          = 100000;
ndata            = nsample;
evalbnds = repmat([-5.5 5.5],4,1);
xx4              = [
   linspace(evalbnds(1,1),evalbnds(1,2),neval);
   linspace(evalbnds(2,1),evalbnds(2,2),neval);
   linspace(evalbnds(3,1),evalbnds(3,2),neval);
   linspace(evalbnds(4,1),evalbnds(4,2),neval);
]';
[xx4mg1, xx4mg2, xx4mg3, xx4mg4] = ndgrid(xx4(:,1),xx4(:,2),xx4(:,3),xx4(:,4));
xx4mg            = [xx4mg1(:) xx4mg2(:) xx4mg3(:) xx4mg4(:)];
xtest4           = randn(nsample,4);
tic;
bss4test         = bsspdfest(xtest4);
tpdf = toc;
tic;
y4mg             = bsseval(xx4mg,bss4test,'pdf');
teval=toc;
y4               = reshape(y4mg,size(xx4mg1));
y4nmg            = mvnpdf(xx4mg,[0 0 0 0],ones(1,4));
y4n              = reshape(y4nmg,size(xx4mg1));

figure;
[H,AX,BigAx,P,PAx]=plotmatrix(xtest4);
 
figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),squeeze(y4(:,:,sliceidx,sliceidx)))
title('Estimated PDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),squeeze(y4n(:,:,sliceidx,sliceidx)))
title('True PDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,:,sliceidx,sliceidx))-squeeze(y4(:,:,sliceidx,sliceidx));
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),tmp)
axis tight
title('PDF Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),squeeze(y4(:,sliceidx,:,sliceidx)))
title('Estimated PDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),squeeze(y4n(:,sliceidx,:,sliceidx)))
title('True PDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,sliceidx,:,sliceidx))-squeeze(y4(:,sliceidx,:,sliceidx));
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),tmp);
axis tight
title('PDF Error');
xlabel('X1');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),squeeze(y4(:,sliceidx,sliceidx,:)))
title('Estimated PDF, dimensions 1 and 4');
xlabel('X1');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),squeeze(y4n(:,sliceidx,sliceidx,:)))
title('True PDF, dimensions 1 and 4');
xlabel('X1');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,sliceidx,sliceidx,:))-squeeze(y4(:,sliceidx,sliceidx,:));
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),tmp)
axis tight
title('PDF Error');
xlabel('X1');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),squeeze(y4(sliceidx,:,:,sliceidx)))
title('Estimated PDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),squeeze(y4n(sliceidx,:,:,sliceidx)))
title('True PDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,:,:,sliceidx))-squeeze(y4(sliceidx,:,:,sliceidx));
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),tmp)
axis tight
title('PDF Error');
xlabel('X2');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),squeeze(y4(sliceidx,:,sliceidx,:)))
title('Estimated PDF, dimensions 2 and 4');
xlabel('X2');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),squeeze(y4n(sliceidx,:,sliceidx,:)))
title('True PDF, dimensions 2 and 4');
xlabel('X2');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,:,sliceidx,:))-squeeze(y4(sliceidx,:,sliceidx,:));
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),tmp)
axis tight
title('PDF Error');
xlabel('X2');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),squeeze(y4(sliceidx,sliceidx,:,:)))
title('Estimated PDF, dimensions 3 and 4');
xlabel('X3');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),squeeze(y4n(sliceidx,sliceidx,:,:)))
title('True PDF, dimensions 3 and 4');
xlabel('X3');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,sliceidx,:,:))-squeeze(y4(sliceidx,sliceidx,:,:));
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),tmp)
axis tight
title('PDF Error');
xlabel('X3');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

tmpx = xx4(1:2,:);
h = diff(tmpx)';

fprintf(fid_out,'Test bsspdfest and bsseval using ND data:\n');
fprintf(fid_out,'   Distribution                 = N(0,I)\n');
fprintf(fid_out,'   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf(fid_out,'   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf(fid_out,'   Evaluation interval, x3      = (%f,%f)\n',evalbnds(3,:));
fprintf(fid_out,'   Evaluation interval, x4      = (%f,%f)\n',evalbnds(4,:));
fprintf(fid_out,'   Number of data points        = %d\n',ndata);
fprintf(fid_out,'   Number of evaluation points  = %d\n',length(xx4mg));
fprintf(fid_out,'   Area under the estimated pdf = %f\n',sum(y4(:))*prod(h));
fprintf(fid_out,'   Area under the N(0,I) pdf    = %f\n',sum(y4n(:))*prod(h));
fprintf(fid_out,'   Integrated absolute error    = %f\n',sum(abs(y4n(:)-y4(:)))*prod(h));
fprintf(fid_out,'   RMSE                         = %f\n',sqrt(sum((y4n(:)-y4(:)).^2)/length(y4n)));
fprintf(fid_out,'   IRMSE                        = %f\n',sqrt(sum((y4n(:)-y4(:)).^2)*prod(h)));
fprintf(fid_out,'\n');

y4mg             = bsseval(xx4mg,bss4test,'cdf');
y4               = reshape(y4mg,size(xx4mg1));
y4nmg            = mvncdf(xx4mg,[0 0 0 0],ones(1,4));
y4n              = reshape(y4nmg,size(xx4mg1));

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),squeeze(y4(:,:,sliceidx,sliceidx)))
title('Estimated CDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),squeeze(y4n(:,:,sliceidx,sliceidx)))
title('True CDF, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,:,sliceidx,sliceidx))-squeeze(y4(:,:,sliceidx,sliceidx));
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),tmp)
axis tight
title('CDF Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),squeeze(y4(:,sliceidx,:,sliceidx)))
title('Estimated CDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),squeeze(y4n(:,sliceidx,:,sliceidx)))
title('True CDF, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,sliceidx,:,sliceidx))-squeeze(y4(:,sliceidx,:,sliceidx));
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),tmp);
axis tight
title('CDF Error');
xlabel('X1');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),squeeze(y4(:,sliceidx,sliceidx,:)))
title('Estimated CDF, dimensions 1 and 4');
xlabel('X1');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),squeeze(y4n(:,sliceidx,sliceidx,:)))
title('True CDF, dimensions 1 and 4');
xlabel('X1');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,sliceidx,sliceidx,:))-squeeze(y4(:,sliceidx,sliceidx,:));
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),tmp)
axis tight
title('CDF Error');
xlabel('X1');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),squeeze(y4(sliceidx,:,:,sliceidx)))
title('Estimated CDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),squeeze(y4n(sliceidx,:,:,sliceidx)))
title('True CDF, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,:,:,sliceidx))-squeeze(y4(sliceidx,:,:,sliceidx));
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),tmp)
axis tight
title('CDF Error');
xlabel('X2');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),squeeze(y4(sliceidx,:,sliceidx,:)))
title('Estimated CDF, dimensions 2 and 4');
xlabel('X2');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),squeeze(y4n(sliceidx,:,sliceidx,:)))
title('True CDF, dimensions 2 and 4');
xlabel('X2');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,:,sliceidx,:))-squeeze(y4(sliceidx,:,sliceidx,:));
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),tmp)
axis tight
title('CDF Error');
xlabel('X2');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),squeeze(y4(sliceidx,sliceidx,:,:)))
title('Estimated CDF, dimensions 3 and 4');
xlabel('X3');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),squeeze(y4n(sliceidx,sliceidx,:,:)))
title('True CDF, dimensions 3 and 4');
xlabel('X3');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,sliceidx,:,:))-squeeze(y4(sliceidx,sliceidx,:,:));
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),tmp)
axis tight
title('CDF Error');
xlabel('X3');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

y4mg             = bsseval(xx4mg,bss4test,'survivor');
y4               = reshape(y4mg,size(xx4mg1));
y4nmg            = 1-mvncdf(xx4mg,[0 0 0 0],ones(1,4));
y4n              = reshape(y4nmg,size(xx4mg1));

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),squeeze(y4(:,:,sliceidx,sliceidx)))
title('Estimated Survivor, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),squeeze(y4n(:,:,sliceidx,sliceidx)))
title('True Survivor, dimensions 1 and 2');
xlabel('X1');
ylabel('X2');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,:,sliceidx,sliceidx))-squeeze(y4(:,:,sliceidx,sliceidx));
surf(squeeze(xx4mg1(:,:,sliceidx,sliceidx)),squeeze(xx4mg2(:,:,sliceidx,sliceidx)),tmp)
axis tight
title('Survivor Error');
xlabel('X1');
ylabel('X2');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),squeeze(y4(:,sliceidx,:,sliceidx)))
title('Estimated Survivor, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),squeeze(y4n(:,sliceidx,:,sliceidx)))
title('True Survivor, dimensions 1 and 3');
xlabel('X1');
ylabel('X3');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,sliceidx,:,sliceidx))-squeeze(y4(:,sliceidx,:,sliceidx));
surf(squeeze(xx4mg1(:,sliceidx,:,sliceidx)),squeeze(xx4mg3(:,sliceidx,:,sliceidx)),tmp);
axis tight
title('Survivor Error');
xlabel('X1');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),squeeze(y4(:,sliceidx,sliceidx,:)))
title('Estimated Survivor, dimensions 1 and 4');
xlabel('X1');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),squeeze(y4n(:,sliceidx,sliceidx,:)))
title('True Survivor, dimensions 1 and 4');
xlabel('X1');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(:,sliceidx,sliceidx,:))-squeeze(y4(:,sliceidx,sliceidx,:));
surf(squeeze(xx4mg1(:,sliceidx,sliceidx,:)),squeeze(xx4mg4(:,sliceidx,sliceidx,:)),tmp)
axis tight
title('Survivor Error');
xlabel('X1');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),squeeze(y4(sliceidx,:,:,sliceidx)))
title('Estimated Survivor, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),squeeze(y4n(sliceidx,:,:,sliceidx)))
title('True Survivor, dimensions 2 and 3');
xlabel('X2');
ylabel('X3');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,:,:,sliceidx))-squeeze(y4(sliceidx,:,:,sliceidx));
surf(squeeze(xx4mg2(sliceidx,:,:,sliceidx)),squeeze(xx4mg3(sliceidx,:,:,sliceidx)),tmp)
axis tight
title('Survivor Error');
xlabel('X2');
ylabel('X3');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),squeeze(y4(sliceidx,:,sliceidx,:)))
title('Estimated Survivor, dimensions 2 and 4');
xlabel('X2');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),squeeze(y4n(sliceidx,:,sliceidx,:)))
title('True Survivor, dimensions 2 and 4');
xlabel('X2');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,:,sliceidx,:))-squeeze(y4(sliceidx,:,sliceidx,:));
surf(squeeze(xx4mg2(sliceidx,:,sliceidx,:)),squeeze(xx4mg4(sliceidx,:,sliceidx,:)),tmp)
axis tight
title('Survivor Error');
xlabel('X2');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(3,1,1)
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),squeeze(y4(sliceidx,sliceidx,:,:)))
title('Estimated Survivor, dimensions 3 and 4');
xlabel('X3');
ylabel('X4');
axis tight
subplot(3,1,2)
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),squeeze(y4n(sliceidx,sliceidx,:,:)))
title('True Survivor, dimensions 3 and 4');
xlabel('X3');
ylabel('X4');
axis tight
subplot(3,1,3)
tmp = squeeze(y4n(sliceidx,sliceidx,:,:))-squeeze(y4(sliceidx,sliceidx,:,:));
surf(squeeze(xx4mg3(sliceidx,sliceidx,:,:)),squeeze(xx4mg4(sliceidx,sliceidx,:,:)),tmp)
axis tight
title('Survivor Error');
xlabel('X3');
ylabel('X4');
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

%
%-------------------------------------------------------------------------
%

%
% 1-D and 2-D bounded domain example figures.
%

%
% 1-d example: uniform distribution
%
rng('default');
xu1 = rand(100000,1);
bssunb1 = bsspdfest(xu1);
bssu1 = bsspdfest(xu1,[0 1]);
xxu1 = linspace(-.1,1.1,5001)';
yut1 = unifpdf(xxu1);
yu = bsseval(xxu1,bssu1);
yunb1 = bsseval(xxu1,bssunb1);

figure;
subplot(2,1,1)
%plot(xxu1,yut1,'k-',xxu1,yu,'k-',xxu1,yunb1,'r-','linewidth',1);
plot(xxu1,[yut1 yu yunb1],'-');
grid;
axis tight;
legend('U(0,1)','Bounded','Unbounded','location','South');
title('Uniform probability density functions');
subplot(2,1,2)
%plot(xxu1,yut1-yu,'k-',xxu1,yut1-yunb1,'r-','linewidth',1);
plot(xxu1,[yut1-yu yut1-yunb1],'-');
grid;
legend('Bounded','Unbounded','location','South');
axis tight;
title('Error')

%
% 1-d example: uniform distribution, m=1, small partition
%
rng('default');
xu1 = rand(100000,1);
bssunb1 = bsspdfest(xu1,[],17,1);
bssu1 = bsspdfest(xu1,[0 1],11,1);
xxu1 = linspace(-.1,1.1,5001)';
yut1 = unifpdf(xxu1);
yu = bsseval(xxu1,bssu1);
yunb1 = bsseval(xxu1,bssunb1);

figure;
subplot(2,1,1)
%plot(xxu1,yut1,'k-',xxu1,yu,'k-',xxu1,yunb1,'r-','linewidth',1);
plot(xxu1,[yut1 yu yunb1],'-');
grid;
axis tight;
legend('U(0,1)','Bounded','Unbounded','location','South');
title('Uniform probability density functions');
subplot(2,1,2)
%plot(xxu1,yut1-yu,'k-',xxu1,yut1-yunb1,'r-','linewidth',1);
plot(xxu1,[yut1-yu yut1-yunb1],'-');
grid;
legend('Bounded','Unbounded','location','North');
axis tight;
title('Error')

%
% 1-d example: exponential distribution
%
rng('default');
xe = exprnd(2,100000,1);
bsse = bsspdfest(xe,[0 22]);
bssenb = bsspdfest(xe,[-1 22]);
xxe = linspace(-1,22,5001)';
yet = exppdf(xxe,2);
ye = bsseval(xxe,bsse);
yenb = bsseval(xxe,bssenb);
figure;
subplot(2,1,1)
%plot(xxe,yet,'k-',xxe,ye,'k-',xxe,yenb,'r-','linewidth',1);
plot(xxe,[yet ye yenb],'-');
grid;
axis tight;
legend('U(0,1)','Bounded','Unbounded','location','NorthEast');
title('Exponential probability density functions, \mu=2');
subplot(2,1,2)
%plot(xxe,yet-ye,'k-',xxe,yet-yenb,'r-','linewidth',1);
plot(xxe,[yet-ye yet-yenb],'-');
grid;
legend('Bounded','Unbounded','location','NorthEast');
axis tight;
title('Error')

%
% 1-d example: truncated normal
%
rng('default');
xtn1 = randn(100000,1);
xtn1 = xtn1(xtn1>-2.0&xtn1<1.5);
bsstn1 = bsspdfest(xtn1,[-2.0 1.5]);
bsstnnb1 = bsspdfest(xtn1);
xxtn1 = linspace(-2.25,1.75,5001)';
inbounds = xxtn1>-2.0 & xxtn1<1.5;
xxtnt1 = xxtn1(inbounds);
ytnt1 = zeros(size(xxtn1));
ytnt1(inbounds) = normpdf(xxtnt1);
ytnt1(inbounds) = ytnt1(inbounds)./trapz(xxtnt1,ytnt1(inbounds));
ytn1 = bsseval(xxtn1,bsstn1);
ytnnb1 = bsseval(xxtn1,bsstnnb1);

figure;
subplot(2,1,1)
%plot(xxe,ytnt1,'k-',xxtn1,ytn,'k-',xxtn1,ytnnb1,'r-','linewidth',1);
plot(xxtn1,[ytnt1 ytn1 ytnnb1],'-');
grid;
axis tight;
legend('N(0,1,-2.0,1.5)','Bounded','Unbounded','location','South');
title('Truncated N(0,1,-2.0,1.5) probability density functions)');
subplot(2,1,2)
%plot(xxtn1,ytnt1-ytn1,'k-',xxtn1,ytnt1-ytnnb1,'r-','linewidth',1);
plot(xxtn1,[ytnt1-ytn1 ytnt1-ytnnb1],'-');
grid;
legend('Bounded','Unbounded','location','South');
axis tight;
title('Error')

%
% 1-d example: truncated normal new and old methods
%
rng('default');
xtn1 = randn(100000,1);
xtn1 = xtn1(xtn1>-2.0&xtn1<1.5);
bsstnd1 = bsspdfest(xtn1,[-2.0 1.5]);
bsstn1 = bsspdfest(xtn1,[-2.0 1.5],11);
xxtn1 = linspace(-2.25,1.75,5001)';
inbounds = xxtn1>-2.0 & xxtn1<1.5;
xxtnt1 = xxtn1(inbounds);
ytnt1 = zeros(size(xxtn1));
ytnt1(inbounds) = normpdf(xxtnt1);
ytnt1(inbounds) = ytnt1(inbounds)./trapz(xxtnt1,ytnt1(inbounds));
ytnd1 = bsseval(xxtn1,bsstnd1);
ytn1 = bsseval(xxtn1,bsstn1);

figure;
subplot(2,1,1)
%plot(xxe,ytnt1,'k-',xxtn1,ytn,'k-',xxtn1,ytnnb1,'r-','linewidth',1);
plot(xxtn1,[ytnt1 ytnd1],'-');
grid;
axis tight;
legend('N(0,1,-2.0,1.5)','Reflection','location','South');
title('Truncated N(0,1,-2.0,1.5) probability density functions)');
subplot(2,1,2)
%plot(xxtn1,ytnt1-ytn1,'k-',xxtn1,ytnt1-ytnnb1,'r-','linewidth',1);
plot(xxtn1,ytnt1-ytnd1,'-');
grid;
legend('Reflection','location','North');
axis tight;
title('Error')

figure;
subplot(2,1,1)
%plot(xxe,ytnt1,'k-',xxtn1,ytn,'k-',xxtn1,ytnnb1,'r-','linewidth',1);
plot(xxtn1,[ytnt1 ytn1],'-');
grid;
axis tight;
legend('N(0,1,-2.0,1.5)','Reflection','location','South');
title('Truncated N(0,1,-2.0,1.5) probability density functions)');
subplot(2,1,2)
%plot(xxtn1,ytnt1-ytn1,'k-',xxtn1,ytnt1-ytnnb1,'r-','linewidth',1);
plot(xxtn1,ytnt1-ytn1,'-');
grid;
legend('Reflection','location','South');
axis tight;
title('Error')

%
% 2-D Uniform distribution
%
rng('default')
xu2 = rand(200000,2);
bssu2 = bsspdfest(xu2,[0 1;0 1],[23 23]');
bssunb2 = bsspdfest(xu2,[],[31 31]');
xx1u = linspace(-0.1,1.1,101)';
xx2u = linspace(-0.1,1.1,101)';
[xx1ug,xx2ug]= ndgrid(xx1u,xx2u);
xxu2 = [xx1ug(:) xx2ug(:)];
yut2 = double(reshape( (xx1ug(:)>=0 & xx1ug(:)<1)& (xx2ug(:)>=0 & xx2ug(:)<1),size(xx1ug)));
yu2 = reshape(bsseval(xxu2,bssu2),size(xx1ug));
yunb2 = reshape(bsseval(xxu2,bssunb2),size(xx1ug));

figure;
plot(xu2(:,1),xu2(:,2),'.','MarkerSize',1);
axis tight;
title('Data');
xlabel('X1');
ylabel('X2');

figure;
subplot(2,1,1)
surf(xx1ug,xx2ug,yu2)
axis tight;
title('Estimated PDF, bounded');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(2,1,2)
surf(xx1ug,xx2ug,yut2-yu2)
axis tight;
title('PDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
subplot(3,1,1)
surf(xx1ug,xx2ug,yu2)
axis tight;
title('Estimated PDF, bounded');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx1ug,xx2ug,yut2)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx1ug,xx2ug,yut2-yu2)
axis tight;
title('PDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
subplot(2,1,1)
surf(xx1ug,xx2ug,yunb2)
axis tight;
title('Estimated PDF, unbounded');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(2,1,2)
surf(xx1ug,xx2ug,yut2-yunb2)
axis tight;
title('PDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])

figure;
subplot(3,1,1)
surf(xx1ug,xx2ug,yunb2)
axis tight;
title('Estimated PDF, unbounded');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx1ug,xx2ug,yut2)
axis tight;
title('True PDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx1ug,xx2ug,yut2-yunb2)
axis tight;
title('PDF Error');
xlabel('X1');
ylabel('X2');
view([-14 52])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

%
% 2-d example: truncated normal
%
rng('default');
xtn2 = randn(200000,2);
xtn2 = xtn2((xtn2(:,1)>-1.5)&(xtn2(:,2)>-1.5),:);
bsstn2 = bsspdfest(xtn2,[-1.5 5;-1.5 5]);
bsstnnb2 = bsspdfest(xtn2);
xx1tn2 = linspace(-2,5,101)';
xx2tn2 = linspace(-2,5,101)';
[xx1tng,xx2tng]= ndgrid(xx1tn2,xx2tn2);
xxtn2 = [xx1tng(:) xx2tng(:)];
inbounds = (xxtn2(:,1)>-1.5) & (xxtn2(:,2)>-1.5);
xxtnt2 = xxtn2(inbounds,:);
ytnt2 = zeros(size(xxtn2,1),1);
ytnt2(inbounds) = mvnpdf(xxtnt2);
ytnt2 = ytnt2./(sum(ytnt2(:)*prod([xx1tn2(2)-xx1tn2(1) xx2tn2(2)-xx2tn2(1)])));
ytnt2 = reshape(ytnt2,size(xx1tng));
ytn2 = reshape(bsseval(xxtn2,bsstn2),size(xx1tng));
ytnnb2 = reshape(bsseval(xxtn2,bsstnnb2),size(xx1tng));

figure;
plot(xtn2(:,1),xtn2(:,2),'.','MarkerSize',1);
%axis tight;
title('Data');
xlabel('X1');
ylabel('X2');

figure;
surf(xx1tng,xx2tng,ytn2)
%axis tight;
title('Truncated normal estimated PDF, bounded');
xlabel('X1');
ylabel('X2');
view([8 15])

figure;
surf(xx1tng,xx2tng,ytnt2-ytn2)
%axis tight;
title('Truncated normal PDF Error, bounded');
xlabel('X1');
ylabel('X2');
view([8 15])

figure;
subplot(2,1,1)
surf(xx1tng,xx2tng,ytn2)
%axis tight;
title('Truncated normal estimated PDF, bounded');
xlabel('X1');
ylabel('X2');
view([8 15])
subplot(2,1,2)
surf(xx1tng,xx2tng,ytnt2-ytn2)
%axis tight;
title('Truncated normal PDF Error, bounded');
xlabel('X1');
ylabel('X2');
view([8 15])

figure;
subplot(3,1,1)
surf(xx1tng,xx2tng,ytn2)
%axis tight;
title('Truncated normal estimated PDF, bounded');
xlabel('X1');
ylabel('X2');
view([8 15])
subplot(3,1,2)
surf(xx1tng,xx2tng,ytnt2)
%axis tight;
title('Truncated normal true PDF');
xlabel('X1');
ylabel('X2');
view([8 15])
subplot(3,1,3)
surf(xx1tng,xx2tng,ytnt2-ytn2)
%axis tight;
title('Truncated normal PDF Error, bounded');
xlabel('X1');
ylabel('X2');
view([8 15])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);

figure;
surf(xx1tng,xx2tng,ytnnb2)
%axis tight;
title('Truncated normal estimated PDF, unbounded');
xlabel('X1');
ylabel('X2');
view([8 15])

figure;
surf(xx1tng,xx2tng,ytnt2-ytnnb2)
%axis tight;
title('Truncated normal PDF Error, unbounded');
xlabel('X1');
ylabel('X2');
view([8 15])

figure;
subplot(2,1,1)
surf(xx1tng,xx2tng,ytnnb2)
%axis tight;
title('Truncated normal estimated PDF, unbounded');
xlabel('X1');
ylabel('X2');
view([8 15])
subplot(2,1,2)
surf(xx1tng,xx2tng,ytnt2-ytnnb2)
%axis tight;
title('Truncated normal PDF Error, unbounded');
xlabel('X1');
ylabel('X2');
view([8 15])

figure;
subplot(3,1,1)
surf(xx1tng,xx2tng,ytnnb2)
%axis tight;
title('Truncated normal estimated PDF, unbounded');
xlabel('X1');
ylabel('X2');
view([8 15])
subplot(3,1,2)
surf(xx1tng,xx2tng,ytnt2)
%axis tight;
title('Truncated normal true PDF');
xlabel('X1');
ylabel('X2');
view([8 15])
subplot(3,1,3)
surf(xx1tng,xx2tng,ytnt2-ytnnb2)
%axis tight;
title('Truncated normal PDF Error, unbounded');
xlabel('X1');
ylabel('X2');
view([8 15])
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 5.0 7.0]);
