disp('Running MATLAB script file gpsc1t.m') ;
%
%    Test program for gpsc1.m
%    General Purpose Significant Concav(vex)ity

ipart = 1 ;     %  0  simple tests, using itestpar below
                %  1  real data sets (density estimation)
                %  2  normal mixture examples (density estimation)
                %  3  real data sets (regression)
                %  4  NR book examples (regression)
                %  5  B-spline test examples (regression)
                %  6  ER Wavelet Examples (regression)

itestpar = 11 ;  %  test various parameter settings (when ipart = 0)
                %  1 kde, entered xgrid
                %  2 kde, entered xgrid and hgrid
                %  3 kde, entered hgrid
                %  4 kde, wierd range no truncation
                %  5 kde, wierd range with truncation
                %  6 kde, like 1, but alpha = .2
                %  7 kde, like 1, but alpha = .01
                %  8 kde, like 1, but ptwise
                %  9 llr, entered xgrid
                %  10 nwr, entered xgrid
                %  11 llr, just 1 row (for additive models)


autop = 0 ;     %  0 - create no postscript files, and pause at each pic
                %  1 - create postscript files, and pause at each pic
                %  2 - create postscript files, and no pause

format compact ;


if ipart == 0 ;   %  then do simple tests

  rand('seed',29384729) ;
  if itestpar == 1 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,2,31]) ;
  elseif itestpar == 2 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,2,31],[.1,5,10]) ;
  elseif itestpar == 3 ;
    data = rand(20,1) ;
    gpsc1(data,0,[.1,5,5]) ;
  elseif itestpar == 4 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,.5,31],0,0) ;
  elseif itestpar == 5 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,.5,31],0,1) ;
  elseif itestpar == 6 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,2,31],0,1,.2) ;
  elseif itestpar == 7 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,2,31],0,1,.01) ;
  elseif itestpar == 8 ;
    data = rand(20,1) ;
    gpsc1(data,[-1,2,31],0,1,.05,0) ;
  elseif itestpar == 9 ;
    data = rand(50,1) ;
    flag = data < .5 ;
    data = [data, flag .* data + (1 - flag) .* (1 - data)] ;
    gpsc1(data,[-0.5,1.5,21],0,1,.05,1,1) ;
  elseif itestpar == 10 ;
    data = rand(50,1) ;
    flag = data < .5 ;
    data = [data, flag .* data + (1 - flag) .* (1 - data)] ;
    gpsc1(data,[-0.5,1.5,21],0,1,.05,1,0) ;
  elseif itestpar == 11 ;
    data = rand(50,1) ;
    data = [data, (data < .7) & (data > .3)] ;
    vecout = gpsc1(data,[0,1,40],[.1,.1,1]) ;
    vecout
  end ;
%function [mapout,xgrid] = gpsc1(data,vxgp,vhgp,eptflag,alpha,simflag,llflag) 

elseif ipart == 1 ;    %  then do real density estimation data sets

  idat = 4 ;     %  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 = ['\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']) ;


    figure(1) ;
      ngrid = 401 ;
      inhp = 11 ;
      maxd = max(data) ;
      mind = min(data) ;
      range = maxd - mind ;
      binw = range / (ngrid - 1) ;
      hmin = 2 * binw ;
      hmax = range ;
      vh = logspace(log10(hmin),log10(hmax),inhp) ;
        [msmoo, xgrid] = gpkde(data,vh,[mind; maxd; ngrid]) ;
      plot(xgrid,msmoo,'c') ;
        title([dtitstr ' Data, ' date]) ;


    figure(2) ;
      gpsz1(data) ;
        title('SiZer') ;
        %  Set up colorful color map
        cocomap = [0,    0,   1; ...
                  .35, .35, .35; ...
                  .5,    0,  .5;
                   1,    0,   0] ; ...
        colormap(cocomap) ;


    figure(3) ;
      gpsc1(data) ;
        title('SiCon') ;
        %  Set up colorful color map
        cocomap = [1,   .5,   0; ...
                  .35, .35, .35; ...
                   0,    1,   0;
                   0,    1,   1] ; ...
        colormap(cocomap) ;


    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                   num2str(ipart) num2str(iidat) '.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    end ;

  end ;

elseif ipart == 2 ;    %  then do normal mixture density 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 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)]) ;

      data = nmdata(nobs,parmat) ;


      figure(1) ;
        ngrid = 401 ;
        inhp = 11 ;
        maxd = max(data) ;
        mind = min(data) ;
        range = maxd - mind ;
        binw = range / (ngrid - 1) ;
        hmin = 2 * binw ;
        hmax = range ;
        vh = logspace(log10(hmin),log10(hmax),inhp) ;
          [msmoo, xgrid] = gpkde(data,vh,[mind; maxd; ngrid]) ;
        plot(xgrid,msmoo,'c') ;
          title([dtitstr ' Data, ' date]) ;


      figure(2) ;
        gpsz1(data) ;
          title('SiZer') ;
          %  Set up colorful color map
          cocomap = [0,    0,   1; ...
                    .35, .35, .35; ...
                    .5,    0,  .5;
                     1,    0,   0] ; ...
          colormap(cocomap) ;


      figure(3) ;
        gpsc1(data) ;
          title('SiCon') ;
          %  Set up colorful color map
          cocomap = [1,   .5,   0; ...
                    .35, .35, .35; ...
                     0,    1,   0;
                     0,    1,   1] ; ...
          colormap(cocomap) ;



%      gpsc1(data) ;
%      title(['#' num2str(idist) ' ' titstr ...
%                  ', n = ' num2str(nobs) ', ' date]) ;

      if autop >= 1 ;
          orient landscape ;
        eval(['print -dps \matlab\steve\ps\egf' ...
                   num2str(ipart) num2str(idist) num2str(inobs) '.ps']) ;
      end ;

      if autop < 2 ;
        disp(['Any key to continue']) ;
        pause ;
      end ;

    end ;

  end ;

elseif ipart == 3 ;    %  then do real regression 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
                 %  8 - Fossil Data
                 %  9 - Share Yield Data

  if idat == 0 ;
    idats = 1 ;
    idate = 9 ;
  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 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' ;
    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 = '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 = ['\matlab\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 ;


subplot(2,1,1) ;
        ngrid = 401 ;
        inhp = 11 ;
        maxd = max(data(:,1)) ;
        mind = min(data(:,1)) ;
        range = maxd - mind ;
        binw = range / (ngrid - 1) ;
        hmin = 2 * binw ;
        hmax = range ;
        vh = logspace(log10(hmin),log10(hmax),inhp) ;
          [msmoo, xgrid] = gpnpr(data,vh,[mind; maxd; ngrid]) ;
        plot(xgrid,msmoo,'c') ;

subplot(2,1,2) ;
    gpsc1(data,0,0,0,0.2) ;
      title([dtitstr ', ' date]) ;
      xlabel(xstr) ;
          %  Set up colorful color map
          cocomap = [1,   .5,   0; ...
                    .35, .35, .35; ...
                     0,    1,   0;
                     0,    1,   1] ; ...
          colormap(cocomap) ;


    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) num2str(iidat) '.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    end ;

  end ;

elseif ipart == 4 ;    %  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('\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)]) ;
%      subplot(2,3,ieg) ;
subplot(2,1,1) ;
plot(data(:,1),data(:,2),'c-',data(:,1),data(:,ieg+2),'yo') ;
subplot(2,1,2) ;
        gpsc1([data(:,1), data(:,ieg+2)]) ;
          title(['NR Book e.g. 1, ' date]) ;
          %  Set up colorful color map
          cocomap = [1,   .5,   0; ...
                    .35, .35, .35; ...
                     0,    1,   0;
                     0,    1,   1] ; ...
          colormap(cocomap) ;
gppause ;

    end ;

    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) '1.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    end ;

  end ;


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

    fid = fopen('\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)]) ;
      subplot(2,3,ieg) ;
        gpsc1([data(:,ieg), data(:,8)]) ;
          title(['NR Book e.g. 2, ' date]) ;
    end ;

    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) '2.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    end ;

  end ;

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

    fid = fopen('\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)]) ;
      subplot(2,3,ieg) ;
        gpsc1([data(:,1), data(:,ieg+2)]) ;
          title(['NR Book e.g. 3, ' date]) ;
    end ;

    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) '3.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    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('\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 ;
    gpsc1([data(:,1), data(:,3)]) ;
      title(['NR Book, e.g. 5a, ' date]) ;

    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) '5a.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    end ;

    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 ;
    gpsc1([data(:,1), data(:,3)]) ;
      title(['NR Book, e.g. 5b, ' date]) ;
          hold on ;
            plot(data(:,1),data(:,3),'.w') ;
            plot(data(:,1),data(:,2),'--b') ;
          hold off ;

    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) '5b.ps']) ;
    end ;

  end ;

elseif ipart == 5 ;    %  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) ;

    gpsc1([data(:,1), data(:,3)]) ;
      title(['B-spline Setting ' num2str(ifile) ', ' date]) ;

    if autop >= 1 ;
        orient landscape ;
      eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) num2str(ifile) '.ps']) ;
    end ;

    if autop < 2 ;
      disp(['Any key to continue']) ;
      pause ;
    end ;

  end ;

elseif ipart == 6 ;    %  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] ; 

      gpsc1(data) ;
          title([tarstr ', ' sigstr ', ' date]) ;

      if autop >= 1 ;
          orient landscape ;
        eval(['print -dps \matlab\steve\ps\egf' ...
                 num2str(ipart) num2str(idist) num2str(isig) '.ps']) ;
      end ;

      if autop < 2 ;
        disp(['Any key to continue']) ;
        pause ;
      end ;

    end ;

  end ; 

end ;

