disp('Running MATLAB script file bwsjpibt.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION bwsjpib,
%    Sheather Jones Plug In bandwidth (binned)
%    Compares results to those from GAUSS program NMSJPI.TST

itest = 4 ;   %  1,2  (test against old GAUSS implementation, etc.)
              %  3    (look carefully at default h range, for real data)            
              %  4    (look carefully at default h range, for nm data)            

%    Last two tests suggest that bandwidth range is an adequate default,
%    since only time was out of range was for:
%        Ray's dust data (VERY challenging, since so "discrete")
%        #16 Distant Bimodal, n = 1000,10000 (very challenging
%              since "sd" is so noninformative).

format compact ;


if itest <= 2 ;
  format long ;

  %  Read in test data:
  load d:\gauss\steve\nmstuff\tests\nmsjpi.dat ;
          %  assumes this ascii file is in the current directory
          %  (usually lives in \gauss\steve\nmstuff\tests)
  data = nmsjpi ;
          %  recall thsi load uses the filename for the variable


  %  Old value, copied from NMSJPI.OUT
  osjpi = 0.20792378354257 ;


  %  Get new value
  if itest == 1 ;       %  Just do quick default SJPI
    nsjpi = bwsjpib(data) ;
    %  Also throw in the binned outside version:
    bcts = gplbinr(data,0) ;
    pbsjpi = bwsjpib(bcts,[min(data),max(data)],0,-1) ;
  elseif itest == 2 ;   %  Pay more attention to getting xgrid right
    nsjpi = bwsjpib(data,[-4,4,401],0,0,0) ;
    %  Also throw in the binned outside version:
    bcts = gplbinr(data,[-4,4,401]) ;
    pbsjpi = bwsjpib(bcts,[-4,4,401],0,-1) ;
  end ;


  %  Compare:
  disp('Old SJPI was:') ;
  osjpi

  disp('New SJPI was:') ;
  nsjpi

  disp('Prebinned SJPI was:') ;
  nsjpi


elseif itest == 3 ;    %  Then do real data examples

  idat = 0 ;     %  0 - to loop through all
                 %  1 - Income Data
                 %  2 - Geyser Data
                 %  3 - Stamp Data
                 %  4 - Chondrite Data
                 %  5 - Buffalo Snowfall Data
                 %  6 - Suicide Data
                 %  7 - Coal Seam Data
                 %  8 - Marine Data
                 %  9 - Dust Data
                 %  10 - Galactic Cluster Data

  if idat == 0 ;
    idats = 1 ;
    idate = 10 ;
  else ;
    idats = idat ;
    idate = idat ;
  end ;

  for iidat = idats:idate ;

    if iidat == 1 ;
      dtitstr = 'Incomes' ;
      infstr = 'inc75.all' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 2 ;
      dtitstr = 'Geyser' ;
      infstr = 'geyser.dat' ;
      iincol = [2; 1] ;    %  data come as 2 columns, want col. 1
    elseif iidat == 3 ;
      dtitstr = 'Stamps' ;
      infstr = 'stamp' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 4 ;
      dtitstr = 'Chondrite' ;
      infstr = 'chrondit.asc' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 5 ;
      dtitstr = 'Snowfalls' ;
      infstr = 'buffalo.sno' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 6 ;
      dtitstr = 'Suicides' ;
      infstr = 'suicide.dat' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 7 ;
      dtitstr = 'Coal Seam' ;
      infstr = 'fishcoal.iva' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 8 ;
      dtitstr = 'Marine' ;
      infstr = 'freddat.asc' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 9 ;
      dtitstr = 'Dust' ;
      infstr = 'raydust.dat' ;
      iincol = 1 ;    %  want all data as 1 column
    elseif iidat == 10 ;
      dtitstr = 'Galaxies' ;
      infstr = 'galaxy.dat' ;
      iincol = 1 ;    %  want all data as 1 column
    end ;

    disp(['    Loading ' dtitstr ' Data']) ;
      infstr = ['d:\gauss\steve\data\' infstr] ;
      fid = fopen(infstr,'r') ;
        if iincol(1) == 1 ;
          data = fscanf(fid,'%g') ;
        elseif iincol(1) == 2 ;
          data = fscanf(fid,'%g') ;
          data = reshape(data,2,length(data)/2)' ;
                   %  since reshape fills down cols first
          data = data(:,iincol(2)) ;
        end ;
      fclose(fid) ;
%    disp(['    Read in ' num2str(length(data)) ' data points']) ;

    hsjpi = bwsjpib(data) ;
    gpkde(data,hsjpi) ;
      title([dtitstr ' Data, h_SJPI = ' num2str(hsjpi)]) ;

    disp(['Any key to continue']) ;
    pause ;


  end ;

elseif itest == 4 ;    %  Then do normal mixture examples

  %  Load Normal mixtures parameters from nnpar.mat
  %  !!!   Careful:  file nmpar.mat needs to be in the current directory  !!!
  %  Then all parameter matrices are in current directory
  load nmpar.mat
  
  %  Loop through distributions
  idiste = 16 ;
  for idist = 1:idiste ;
    eval(['parmat = nmpar' num2str(idist) ' ;']) ;
    eval(['titstr = nmtit' num2str(idist) ' ;']) ;
    disp(['    Working on ' titstr]) ;

    clf ;
    xgrid = linspace(-3,3,401)' ;
    f = nmf(xgrid,parmat) ;
    plot(xgrid,f,'r') ;
      title(['#' num2str(idist) ' ' titstr]) ;

    hold on ;
    for inobs = 2:4 ;
      nobs = 10^inobs ;
      disp(['      Working on n = ' num2str(nobs)]) ;

      %  Generate psuedo data and plot 
      data = nmdata(nobs,parmat) ;
      hsjpi = bwsjpib(data) ;
      gpkde(data,hsjpi) ;
    end ;
    hold off ;

    disp('Any key to continue') ;
    pause

  end ;


end ;


