disp('Running MATLAB script file gpanal2t.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION gpanal2,
%    General Purpose data ANALysis,
%    This does family, SiZer, and SiCon, all on one plot


itest = 1 ;    %  1  real data tests (density estimation)
                %  2  normal mixture tests (density estimation)
                %  11 real data tests (regression)
                %  12 NR book examples (regression)
                %  13 B-spline test examples (regression)
                %  14 ER Wavelet Examples (regression)
                %  21 parameter test 1 (gray color)
                %  22 parameter test 2 (restricted range)
                %  23 parameter test 3 (21 curves)
                %  24 parameter test 4 (print gray .ps file)
                %  25 parameter test 5 (print color .ps file)
                %  26 parameter test 6 (restricted grid size)


format compact ;



if itest == 1 ;    %  then do real data tests

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

    disp(['    Loading ' dtitstr ' Data']) ;
      infstr = ['\matlab\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']) ;

    if length(vrange) == 2 ;
      gpanal2(data,1,vrange) ;
    else ;
      gpanal2(data) ;
    end ;
      title([dtitstr ' Data']) ;

      orient tall ;
    eval(['print -dpsc \matlab\steve\ps\g2a' ...
                               num2str(iidat)  '.ps']) ;

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

  end ;

elseif itest == 2 ;    %  then do normal mixture tests

  %  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 \matlab\steve\nmstuff\nmpar.mat
  
  %  Loop through distributions
  idiste = 16 ;
  for idist = 1:idiste ;
    eval(['parmat = nmpar' num2str(idist) ' ;']) ;
    eval(['titstr = nmtit' num2str(idist) ' ;']) ;
    disp(['    Working on #' num2str(idist) titstr]) ;

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

      %  Generate psuedo data and do family plot 
      data = nmdata(nobs,parmat) ;
      gpanal2(data) ;
      title(['#' num2str(idist) ' ' titstr ', n = ' num2str(nobs)]) ;

        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2b' ...
                                  num2str(idist) num2str(inobs)  '.ps']) ;

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

    end ;

  end ;


elseif itest == 11 ;    %  then do real data regression

  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
                 %  8 - Fossil Data
                 %  9 - Share Yield Data

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

  for iidat = idats:idate ;

    clf ;
    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)' ;
    elseif iidat == 8 ;
      dtitstr = 'Fossils' ;
      %  from Tim Bralower
      infstr = 'fossils.dat' ;
      frmtstr = '%g %g %g' ;
      insize = [3, inf] ;
      vicol = [1, 2] ;
      datatran = [] ;
      xstr = 'Age (mil. years)' ;
      ystr = 'Strontium Ratio' ;
    elseif iidat == 9 ;
      dtitstr = 'Share Yield' ;
      %  from Probal Chaudhuri
      %  Index of Stock Profits on "ordinary shares", Indian Stock Market
      %  on scale of "interest rates"
      %  from Reserve Bank of India Bulletin
      %  time series of monthly data, April 1982 - March 1991
      infstr = 'sharyiel.dat' ;
      frmtstr = '%g' ;
      insize = [1, inf] ;
      vicol = [1] ;
      datatran = [] ;
      xstr = 'month' ;
      ystr = 'Share yield' ;
    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 ;

    gpanal2(data) ;
      title(dtitstr) ;
      xlabel(xstr) ;
      ylabel(ystr) ;

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

      orient tall ;
    eval(['print -dpsc \matlab\steve\ps\g2c' ...
                               num2str(iidat)  '.ps']) ;



  end ;

elseif itest == 12 ;    %  then do NR Book regression 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 ;
    for ieg = 1:6 ;
      disp(['        Working on part ' num2str(ieg)]) ;
      gpanal2([data(:,1), data(:,ieg+2)]) ;
        title('NR Book e.g. 1') ;

        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2d' ...
                                       '1' num2str(ieg)  '.ps']) ;

    end ;

  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 ;
    for ieg = 1:6 ;
      disp(['        Working on part ' num2str(ieg)]) ;
      gpanal2([data(:,ieg), data(:,8)]) ;
        title('NR Book e.g. 2') ;

        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2e' ...
                                       '2' num2str(ieg)  '.ps']) ;
    end ;

  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 ;
    for ieg = 1:5 ;
      disp(['        Working on part ' num2str(ieg)]) ;
      gpanal2([data(:,1), data(:,ieg+2)]) ;
        title('NR Book e.g. 3') ;

        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2f' ...
                                       '3' num2str(ieg)  '.ps']) ;
    end ;

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

    clf ;
    gpanal2([data(:,1), data(:,3)]) ;
      title('NR Book, e.g. 5a') ;

        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2g' ...
                                       '5a.ps']) ;

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

    clf ;
    gpanal2([data(:,1), data(:,3)]) ;
      title('NR Book, e.g. 5b') ;
        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2h' ...
                                       '5b.ps']) ;

  end ;

elseif itest == 13 ;    %  then do B-spline regression examples 

  clf
  for ifile = 1:9 ;
    disp(['        Working on part ' num2str(ifile)]) ;

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

    gpanal2([data(:,1), data(:,3)]) ;
      title('B-spline Targets') ;

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

      orient tall ;
    eval(['print -dpsc \matlab\steve\ps\gp2i' ...
                               num2str(ifile)  '.ps']) ;

  end ;

elseif itest == 14 ;    %  then do ER wavelet regression examples 

  clf ;
  nobs = 1024 ;
  xgrid = linspace(1/(2*nobs), 1 - (1/(2*nobs)), nobs)' ;

  for idist = 1:10 ;

    %  Following lines are modifications of Sudeshna Adak's Gen_sign.m
    if idist == 1 ;
      tarstr = 'Step' ;
      m = 0.2 + 0.6*(xgrid > 1/3 & xgrid <= 0.75); 
    elseif idist == 2 ;
      tarstr = 'Wave' ;
      m = 0.5 + (0.2.*cos(4*pi*xgrid)) + (0.1.*cos(24*pi*xgrid));
    elseif idist == 3 ;
      tarstr = 'Blip' ;
       m = (0.32 + (0.6.*xgrid) + ...
                      0.3*exp(-100*((xgrid-0.3).^2))).*(xgrid <= 0.8) + ...
           (-0.28 + (0.6.*xgrid) + ... 
                      0.3*exp(-100*((xgrid-1.3).^2))).*(xgrid > 0.8);
    elseif idist == 4 ;
      tarstr = 'Blocks' ;
      pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
      hgt = [4 (-5) 3 (-4) 5 (-4.2) 2.1 4.3  (-3.1) 2.1 (-4.2)];
      m = 2*ones(size(xgrid));
      for j=1:length(pos)
        m = m + (1 + sign(xgrid-pos(j))).*(hgt(j)/2) ;
      end
      m = (0.6/9.2)*m + 0.2;
    elseif idist == 5 ;
      tarstr = 'Bumps' ;
      pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
      hgt = [ 4  5   3   4  5  4.2 2.1 4.3  3.1 5.1 4.2];
      wth = [.005 .005 .006 .01 .01 .03 .01 .01  .005 .008 .005];
      m = zeros(size(xgrid));
      for j =1:length(pos)
        m = m + hgt(j)./( 1 + abs((xgrid - pos(j))./wth(j))).^4;
      end 
      m = ((0.6/5.3437952)*m) + 0.2;
    elseif idist == 6 ;
      tarstr = 'HeaviSine' ;
      m = 4.*sin(4*pi.*xgrid) - sign(xgrid - .3) - sign(.72 - xgrid) + 5;
      m = (0.6/9)*m + 0.2;
    elseif idist == 7 ;
      tarstr = 'Doppler' ;
      m = sqrt(xgrid.*(1-xgrid)).*sin((2*pi*1.05) ./(xgrid+.05)) + 0.5;
      m = 0.6*m + 0.2;
    elseif idist == 8 ;
      tarstr = 'Angles' ;
      m = ((2*xgrid + 0.5).*(xgrid <= 0.15)) + ...
        ((-12*(xgrid-0.15) + 0.8).*(xgrid > 0.15 & xgrid <= 0.2)) + ...
        0.2*(xgrid > 0.2 & xgrid <= 0.5) + ...
        ((6*(xgrid - 0.5) + 0.2).*(xgrid > 0.5 & xgrid <= 0.6)) + ...
        ((-10*(xgrid - 0.6) + 0.8).*(xgrid > 0.6 & xgrid <= 0.65)) + ...
        ((-0.5*(xgrid - 0.65) + 0.3).*(xgrid > 0.65 & xgrid <= 0.85)) + ...
        ((2*(xgrid - 0.85) + 0.2).*(xgrid > 0.85));
    elseif idist == 9 ;
      tarstr = 'Parabolas' ;
      pos = [0.1 0.2 0.3 0.35 0.37 0.41 0.43 0.5 0.7 0.9];
      hgt = [(-30) 60 (-30) 500 (-1000) 1000 (-500) 7.5 (-15) 7.5];
      m = zeros(size(xgrid));
      for j =1:length(pos)
        m = m + hgt(j).*((xgrid-pos(j)).^2).*(xgrid > pos(j));
      end
      m = m + 0.8;
    elseif idist == 10 ;
      tarstr = 'Time Shifted Sine' ;
      u = xgrid;
      for j =1:4,
        u = 0.5*(1-cos(pi*u));
      end
      m = 0.3*sin(3*pi*(u+xgrid)) + 0.5;
    end ;

    for isig = 1:2 ;

      if isig == 1 ;
        sigstr = 'Low Noise' ;
        sig = .02 ;
      elseif isig == 2 ;
        sigstr = 'High Noise' ;
        sig = .1 ;
      end ;

        randn('seed',92394873) ;
      ydat = m + sig * randn(nobs,1) ;
      data = [xgrid, ydat] ; 

      gpanal2(data) ;

      subplot(3,1,1) ;
      hold on ;
        plot(xgrid,m,'y--') ;
      hold off ;


        orient tall ;
      eval(['print -dpsc \matlab\steve\ps\g2j' ...
                                 num2str(idist) num2str(isig) '.ps']) ;

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

    end ;

  end ; 


elseif itest >= 21 ;    %  then do parameter tests

  dtitstr = 'Incomes' ;
  infstr = 'inc75.all' ;
  iincol = 1 ;    %  want all data as 1 column

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


  if itest == 21 ;
    gpanal2(data,0) ;
      title([dtitstr ' Grey color']) ;
  elseif itest == 22 ;
    gpanal2(data,1,[0,3]) ;
      title([dtitstr ' Restricted range']) ;
  elseif itest == 23 ;
    gpanal2(data,1,[0,3],21) ;
      title([dtitstr ' 21 curves']) ;
  elseif itest == 24 ;  
    gpanal2(data,0,[0,3],0,'gpantest') ;
      title([dtitstr ' make ps file, grey']) ;
  elseif itest == 25 ;  
    gpanal2(data,1,[0,3],0,'gpantest') ;
      title([dtitstr ' make ps file, color']) ;
  elseif itest == 26 ;
    gpanal2(data,1,[0,3,101]) ;
      title([dtitstr ' Restricted grid size']) ;
  end ;


end ;

