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


itest = 203 ;     %  1   Toy Example (demonstrating Q-Q idea)

                %  11  6 plots for N(0,1) data
                %  12  6 plots for N(-2,(.5)^2) data

                %  21  6 plots for Pareto(10) data
                %  22  6 plots for Pareto(2.5) data
                %  23  6 plots for Pareto(1.5) data

                %  31  6 plots for Weibull(1) data
                %  32  6 plots for Weibull(.1) data
                %  33  6 plots for Weibull(10) data

                %  101,...,140   parameter tests

                %  > 200  fiddle with large sample sizes


format compact ;



if itest == 1 ;    %  Make Toy example

  n = 5 ;
  seed = 34957345 ;
  randn('seed',seed) ;
  data = randn(n,1) ;

  left = -3 ;
  right = 3 ;
  madat = max(abs(data)) ;


  ngrid = 200 ;
  qgrid = linspace(left,right,ngrid) ;
  pgrid = phi(qgrid) ;


  sdata = sort(data) ;
  q2ep = reshape([sdata';sdata'],2*n,1) ;
  q2ep = [left; q2ep; right] ;
  p2ep = linspace(0,1,n+1) ;
  p2ep = reshape([p2ep;p2ep],2*(n+1),1) ;


  p2 = linspace(1/(2*n),1-1/(2*n),n)' ;
  q2e = sdata ;
  q2t = phiinv(p2) ;

  subplot(2,1,1) ;
    cdfh = plot(qgrid,pgrid,'b-', ...
                q2ep, p2ep,'r-', ...
                [q2t';q2t'],[zeros(1,n);p2'],'b-', ...
                [q2e';q2e'],[zeros(1,n);p2'],'r-', ...
                [left*ones(1,n);right*ones(1,n)],[p2'; p2'],'w--') ;
      set(cdfh(1),'LineWidth',3) ;
      set(cdfh(2),'LineWidth',3) ;
      title('Empirical and Theoretical C.D.F.s', ...
                 'FontSize',18) ;
      xlabel('quantile') ;
      ylabel('probability') ;


  subplot(2,1,2) ;
    qqh = plot(q2t,q2e,'m-', ...
               [-madat;madat],[-madat;madat],'g-', ...
               [-madat*ones(1,n);q2t'],[q2e';q2e'],'r-', ...
               [q2t';q2t'],[-madat*ones(1,n);q2e'],'b-') ;
      set(qqh(1),'LineWidth',3) ;
      set(qqh(2),'LineWidth',2) ;
      title('Q-Q Plot','FontSize',18) ;
      axis equal ;
      axis([-madat,madat,-madat,madat]) ;
      xlabel('theoretical quantile') ;
      ylabel('empirical quantile') ;




  orient tall ;

  prstr = ['print -dpsc \matlab\steve\ps\qqt' ...
                               num2str(itest) ...
                                        '.ps ;'] ;
  eval(prstr) ;




elseif  itest > 10  &  itest < 40  ;    %  work with 6 plot format



  iview = 2 ;    %  1 - top 3 using estimated parameters
                 %          bottom 3 using given parameters
                 %  2 - all estimated, direct on top
                 %          log dist's on bottom


  nobs = 500 ;


  if itest == 11 ;    %  N(0,1) data
    datstr = 'N(0,1)' ;
    mu = 0 ;
    sigma = 1 ;
    alpha = 1 ;

    randn('seed',23498575) ;
    data = mu + sigma * randn(nobs,1) ;

  elseif itest == 12 ;    %  N(10,(.1)^2) data
    datstr = 'N(10,(.1)^2)' ;
    mu = 10 ;
    sigma = 0.1 ;
    alpha = 1 ;

    randn('seed',23489553) ;
    data = mu + sigma * randn(nobs,1) ;

  elseif itest == 21 ;    %  Pareto(10) data
    datstr = 'Pareto(10)' ;
    alpha = 10 ;
    sigma = 1 ;
    mu = 2 * sigma ;

    rand('seed',67942384) ;
    data = sigma * (rand(nobs,1).^(-1/alpha) - 1) ;

  elseif itest == 22 ;    %  Pareto(2.5) data
    datstr = 'Pareto(2.5)' ;
    alpha = 2.5 ;
    sigma = 1 ;
    mu = 2 * sigma ;

    rand('seed',59328745) ;
    data = sigma * (rand(nobs,1).^(-1/alpha) - 1) ;

  elseif itest == 23 ;    %  Pareto(1.5) data
    datstr = 'Pareto(1.5)' ;
    alpha = 1.5 ;
    sigma = 1 ;
    mu = 2 * sigma ;

    rand('seed',83648563) ;
    data = sigma * (rand(nobs,1).^(-1/alpha) - 1) ;

  elseif itest == 31 ;    %  Weibull(1) data
    datstr = 'Weibull(1)' ;
    alpha = 1 ;
    sigma = 1 ;
    mu = 2 * sigma ;

    rand('seed',20987455) ;
    data = sigma * (-log(1 - rand(500,1))).^(1/alpha) ;

  elseif itest == 32 ;    %  Weibull(.1) data
    datstr = 'Weibull(.1)' ;
    alpha = .1 ;
    sigma = 1 ;
    mu = 2 * sigma ;

    rand('seed',93874553) ;
    data = sigma * (-log(1 - rand(500,1))).^(1/alpha) ;

  elseif itest == 33 ;    %  Weibull(10) data
    datstr = 'Weibull(10)' ;
    alpha = 10 ;
    sigma = 1 ;
    mu = 2 * sigma ;

    rand('seed',73403567) ;
    data = sigma * (-log(1 - rand(500,1))).^(1/alpha) ;

  end ;



  % Start graphics
  %
  figure(1) ;
  clf ;

  tstr = ['Q-Q, n = ' num2str(nobs) ', from ' datstr] ;

  paramstruct = struct('iscreenwrite',1, ...
                         'titlestr',tstr, ...
                         'ishowpar',1, ...
                         'idist',1 ...
                                          ) ;




  subplot(2,3,1) ;  %  Gaussian Dist, estimated mu and sigma
  paramstruct.idist = 1 ;
    qq(data,paramstruct) ;



  subplot(2,3,2) ;  %  Pareto Dist, estimated alpha and sigma
  paramstruct.idist = 2 ;
    qq(data,paramstruct) ;



  subplot(2,3,3) ;  %  Weibull Dist, estimated alpha and sigma
  paramstruct.idist = 3 ;
    qq(data,paramstruct) ;



  if iview == 1 ;


    %  Set up nonestimated frames
    %
    paramstruct.mu = mu ;
    paramstruct.sigma = sigma ;
    paramstruct.alpha = alpha ;


    subplot(2,3,4) ;  %  Gaussian Dist, given mu and sigma
    paramstruct.idist = 1 ;
      qq(data,paramstruct) ;


    subplot(2,3,5) ;  %  Pareto Dist, given alpha and sigma
    paramstruct.idist = 2 ;
      qq(data,paramstruct) ;


    subplot(2,3,6) ;  %  Weibull Dist, given alpha and sigma
    paramstruct.idist = 3 ;
      qq(data,paramstruct) ;



  elseif iview == 2 ;


    subplot(2,3,4) ;  %  lognormal Dist
    paramstruct.idist = 11 ;
      qq(data,paramstruct) ;


    subplot(2,3,5) ;  %  log - Pareto Dist
    paramstruct.idist = 12 ;
      qq(data,paramstruct) ;


    subplot(2,3,6) ;  %  log - Weibull Dist
    paramstruct.idist = 13 ;
      qq(data,paramstruct) ;



  end ;




  orient landscape ;

  prstr = ['print -dpsc \matlab\steve\ps\qqt' ...
                               num2str(itest) ...
                               'v'num2str(iview) ...
                                        '.ps ;'] ;
  eval(prstr) ;







elseif  itest > 50  &  itest < 60  ;    %  work with Weibull(10,1) data

  figure(1) ;
  clf ;

  rand('seed',34572093) ;
  data = (-log(rand(500,1))).^(1/10) ;

  if itest == 51 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 52 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'idist',2 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 53 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'idist',3 ...
                                          ) ;
    qq(data,paramstruct) ;

  end ;









elseif  100 < itest  &  itest < 200 ;    %  then do parameter tests

  figure(1) ;
  clf ;

  randn('state',sum(100*clock))
  data = randn(500,1) ;


  if itest == 101 ;    %  Try pure defaults
    qq(data) ;

  elseif itest == 102 ;    %  Try dummy paramstruct
    paramstruct = 3 ;
    qq(data,paramstruct) ;

  elseif itest == 103 ;
    paramstruct = struct('iscreenwrite',1 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 104 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'nsim',10 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 105 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'nsim',500 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 106 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'icolor',1 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 107 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'icolor',0 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 108 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'nsim',0, ...
                         'icolor',0 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 109 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'savestr','test1' ...
                                          ) ;
    qq(data,paramstruct) ;
    dir ;

  elseif itest == 110 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'savestr','\temp\test1' ...
                                          ) ;
    qq(data,paramstruct) ;
    dir \temp\*.ps ;

  elseif itest == 111 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'savestr','\temp\test1' ...
                                          ) ;
    qq(data,paramstruct) ;
    disp('now look at file \temp\test1.ps')

  elseif itest == 112 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'icolor',0, ...
                         'savestr','\temp\test1' ...
                                          ) ;
    qq(data,paramstruct) ;
    disp('now look at file \temp\test1.ps')

  elseif itest == 113 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'xlabelstr','Test_X' ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 114 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'ylabelstr','Test_Y' ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 115 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'titlestr','Try a title' ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 116 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'left',-1 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 117 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'right',5 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 118 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'bottom',-5 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 119 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'top',1 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 120 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'ioverlay',0 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 121 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'nsim',0, ...
                         'ioverlay',0 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 122 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'nsim',0 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 123 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'nsim',0, ...
                         'ioverlay',2 ...
                                          ) ;
    qq(data,paramstruct) ;

  elseif itest == 124 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'mu',0, ...
                         'sigma',1 ...
                                          ) ;
    qq(1.5*data+1,paramstruct) ;

  elseif itest == 125 ;
    paramstruct = struct('iscreenwrite',1, ...
                         'ioverlay',2, ...
                         'mu',0, ...
                         'sigma',1 ...
                                          ) ;
    qq(1.5*data+1,paramstruct) ;

  elseif itest == 126 ;
    vqalign = [.25; .75] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data),paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;

  elseif itest == 127 ;
    vqalign = [.60; .80] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data),paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;

  elseif itest == 128 ;
    vqalign = [.80; .95] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data),paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;

  elseif itest == 129 ;
    vqalign = [.60; .80] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'idist',3, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data),paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;

  elseif itest == 130 ;
    vqalign = [.80; .95] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'idist',3, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data),paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;


  elseif itest == 131 ;
    vqalign = [.5; .8] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'idist',2, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data).^(-.1)-1,paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;


  elseif itest == 132 ;
    vqalign = [.5; .8] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'idist',2, ...
                  'maxstep',10, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data).^(-.1)-1,paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;


  elseif itest == 133 ;
    vqalign = [.5; .8] ;
    paramstruct = struct('iscreenwrite',1, ...
                  'idist',2, ...
                  'relaccthreshold',.01, ...
                  'vqalign',vqalign ...
                                          ) ;
    qq(phi(data).^(-.1)-1,paramstruct) ;
    disp(['  check these meets at percentiles: ' ...
                 num2str(vqalign(1)) ', ' ...
                 num2str(vqalign(2))]) ;


  elseif itest == 134 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;


  elseif itest == 135 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'parfontsize',18, ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;


  elseif itest == 136 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'parfontsize',10, ...
                  'titlefontsize',24, ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;


  elseif itest == 137 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'idist',1 ...
                                          ) ;
    qq(exp(data),paramstruct) ;


  elseif itest == 138 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'idist',11 ...
                                          ) ;
    qq(exp(data),paramstruct) ;


  elseif itest == 139 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'idist',12 ...
                                          ) ;
    qq(phi(data).^(-.1)-1,paramstruct) ;


  elseif itest == 140 ;
    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'idist',13 ...
                                          ) ;
    qq(-log(phi(data)),paramstruct) ;




  end ;    %  of itest if block



elseif itest > 200 ;    %  the fiddle with large sample size



  if itest == 201 ;

    figure(1) ;
    clf ;

    randn('state',sum(100*clock))
    data = randn(50000,1) ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;



  elseif itest == 202 ;

    figure(1) ;
    clf ;

    randn('state',sum(100*clock))
    data = randn(50000,1) ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',50000, ...
                  'titlestr','nsimplotval = 50000', ...
                  'savestr','qqtest50000', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;


    figure(2) ;
    clf ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',3000, ...
                  'titlestr','nsimplotval = 3000', ...
                  'savestr','qqtest3000', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;


    figure(3) ;
    clf ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',900, ...
                  'titlestr','nsimplotval = 900', ...
                  'savestr','qqtest900', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;

    figure(4) ;
    clf ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',300, ...
                  'titlestr','nsimplotval = 300', ...
                  'savestr','qqtest300', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;


    figure(5) ;
    clf ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',30, ...
                  'titlestr','nsimplotval = 30', ...
                  'savestr','qqtest30', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;



  elseif itest == 203 ;

    figure(1) ;
    clf ;

    randn('state',sum(100*clock))
    data = randn(1000,1) ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',1000, ...
                  'titlestr','n = 1000, nsimplotval = 1000', ...
                  'savestr','qqtestn1000p1000', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;



    figure(2) ;
    clf ;

    paramstruct = struct('iscreenwrite',1, ...
                  'ishowpar',1, ...
                  'nsimplotval',900, ...
                  'titlestr','n = 1000, nsimplotval = 900', ...
                  'savestr','qqtestn1000p900', ...
                  'idist',1 ...
                                          ) ;
    qq(data,paramstruct) ;



  end ;


end ;  %  of itest if block



