disp('Running MATLAB script file CHSZeg1.m') ;
%
%    Censored Hazard SiZer, Example 1. 
%    
%    Creates examples for the Censored Hazard SiZer study
%



ipart = 5 ;      %  1 - Bimodal Gaussian with rounding
                 %  2 - Bimodal Gaussian with censoring
                 %  3 - Device Put Data from Jiang, Jian-Cheng
                 %  4 - Stanford Heart Transplant
                 %  5 - Stanford Heart Transplant, log scale



if ipart == 1 ;

  n = 500 ;
  left = -2.5 ;
  right = 2.5 ;
  rand('seed',93248548) ;
  randn('seed',36982643) ;

  vdata = randn(n,1) / 2 ;
          %  column vector of i.i.d. Gaussian(0,1/4)s
  vdata = vdata + 1 * (1 - 2 * (rand(n,1) > .5)) ;
          %  add random +- 1

  vdata = round(vdata * 10) / 10 ;
          %  Round to nearest 0.01

    foltit = ['n = ' num2str(n) ...
              ', .5N(-1,1/4) + .5N(1,1/4) data, rounded to 0.1'] ;
  paramstruct = struct('iout',3,...
                       'famoltitle',foltit,...
                       'sizertitle','Regular SiZer Map',...
                       'titlefontsize',15,...
                       'imovie',0,...
                       'minx',left,...
                       'maxx',right,...
                       'iscreenwrite',1) ;

  figure(1) ;
  clf ;
  sizer2(vdata,paramstruct) ;

    orient tall ;
    pstr = ['print -dpsc \matlab\steve\ps\CHSZeg1' ...
                           'p' num2str(ipart) '.ps'] ;
  eval(pstr) ;




elseif ipart == 2 ;

  n = 500 ;
  left = -2.5 ;
  right = 2 ;
  rand('seed',43570293) ;
  randn('seed',19237419) ;


  %  Generate data
  %
  vt = randn(n,1) / 2 ;
          %  column vector of i.i.d. Gaussian(0,1/4)s
  vt = vt + 1 * (1 - 2 * (rand(n,1) > .5)) ;
          %  add random +- 1

  vc = randn(n,1) / 2 ;
          %  column vector of i.i.d. Gaussian(0,1/4)s
  vc = vc + 1 * (1 - 2 * (rand(n,1) > .5)) ;
          %  add random +- 1

  vdata = min(vt,vc) ;
  vdel = vdata == vt ;

  vucdata = vdata(vdel) ;
  nuc = length(vucdata) ;
  vcdata = vdata(~vdel) ;


  figure(1) ;
  clf ;


  subplot(2,2,1) ;    %  upper left, uncensored family
    foltit = ['n = ' num2str(nuc+1) ...
              ' only uncensored data'] ;
    paramstruct = struct('iout',4,...
                         'famoltitle',foltit,...
                         'titlefontsize',18,...
                         'imovie',0,...
                         'minx',left,...
                         'maxx',right,...
                         'ihhighlight',0,...
                         'ndataoverlay',0,...
                         'iscreenwrite',1) ;
  sizer2(vucdata,paramstruct) ;
  hold on ;    %  add jitter plot
    vax = axis ;
    vjit1 = rand(nuc,1) ;
    randht = vax(3) + (0.6 + 0.3 * vjit1) * ...
                              (vax(4) - vax(3)) ;
    plot(vucdata,randht,'g.') ;
  hold off ;


  subplot(2,2,3) ;    %  upper left, uncensored SiZer
    paramstruct = struct('iout',6,...
                         'sizertitle','Uncensored SiZer Map',...
                         'titlefontsize',18,...
                         'imovie',0,...
                         'minx',left,...
                         'maxx',right,...
                         'ihhighlight',0,...
                         'ndataoverlay',0,...
                         'iscreenwrite',1) ;
  sizer2(vucdata,paramstruct) ;


  subplot(2,2,2) ;    %  upper right, censored family
    foltit = ['n = ' num2str(n) ...
              ' full, censored data'] ;
    paramstruct = struct('iout',4,...
                         'famoltitle',foltit,...
                         'titlefontsize',18,...
                         'icensor',1,...
                         'imovie',0,...
                         'minx',left,...
                         'maxx',right,...
                         'ndataoverlay',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vdel],paramstruct) ;
  hold on ;    %  add jitter plot
    vax = axis ;
    randht = vax(3) + (0.6 + 0.3 *vjit1) * ...
                              (vax(4) - vax(3)) ;
          %  Note, this uses same vjit1 as defined above
    plot(vucdata,randht,'g.') ;

    vjit2 = rand(n-nuc,1) ;

    randht = vax(3) + (0.6 + 0.3 * vjit2) * ...
                              (vax(4) - vax(3)) ;
          %  Note, this uses same vjit1 as defined above
    plot(vcdata,randht,'y.') ;
  hold off ;


  subplot(2,2,4) ;    %  lower right, censored SiZer
    paramstruct = struct('iout',6,...
                         'sizertitle','Censored SiZer Map',...
                         'titlefontsize',18,...
                         'icensor',1,...
                         'imovie',0,...
                         'minx',left,...
                         'maxx',right,...
                         'ndataoverlay',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vdel],paramstruct) ;


    orient landscape ;
    pstr = ['print -dpsc \matlab\steve\ps\CHSZeg1' ...
                           'p' num2str(ipart) '.ps'] ;
  eval(pstr) ;




elseif ipart == 3 ;

  fid = fopen('/matlab/steve/data/device.put.dat','rt') ;
          %  'rt' is for "read only" and "text"

    [vdata,cnt] = fscanf(fid,'%f') ;

  fclose(fid) ;

  titstr = 'Device Lifetimes' ;

  subplot(2,2,1) ;
    paramstruct = struct('iout',4,...
                         'famoltitle',['Density Estimates, ' titstr],...
                         'titlefontsize',18,...
                         'ihazard',0,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2(vdata,paramstruct) ;

  subplot(2,2,3) ;
    paramstruct = struct('iout',6,...
                         'sizertitle','Density SiZer Map',...
                         'titlefontsize',18,...
                         'ihazard',0,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2(vdata,paramstruct) ;

  subplot(2,2,2) ;
    paramstruct = struct('iout',4,...
                         'famoltitle',['Hazard Estimates, ' titstr],...
                         'titlefontsize',18,...
                         'ihazard',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2(vdata,paramstruct) ;

  subplot(2,2,4) ;
    paramstruct = struct('iout',6,...
                         'sizertitle','Hazard SiZer Map',...
                         'titlefontsize',18,...
                         'ihazard',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2(vdata,paramstruct) ;


    orient landscape ;
    pstr = ['print -dpsc \matlab\steve\ps\CHSZeg1' ...
                           'p' num2str(ipart) '.ps'] ;
  eval(pstr) ;




elseif ipart == 4 ;

  fid = fopen('/matlab/steve/data/stanford.dat','rt') ;
          %  'rt' is for "read only" and "text"

    [vdata,cnt] = fscanf(fid,'%f') ;

  fclose(fid) ;

  nobs = cnt / 2 ;
  vcensor = vdata(1:nobs) ;
  vdata = vdata((nobs+1):cnt) ;

  titstr = 'Heart Transplants' ;

  subplot(2,2,1) ;
    paramstruct = struct('iout',4,...
                         'famoltitle',['Density Estimates, ' titstr],...
                         'titlefontsize',18,...
                         'xlabelstr','Days to death',...
                         'labelfontsize',15,...
                         'ihazard',0,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;

  subplot(2,2,3) ;
    paramstruct = struct('iout',6,...
                         'sizertitle','Density SiZer Map',...
                         'titlefontsize',18,...
                         'xlabelstr','Days to death',...
                         'labelfontsize',15,...
                         'ihazard',0,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;

  subplot(2,2,2) ;
    paramstruct = struct('iout',4,...
                         'famoltitle',['Hazard Estimates, ' titstr],...
                         'titlefontsize',18,...
                         'xlabelstr','Days to death',...
                         'labelfontsize',15,...
                         'ihazard',1,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;

  subplot(2,2,4) ;
    paramstruct = struct('iout',6,...
                         'sizertitle','Hazard SiZer Map',...
                         'titlefontsize',18,...
                         'xlabelstr','Days to death',...
                         'labelfontsize',15,...
                         'ihazard',1,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;


    orient landscape ;
    pstr = ['print -dpsc \matlab\steve\ps\CHSZeg1' ...
                           'p' num2str(ipart) '.ps'] ;
  eval(pstr) ;




elseif ipart == 5 ;

  fid = fopen('/matlab/steve/data/stanford.dat','rt') ;
          %  'rt' is for "read only" and "text"

    [vdata,cnt] = fscanf(fid,'%f') ;

  fclose(fid) ;

  nobs = cnt / 2 ;
  vcensor = vdata(1:nobs) ;
  vdata = vdata((nobs+1):cnt) ;
  vdata = log10(vdata) ;

  titstr = 'Heart Transplants' ;

  subplot(2,2,1) ;
    paramstruct = struct('iout',4,...
                         'famoltitle',['Density Estimates, ' titstr],...
                         'titlefontsize',18,...
                         'xlabelstr','log_{10}(Days to death)',...
                         'labelfontsize',15,...
                         'ihazard',0,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;

  subplot(2,2,3) ;
    paramstruct = struct('iout',6,...
                         'sizertitle','Density SiZer Map',...
                         'titlefontsize',18,...
                         'xlabelstr','log_{10}(Days to death)',...
                         'labelfontsize',15,...
                         'ihazard',0,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;

  subplot(2,2,2) ;
    paramstruct = struct('iout',4,...
                         'famoltitle',['Hazard Estimates, ' titstr],...
                         'titlefontsize',18,...
                         'xlabelstr','log_{10}(Days to death)',...
                         'labelfontsize',15,...
                         'ihazard',1,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;

  subplot(2,2,4) ;
    paramstruct = struct('iout',6,...
                         'sizertitle','Hazard SiZer Map',...
                         'titlefontsize',18,...
                         'xlabelstr','log_{10}(Days to death)',...
                         'labelfontsize',15,...
                         'ihazard',1,...
                         'icensor',1,...
                         'imovie',0,...
                         'iscreenwrite',1) ;
  sizer2([vdata,vcensor],paramstruct) ;


    orient landscape ;
    pstr = ['print -dpsc \matlab\steve\ps\CHSZeg1' ...
                           'p' num2str(ipart) '.ps'] ;
  eval(pstr) ;




end ;    %  of ipart if block



