disp('Running MATLAB script file gpkdet.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION gpkde,
%    General Purpose Kernel Density Estimation

itest = 9 ;     %  1,2,3,4,5,6,7,8,9

format compact ;
format short ;


if itest == 1 ;   %  Simplest visual check of implementations
  indat = [1; 2; 2; 3; 3.5] ;
  invh = .2 ;
  vxg = [0; 4] ;
  titstr = 'matrix implementation'
  subplot(2,2,1) ;
    tic ;
    gpkde(indat,invh,vxg,1) ;
    toc
      title(titstr) ;
  titstr = 'looped implementation'
  subplot(2,2,2) ;
    tic ;
    gpkde(indat,invh,vxg,2) ;
    toc
      title(titstr) ;
  titstr = 'binned implementation' 
  subplot(2,2,3) ;
    tic ;
    gpkde(indat,invh,vxg,0) ;
    toc
      title(titstr) ;
 
elseif itest == 2 ;   %  Check of numerical integrals
  format long ;
  indat = 2 * rand(100,1) ;
  invh = .2 ;
  vxg = [-1, 3, 9] ;
  disp('Check these are 1 (9 grid points):') ;
    disp('    Matrix:') ;  
    vkde = gpkde(indat,invh,vxg,1) ;
  int = sum(vkde * .5) 
    disp('    Looped:') ;  
    vkde = gpkde(indat,invh,vxg,2) ;
  int = sum(vkde * .5) 
    disp('    Binned:') ;  
    vkde = gpkde(indat,invh,vxg,0) ;
  int = sum(vkde * .5) 

  vxg = [-1, 3, 17] ;
  disp('Check these are 1 (17 grid points):') ;
    disp('    Matrix:') ;  
    vkde = gpkde(indat,invh,vxg,1) ;
  int = sum(vkde * .25) 
    disp('    Looped:') ;  
    vkde = gpkde(indat,invh,vxg,2) ;
  int = sum(vkde * .25) 
    disp('    Binned:') ;  
    vkde = gpkde(indat,invh,vxg,0) ;
  int = sum(vkde * .25) 

  vxg = [-1, 3] ;
  disp('Check these are 1 (401 grid points):') ;
    disp('    Matrix:') ;  
    vkde = gpkde(indat,invh,vxg,1) ;
  int = sum(vkde * .01) 
    disp('    Looped:') ;  
    vkde = gpkde(indat,invh,vxg,2) ;
  int = sum(vkde * .01) 
    disp('    Binned:') ;  
    vkde = gpkde(indat,invh,vxg,0) ;
  int = sum(vkde * .01) 

elseif itest == 3 ;   %  Check multiple bandwidths
  indat = [(randn(25,1) - 1.5); (randn(25,1) + 1.5)] ;
  invh = [.2 .6 1.8] ;
  vxg = [-4; 4] ;
  titstr = 'matrix implementation'
  subplot(2,2,1) ;
    tic ;
    gpkde(indat,invh,vxg,1) ;
    toc
      title(titstr) ;
  titstr = 'looped implementation'
  subplot(2,2,2) ;
    tic ;
    gpkde(indat,invh,vxg,2) ;
    toc
      title(titstr) ;
  titstr = 'binned implementation' 
  subplot(2,2,3) ;
    tic ;
    gpkde(indat,invh,vxg,0) ;
    toc
      title(titstr) ;
  titstr = 'binned impl., all defaults' 
  subplot(2,2,4) ;
    tic ;
    gpkde(indat,invh) ;
    toc
      title(titstr) ;
  
elseif itest == 4 ;   %  Check implementations against each other
  format long ;
  indat = randn(100,1) ;
  invh = [.2 .6 1.8] ;
  vxg = [-3, 3, 13] ;
  disp('Check these are 0 (13 grid points):') ;
  mkde = gpkde(indat,invh,vxg,1) ;
  lkde = gpkde(indat,invh,vxg,2) ;
  bkde = gpkde(indat,invh,vxg,0) ;
  disp('    Matrix - Looped:') ;  
  max(abs(mkde - lkde)) 
  disp('    Matrix - Binned:') ;  
  max(abs(mkde - bkde)) 

  vxg = [-3, 3] ;
  disp('Check these are 0 (401 grid points):') ;
  mkde = gpkde(indat,invh,vxg,1) ;
  lkde = gpkde(indat,invh,vxg,2) ;
  bkde = gpkde(indat,invh,vxg,0) ;
  disp('    Matrix - Looped:') ;  
  max(abs(mkde - lkde)) 
  disp('    Matrix - Binned:') ;  
  max(abs(mkde - bkde)) 

elseif itest == 5 ;   %  Check grid choice defaults
  indat = [(randn(250,1) - 1.5); (randn(250,1) + 1.5)] ;
  invh = [.2 .6 1.8] ;
  subplot(2,2,1) ;
    gpkde(indat,invh) ;
  subplot(2,2,2) ;
    gpkde(indat,invh,[-5,5]) ;
  subplot(2,2,3) ;
    gpkde(indat,invh,[-5,5,401]) ;
  subplot(2,2,4) ;
    gpkde(indat,invh,[-5,5,21]) ;

elseif itest == 6 ;   %  Check endpt handling
  indat = randn(1000,1) ;
  invh = [.1 .3 .9] ;
  subplot(2,2,1) ;
    gpkde(indat,invh,[-1,1],0) ;
  subplot(2,2,2) ;
    gpkde(indat,invh,[-1,1,17],0,0) ;
  subplot(2,2,3) ;
    gpkde(indat,invh,[-1,1],0,1) ;
  subplot(2,2,4) ;
    gpkde(indat,invh,[-1,1,17],0,1) ;

elseif itest == 7 ;   %  Check already binned data
  format long ;
  indat = randn(1000,1) ;
  invh = [.1 .3 .9] ;
  kde1 = gpkde(indat,invh,0,0) ;
  bcts = gplbinr(indat) ;
  kde2 = gpkde(bcts,invh,[min(indat) max(indat)],-1) ;
  disp('Check this is 0:') ;
  max(max(abs(kde1 - kde2))) 

elseif itest == 8 ;   %  Check data based bandwidth selection
  indat = [(randn(250,1) - 1.5); (randn(250,1) + 1.5)] ;
  subplot(2,2,1) ;
    gpkde(indat,-1) ;
      title('Simple Normal Reference') ;
  subplot(2,2,2) ;
    gpkde(indat,-2) ;
      title('ROT2') ;
  subplot(2,2,3) ;
    gpkde(indat,-3) ;
      title('Oversmoothed') ;
  subplot(2,2,4) ;
    gpkde(indat) ;
      title('SJPI') ;
    
elseif itest == 9 ;   %  Check kernel plot
  indat = [(randn(250,1) - 1.5); (randn(250,1) + 1.5)] ;

  [okde,oxgrid,omker] = gpkde(indat,[.2; .6; 1.8],[-5; 5],0,0) ;
  subplot(2,1,1) ;
    plot(oxgrid,[okde omker]) ;

  [okde,oxgrid,omker] = gpkde(indat,.6) ;
  subplot(2,1,2) ;
    plot(oxgrid,[okde omker]) ;

end ;
% gpkde(data,vh,vxgrid,imptyp,eptflag) 


