%
% Examples and tests for bsspdfest and bsseval without the statistics and
% machine learning toolbox.
%

%
% Copyright 2014-2023 Biometrics Northwest LLC
%

%
% Modified: 2023-05-15
% Author  : Kevin R. Gehringer
%
%    Reduced the partition size and the number of evaluation points for the
%    N-D test to make it complete more rapidly. Now it takes about an hour
%    rather than 3 or more hours.
%
%    Added a diary file for the results.
%
%    modified the tests so that they all agreed with the script
%    test_bsspdfest.m
%
% Modified: 2023-03-23
% Author  : Kevin R. Gehringer
%
%    Modified for the interface changes from positional arguments to varargin.
%
% Modified: 2022-12-20
% Author  : Kevin R. Gehringer
%
%    Added tests for 3-D and N-D distributions.
%
% Modified: 2016-06-25
% Author  : Kevin R. Gehringer
%
%    Added different examples and examples with various probability
%    functions.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

%
% In each example below the following abbreviations are used:
%
% bss:
%     Estimate the B-spline series and evaluate the probability function.
%     This is comparable to what MATLAB's ksdensity and mvksdensity
%     functions do.
%
% bssci:
%     Create a gridded interpolant for very fast approximate evaluation
%     of a B-spline series.
%
% bssi:
%    Use the gridded interpolant to evaluate the desired probability
%    function.
%
% normpdftestbss:
%     Compute the true one dimensional normal probability density function
%
% mvnpdftestbss:
%     Compute the true multi-dimensional probability density function
%

%
%-----------------------------------------------------------------------
% One dimensional test:
%
%    PDF: Mixture of 3 normal distributions with 2 gaps and zero values
%        between them.
%
%        (1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)
%
%-----------------------------------------------------------------------
%

diary_file = 'test_bsspdfest_results_nostats.txt';
if ( exist(diary_file,'file') )
    delete(diary_file);
end
diary(diary_file);

fprintf('PDF: (1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)');
rng('default');
neval   = 1001;
nsample = 500;
ndata   = 3*nsample;
fprintf(', %d data points\n',ndata);
xtest1 = [0.25*randn(nsample,1); 5+0.5*randn(nsample,1); 11+0.25*randn(nsample,1)];
evalbnds = [-6 15];
xx1    = linspace(evalbnds(1),evalbnds(2),neval)';

fprintf('bss......... ')
tic;
bsstest1 = bsspdfest(xtest1);
y1    = bsseval(xx1,bsstest1,'function','pdf');
toc;

fprintf('bssci....... ')
tic;
[~,bssi1]  = bsseval([],bsstest1,'function','pdf','createinterpolant');
toc;

fprintf('bssi........ ')
tic;
yi1   = bsseval(xx1,bssi1,'function','pdf');
toc;

fprintf('normpdf..... ');
tic;
yn1 = (1/3)*normpdftestbss(xx1,0,0.25) + (1/3)*normpdftestbss(xx1,5,0.5) + (1/3)*normpdftestbss(xx1,11,0.25);
toc;

figure;
subplot(2,1,1);
plot(xx1,[yn1 y1 yi1]);
axis tight;
title('PDF:(1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)')
legend('true','bss','bssi','location','northeastoutside');
subplot(2,1,2);
plot(xx1,[yn1-y1 yn1-yi1]);
axis tight;
title('PDF error:(1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)')
legend('bss','bssi','location','northeastoutside');
fprintf('\n');

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

%
%-----------------------------------------------------------------------
% One dimensional test:
%
%    ICDF: Mixture of 3 normal distributions with 2 gaps and zero values
%        between them.
%
%        (1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)
%
%-----------------------------------------------------------------------
%

fprintf('ICDF: (1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)');
rng('default');
neval   = 1001;
nsample = 500;
ndata   = 3*nsample;
fprintf(', %d data points\n',ndata);
xtest1 = [0.25*randn(nsample,1); 5+0.5*randn(nsample,1); 11+0.25*randn(nsample,1)];
xxp_1   = linspace(0,1,neval)';

fprintf('bss......... ')
tic;
bsstest1 = bsspdfest(xtest1);
yicdf_1   = bsseval(xxp_1,bsstest1,'function','icdf');
toc;

% fprintf('bssci....... ')
% tic;
% [~,bsstesti_1]   = bsseval([],bsstest1,'function','icdf');
% toc;
% 
% fprintf('bssi........ ')
% tic;
% yicdfi_1   = bsseval(xxp_1,bsstesti_1,'function','icdf');
% toc;

fprintf('normpdf..... ');
tic;
xx1 = linspace(-3,14,2501)';
yn1 = (1/3)*normpdftestbss(xx1,0,0.25) + (1/3)*normpdftestbss(xx1,5,0.5) + (1/3)*normpdftestbss(xx1,11,0.25);
yncdf_1 = cumtrapz(xx1,yn1);
yncdf_1(yncdf_1>1)=1;
idx0 = find(yncdf_1<sqrt(eps),1,'last');
idx1 = find(abs(1-yncdf_1)<sqrt(eps),1,'first');
while ( abs(yncdf_1(idx0+1)-yncdf_1(idx0)) < sqrt(eps) )
    idx0 = idx0+1;
end
while ( abs(yncdf_1(idx1)-yncdf_1(idx1-1)) < sqrt(eps) )
    idx1 = idx1-1;
end
idxbad   = find(diff(yncdf_1)==0);
idxbad   = idxbad((idxbad>=idx0 & idxbad <=idx1));
idxgood  = (idx0:idx1)';
idxgood  = setdiff(idxgood,idxbad);
ynicdf_1 = interp1(yncdf_1(idxgood),xx1(idxgood),xxp_1,'linear','extrap');
toc;

figure;
subplot(2,1,1);
%plot(xxp_1,[ynicdf_1 yicdf_1 yicdfi_1]);
plot(xxp_1,[ynicdf_1 yicdf_1]);
axis tight;
title('ICDF:(1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)')
%legend('true','bss','bssi','location','northeastoutside');
legend('true','bss','location','northeastoutside');
subplot(2,1,2);
%plot(xxp_1,[ynicdf_1-yicdf_1 ynicdf_1-yicdfi_1]);
plot(xxp_1,ynicdf_1-yicdf_1);
axis tight;
title('ICDF error:(1/3)*N(0,0.25) + (1/3)*(5,0.5) + (1/3)*N(11,0.25)')
%legend('bss','bssi','location','northeastoutside');
legend('bss','location','northeastoutside');
fprintf('\n');

%
%-----------------------------------------------------------------------
% Two dimensional test:
%
%    PDF: Mixture of 5 two dimensional normal distributions
%
%
%-----------------------------------------------------------------------
%
fprintf('PDF: Mixture of 5 2-d normal distributions');
rng('default');
neval    = 121;
nsample  = 5000;
ndata    = 5*nsample;
fprintf(', %d data points',ndata);
fprintf(', %d per component\n',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  = chol(xcov0);
chxcov3  = chol(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;
    ];
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(:)];

fprintf('bss......... ')
tic;
bss2test = bsspdfest(xtest2);
y2mg     = bsseval(xx2mg,bss2test,'function','pdf');
y2    = reshape(y2mg,size(xx2mg1));
toc;

fprintf('bssci....... ')
tic;
[~,bss2testi] = bsseval([],bss2test,'function','pdf','createinterpolant');
toc;

fprintf('bssi........ ')
tic;
y2mg = bsseval(xx2mg,bss2testi,'function','pdf');
y2pdfi  = reshape(y2mg,size(xx2mg1));
toc;

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

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

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2)
axis tight;
title('Estimated PDF: bss');
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: bss');
xlabel('X1');
ylabel('X2');

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

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2pdfi)
axis tight;
title('Estimated PDF: bssi');
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-y2pdfi)
axis tight;
title('PDF Error: bssi');
xlabel('X1');
ylabel('X2');

figure;
surf(xx2mg1,xx2mg2,y2n-y2pdfi)
axis tight;
title('PDF Error: bssi');
xlabel('X1');
ylabel('X2');

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

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

%
%-----------------------------------------------------------------------
% Two dimensional test:
%
%    CDF: Mixture of 5 two dimensional normal distributions
%
%
%-----------------------------------------------------------------------
%
fprintf('CDF: Mixture of 5 2-d normal distributions');
rng('default');
neval    = 121;
nsample  = 5000;
ndata    = 5*nsample;
fprintf(', %d data points',ndata);
fprintf(', %d per component\n',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  = chol(xcov0);
chxcov3  = chol(xcov3);
evalbnds = [-9 9;-7 7];

%
% Generate the data and estimate cdf
%
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;
    ];
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(:)];

fprintf('bss......... ')
tic;
bss2test = bsspdfest(xtest2);
y2mg     = bsseval(xx2mg,bss2test,'function','cdf');
y2cdf    = reshape(y2mg,size(xx2mg1));
toc;

fprintf('bssci....... ')
tic;
[~,bss2testi] = bsseval([],bss2test,'function','cdf','createinterpolant');
toc;

fprintf('bssi........ ')
tic;
y2mg   = bsseval(xx2mg,bss2testi,'function','cdf');
y2cdfi = reshape(y2mg,size(xx2mg1));
toc;

%
% Compute the true density function
%
fprintf('mvncdf...... ');
tic;
y2nmg        = zeros(length(xx2mg),1);
for i = 1:5
   if (i == 1)
      y2nmg = y2nmg + (1/5)*mvnpdftestbss(xx2mg,xmu(i,:),xcov0);
   else
      y2nmg = y2nmg + (1/5)*mvnpdftestbss(xx2mg,xmu(i,:),xcov3);
   end
end
y2n = reshape(y2nmg,size(xx2mg1));
y2n = prod(diff(xx2(1:2,:)))*cumsum(cumsum(y2n,2));
toc;
fprintf('\n');

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

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2cdf)
axis tight;
title('Estimated CDF: bss');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2cdf)
axis tight;
title('CDF Error: bss');
xlabel('X1');
ylabel('X2');

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

figure;
subplot(3,1,1)
surf(xx2mg1,xx2mg2,y2cdfi)
axis tight;
title('Estimated CDF: bssi');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,2)
surf(xx2mg1,xx2mg2,y2n)
axis tight;
title('True CDF');
xlabel('X1');
ylabel('X2');
view([-14 52])
subplot(3,1,3)
surf(xx2mg1,xx2mg2,y2n-y2cdfi)
axis tight;
title('CDF Error: bssi');
xlabel('X1');
ylabel('X2');

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

%
%-----------------------------------------------------------------------
% Three dimensional test:
%
%    PDF: Standard three dimensional normal distribution
%
%
%-----------------------------------------------------------------------
%

%
% Generate data from a standard 3D normal distribution and estimate
% the probability density function using the 3D functions:
% bsspdfest and bsseval.
%
% N(0,I)
%

fprintf('PDF: Standard 3-D normal distribution');
rng('default');
neval    = 121;
sliceidx = 1+floor(neval/2);
nsample  = 25000;
ndata    = nsample;
fprintf(', %d data points\n',ndata);
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);

fprintf('bss......... ')
tic;
bss3test         = bsspdfest(xtest3);
y3mg             = bsseval(xx3mg,bss3test,'function','pdf');
y3               = reshape(y3mg,size(xx3mg1));
toc;

fprintf('bssci....... ')
tic;
[~,bss3testi] = bsseval([],bss3test,'function','pdf','createinterpolant');
toc;

fprintf('bssi........ ')
tic;
y3mg   = bsseval(xx3mg,bss3testi,'function','pdf');
y3pdfi = reshape(y3mg,size(xx3mg1));
toc;

fprintf('mvnpdf...... ');
tic;
y3nmg = mvnpdftestbss(xx3mg,[0 0 0],eye(3));
y3n   = reshape(y3nmg,size(xx3mg1));
toc;

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('\n');
fprintf('Test bsspdfest and bsseval using 3D data:\n');
fprintf('   Distribution                 = N(0,I)\n');
fprintf('   Evaluation interval, x1      = (%f,%f)\n',evalbnds(1,:));
fprintf('   Evaluation interval, x2      = (%f,%f)\n',evalbnds(2,:));
fprintf('   Evaluation interval, x3      = (%f,%f)\n',evalbnds(3,:));
fprintf('   Number of data points        = %d\n',ndata);
fprintf('   Number of evaluation points  = %d\n',length(xx3mg));
fprintf('   Area under the true pdf      = %f\n',sum(y3n(:))*prod(h));
fprintf('   Estimator                    = %12s %12s\n','bsspdfest','mvksdensity');
fprintf('   Area under the estimated pdf = %12.6f %12.6f\n',[sum(y3(:))*prod(h) nan]);
fprintf('   Integrated absolute error    = %12.6f %12.6f\n',[sum(abs(y3n(:)-y3(:)))*prod(h) nan]);
fprintf('   RMSE                         = %12.6f %12.6f\n',[sqrt(sum((y3n(:)-y3(:)).^2)/length(y3n)) nan]);
fprintf('   IRMSE                        = %12.6f %12.6f\n',[sqrt(sum((y3n(:)-y3(:)).^2)*prod(h)) nan]);
fprintf('\n');

%
%-----------------------------------------------------------------------
% N dimensional test:
%
%    PDF: Standard N (four) dimensional normal distribution
%
%
%-----------------------------------------------------------------------
%

%
% Generate data from a standard 4D normal distribution and estimate
% the probability density function using the functions:
% bsspdfest and bsseval.
%
% N(0,I)
%

fprintf('PDF: Standard N-D (4-D) normal distribution');
rng('default');
neval            = 33;
sliceidx = 1+floor(neval/2);
nsample          = 50000;
ndata            = nsample;
fprintf(', %d data points\n',ndata);
fprintf('   This test can take over an hour to complete.\n');
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);

fprintf('bss......... ')
tic;
bss4test         = bsspdfest(xtest4,'partitionsize',21);
y4mg             = bsseval(xx4mg,bss4test,'function','pdf');
toc;
y4               = reshape(y4mg,size(xx4mg1));

fprintf('bssci....... ')
tic;
[~,bss4testi] = bsseval([],bss4test,'function','pdf','createinterpolant');
toc;

fprintf('bssi........ ')
tic;
y4mg   = bsseval(xx4mg,bss4testi,'function','pdf');
y4pdfi = reshape(y4mg,size(xx4mg1));
toc;

fprintf('mvnpdf...... ');
tic;
y4nmg = mvnpdftestbss(xx4mg,[0 0 0 0],eye(4));
y4n   = reshape(y4nmg,size(xx4mg1));
toc;

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)';

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

diary off
