disp('Running MATLAB script file bwrswbt.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION bwrswb,
%    Ruppert Sheather Wand bandwidth (binned, for regression)

itest = 2 ;   %  1    (simple test)
               %  1.1  (simple test with linear data)
               %  2    (look at real data examples)            
               %  3    (NR book examples)
               %  4    (B-spline test examples)
               %  10   Do detailed output

% start new stuff here

format compact ;


if itest == 1 ;

  indat = rand(500,1) .^ 2 ;
  indat = [indat, (indat - .4).^4 + .01 * randn(500,1)] ;

  qrswrot = bwrswb(indat,-1) 
  subplot(2,2,1) ;
    gpnpr(indat,qrswrot) ;
      title('RSW - ROT, quick defaults') ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.w') ;
      hold off ;

  arswrot = bwrswb(indat,-1,[0,1,401],0) 
  subplot(2,2,2) ;
    gpnpr(indat,arswrot) ;
      title('RSW - ROT, args') ;

  qrswdpi = bwrswb(indat) 
  subplot(2,2,3) ;
    gpnpr(indat,qrswdpi) ;
      title('RSW - DPI, quick defaults') ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.w') ;
      hold off ;

  arswdpi = bwrswb(indat,0,[0,1,401],0) 
  subplot(2,2,4) ;
    gpnpr(indat,arswdpi) ;
      title('RSW - DPI, args') ;

  disp('See how close these are:') ;
  dif1 = abs(qrswrot - arswrot) 
  dif2 = abs(qrswdpi - arswdpi) 

  disp('See how close these are:') ;
  dif3 = abs(qrswrot - qrswdpi) 

  %  output stuff to file for emailing
  fid = fopen('/matlab/steve/bwrswbt.out','wt') ;

  titlstr1 = 'Output from the MATLAB Script bwrswbt.m,   ' ;
cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;

  forstr = '   h_ROT = %12.10f\n\n' ;
cntbytes = fprintf(fid,forstr,arswrot) ;

  forstr = '   h_DPI = %12.10f\n\n' ;
cntbytes = fprintf(fid,forstr,arswdpi) ;

  hdrstr = 'For the raw data:   ' ;
cntbytes = fprintf(fid,'%1s\n',hdrstr) ;

  forstr = '  %12.10f  %12.10f\n' ;
cntbytes = fprintf(fid,forstr,indat') ;

  fclose('all') ;


elseif itest == 1.1 ;    %  Then do linear data example

  indat = 1:6 ;
  indat = [indat', indat'] ;

  hout = bwrswb(indat,0)


elseif itest == 2 ;    %  Then do real indat examples

  idat = 0 ;     %  0 - to loop through all
                 %  1 - Canadian Earning Data
                 %  2 - Canadian Lynx Data
                 %  3 - Cars Data
                 %  4 - Motorcycle Data
                 %  5 - NIcFoo Data
                 %  6 - Swamp Data
                 %  7 - Sunspot Data

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

  for iidat = idats:idate ;

    if iidat == 1 ;
      dtitstr = 'Earning Power' ;
           %  Data from Aman Ullah, from Canadian survey of 
           %  earning power as a function of age.  
           %  Used in Chu and Marron Stat. Sci. paper.
      infstr = 'canearn.dat' ;
      frmtstr = '%g %g' ;
      insize = [2, inf] ;
      vicol = [1, 2] ;
      datatran = [] ;
      xstr = 'Age' ;
      ystr = 'Income' ;
    elseif iidat == 2 ;
      dtitstr = 'Canadian Lynx' ;
           %  Famous time series data, from lynx trappings
      infstr = 'canlynx.dat' ;
      frmtstr = '%g' ;
      insize = [1,inf] ;
      vicol = [1] ;
      datatran = 'data(:,2) = log10(data(:,2) + 1) ;' ;
      xstr = 'Year' ;
      ystr = 'log10(# Lynx trapped + 1)' ;
    elseif iidat == 3 ;
      dtitstr = 'Cars' ;
           %  Subset of famous cars data, found by Mat Wand as
           %  good example of need for monotonic smoothing
      infstr = 'cars.mat' ;
      frmtstr = '%g %g' ;
      insize = [2, inf] ;
      vicol = [1, 2] ;
      datatran = [] ;
      xstr = 'Displacement' ;
      ystr = 'Mileage' ;
    elseif iidat == 4 ;
      dtitstr = 'Motorcycle' ;
           %  Full version of Haerdle's motorcycle crash data
      infstr = 'motodata.txt' ;
      frmtstr = '%g %g %s' ;
      insize = [5, inf] ;
      vicol = [1, 2] ;
      datatran = 'data = data(1:800,:) ;' ;
         %  First 800 are first trace
      xstr = 'Time' ;
      ystr = 'Acceleration' ;
    elseif iidat == 5 ;
      dtitstr = 'Food Expenditures' ;
      infstr = 'nicfoo.dat' ;
      frmtstr = '%g %g' ;
      insize = [2, inf] ;
      vicol = [1, 2] ;
      datatran = [] ;
      xstr = 'Net Income' ;
      ystr = 'Exp. for Food' ;
    elseif iidat == 6 ;
      dtitstr = 'Swamp Phosphorus' ;
      infstr = 's7xform.sho' ;
      frmtstr = '%g %g' ;
      insize = [2, inf] ;
      vicol = [1, 2] ;
%above lines are for data with "big gap", 
%which causes problem for hrsw
%following is alternate version, without this problem
%(but uninteresting data)
%      infstr = 's7ts.dat' ;
%      frmtstr = '%g %g %g %g %g %g' ;
%      insize = [6, inf] ;
%      vicol = [4, 5] ;
      datatran = 'data(:,2) = log10(data(:,2)) ;' ;
      xstr = 'Day' ;
      ystr = 'log10(Phosphorus Conc.)' ;
    elseif iidat == 7 ;
      dtitstr = 'Sunspots' ;
      infstr = 'sunspot.dat' ;
      frmtstr = '%g' ;
      insize = [1, inf] ;
      vicol = [1] ;
      datatran = 'data(:,2) = sqrt(data(:,2)) ;' ;
      xstr = 'Year' ;
      ystr = 'sqrt(# Sunspots)' ;
    end ;


    disp(['    Loading ' dtitstr ' Data']) ;
      infstr = ['d:\gauss\steve\data\' infstr] ;
      fid = fopen(infstr,'r') ;
        data = fscanf(fid,frmtstr,insize) ;
        data = data' ;
          %  since data rows are read in as columns
        data = data(:,vicol) ;
        if size(data,2) == 1 ;  %  then have only one column
                                %  so add dummy column for x's
          data = [(1:length(data))', data] ;
        end ;
      fclose(fid) ;
    disp(['    Read in ' num2str(size(data,1)) ' data points']) ;


    if length(datatran) ~= 0 ;    %  Then do some data transformation
      eval(datatran) ;
    end ;


    hrot = bwrswb(data,-1) ;
    hdpi = bwrswb(data) ;
    gpnpr(data,hdpi) ;
      title([dtitstr ' Data, h_ROT = ' num2str(hrot) ...
                                ', h_DPI = ' num2str(hdpi)]) ;
      xlabel(xstr) ;
      ylabel(ystr) ;
      hold on ;
        [npr, xgrid] = gpnpr(data,hrot) ;
        plot(xgrid,npr,'--r') ;
        plot(data(:,1),data(:,2),'.w') ;
      hold off ;

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


  end ;

elseif itest == 3 ;    %  Then do NR book examples

  iset = 0 ;    %  0 - Loop through all
                %  1 - Range of Noise Levels
                %  2 - Shifting Design Points
                %  3 - Heteroscedasticity
                %  4 - (Not ready yet)
                %  5 - 2 extreme data sets


  if iset == 1 | iset == 0 ;
    disp('Running example 1') ;

    fid = fopen('d:\gauss\steve\nrbook\nrdat1.dat','r') ;
      data = fscanf(fid,'%g %g %g %g %g %g %g %g',[8, inf]) ;
      data = data' ;
          %  since data rows are read in as columns
    fclose(fid) ;

    clf ;
    mest = [] ;
    for ieg = 1:6 ;
      hdpi = bwrswb([data(:,1), data(:,ieg+2)]) ;
      [est, xgrid] = gpnpr([data(:,1), data(:,ieg+2)],hdpi) ;
      subplot(2,3,ieg) ;
        plot(xgrid,est) ;
          title(['h_DPI = ' num2str(hdpi)]) ;
          hold on ;
            plot(data(:,1),data(:,ieg+2),'.w') ;
            plot(data(:,1),data(:,2),'--b') ;
          hold off ;
      mest = [mest, est] ;
    end ;
    disp('Any key to continue') ;
    pause ;

    clf ;
    plot(xgrid,mest) ;
      title('NR Book e.g. 1') ;

    disp('Any key to continue') ;
    pause ;
  end ;


  if iset == 2 | iset == 0 ;
    disp('Running example 2') ;

    fid = fopen('d:\gauss\steve\nrbook\nrdat2.dat','r') ;
      data = fscanf(fid,'%g %g %g %g %g %g %g %g',[8, inf]) ;
      data = data' ;
          %  since data rows are read in as columns
    fclose(fid) ;

    clf ;
    mest = [] ;
    for ieg = 1:6 ;
      hdpi = bwrswb([data(:,ieg), data(:,8)]) ;
      [est, xgrid] = gpnpr([data(:,ieg), data(:,8)],hdpi) ;
      subplot(2,3,ieg) ;
        plot(xgrid,est) ;
          title(['h_DPI = ' num2str(hdpi)]) ;
          hold on ;
            plot(data(:,ieg),data(:,8),'.w') ;
            plot(data(:,ieg),data(:,7),'--b') ;
          hold off ;
      mest = [mest, est] ;
    end ;
    disp('Any key to continue') ;
    pause ;

    clf ;
    plot(xgrid,mest) ;
      title('NR Book e.g. 2') ;

    disp('Any key to continue') ;
    pause ;
  end ;

  if iset == 3 | iset == 0 ;
    disp('Running example 3') ;

    fid = fopen('d:\gauss\steve\nrbook\nrdat3.dat','r') ;
      data = fscanf(fid,'%g %g %g %g %g %g %g',[7, inf]) ;
      data = data' ;
          %  since data rows are read in as columns
    fclose(fid) ;

    clf ;
    mest = [] ;
    for ieg = 1:5 ;
      hdpi = bwrswb([data(:,1), data(:,ieg+2)]) ;
      [est, xgrid] = gpnpr([data(:,1), data(:,ieg+2)],hdpi) ;
      subplot(2,3,ieg) ;
        plot(xgrid,est) ;
          title(['h_DPI = ' num2str(hdpi)]) ;
          hold on ;
            plot(data(:,1),data(:,ieg+2),'.w') ;
            plot(data(:,1),data(:,2),'--b') ;
          hold off ;
      mest = [mest, est] ;
    end ;
    disp('Any key to continue') ;
    pause ;

    clf ;
    plot(xgrid,mest) ;
      title('NR Book e.g. 3') ;

    disp('Any key to continue') ;
    pause ;
  end ;

  if iset == 4 | iset == 0 ;
    disp('Example 4 was not finished') ;
  end ;

  if iset == 5 | iset == 0 ;
    disp('Running example 5') ;

    fid = fopen('d:\gauss\steve\nrbook\nrdat5a.dat','r') ;
      data = fscanf(fid,'%g %g %g',[3, inf]) ;
      data = data' ;
          %  since data rows are read in as columns
    fclose(fid) ;

    hdpi = bwrswb([data(:,1), data(:,3)]) ;
    gpnpr([data(:,1), data(:,3)],hdpi) ;
      title(['NR Book, e.g. 5a, h_DPI = ' num2str(hdpi)]) ;
          hold on ;
            plot(data(:,1),data(:,3),'.w') ;
            plot(data(:,1),data(:,2),'--b') ;
          hold off ;
    disp('Any key to continue') ;
    pause ;

    fid = fopen('d:\gauss\steve\nrbook\nrdat5b.dat','r') ;
      data = fscanf(fid,'%g %g %g',[3, inf]) ;
      data = data' ;
          %  since data rows are read in as columns
    fclose(fid) ;

    hdpi = bwrswb([data(:,1), data(:,3)]) ;
    gpnpr([data(:,1), data(:,3)],hdpi) ;
      title(['NR Book, e.g. 5b, h_DPI = ' num2str(hdpi)]) ;
          hold on ;
            plot(data(:,1),data(:,3),'.w') ;
            plot(data(:,1),data(:,2),'--b') ;
          hold off ;

  end ;

elseif itest == 4 ;    %  Then do B-spline examples

  for ifile = 1:9 ;

    infstr = ['d:\gauss\steve\bspline\bs3' num2str(ifile) '.dat'] ;
    fid = fopen(infstr,'r') ;
      data = fscanf(fid,'%g %g %g',[3, inf]) ;
      data = data' ;
          %  since data rows are read in as columns
    fclose(fid) ;

    hdpi = bwrswb([data(:,1), data(:,3)]) ;
    subplot(3,3,ifile) ;
      gpnpr([data(:,1), data(:,3)],hdpi) ;
        title(['h_DPI = ' num2str(hdpi)]) ;
        hold on ;
          plot(data(:,1),data(:,3),'.w') ;
          plot(data(:,1),data(:,2),'--b') ;
        hold off ;
        axis([0,1,0,1]) ;
          %  Comment this last out to see all of data
  end ;


elseif itest == 10 ;     %  work with modified version for detailed output
                     %  using "cars" data

    dtitstr = 'Cars' ;
           %  Subset of famous cars data, found by Matt Wand as
           %  good example of need for monotonic smoothing
    infstr = 'cars.mat' ;
    frmtstr = '%g %g' ;
    insize = [2, inf] ;
    vicol = [1, 2] ;
    datatran = [] ;
    xstr = 'Displacement' ;
    ystr = 'Mileage' ;


  disp(['    Loading ' dtitstr ' Data']) ;
    infstr = ['d:\gauss\steve\data\' infstr] ;
    fid = fopen(infstr,'r') ;
      data = fscanf(fid,frmtstr,insize) ;
      data = data' ;
          %  since data rows are read in as columns
      data = data(:,vicol) ;
      if size(data,2) == 1 ;  %  then have only one column
                              %  so add dummy column for x's
        data = [(1:length(data))', data] ;
      end ;
    fclose(fid) ;
  disp(['    Read in ' num2str(size(data,1)) ' data points']) ;

  diary bwrswbt.ou3 ;

  comment = 'Intermediate values for the RSW bandwidth, and the cars data:'

  hrsw = bwrswb2(data) ;

  hrsw

  data

  diary off ;


end ;


