disp('Running MATLAB script file UNC2001RespSiz1.m') ;
%
%    For making analysis of Response Data, based on UNC traces
%    gathered in April 2001, by Don Smith and Co.
%    and filtered by Felix Hernandez Campos.
%
%    Because of the size of the data, they are summarized as:
%         Summary statistics:
%              rsp01power.xls
%              rsp01morestats.xls  (with trimmed down version in 
%                                           rsp01morestatsraw.xls
%                                   and just data saved in
%                                           rsp01morestatsraw.txt)
%         Quantiles of body of distribution:
%              responsesizequantiles.xls
%         Individual responses:
%              20010420_1300.raw
%                     to
%              20010429_800.raw
%                (21 time blocks)
%
%    For further details, see the files:
%         responsesizetailsnotes.txt
%         daydatekey.txt


ipart = 84 ;    %  0  -  save Data Summaries and Quantiles in .mat files
               %  1  -  make simple summary graphics
               %  2  -  make tail duration vs. size scatterplots (biggest 2,000)
               %  3  -  Do Individual tail index analyses 
               %                  (and save avg. Hill estimates)
               %  4  -  Do rescaled duration vs. size scatterplots
               %  5  -  Square root rescaling of plot in 4
               %  6  -  Do Ratio tail index estimates
               %  8  -  make tail rate vs. size scatterplots (biggest 2,000)
               %  9  -  make tail inverse rate vs. time scatterplots (biggest 2,000)
               %  11  -  Box-Cox transformed family of plots like ipart == 2,
               %              tail duration vs. size scatterplots (biggest 2,000)
               %  21  -  Log-Log Complementary CDF
               %  22  -  slopes of Log-Log Complementary CDF
    %  23  -  improved slopes of Log-Log Complementary CDF
               %  31  -  Asymptotic Indep. Angle Plot (Time vs. Size)
               %  32  -  Tail Index Analysis for Ratio  (Size / Time)
               %  33  -  Asymptotic Indep. Angle Plot (Ratio vs. Size)
               %  34  -  Asymptotic Indep. Angle Plot (Time vs. Ratio)
               %  41  -  Full Data Log-Log Complementary CDF
               %  42  -  Full Data slopes of Log-Log Complementary CDF
               %  43  -  Full Data Pareto Q-Q 
               %  44  -  Full Data log-normal Q-Q 
               %  51  -  Full Data Log-Log CCDF, 
               %             plus Double-Pareto log-normal (fit by eye)
               %  52  -  Full Data Log-Log CCDF, 
               %             plus log-normal 2-mixture (fit by eye)
               %  53  -  Full Data Log-Log CCDF, 
               %             plus Double-Pareto log-normal 2-mixture (fit by eye)
               %  54  -  Full Data Log-Log CCDF, 
               %             plus log-normal 3-mixture (fit by eye)
               %  55  -  Full Data Log-Log CCDF, 
               %             plus Double-Pareto log-normal 3-mixture (fit by eye)
               %  61  -  Same as 51, but add lines showing variability
               %  62  -  Same as 52, but add lines showing variability
               %  63  -  Same as 53, but add lines showing variability
               %  64  -  Same as 54, but add lines showing variability
               %  65  -  Same as 55, but add lines showing variability
               %  71  -  Same as 41, but add 100 simulated Pareto's
               %                     (from part 43) to upper left
               %  72  -  Same as 41, but add 100 sim'd Double-Pareto log-normal's
               %                     (from part 51) to upper left
               %  73  -  Same as 41, but add 100 sim'd D-P log-n 3-mix's
               %                     (from part 54) to upper left
               %  74  -  Same as 41, but add 100 sim'd log-gam's
               %                     (from part 54) to upper left
               %  81  -  extremal index analysis > 100k, response size
               %  82  -  extremal index analysis > 100k, response times
               %  83  -  extremal index analysis > 100k, response rates
               %  84  -  extremal index analysis > 100k, response inverse rates
               %                  (these are >100 k versions of ipart = 21-28
               %                     from UNC2001RespSiz2.m)



%vitimeb = 1:21 ;    %  vector of indices of time blocks
%vitimeb = 12 ;    %  vector of indices of time blocks
%vitimeb = 19 ;    %  vector of indices of time blocks
vitimeb = [11 19] ;    %  vector of indices of time blocks
%vitimeb = 10 ;    %  vector of indices of time blocks
                    %  1 - Monday, April 23, 8:00-12:00
                    %  2 - Monday, April 23, 13:00-17:00
                    %  3 - Monday, April 23, 19:30-23:30
                    %  4 - Tuesday, April 24, 8:00-12:00
                    %  5 - Tuesday, April 24, 13:00-17:00
                    %  6 - Tuesday, April 24, 19:30-23:30
                    %  7 - Wednesday, April 25, 8:00-12:00
                    %  8 - Wednesday, April 25, 13:00-17:00
                    %  9 - Wednesday, April 25, 19:30-23:30
                    %  10 - Thursday, April 26, 8:00-12:00
                    %  11 - Thursday, April 26, 13:00-17:00
                    %  12 - Thursday, April 26, 19:30-23:30
                    %  13 - Friday, April 20, 8:00-12:00
                    %  14 - Friday, April 20, 13:00-17:00
                    %  15 - Friday, April 20, 19:30-23:30
                    %  16 - Saturday, April 21, 8:00-12:00
                    %  17 - Saturday, April 21, 13:00-17:00
                    %  18 - Saturday, April 21, 19:30-23:30
                    %  19 - Sunday, April 29, 8:00-12:00
                    %  20 - Sunday, April 29, 13:00-17:00
                    %  21 - Sunday, April 29, 19:30-23:30






if ipart == 0 ;    %  then store summary statistics as .mat file


  %  Load Data Summaries
  %
  fid = fopen('DataUNC2001\rsp01morestatsraw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [summdata,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  summdata = reshape(summdata,5,21)' ;
          %  use this for 21x5 matrix, since reshape works column wise

%  summdata
      %  used to check read in


  if vitimeb ~= 1:21 ;
    disp('!!!   Warning from UNC2001RespSiz1.m:') ;
    disp('!!!   Wrong value of vitimeb,') ;
    disp('!!!   Resetting to correct value') ;
    vitimeb = 1:21 ;
  end ;


  vnumresp = [] ;
  vtotsize = [] ;
  vmeansize = [] ;
  vsdsize = [] ;
  virowraw = [] ;
  daystrcell = cell(21,1) ;
  timestrcell = cell(21,1) ;
  filenamecell = cell(21,1) ;
  for itimeb = vitimeb ;    %  loop through time blocks
    if itimeb == 1 ;    %  Monday, April 23, 8:00-12:00
      daystr = 'Monday' ;
      timestr = 'Morning' ;
      irowraw = 9 ;
      filename = '20010423_800.raw' ;
    elseif itimeb == 2 ;    %  Monday, April 23, 13:00-17:00
      daystr = 'Monday' ;
      timestr = 'Afternoon' ;
      irowraw = 7 ;
      filename = '20010423_1300.raw' ;
    elseif itimeb == 3 ;    %  Monday, April 23, 19:30-23:30
      daystr = 'Monday' ;
      timestr = 'Evening' ;
      irowraw = 8 ;
      filename = '20010423_1930.raw' ;
    elseif itimeb == 4 ;    %  Tuesday, April 24, 8:00-12:00
      daystr = 'Tuesday' ;
      timestr = 'Morning' ;
      irowraw = 12 ;
      filename = '20010424_800.raw' ;
    elseif itimeb == 5 ;    %  Tuesday, April 24, 13:00-17:00
      daystr = 'Tuesday' ;
      timestr = 'Afternoon' ;
      irowraw = 10 ;
      filename = '20010424_1300.raw' ;
    elseif itimeb == 6 ;    %  Tuesday, April 24, 19:30-23:30
      daystr = 'Tuesday' ;
      timestr = 'Evening' ;
      irowraw = 11 ;
      filename = '20010424_1930.raw' ;
    elseif itimeb == 7 ;    %  Wednesday, April 25, 8:00-12:00
      daystr = 'Wednesday' ;
      timestr = 'Morning' ;
      irowraw = 15 ;
      filename = '20010425_800.raw' ;
    elseif itimeb == 8 ;    %  Wednesday, April 25, 13:00-17:00
      daystr = 'Wednesday' ;
      timestr = 'Afternoon' ;
      irowraw = 13 ;
      filename = '20010425_1300.raw' ;
    elseif itimeb == 9 ;    %  Wednesday, April 25, 19:30-23:30
      daystr = 'Wednesday' ;
      timestr = 'Evening' ;
      irowraw = 14 ;
      filename = '20010425_1930.raw' ;
    elseif itimeb == 10 ;    %  Thursday, April 26, 8:00-12:00
      daystr = 'Thursday' ;
      timestr = 'Morning' ;
      irowraw = 18 ;
      filename = '20010426_800.raw' ;
    elseif itimeb == 11 ;    %  Thursday, April 26, 13:00-17:00
      daystr = 'Thursday' ;
      timestr = 'Afternoon' ;
      irowraw = 16 ;
      filename = '20010426_1300.raw' ;
    elseif itimeb == 12 ;    %  Thursday, April 26, 19:30-23:30
      daystr = 'Thursday' ;
      timestr = 'Evening' ;
      irowraw = 17 ;
      filename = '20010426_1930.raw' ;
    elseif itimeb == 13 ;    %  Friday, April 20, 8:00-12:00
      daystr = 'Friday' ;
      timestr = 'Morning' ;
      irowraw = 3 ;
      filename = '20010420_800.raw' ;
    elseif itimeb == 14 ;    %  Friday, April 20, 13:00-17:00
      daystr = 'Friday' ;
      timestr = 'Afternoon' ;
      irowraw = 1 ;
      filename = '20010420_1300.raw' ;
    elseif itimeb == 15 ;    %  Friday, April 20, 19:30-23:30
      daystr = 'Friday' ;
      timestr = 'Evening' ;
      irowraw = 2 ;
      filename = '20010420_1930.raw' ;
    elseif itimeb == 16 ;    %  Saturday, April 21, 8:00-12:00
      daystr = 'Saturday' ;
      timestr = 'Morning' ;
      irowraw = 6 ;
      filename = '20010421_800.raw' ;
    elseif itimeb == 17 ;    %  Saturday, April 21, 13:00-17:00
      daystr = 'Saturday' ;
      timestr = 'Afternoon' ;
      irowraw = 4 ;
      filename = '20010421_1300.raw' ;
    elseif itimeb == 18 ;    %  Saturday, April 21, 19:30-23:30
      daystr = 'Saturday' ;
      timestr = 'Evening' ;
      irowraw = 5 ;
      filename = '20010421_1930.raw' ;
    elseif itimeb == 19 ;    %  Sunday, April 29, 8:00-12:00
      daystr = 'Sunday' ;
      timestr = 'Morning' ;
      irowraw = 21 ;
      filename = '20010429_800.raw' ;
    elseif itimeb == 20 ;    %  Sunday, April 29, 13:00-17:00
      daystr = 'Sunday' ;
      timestr = 'Afternoon' ;
      irowraw = 19 ;
      filename = '20010429_1300.raw' ;
    elseif itimeb == 21 ;    %  Sunday, April 29, 19:30-23:30
      daystr = 'Sunday' ;
      timestr = 'Evening' ;
      irowraw = 20 ;
      filename = '20010429_1930.raw' ;
    end ;


    %  create main outputs
    %
    vnumresp = [vnumresp; summdata(irowraw,1)] ;
    vtotsize = [vtotsize; summdata(irowraw,2)] ;
    vmeansize = [vmeansize; summdata(irowraw,3)] ;
    vsdsize = [vsdsize; summdata(irowraw,5)] ;
    virowraw = [virowraw ; irowraw] ;

    daystrcell{itimeb} = daystr ;
    timestrcell{itimeb} = timestr ;
    filenamecell{itimeb} = filename ;

    disp(['  ' daystr ' ' timestr ' ' filename])  ;



  end ;    %  of itimeb loop
  disp('  ') ;



  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise

  %mq
      %  used to check read in



  %  Check everything by looping through time blocks again to check maxima
  %
  vmaxsize = [] ;
  for itimeb = vitimeb ;    %  loop through time blocks

    maxsize = mq(1000,virowraw(itimeb)) ;
    vmaxsize = [vmaxsize; maxsize] ;

    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;


    %size(mtailresp)
    %pauseSM
        %  lines used for checking


    %  do formal check of max from quantiles, against max from tails
    %
    maxtailresp = max(mtailresp(:,1)) ;
    disp(['      ' num2str(itimeb) '.  Check this is zero: ' ...
                num2str(abs(maxtailresp - maxsize))]) ;



  end ;    %  of itimeb loop
  disp('  ') ;




  %  save Data Summaries in a .mat file
  %
  save('DataUNC2001\UNC2001RespSiz1DatSumm.mat','vnumresp', ...
               'vtotsize','vmeansize','vsdsize','virowraw','vmaxsize', ...
               'daystrcell','timestrcell','filenamecell') ;



  %  save Quantiles in a .mat file
  %
  save('DataUNC2001\UNC2001RespSiz1Quant.mat','mq') ;

  disp('   UNC2001RespSiz1.m finished data saves') ;



elseif ipart == 1 ;    %  then make simple summary graphics


  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  visubpart = 1:6 ;
%  visubpart = 1 ;    %  1 - Number of Responses Summary Plot
                      %  2 - Total Size Summary Plot
                      %  3 - Mean Size Summary Plot
                      %  4 - S.D. Size Summary Plot
                      %  5 - Coeff. Var. Size Summary Plot
                      %  6 - Max Size Summary Plot


  for isp = visubpart ;

    if isp == 1 ;    %  then make Number of Responses Plot
      vdat = vnumresp ;
      titstr = 'Number of HTTP Responses' ;
      pstr = 'UNC2001RS1NumRes.ps' ;
    elseif isp == 2 ;    %  then make Total Size Summary Plot
      vdat = vtotsize ;
      titstr = 'Total Response Size' ;
      pstr = 'UNC2001RS1TotSiz.ps' ;
    elseif isp == 3 ;    %  then make Mean Size Summary Plot
      vdat = vmeansize ;
      titstr = 'Mean Response Size' ;
      pstr = 'UNC2001RS1MeanSiz.ps' ;
    elseif isp == 4 ;    %  then make S.D. Size Summary Plot
      vdat = vsdsize ;
      titstr = 'S.D. Response Size' ;
      pstr = 'UNC2001RS1SdSiz.ps' ;
    elseif isp == 5 ;    %  then make Coeff. Var. Summary Plot
      vdat = vsdsize ./ vmeansize ;
      titstr = 'Coefficient of Variation' ;
      pstr = 'UNC2001RS1CoeVar.ps' ;
    elseif isp == 6 ;    %  then make Max Size Summary Plot
      vdat = vmaxsize ;
      titstr = 'Max Response Size' ;
      pstr = 'UNC2001RS1MaxSiz.ps' ;
    end ;


    figure(1) ;
    clf ;
    colormap([0 0 1; 1 0 1; 1 0 0]) ;
        %  set colormap to only have colors Blue, Magenta, Red


    subplot(2,2,1) ;    %  Mornings only
      vind = 3 * (1:7) - 2 ;
      vdatsub = vdat(vind) ;
      top = 1.1 * max(vdatsub) ;
      left = 0.5 ;
      right = 7.5 ;

      ph = bar((1:7)',vdatsub) ;
        shading flat ;
        axis([left,right,0,top]) ;
        title([titstr ', Morning'],'FontSize',15) ;
        set(ph,'FaceColor','b') ;

      for i = 1:7 ;
        daystr = daystrcell{vind(i)} ;
        tx = i ;
        ty = 0 + 0.5 * top ;
        th = text(tx,ty,daystr,'Rotation',90,'FontSize',15) ;
      end ;


    subplot(2,2,2) ;    %  Afternoons only
      vind = 3 * (1:7) - 1 ;
      vdatsub = vdat(vind) ;
      top = 1.1 * max(vdatsub) ;
      left = 0.5 ;
      right = 7.5 ;

      ph = bar((1:7)',vdatsub) ;
        shading flat ;
        axis([left,right,0,top]) ;
        title([titstr ', Afternoon'],'FontSize',15) ;
        set(ph,'FaceColor','m') ;

      for i = 1:7 ;
        daystr = daystrcell{vind(i)} ;
        tx = i ;
        ty = 0 + 0.5 * top ;
        th = text(tx,ty,daystr,'Rotation',90,'FontSize',15) ;
      end ;


    subplot(2,2,3) ;    %  Evenings only
      vind = 3 * (1:7) ;
      vdatsub = vdat(vind) ;
      top = 1.1 * max(vdatsub) ;
      left = 0.5 ;
      right = 7.5 ;

      ph = bar((1:7)',vdatsub) ;
        shading flat ;
        axis([left,right,0,top]) ;
        title([titstr ', Evening'],'FontSize',15) ;
        set(ph,'FaceColor','r') ;

      for i = 1:7 ;
        daystr = daystrcell{vind(i)} ;
        tx = i ;
        ty = 0 + 0.5 * top ;
        th = text(tx,ty,daystr,'Rotation',90,'FontSize',15) ;
      end ;


    subplot(2,2,4) ;    %  Everything
      vdatsub = reshape(vdat,3,7)' ;
      top = 1.1 * max(max(vdatsub)) ;
      left = 0.5 ;
      right = 7.5 ;

      ph = bar((1:7)',vdatsub,'grouped') ;
        shading flat ;
        axis([left,right,0,top]) ;
        title([titstr ', All Times'],'FontSize',15) ;

      for i = 1:7 ;
        daystr = daystrcell{vind(i)} ;
        tx = i ;
        ty = 0 + 0.5 * top ;
        th = text(tx,ty,daystr,'Rotation',90,'FontSize',15) ;
      end ;


      orient landscape ;
    print('-dpsc2',pstr) ;


  end ;




elseif  ipart == 2  | ...
        ipart == 8  | ...
        ipart == 9  ;    %  then make tail scatterplots (biggest 2,000):
                         %    -  duration vs. size 
                         %    -  rate vs. size
                         %    -  inverse rate vs. time

  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;
  
  
  %  loop through time blocks, read big data files, and make plots
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailresp = mtailresp(:,[1 4]) ;
        %  keep only Response Sizes and durations

    [temp,vind] = sort(mtailresp(:,1)) ;
    mtailresp = flipud(mtailresp(vind,:)) ;
        %  sort on first column, put in decreasing order

    vsize = mtailresp(1:ntop,1) ;
    vtime = mtailresp(1:ntop,2) ;
        %  work with only top values

    if ipart == 2 ;

      vx = vsize ;
      vy = vtime ;
      xlabelstr = 'Response Size (bytes)' ;
      ylabelstr = 'Response Duration (sec)' ;
      pstr = ['UNC2001RS1SPTimeSize' num2str(itimeb) '.ps'] ;

    elseif ipart == 8 ;

      vx = vsize ;
      vy = vsize ./ vtime ;
      xlabelstr = 'Response Size (bytes)' ;
      ylabelstr = 'Rate (bytes / sec)' ;
      pstr = ['UNC2001RS1SPRateSize' num2str(itimeb) '.ps'] ;

    elseif ipart == 9 ;

      vx = vtime ;
      vy = vtime ./ vsize ;
      xlabelstr = 'Response Duration (sec)' ;
      ylabelstr = 'Inverse Rate (sec / byte)' ;
      pstr = ['UNC2001RS1SPIrateTime' num2str(itimeb) '.ps'] ;

    end ;


    %  make main plot
    %
    left = min(vx) ;
    right = max(vx) ;
    xrange = right - left ;
    left = left - 0.01 * xrange ;
    right = right + 0.01 * xrange ;

    bottom = min(vy) ;
    top = max(vy) ;
    yrange = top - bottom ;
    bottom = bottom - 0.01 * yrange ;
    top = top + 0.01 * yrange ;

    plot(vx,vy,'b+') ;
      title(['HTTP Response Size Scatterplot, ' daystrcell{itimeb} ...
                             ', ' timestrcell{itimeb}],'FontSize',18) ;
      axis([left,right,bottom,top]) ;
      xlabel(xlabelstr,'FontSize',18) ;
      ylabel(ylabelstr,'FontSize',18) ;

      tx = left + 0.2 * xrange ;
      ty = bottom + 0.9 * yrange ;
      text(tx,ty,['Largest ' num2str(ntop) ' Sizes'],'FontSize',18) ;

      ty = bottom + 0.8 * yrange ;
      text(tx,ty,['Out of ' num2str(vnumresp(itimeb)) ' Total'],'FontSize',18) ;


      orient landscape ;
    print('-dpsc2',pstr) ;



  end ;




elseif ipart == 3 ;    %  then do average Hill Estimates (biggest 2,000)

  if vitimeb ~= 1:21 ;
    disp('!!!   Warning from UNC2001RespSiz1.m:') ;
    disp('!!!   Wrong value of vitimeb,') ;
    disp('!!!   Resetting to correct value') ;
    vitimeb = 1:21 ;
  end ;
      %  need this since are saving Hill estimates to a file


  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  set up preliminaries
  %
  minlog10k = log10(100) ;
  maxlog10k = log10(3400) ;
  vlog10k = linspace(minlog10k,maxlog10k,50) ;
  bottom = 0.3 ;
  top = 2.7 ;

  
  %  loop through time blocks, read big data files, and make plots
  %
  vahillsize = [] ;
  vahilltime = [] ;
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailsize = mtailresp(:,1) ;
    mtailtime = mtailresp(:,4) ;
        %  keep only Response Sizes and durations
    
    mtailsize = flipud(sort(mtailsize)) ;
    mtailtime = flipud(sort(mtailtime)) ;
        %  sort, put in decreasing order


    %  loop through thresholds, k
    %
    valphasize = [] ;
    valphatime = [] ;
    for log10k = vlog10k ;

      k = round(10^log10k) ;

      sizetrunc = mtailsize(1:k) ;
      timetrunc = mtailtime(1:k) ;

      sizetrunc = log(sizetrunc / mtailsize(k+1)) ;
      timetrunc = log(timetrunc / mtailtime(k+1)) ;

      alphasize = 1 / mean(sizetrunc) ;
      alphatime = 1 / mean(timetrunc) ;

      valphasize = [valphasize; alphasize] ;
      valphatime = [valphatime; alphatime] ;

    end ;


    %  make Hill estimator plot
    %
    flagLT2000 = vlog10k < log10(2000) ;
    subplot(2,1,1) ;
    plot(vlog10k',valphasize,'r-','LineWidth',2) ;
      title(['Hill estimates for HTTP Response Size, ' daystrcell{itimeb} ...
                             ', ' timestrcell{itimeb}],'FontSize',18) ;
      xlabel('log_{10}(k)','FontSize',18) ;
      ylabel('\alpha hat','FontSize',18) ;
      axis([minlog10k,maxlog10k,bottom,top]) ;

      avghillsize = mean(valphasize(flagLT2000)) ;
      vahillsize = [vahillsize; avghillsize] ;

      hold on ;
        plot([minlog10k;maxlog10k],[1;1],'k--') ;
        plot([minlog10k;maxlog10k],[2;2],'k--') ;
        plot([log10(2000);log10(2000)],[bottom;top],'k-') ;

        tx = log10(2000) + 0.01 * (maxlog10k - minlog10k) ;
        ty = bottom + 0.9 * (top - bottom) ;
        text(tx,ty,'k = 2000','FontSize',15) ;

        plot([minlog10k;maxlog10k],[avghillsize; avghillsize],'r-') ;
        tx = minlog10k + 0.3 * (maxlog10k - minlog10k) ;
        ty = avghillsize + 0.05 * (top - bottom) ;
        text(tx,ty,['Avg. Hill estimate = ' num2str(avghillsize)], ...
                          'Color','r','FontSize',15) ;

      hold off ;

    subplot(2,1,2) ;
    plot(vlog10k',valphatime,'b-','LineWidth',2) ;
      title(['Hill estimates for HTTP Response Durations, ' daystrcell{itimeb} ...
                             ', ' timestrcell{itimeb}],'FontSize',18) ;
      xlabel('log_{10}(k)','FontSize',18) ;
      ylabel('\alpha hat','FontSize',18) ;
      axis([minlog10k,maxlog10k,bottom,top]) ;

      avghilltime = mean(valphatime(flagLT2000)) ;
      vahilltime = [vahilltime; avghilltime] ;

      hold on ;
        plot([minlog10k;maxlog10k],[1;1],'k--') ;
        plot([minlog10k;maxlog10k],[2;2],'k--') ;
        plot([log10(2000);log10(2000)],[bottom;top],'k-') ;

        tx = log10(2000) + 0.01 * (maxlog10k - minlog10k) ;
        ty = bottom + 0.9 * (top - bottom) ;
        text(tx,ty,'k = 2000','FontSize',15) ;

        plot([minlog10k;maxlog10k],[avghilltime; avghilltime],'b-') ;
        tx = minlog10k + 0.3 * (maxlog10k - minlog10k) ;
        ty = avghilltime + 0.05 * (top - bottom) ;
        text(tx,ty,['Avg. Hill estimate = ' num2str(avghilltime)], ...
                          'Color','b','FontSize',15) ;

      hold off ;


      orient tall ;
    print('-dpsc2',['UNC2001RS1HillEst' num2str(itimeb) '.ps']) ;


  end ;  



  %  save Hill Estimates in a .mat file
  %
  save('DataUNC2001\UNC2001RespSiz1AvgHillEst.mat','vahillsize','vahilltime') ;




elseif ipart == 4 ;    %  then make tail duration vs. size scatterplots (biggest 2,000)

  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Avg. Hill Estimates
  %
  load('DataUNC2001\UNC2001RespSiz1AvgHillEst.mat') ;

  
  %  loop through time blocks, read big data files, and make plots
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailsize = mtailresp(:,1) ;
    mtailtime = mtailresp(:,4) ;
        %  keep only Response Sizes and durations

    mtailsize = mtailsize.^vahillsize(itimeb) ;
    mtailtime = mtailtime.^vahilltime(itimeb) ;
        %  do Hill Estimate rescaling

    mtailsummary = mtailsize.^2 + mtailtime.^2 ;
    [temp,vind] = sort(mtailsummary) ;
    mtailsize = flipud(mtailsize(vind,:)) ;
    mtailtime = flipud(mtailtime(vind,:)) ;
        %  sort in tail summary order, then put in decreasing order

    vx = mtailsize(1:ntop) ;
    vy = mtailtime(1:ntop) ;
        %  work with only ntop values


    %  make main plot
    %
    left = min(vx) ;
    right = max(vx) ;
    xrange = right - left ;
    left = left - 0.01 * xrange ;
    right = right + 0.01 * xrange ;

    bottom = min(vy) ;
    top = max(vy) ;
    yrange = top - bottom ;
    bottom = bottom - 0.01 * yrange ;
    top = top + 0.01 * yrange ;

    plot(vx,vy,'m+') ;
      title(['HTTP Resp. S.P. (marginal avg.  \alpha rescaled), ' daystrcell{itimeb} ...
                             ', ' timestrcell{itimeb}],'FontSize',18) ;
      axis([left,right,bottom,top]) ;
      xlabel(['Response Size (bytes)^{' num2str(vahillsize(itimeb)) '}'], ...
                                   'FontSize',18,'Color','r') ;
      ylabel(['Response Duration (sec)^{' num2str(vahilltime(itimeb)) '}'], ...
                                   'FontSize',18,'Color','b') ;

      tx = left + 0.2 * xrange ;
      ty = bottom + 0.9 * yrange ;
      text(tx,ty,['Largest ' num2str(ntop) ' (w.r.t. L^2)'],'FontSize',18) ;

      ty = bottom + 0.8 * yrange ;
      text(tx,ty,['Out of ' num2str(vnumresp(itimeb)) ' Total'],'FontSize',18) ;


      orient landscape ;
    print('-dpsc2',['UNC2001RS1SPTSavgHill' num2str(itimeb) '.ps']) ;


  end ;




elseif ipart == 5 ;    %  then make tail duration vs. size scatterplots
                       % (biggest 2,000) on square root scale

  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Avg. Hill Estimates
  %
  load('DataUNC2001\UNC2001RespSiz1AvgHillEst.mat') ;

  
  %  loop through time blocks, read big data files, and make plots
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailsize = mtailresp(:,1) ;
    mtailtime = mtailresp(:,4) ;
        %  keep only Response Sizes and durations

    mtailsize = mtailsize.^vahillsize(itimeb) ;
    mtailtime = mtailtime.^vahilltime(itimeb) ;
        %  do Hill Estimate rescaling

    mtailsummary = mtailsize.^2 + mtailtime.^2 ;
    [temp,vind] = sort(mtailsummary) ;
    mtailsize = flipud(mtailsize(vind,:)) ;
    mtailtime = flipud(mtailtime(vind,:)) ;
        %  sort in tail summary order, then put in decreasing order

    vx = mtailsize(1:ntop) ;
    vy = mtailtime(1:ntop) ;
        %  work with only ntop values


    %  make main plot
    %
    left = sqrt(min(vx)) ;
    right = sqrt(max(vx)) ;
    xrange = right - left ;
    left = left - 0.01 * xrange ;
    right = right + 0.01 * xrange ;

    bottom = sqrt(min(vy)) ;
    top = sqrt(max(vy)) ;
    yrange = top - bottom ;
    bottom = bottom - 0.01 * yrange ;
    top = top + 0.01 * yrange ;

    plot(sqrt(vx),sqrt(vy),'m+') ;
      title(['HTTP Resp. S.P. (marginal avg.  \alpha rescaled), ' daystrcell{itimeb} ...
                             ', ' timestrcell{itimeb}],'FontSize',18) ;
      axis([left,right,bottom,top]) ;
      xlabel(['sqrt(Response Size (bytes)^{' num2str(vahillsize(itimeb)) '})'], ...
                                   'FontSize',18,'Color','r') ;
      ylabel(['sqrt(Response Duration (sec)^{' num2str(vahilltime(itimeb)) '})'], ...
                                   'FontSize',18,'Color','b') ;

      tx = left + 0.2 * xrange ;
      ty = bottom + 0.9 * yrange ;
      text(tx,ty,['Largest ' num2str(ntop) ' (w.r.t. L^2)'],'FontSize',18) ;

      ty = bottom + 0.8 * yrange ;
      text(tx,ty,['Out of ' num2str(vnumresp(itimeb)) ' Total'],'FontSize',18) ;


      orient landscape ;
    print('-dpsc2',['UNC2001RS1SPTSavgHillSqrt' num2str(itimeb) '.ps']) ;


  end ;




elseif ipart == 6 ;    %  then do ratio Hill Estimates (biggest 2,000)

  if vitimeb ~= 1:21 ;
    disp('!!!   Warning from UNC2001RespSiz1.m:') ;
    disp('!!!   Wrong value of vitimeb,') ;
    disp('!!!   Resetting to correct value') ;
    vitimeb = 1:21 ;
  end ;
      %  need this since are saving Hill estimates to a file


  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  set up preliminaries
  %
  minlog10k = log10(100) ;
  maxlog10k = log10(3400) ;
  vlog10k = linspace(minlog10k,maxlog10k,50) ;
  bottom = 0.5 ;
  top = 2.5 ;

  
  %  loop through time blocks, read big data files, and make plots
  %
  vahillratio = [] ;
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailsize = mtailresp(:,1) ;
    mtailtime = mtailresp(:,4) ;
        %  keep only Response Sizes and durations
    
    mtailsize = flipud(sort(mtailsize)) ;
    mtailtime = flipud(sort(mtailtime)) ;
        %  sort, put in decreasing order


    %  loop through thresholds, k
    %
    valpharatio = [] ;
    for log10k = vlog10k ;

      k = round(10^log10k) ;

      sizetrunc = mtailsize(1:k) ;
      timetrunc = mtailtime(1:k) ;

      sizetrunc = log(sizetrunc / mtailsize(k+1)) ;
      timetrunc = log(timetrunc / mtailtime(k+1)) ;

      alpharatio = mean(timetrunc) / mean(sizetrunc) ;
          %  estimates ratio of alphasize / alphatime

      valpharatio = [valpharatio; alpharatio] ;

    end ;


    %  make Hill ratio plot
    %
    flagLT2000 = vlog10k < log10(2000) ;

    plot(vlog10k',valpharatio,'m-','LineWidth',2) ;
      title(['Ratio Hill estimates for HTTP Response Durations, ' daystrcell{itimeb} ...
                             ', ' timestrcell{itimeb}],'FontSize',18) ;
      xlabel('log_{10}(k)','FontSize',18) ;
      ylabel('\alpha_{size} / \alpha_{time}','FontSize',18) ;
      axis([minlog10k,maxlog10k,bottom,top]) ;

      avghillratio = mean(valpharatio(flagLT2000)) ;
      vahillratio = [vahillratio; avghillratio] ;

      hold on ;
        plot([minlog10k;maxlog10k],[1;1],'k--') ;
        plot([log10(2000);log10(2000)],[bottom;top],'k-') ;

        tx = log10(2000) + 0.01 * (maxlog10k - minlog10k) ;
        ty = bottom + 0.9 * (top - bottom) ;
        text(tx,ty,'k = 2000','FontSize',15) ;

        plot([minlog10k;maxlog10k],[avghillratio; avghillratio],'m-') ;
        tx = minlog10k + 0.3 * (maxlog10k - minlog10k) ;
        ty = avghillratio + 0.05 * (top - bottom) ;
        text(tx,ty,['Avg. Ratio Hill estimate = ' num2str(avghillratio)], ...
                          'Color','m','FontSize',15) ;

      hold off ;


      orient landscape ;
    print('-dpsc2',['UNC2001RS1HillRat' num2str(itimeb) '.ps']) ;


  end ;  



%  %  save Hill Estimates in a .mat file
%  %
%  save('DataUNC2001\UNC2001RespSiz1AvgHillEst.mat','vahillsize','vahilltime') ;
%
%    Change this if actually use it....




elseif ipart == 11 ;    %  then make Box-Cox transformed version of 
                        %      tail duration vs. size scatterplots (biggest 2,000)

  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;
  
  
  %  loop through time blocks, read big data files, and make plots
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailresp = mtailresp(:,[1 4]) ;
        %  keep only Response Sizes and durations

    [temp,vind] = sort(mtailresp(:,1)) ;
    mtailresp = flipud(mtailresp(vind,:)) ;
        %  sort on first column, put in decreasing order

    vx = mtailresp(1:ntop,1) ;
    vy = mtailresp(1:ntop,2) ;
        %  work with only top values


    %  loop though transformations
    %
    vpower = 0:0.1:1 ;
    for power = vpower ;    %  loop through powers

      if power == 0 ;
        tvx = log(vx) ;
        tvy = log(vy) ;
        tstr = 'log trans. ' ;
      else ;
        tvx = (vx.^power - 1) / power ;
        tvy = (vy.^power - 1) / power ;
        tstr = [num2str(power) ' - power trans. '] ;
      end ;


      %  make main plot
      %
      left = min(tvx) ;
      right = max(tvx) ;
      xrange = right - left ;
      left = left - 0.01 * xrange ;
      right = right + 0.01 * xrange ;

      bottom = min(tvy) ;
      top = max(tvy) ;
      yrange = top - bottom ;
      bottom = bottom - 0.01 * yrange ;
      top = top + 0.01 * yrange ;

      plot(tvx,tvy,'b+') ;
        title(['Box - Cox Transformed HTTP Resp. Scatterplot, ' daystrcell{itimeb} ...
                               ', ' timestrcell{itimeb}],'FontSize',18) ;
        axis([left,right,bottom,top]) ;
        xlabel([tstr 'Size (bytes)'],'FontSize',18) ;
        ylabel([tstr 'Duration (sec)'],'FontSize',18) ;

        tx = left + 0.2 * xrange ;
        ty = bottom + 0.9 * yrange ;
        text(tx,ty,['Largest ' num2str(ntop) ' Sizes'],'FontSize',18) ;

        ty = bottom + 0.8 * yrange ;
        text(tx,ty,['Out of ' num2str(vnumresp(itimeb)) ' Total'],'FontSize',18) ;


      if itimeb == 12 ;
        boxbottom = 1600 ;
        boxleft = 0.6 * 10^7 ;
        
        if power == 0 ;
          tboxbottom = log(boxbottom) ;
          tboxleft = log(boxleft) ;
        else ;
          tboxbottom = (boxbottom.^power - 1) / power ;
          tboxleft = (boxleft.^power - 1) / power ;
        end ;

        hold on ;
          plot([tboxleft;right],[tboxbottom;tboxbottom],'r-') ;
          plot([tboxleft;tboxleft],[tboxbottom;top],'r-') ;
        hold off ;

      end ;


        orient landscape ;
      print('-dpsc2',['UNC2001RS1SPTimeSizeT' num2str(itimeb) ...
                               'L' num2str(10*power) '.ps']) ;


    end ;    %  of ipower loop



  end ;




elseif ipart == 21 ;    %  then plot Log-Log Complementary CDF

  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  Check everything by looping through time blocks again to check maxima
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    vq = mq(:,virowraw(itimeb)) ;
        %  unpack vector of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vq) ;
    vccdf = flipud((1:nq)' / nq) ;


    %  1st lay down all in thin blue
    %
    left = 1 ;
    right = 9 ;
    bottom = -3.01 ;
    top = 0.01 ;
    plot(log10(mq),log10(vccdf),'b-') ;
      axis([left, right, bottom, top]) ;
      title(['log-log CCDF, HTTP R. S., n = ' num2str(n)],'FontSize',18) ;
      xlabel('log_{10}(quantile)','FontSize',18) ;
      ylabel('log_{10}(1 - prob.)','FontSize',18) ;


      %  now add this particular one in thick red
      %
      hold on ;
        plot(log10(vq),log10(vccdf),'r-','LineWidth',2) ;

        vprobplot = [0.2, 0.5, 0.8, 0.9, 0.95 0.99, 0.995, 0.999] ;
        for probplot = vprobplot ;
          ty = log10(1 - probplot + (1/nq)) ;
          [temp,ind2min] = min(abs(vccdf - (1 - probplot + (1/nq)))) ;
          tx = log10(vq(ind2min)) ;
          plot(tx,ty,'k+','MarkerSize',10,'LineWidth',1) ;
          tx = tx + 0.02 * (right - left) ;
          text(tx,ty,[num2str(probplot) ' quantile'], ...
                                       'FontSize',15) ;
        end ;

          tx = left + 0.65 * (right - left) ;
          ty = bottom + 0.9 * (top - bottom) ;
        text(tx,ty,[daystrcell{itimeb} ...
                        ', ' timestrcell{itimeb}], ...
                        'FontSize',18,'Color','r') ;

      hold off ;


      orient landscape ;
    print('-dpsc2',['UNC2001RS1CCDF' num2str(itimeb) '.ps']) ;


  end ;    %  of itimeb loop




elseif ipart == 22 ;    %  then plot slope of Log-Log Complementary CDF

  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  Check everything by looping through time blocks again to check maxima
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    vq = mq(:,virowraw(itimeb)) ;
        %  unpack vector of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vq) ;
    vccdf = flipud((1:nq)' / nq) ;

    sbottom = -4.01 ;
    stop = 0.01 ;

    vcent = mean([log10(vq(2:nq)), log10(vq(1:(nq-1)))],2) ;
    vdenom = (log10(vq(2:nq)) - log10(vq(1:(nq-1)))) ;
    flag0 = vdenom <= 0 ;
    nflag0 = sum(flag0) ;
    if nflag0 > 0 ;
      vdenom(flag0) = ones(nflag0,1) ;
      vslope = (log10(vccdf(2:nq)) - log10(vccdf(1:(nq-1)))) ./ vdenom ;
      vslope(flag0) = sbottom * ones(nflag0,1) ;
    else ;
      vslope = (log10(vccdf(2:nq)) - log10(vccdf(1:(nq-1)))) ./ vdenom ;
    end ;

    subplot(2,1,1) ;
      plot(log10(vq),log10(vccdf),'b-','LineWidth',2) ;
        axis([1, 9, -3.01, 0.01]) ;
        title(['log-log CCDF, HTTP R. S., n = ' num2str(n) ...
                                 ', ' daystrcell{itimeb} ...
                                 ', ' timestrcell{itimeb}],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('log_{10}(1 - prob.)','FontSize',18) ;

    subplot(2,1,2) ;
      plot(vcent,vslope,'r-','LineWidth',2) ;
        axis([1, 9, sbottom, stop]) ;
        title(['Corresponding Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([1; 9],[-1; -1],'k--') ;
          plot([1; 9],[-2; -2],'k--') ;
        hold off ;

      orient landscape ;
    print('-dpsc2',['UNC2001RS1CCDFslope' num2str(itimeb) '.ps']) ;


  end ;    %  of itimeb loop




elseif  ipart == 31  | ...
        ipart == 32  | ...
        ipart == 33  | ...
        ipart == 34  ;     %      Everything based on L2  biggest 2,000
                       %  31  -  Asymptotic Indep. Angle Plot (Time vs. Size)
                       %  32  -  Tail Index Plot for ratios  (Size / Time)
                       %  33  -  Asymptotic Indep. Angle Plot (Ratio vs. Size)
                       %  34  -  Asymptotic Indep. Angle Plot (Time vs. Ratio)

  figure(1) ;
  clf ;
  ntop = 2000 ;
      %  show scatterplot of only largest ntop values


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Avg. Hill Estimates
  %
  load('DataUNC2001\UNC2001RespSiz1AvgHillEst.mat') ;


  if ipart == 32 ;

    if vitimeb ~= 1:21 ;
      disp('!!!   Warning from UNC2001RespSiz1.m:') ;
      disp('!!!   Wrong value of vitimeb,') ;
      disp('!!!   Resetting to correct value') ;
      vitimeb = 1:21 
    end ;
        %  need this since are saving Hill estimates to a file

    %  set up preliminaries
    %
    minlog10k = log10(100) ;
    maxlog10k = log10(ntop-1) ;
    vlog10k = linspace(minlog10k,maxlog10k,50) ;
    bottom = 0 ;
    top = 2.5 ;

    vahillratio = [] ;

  elseif  ipart == 33  | ...
          ipart == 34  ;

    %  Then load ratio Hill Estimates
    %
    load('DataUNC2001\UNC2001RespSiz1AHillEstRat.mat') ;
        %  Hill estimates for Ratio are now in vahillratio


  end ;


  
  %  loop through time blocks, read big data files, and make plots
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    disp(['  For ' daystrcell{itimeb} ', ' timestrcell{itimeb}]) ;

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    vtailsize = mtailresp(:,1) ;
    vtailtime = mtailresp(:,4) ;
        %  keep only Response Sizes and durations

    vtailsize = vtailsize.^vahillsize(itimeb) ;
    vtailtime = vtailtime.^vahilltime(itimeb) ;
        %  do Hill Estimate rescaling

    vtailsize = vtailsize ./ median(vtailsize) ;
    vtailtime = vtailtime ./ median(vtailtime) ;
        %  put on reasonable scales for L2 analysis

    vtailsummary = vtailsize.^2 + vtailtime.^2 ;
    [temp,vind] = sort(vtailsummary) ;
    vtailsize = flipud(vtailsize(vind)) ;
    vtailtime = flipud(vtailtime(vind)) ;
        %  sort in tail summary order, then put in decreasing order

    vx = vtailsize(1:ntop) ;
    vy = vtailtime(1:ntop) ;
        %  work with only ntop values
        
    vx = vx ./ median(vx) ;
    vy = vy ./ median(vy) ;
        %  put on reasonable scales for L2 analysis


    vratio = vx ./ vy ;

    n = vnumresp(itimeb) ;
        %  unpack full sample size



    if  ipart == 31  | ...
        ipart == 33  | ...
        ipart == 34  ;     %  Make Asy. Indep. Angle Plot


      if ipart == 31 ;
        vangle = acos(vx ./ sqrt(vx.^2 + vy.^2)) ;
        xtitstr = 'Size' ;
        ytitstr = 'Time' ;
        pstr = ['UNC2001RS1AIAPtsT' num2str(itimeb) '.ps'] ;
      elseif ipart == 33 ;
        xtitstr = 'Size' ;
        ytitstr = 'Rate' ;
        vratio = vratio.^vahillratio(itimeb) ;
        vratio = vratio / median(vratio) ;
        vangle = acos(vx ./ sqrt(vx.^2 + vratio.^2)) ;
        pstr = ['UNC2001RS1AIAPrsT' num2str(itimeb) '.ps'] ;
      elseif ipart == 34 ;
        xtitstr = 'Rate' ;
        ytitstr = 'Time' ;
        vratio = vratio.^vahillratio(itimeb) ;
        vratio = vratio / median(vratio) ;
        vangle = acos(vratio ./ sqrt(vratio.^2 + vy.^2)) ;
        pstr = ['UNC2001RS1AIAPtrT' num2str(itimeb) '.ps'] ;
      end ;

                       %  31  -  Asymptotic Indep. Angle Plot (Time vs. Size)
                       %  32  -  Tail Index Analysis for Ratio  (Size / Time)
                       %  33  -  Asymptotic Indep. Angle Plot (Ratio vs. Size)
                       %  34  -  Asymptotic Indep. Angle Plot (Time vs. Ratio)



      %  make main plot
      %
      ftitstr = ['HTTP Resp. Asy. Indep. Angles, ' ytitstr ' vs. ' xtitstr ...
                        ', top ' num2str(ntop)] ;
      stitstr = ['SiZer, ' daystrcell{itimeb} ', ' timestrcell{itimeb} ...
                        ', n = ' num2str(n) ...
                        ' (marg. avg.  \alpha rescaled / median)'] ;
      paramstruct = struct('imovie',0,...
                           'famoltitle',ftitstr, ...
                           'sizertitle',stitstr, ...
                           'titlefontsize',18, ...
                           'xlabelstr','Angle * (2/\pi)', ...
                           'labelfontsize',15, ...
                           'ndataoverlay',2, ...
                           'ibdryadj',1, ...
                           'minx',0, ...
                           'maxx',1, ...
                           'fhmin',10^(-2.5), ...
                           'fhmax',10^(-0.5), ...
                           'shmin',10^(-2.5), ...
                           'shmax',10^(-0.5), ...
                           'hhighlight',0, ...
                           'iscreenwrite',1) ;
      sizerSM((2/pi) * vangle,paramstruct) ;
      



        orient landscape ;
      print('-dpsc2',pstr) ;




    elseif  ipart == 32 ;

    vratiosort = sort(vratio) ;
    vratiosort = flipud(vratiosort) ;


      %  loop through thresholds, k
      %
      valpharatio = [] ;
      for log10k = vlog10k ;

        k = round(10^log10k) ;

        ratiotrunc = vratiosort(1:k) ;
        ratiotrunc = log(ratiotrunc / vratiosort(k+1)) ;
        alpharatio = 1 / mean(ratiotrunc) ;

        valpharatio = [valpharatio; alpharatio] ;

      end ;


      %  make Hill estimator plot
      %
      flagLT2000 = vlog10k < log10(2000) ;
      plot(vlog10k',valpharatio,'m-','LineWidth',2) ;
        title(['Hill estimates, HTTP Ratio (Size/Time), ' daystrcell{itimeb} ...
                               ', ' timestrcell{itimeb}],'FontSize',18) ;
        xlabel('log_{10}(k)','FontSize',18) ;
        ylabel('\alpha hat','FontSize',18) ;
        axis([minlog10k,maxlog10k,bottom,top]) ;

        avghillratio = mean(valpharatio(flagLT2000)) ;
        vahillratio = [vahillratio; avghillratio] ;

        hold on ;
          plot([minlog10k;maxlog10k],[1;1],'k--') ;
          plot([minlog10k;maxlog10k],[2;2],'k--') ;

          plot([minlog10k;maxlog10k],[avghillratio; avghillratio],'m-') ;
          tx = minlog10k + 0.3 * (maxlog10k - minlog10k) ;
          ty = avghillratio + 0.05 * (top - bottom) ;
          text(tx,ty,['Avg. Hill estimate = ' num2str(avghillratio)], ...
                            'Color','m','FontSize',15) ;

        hold off ;



        orient landscape ;
      print('-dpsc2',['UNC2001RS1HillEstRat' num2str(itimeb) '.ps']) ;



    end ;



    if ipart == 32 ;

      %  save Hill Estimates in a .mat file
      %
      save('DataUNC2001\UNC2001RespSiz1AHillEstRat.mat','vahillratio') ;

    end ;



  end ;




elseif ipart == 41 ;    %  then plot Log-Log Complementary CDF (full data version)


  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  loop through time blocks to lay down initial blue envelope
  %
  subplot(1,1,1) ;
    left = 1 ;
    right = 9 ;
    bottom = -7 ;
    top = 0.01 ;
  plot([5; 5], [bottom; top],'k--') ;
    axis([left, right, bottom, top]) ;
  vph = [] ;
  vn = [] ;
  for itimeb = vitimeb ;    %  loop through time blocks

    disp(['       plotting for time block ' num2str(itimeb)]) ;

    %  first get part that appeared earlier
    %  (i.e. for given quantiles)
    %
    vq = mq(:,virowraw(itimeb)) ;
        %  unpack vector of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vq) ;
    vccdf = flipud((1/n) + (0:(nq-1))' / nq) ;


    %  second load large data values, and integrate
    %
    filename = filenamecell{itimeb} ;
    mtailresp = load(['DataUNC2001\' filename]) ;
    vtailsize = mtailresp(:,1) ;
    ntail = length(vtailsize) ;

    vq = [vq; sort(vtailsize)] ;
    vccdf = [vccdf; flipud((1:ntail)' / n)] ;
    
    [vq, vindsort] = sort(vq) ;
    vccdf = vccdf(vindsort) ;


    %  lay down this one in thin blue
    %
    hold on ;
      ph = plot(log10(vq),log10(vccdf),'b-') ;
    hold off ;

    vph = [vph; ph] ;
    vn = [vn; n] ;


  end ;    %  of itimeb loop



  %  do part common to all plots
  %
  hold on ;

    tith = title('title string','FontSize',18) ;
    xlabel('log_{10}(quantile)','FontSize',18) ;
    ylabel('log_{10}(1 - prob.)','FontSize',18) ;

      tx = left + 0.65 * (right - left) ;
      ty = bottom + 0.9 * (top - bottom) ;
    texth = text(tx,ty,'text string', ...
                  'FontSize',18,'Color','r') ;

    vprobplot = [0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999, 0.99999] ;
    npp = length(vprobplot) ;
    vpph = [] ;
    vppth = [] ;
    for probplot = vprobplot ;
      ty = log10(1 - probplot) ;
      tx = 5 ;
      pph = plot(tx,ty,'k+','MarkerSize',10,'LineWidth',1) ;
      tx = tx + 0.02 * (right - left) ;
      ppth = text(tx,ty,[num2str(probplot) ' quantile'], ...
                                   'FontSize',15) ;
      vpph = [vpph; pph] ;
      vppth = [vppth; ppth] ;
    end ;

  hold off ;



  %  loop through time blocks again, highlight individuals and print
  %
  vachil = get(gca,'Children') ;
  for itimeb = vitimeb ;    %  loop through time blocks


    %  set appropriate values for this time block
    %
      titstr = ['log-log CCDF, HTTP R. S., n = ' num2str(vn(itimeb))] ;
    set(tith,'String',titstr) ;

    set(texth,'String',[daystrcell{itimeb} ...
                        ', ' timestrcell{itimeb}]) ;

    vxdat = get(vph(itimeb),'XData') ;
    vydat = get(vph(itimeb),'YData') ;

    hold on ;
      newph = plot(vxdat,vydat,'r-','LineWidth',2) ;
          %  replot, to highlith
    hold off ;

    for ipp = 1:npp ;
      ppy = get(vpph(ipp),'YData') ;
      [temp,ind2min] = min(abs(vydat - ppy)) ;
      ppx = vxdat(ind2min) ;
      set(vpph(ipp),'XData',ppx) ;
        vpos = get(vppth(ipp),'Position') ;
        vpos(1) = ppx + 0.02 * (right - left) ;
      set(vppth(ipp),'Position',vpos) ;
    end ;



      orient landscape ;
    print('-dpsc2',['UNC2001RS1allCCDFfull' num2str(itimeb) '.ps']) ;


    %  Reset, i.e. eliminate new highlighted curve
    %
    delete(newph) ;


  end ;    %  of itimeb loop




elseif ipart == 42 ;    %  then plot slope of Log-Log Complementary CDF (full data version)


  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  loop through time blocks
  %
  for itimeb = vitimeb ;    %  loop through time blocks


    %  first get part that appeared earlier
    %  (i.e. for given quantiles)
    %
    vq = mq(:,virowraw(itimeb)) ;
        %  unpack vector of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vq) ;
    vccdf = flipud((1/n) + (0:(nq-1))' / nq) ;


    %  load large data values, and integrate
    %
    filename = filenamecell{itimeb} ;
    mtailresp = load(['DataUNC2001\' filename]) ;
    vtailsize = mtailresp(:,1) ;
    ntail = length(vtailsize) ;

    vq = [vq; sort(vtailsize)] ;
    vccdf = [vccdf; flipud((1:ntail)' / n)] ;
    
    [vq, vindsort] = sort(vq) ;
    vccdf = vccdf(vindsort) ;
    nq = length(vq) ;
    
    flagsame = vq(1:(nq-1)) == vq(2:nq) ;
        %  one where next entry is the same
    vq = vq(~flagsame) ;
    vccdf = vccdf(~flagsame) ;
        %  keep only entries that are different (always keeping last one)
    nq = length(vq) ;


    %  compute direct difference quotients
    %
    sbottom = -4.01 ;
    stop = 0.01 ;

    vcent = mean([log10(vq(2:nq)), log10(vq(1:(nq-1)))],2) ;
    vdenom = (log10(vq(2:nq)) - log10(vq(1:(nq-1)))) ;
    flag0 = vdenom <= 0 ;
    nflag0 = sum(flag0) ;
    if nflag0 > 0 ;
      vdenom(flag0) = ones(nflag0,1) ;
      vslope = (log10(vccdf(2:nq)) - log10(vccdf(1:(nq-1)))) ./ vdenom ;
      vslope(flag0) = sbottom * ones(nflag0,1) ;
    else ;
      vslope = (log10(vccdf(2:nq)) - log10(vccdf(1:(nq-1)))) ./ vdenom ;
    end ;

    left = 1 ;
    right = 9 ;
    bottom = -7 ;
    top = 0.01 ;


    %  Make gridded difference quotients
    %
    del1 = 0.3 ;
    del2 = 0.03 ;
        %  grid widths
    vgrid1 = (left:del1:right)' ;
    vgrid2 = (left:del2:right)' ;
        %  grid points covering full range
    minq = log10(min(vq)) ;
    maxq = log10(max(vq)) ;
        %  endpoints where have data beyond
    vgrid1 = vgrid1(minq < vgrid1) ;
    vgrid2 = vgrid2(minq < vgrid2) ;
    vgrid1 = vgrid1(vgrid1 < maxq) ;
    vgrid2 = vgrid2(vgrid2 < maxq) ;
        %  reduce to grid within the data
    vheightgrid1 = interp1(log10(vq),log10(vccdf),vgrid1) ;
    vheightgrid2 = interp1(log10(vq),log10(vccdf),vgrid2) ;
        %  interpolate to grid
    ngrid1 = length(vheightgrid1) ;
    ngrid2 = length(vheightgrid2) ;
    vcent1 = mean([vgrid1(2:ngrid1), vgrid1(1:(ngrid1-1))],2) ;
    vcent2 = mean([vgrid2(2:ngrid2), vgrid2(1:(ngrid2-1))],2) ;
        %  center points for plotting
    vslope1 = (vheightgrid1(2:ngrid1) - vheightgrid1(1:(ngrid1-1))) ./ del1 ;
    vslope2 = (vheightgrid2(2:ngrid2) - vheightgrid2(1:(ngrid2-1))) ./ del2 ;
        %  get difference quotients (slopes)



    %  Make graphics
    %
    subplot(2,2,1) ;
      plot(log10(vq),log10(vccdf),'b-','LineWidth',2) ;
        axis([left, right, bottom, top]) ;
        title(['log-log CCDF, HTTP R. S., n = ' num2str(n)],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('log_{10}(1 - prob.)','FontSize',18) ;
          tx = left + 0.1 * (right - left) ;
          ty = bottom + 0.1 * (top - bottom) ;
        text(tx,ty,[daystrcell{itimeb} ...
                    ', ' timestrcell{itimeb}],'FontSize',18) ;

        vprobplot = [0.5, 0.9, 0.999, 0.99999] ;
        vlq = [] ;
        for probplot = vprobplot ;
          ty = log10(1 - probplot) ;
          [temp,ind2min] = min(abs(vccdf - (1 - probplot))) ;
          tx = log10(vq(ind2min)) ;
          vlq = [vlq; tx] ;
          hold on ;
            plot([tx; tx],[bottom,top],'g--') ;
            tx = tx + 0.02 * (right - left) ;
            text(tx,ty,[num2str(probplot) ' quantile'], ...
                                  'FontSize',15,'Color','g') ;
          hold off ;
        end ;
        nlq = length(vlq) ;

    subplot(2,2,3) ;
      plot(vcent,vslope,'r-','LineWidth',2) ;
        axis([left,right, sbottom, stop]) ;
        title(['Corresponding Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([left; right],[-1; -1],'k--') ;
          plot([left; right],[-2; -2],'k--') ;
          plot([vlq'; vlq'],[sbottom * ones(1,nlq); stop * ones(1,nlq)],'g--') ;
        hold off ;

    subplot(2,2,2) ;
      plot(vcent1,vslope1,'m-','LineWidth',2) ;
        axis([left,right, sbottom, stop]) ;
        title(['\Delta = ' num2str(del1) ' Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([left; right],[-1; -1],'k--') ;
          plot([left; right],[-2; -2],'k--') ;
          plot([vlq'; vlq'],[sbottom * ones(1,nlq); stop * ones(1,nlq)],'g--') ;
        hold off ;

    subplot(2,2,4) ;
      plot(vcent2,vslope2,'m-','LineWidth',2) ;
        axis([left,right, sbottom, stop]) ;
        title(['\Delta = ' num2str(del2) ' Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([left; right],[-1; -1],'k--') ;
          plot([left; right],[-2; -2],'k--') ;
          plot([vlq'; vlq'],[sbottom * ones(1,nlq); stop * ones(1,nlq)],'g--') ;
        hold off ;

      orient landscape ;
    print('-dpsc2',['UNC2001RS1allCCDFSfull' num2str(itimeb) '.ps']) ;


  end ;    %  of itimeb loop




elseif  ipart == 43  | ...
        ipart == 44  ;    %  then do Q-Q plot (full data version)


  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  loop through time blocks
  %
  for itimeb = vitimeb ;    %  loop through time blocks


    %  first get part that appeared earlier
    %  (i.e. for given quantiles)
    %
    vqdata = mq(:,virowraw(itimeb)) ;
        %  unpack vector of small fixed set of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vqdata) ;
    vccdf = flipud((1/n) + (0:(nq-1))' / nq) ;


    %  load large data values, and integrate
    %
    filename = filenamecell{itimeb} ;
    mtailresp = load(['DataUNC2001\' filename]) ;
    vtailsize = mtailresp(:,1) ;
    ntail = length(vtailsize) ;

    qdata = [vqdata; sort(vtailsize)] ;
    vccdf = [vccdf; flipud((1:ntail)' / n)] ;
    
    [qdata, vindsort] = sort(qdata) ;
    vccdf = vccdf(vindsort) ;
    nq = length(qdata) ;
    

    %  Use lines copied and modified from qqSM.m
    %
    if ipart == 43 ;    %  Pareto
%      vqalign = [0.99; 0.99999] ;
%      vqalign = [0.9; 0.9999] ;
%      vqalign = [0.9; 0.999] ;
%      vqalign = [0.9; 0.99] ;
%      vqalign = [0.5; 0.99] ;
      vqalign = [0.8; 0.99] ;
      titlestr = ['HTTP R. S., Pareto Q-Q, ' num2str(vqalign(1)) ...
                       ' & ' num2str(vqalign(2)) ... 
                       ' quant. matched, n = ' num2str(n)] ;
      diststr = 'Pareto-log' ;
    elseif ipart == 44 ;    %  log-normal
      vqalign = [0.8; 0.99] ;
      titlestr = ['HTTP R. S., log-normal Q-Q, ' num2str(vqalign(1)) ...
                       ' & ' num2str(vqalign(2)) ... 
                       ' quant. matched, n = ' num2str(n)] ;
      diststr = 'log-normal' ;
    end ;
    maxstep = 1000 ;
    relaccthreshold = 10^(-6) ;
    colorcell = {'r' 'g' 'b'} ;
            %  1st - data
            %  2nd - line
            %  3rd - simulated versions
    ltypestr = '-' ;
    titlefontsize = 18 ;
    xlabelstr = ', Q' ;
    xlabelstr = [diststr, xlabelstr] ;
    ylabelstr = 'Data Q' ;
    ylabelstr = ['log ' ylabelstr] ;
    labelfontsize = 18 ;
    parfontsize = 18 ;
    vshowq = [0.1, 0.5, 0.9, 0.99, 0.999, 0.9999, 0.99999] ;
    nsim = 100 ;
    nsimplotval = 900 ;


    if ipart == 43 ;    %  estimate Pareto parameters

      p1 = vqalign(1) ;
      p2 = vqalign(2) ;
      vq = cquantSM(vqdata,vqalign,0) ;
           %  0 to indicate presorted data
           %  (only work with small fixed quantiles)
      q1hat = vq(1) ;
      q2hat = vq(2) ;

      numera = log(1 - p1) - log(1 - p2);
      oldsigma = 0 ;
            %  common and starting values

      for istep = 1:maxstep ;

        q1z = q1hat + oldsigma ;
        q2z = q2hat + oldsigma ;

        alpha = numera / (log(q2z) - log(q1z)) ;
        sigma = q1z * (1 - p1)^(1/alpha) ;

        relacc = abs(sigma - oldsigma) / sigma ;

        disp(['     For iteration step ' num2str(istep) ...
                    ', relacc = ' num2str(relacc) ...
                    ', alpha = ' num2str(alpha) ...
                    ', sigma = ' num2str(sigma)]) ;

      
        if relacc < relaccthreshold ;
          errflag = 0 ;
          disp(['     For iteration step ' num2str(istep) ...
                      ', relacc = ' num2str(relacc) ...
                      ', alpha = ' num2str(alpha) ...
                      ', sigma = ' num2str(sigma)]) ;
          break ;
        else ;
          errflag = 1 ;
          oldsigma = sigma ;
        end ;

      end ;

      if errflag ~= 0 ;
        disp('!!!   Warning:  Pareto fit may be unstable   !!!') ;
      end ;

      ididqalign = 1 ;


      %  for these estimated parameters, get theoretical quantiles
      %
      qtheory = (vccdf.^(-1/alpha) - 1) * sigma ;
      pgrid = (1:n)' / (n + 1) ;
      qtheoryall = ((1 - pgrid).^(-1/alpha) - 1) * sigma ;


      %  put on log scale
      %
      qdata = log(qdata) ;
      qtheory = log(qtheory) ;
        q1hat = log(q1hat) ;
        q2hat = log(q2hat) ;
      qtheoryall = log(qtheoryall) ;


      par1str = '\alpha' ;
      par1val = alpha ;
      par2str = '\sigma' ;
      par2val = sigma ;



    elseif ipart == 44 ;    %  estimate log-normal parameters

      qdata = log(qdata) ;
      vqdata = log(vqdata) ;

      p1 = vqalign(1) ;
      p2 = vqalign(2) ;
      vq = cquantSM(vqdata,vqalign,0) ;
           %  0 to indicate presorted data
           %  (only work with small fixed quantiles)
      q1hat = vq(1) ;
      q2hat = vq(2) ;

      sigma = (q1hat - q2hat) / (norminv(p1) - norminv(p2)) ;
      mu = q1hat - sigma * norminv(p1) ;


      %  for these estimated parameters, get theoretical quantiles
      %
      qtheory = mu + sigma * norminv(1 - vccdf) ;
      pgrid = (1:n)' / (n + 1) ;
      qtheoryall = mu + sigma * norminv(pgrid) ;


      par1str = '\mu' ;
      par1val = mu ;
      par2str = '\sigma' ;
      par2val = sigma ;


    end ;



    left = qtheory(1) ; 
    right = qtheory(nq) ; 
    bottom = qdata(1) ;
    top = qdata(nq) ;


    %  Make main graphic in current axes
    %
    plot(qtheory,qdata, ...
        [colorcell{1} '-'], ...
        'LineWidth',3) ;
      axis equal ;
      axis([left,right,bottom,top]) ;

      th = title(titlestr) ;
      set(th,'FontSize',titlefontsize) ;

      xlh = xlabel(xlabelstr) ;    
      ylh = ylabel(ylabelstr) ;    
      set(xlh,'FontSize',labelfontsize) ;
      set(ylh,'FontSize',labelfontsize) ;


    %    overlay 45 degree line
    %
    minmin = min(left,bottom) ;
    maxmax = max(right,top) ;
    hold on ;
      plot([minmin,maxmax],[minmin,maxmax], ...
           [colorcell{2} ltypestr], ...
           'LineWidth',2) ;
    hold off ;


    %  Add simulated realizations (if needed)
    %
    simseed = 34875094 ;
    rand('seed',simseed) ;
        %  set seed

    for isim = 1:nsim ;    %  loop through simulation steps

      disp(['      working on simulated data set ' num2str(isim) ...
                                            ' of ' num2str(nsim)]) ;

      if ipart == 43 ;    %  Pareto

        vsimdat = sigma * ((1 - rand(n,1)).^(-1/alpha) - 1) ;
        vsimdat = log(vsimdat) ;

      elseif ipart == 44 ;    %  log-normal

        vsimdat = mu + sigma * randn(n,1) ;

      end;


      vsimdat = sort(vsimdat) ;
          %  sort this column

      nspvo3 = floor(nsimplotval / 3) ;
            %  one third of total, to put at each end
      vindleft = (1:nspvo3)' ;
            %  left end  (include all)
      vindright = ((n-nspvo3+1):n)' ;
            %  right end  (include all)

      nspvlo = nsimplotval - length(vindleft) ...
                           - length(vindright) ;
            %  number of grid points left over (for use in middle)
      vindmid = round(linspace(nspvo3+1,n-nspvo3,nspvlo)') ;

      vind = [vindleft; vindmid; vindright] ;
            %  vector of indices, full

      qtheoryp = qtheoryall(vind,:) ;
      vsimdatp = vsimdat(vind,:) ;

      hold on  ;
        plot(qtheoryp,vsimdatp,'-','Color',colorcell{3}) ;
      hold off ;

    end ;    %  of loop through simulated data sets


    %  replot important stuff
    %
    hold on ;
      plot(qtheory,qdata, ...
           [colorcell{1} '-'], ...
           'LineWidth',3) ;
      plot([minmin,maxmax],[minmin,maxmax], ...
           [colorcell{2} ltypestr], ...
           'LineWidth',2) ;
    hold off ;


    %  Add text for parameters and time block
    %
      tx = left + 0.75 * (right - left) ;
      ty = bottom + 0.25 * (top - bottom) ;
    th = text(tx,ty,[par1str ' = ' num2str(par1val)]) ;
    set(th,'FontSize',parfontsize) ;

      ty = bottom + 0.1 * (top - bottom) ;
    th = text(tx,ty,[par2str ' = ' num2str(par2val)]) ;
    set(th,'FontSize',parfontsize) ;

      tx = left + 0.1 * (right - left) ;
      ty = bottom + 0.9 * (top - bottom) ;
    th = text(tx,ty,[daystrcell{itimeb} ', ' timestrcell{itimeb}]) ;
    set(th,'FontSize',parfontsize) ;


    %  show fit quantiles
    %
    hold on ;
      sch = plot([q1hat; q2hat], [q1hat; q2hat], 'ko') ;
    hold off ;


    %  show additional quantiles
    %
    for isq = 1:length(vshowq) ;
  
      sq = vshowq(isq) ;
    
      if  0 < sq   &  sq < 1  ;    %  then have valid probability

        [temp,qi] = min(abs((1 - vccdf) - sq)) ;
            %  gets index of this quantile
            %  i.e. where prob is closest to sq prob

        qehat = qdata(qi) ;
        qthat = qtheory(qi) ;

        qhx = qthat + 0.02 * (right - left) ;

        hold on ;
          plot([qthat], [qehat], 'k+') ;
          th = text(qhx,qehat,[num2str(sq) ' quantile']) ;
          if ~isempty(parfontsize) ;
            set(th,'FontSize',parfontsize) ;
          end ;
        hold off ;

      else ;    %  then don't have valid probability, give error message

        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
        disp('!!!   Warning from qqSM:                       !!!') ;
        disp('!!!   Invalid entry (not in (0,1)) in vshowq   !!!') ;
        disp('!!!   So won''t show this quantile              !!!') ;
        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

      end ;    %  of sq if-block  


    end ;    %    of isq loop through quantiles to show




      orient landscape ;
    if ipart == 43 ;
      print('-dpsc2',['UNC2001RS1allQQpar' num2str(itimeb) '.ps']) ;
    elseif ipart == 44 ;
      print('-dpsc2',['UNC2001RS1allQQln' num2str(itimeb) '.ps']) ;
    end ;


  end ;    %  of itimeb loop


elseif  ipart == 51  | ...
        ipart == 52  | ...
        ipart == 53  | ...
        ipart == 54  | ...
        ipart == 55  | ...
        ipart == 61  | ...
        ipart == 62  | ...
        ipart == 63  | ...
        ipart == 64  | ...
        ipart == 65  ;    %  do visual fits to Full Data Log-Log CCDF


  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  set parameters
  %
  if  ipart == 51  | ...
      ipart == 61  ;    %  Double-Pareto log-normal (fit by eye)

    diststr = 'Double Pareto log-normal' ;
    outdistr = 'dpln' ;
    if ipart == 61 ;
      outdistr = [outdistr 'sim'] ;
    end ;

    alpha = 1.2 ;    %  "tail parameter" near origin
    beta = 1.2 ;    %  tail parameter
    nu = 7 ;    %  mean of log-normal
    tau = 1.0 ;    %  sd of log-normal
        %  these chosen by eye, for good fit in ipart = 10, Thurday morning

  elseif  ipart == 52  | ...
          ipart == 62  ;    %  log-normal 2-mixture (fit by eye)

    diststr = 'log-normal 2-mix' ;
    outdistr = 'ln2m' ;
    if ipart == 62 ;
      outdistr = [outdistr 'sim'] ;
    end ;

    mu1 = 7.5 ;
    sig1 = 1.4 ;
    mu2 = 13.05 ;
    sig2 = 1.55 ;
    w1 = .9988 ;
    w2 = 1 - w1 ;
        %  these chosen by eye, for good fit in ipart = 10, Thurday morning

  elseif  ipart == 53  | ...
          ipart == 63 ;    %  Double-Pareto log-normal 2-mixture (fit by eye)

    diststr = 'D-Par log-N 2-mix' ;
    outdistr = 'dpln2m' ;
    if ipart == 63 ;
      outdistr = [outdistr 'sim'] ;
    end ;

    alpha1 = 2 ;    %  "tail parameter" near origin
    beta1 = 2 ;    %  tail parameter
    nu1 = 7.8 ;    %  mean of log-normal
    tau1 = 1.0 ;    %  sd of log-normal
    alpha2 = 1.2 ;    %  "tail parameter" near origin
    beta2 = 1.2 ;    %  tail parameter
    nu2 = 13.5 ;    %  mean of log-normal
    tau2 = 0.8 ;    %  sd of log-normal
    w1 = .9992 ;
    w2 = 1 - w1 ;
        %  these chosen by eye, for good fit in ipart = 10, Thurday morning

  elseif  ipart == 54  | ...
          ipart == 64  ;    %  log-normal 3-mixture (fit by eye)

    diststr = 'log-normal 3-mix' ;
    outdistr = 'ln3m' ;
    if ipart == 64 ;
      outdistr = [outdistr 'sim'] ;
    end ;

    mu1 = 5.7 ;
    sig1 = 0.6 ;
    mu2 = 8.45 ;
    sig2 = 1.2 ;
    mu3 = 13.05 ;
    sig3 = 1.55 ;
    w1 = 0.55 ;
    w2 = 0.4488 ;
    w3 = 1 - w1 - w2 ;
        %  these chosen by eye, for good fit in ipart = 10, Thurday morning

  elseif  ipart == 55  | ...
          ipart == 65  ;    %  Double-Pareto log-normal 3-mixture (fit by eye)

    diststr = 'D-Par log-N 3-mix' ;
    outdistr = 'dpln3m' ;
    if ipart == 65 ;
      outdistr = [outdistr 'sim'] ;
    end ;

    alpha1 = 2 ;    %  "tail parameter" near origin
    beta1 = 2 ;    %  tail parameter
    nu1 = 5.7 ;    %  mean of log-normal
    tau1 = 0.6 ;    %  sd of log-normal
    alpha2 = 10 ;    %  "tail parameter" near origin
    beta2 = 10 ;    %  tail parameter
    nu2 = 8.45 ;    %  mean of log-normal
    tau2 = 1.2 ;    %  sd of log-normal
    alpha3 = 1.3 ;    %  "tail parameter" near origin
    beta3 = 1.3 ;    %  tail parameter
    nu3 = 13.35 ;    %  mean of log-normal
    tau3 = 0.75 ;    %  sd of log-normal
    w1 = 0.6 ;
    w2 = 0.399 ;
    w3 = 1 - w1 - w2 ;
        %  these chosen by eye, for good fit in ipart = 10, Thurday morning
  end ;


  %  loop through time blocks and plot CCDF
  %
  subplot(1,1,1) ;
    left = 1 ;
    right = 9 ;
    bottom = -7 ;
    top = 0.01 ;
  for itimeb = vitimeb ;    %  loop through time blocks

    disp(['       plotting for time block ' num2str(itimeb)]) ;

    %  first get part that appeared earlier
    %  (i.e. for given quantiles)
    %
    vq = mq(:,virowraw(itimeb)) ;
        %  unpack vector of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vq) ;
    vccdf = flipud((1/n) + (0:(nq-1))' / nq) ;


    %  second load large data values, and integrate
    %
    filename = filenamecell{itimeb} ;
    mtailresp = load(['DataUNC2001\' filename]) ;
    vtailsize = mtailresp(:,1) ;
    ntail = length(vtailsize) ;

    vq = [vq; sort(vtailsize)] ;
    vccdf = [vccdf; flipud((1:ntail)' / n)] ;
    
    [vq, vindsort] = sort(vq) ;
    vccdf = vccdf(vindsort) ;


    %  Plot log-log CCDF of Data
    %
    vxdat = log10(vq) ;
    vydat = log10(vccdf) ;
    ph = plot(vxdat,vydat,'b-','LineWidth',2) ;
      axis([left, right, bottom, top]) ;
      titstr = ['log-log CCDF, HTTP R. S., n = ' ...
                         num2str(n) ' & ' diststr] ;
      tith = title(titstr,'FontSize',18) ;
      xlabel('log_{10}(quantile)','FontSize',18) ;
      ylabel('log_{10}(1 - prob.)','FontSize',18) ;


    hold on ;

        tx = left + 0.65 * (right - left) ;
        ty = bottom + 0.9 * (top - bottom) ;
        tstr = [daystrcell{itimeb} ', ' timestrcell{itimeb}] ;
      texth = text(tx,ty,tstr,'FontSize',18,'Color','r') ;

      vprobplot = [0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999, 0.99999] ;
      npp = length(vprobplot) ;
      for probplot = vprobplot ;
        ty = log10(1 - probplot) ;
        [temp,ind2min] = min(abs(vydat - ty)) ;
        tx = vxdat(ind2min) ;
        pph = plot(tx,ty,'k+','MarkerSize',10,'LineWidth',1) ;
        tx = tx + 0.02 * (right - left) ;
        ppth = text(tx,ty,[num2str(probplot) ' quantile'], ...
                                     'FontSize',15) ;
      end ;

    hold off ;



    %  add simulated data
    %
    if  ipart == 61  | ...
        ipart == 62  | ...
        ipart == 63  | ...
        ipart == 64  | ...
        ipart == 65  ;

      %  set up preliminaries
      %
      nsim = 100 ;

      vindkeep = round(linspace(n/(2 * 1000), ...
                          n * (1 - 1/(2 * 1000)),1000)') ;
          %  essentially 1000 equally spaced points
      vindkeep = [vindkeep; ((n - 1000 + 1):n)'] ;
          %  add in last 1000 data points
      vindkeep = sort(vindkeep) ;
          %  put in increasing order

      pgrid = linspace(1 / (2 * n), (1 - 1 / (2 * n)), n)' ;
      pgrid = 1 - pgrid ;
      pgrid = pgrid(vindkeep) ;
      


      for isim = 1:nsim ;

        disp(['      working on simulated data set ' num2str(isim) ...
                                            ' of ' num2str(nsim)]) ;

        %  Generate curve
        %
        if ipart == 61  ;    %  Double-Pareto log-normal (fit by eye)

          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag = vunif < (alpha ./ (alpha + beta)) ;
              %  one with probability proportional to alpha
          vunif = rand(n,1) ;
              %  Unif(0,1)
          vu = flag .* vunif .^ (1 ./ beta) + ...
               (~flag) .* vunif .^ (-1 ./ alpha) ;

          vv = nu + tau .* randn(n,1) ;
          vv = exp(vv) ;
          
          vsimdat = vu .* vv ;

        elseif ipart == 62 ;    %  log-normal 2-mixture (fit by eye)

          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag1 = vunif < w1 ;
              %  one with probability w1
          flag2 = ~flag1 ;
              %  one when not flag1
          vmu = flag1 * mu1 + flag2 * mu2 ;
          vsig = flag1 * sig1 + flag2 * sig2 ;
          vsimdat = vmu + vsig .* randn(n,1) ;
          vsimdat = exp(vsimdat) ;

        elseif ipart == 63  ;    %  Double-Pareto log-normal 2-mixture (fit by eye)

          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag1 = vunif < w1 ;
              %  one with probability w1
          flag2 = ~flag1 ;
              %  one when not flag1
          valpha = flag1 * alpha1 + flag2 * alpha2 ;
          vbeta = flag1 * beta1 + flag2 * beta2 ;
          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag = vunif < (valpha ./ (valpha + vbeta)) ;
              %  one with probability proportional to alpha
          vunif = rand(n,1) ;
              %  Unif(0,1)
          vu = flag .* vunif .^ (1 ./ vbeta) + ...
               (~flag) .* vunif .^ (-1 ./ valpha) ;

          vnu = flag1 * nu1 + flag2 * nu2 ;
          vtau = flag1 * tau1 + flag2 * tau2 ;
          vv = vnu + vtau .* randn(n,1) ;
          vv = exp(vv) ;
          
          vsimdat = vu .* vv ;

        elseif ipart == 64 ;    %  log-normal 3-mixture (fit by eye)

          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag1 = vunif < w1 ;
              %  one with probability w1
          flag3 = (w1 + w2) < vunif ;
              %  one with probability w3 = 1 - (w1 + w2)
          flag2 = ~(flag1 | flag3) ;
              %  one when neither flag1 nor flag3 are
          vmu = flag1 * mu1 + flag2 * mu2 + flag3 * mu3 ;
          vsig = flag1 * sig1 + flag2 * sig2 + flag3 * sig3 ;
          vsimdat = vmu + vsig .* randn(n,1) ;
          vsimdat = exp(vsimdat) ;

        elseif ipart == 65  ;    %  Double-Pareto log-normal 3-mixture (fit by eye)

          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag1 = vunif < w1 ;
              %  one with probability w1
          flag3 = (w1 + w2) < vunif ;
              %  one with probability w3 = 1 - (w1 + w2)
          flag2 = ~(flag1 | flag3) ;
              %  one when neither flag1 nor flag3 are
          valpha = flag1 * alpha1 + flag2 * alpha2 + flag3 * alpha3 ;
          vbeta = flag1 * beta1 + flag2 * beta2 + flag3 * beta3 ;
          vunif = rand(n,1) ;
              %  Unif(0,1)
          flag = vunif < (valpha ./ (valpha + vbeta)) ;
              %  one with probability proportional to alpha
          vunif = rand(n,1) ;
              %  Unif(0,1)
          vu = flag .* vunif .^ (1 ./ vbeta) + ...
               (~flag) .* vunif .^ (-1 ./ valpha) ;

          vnu = flag1 * nu1 + flag2 * nu2 + flag3 * nu3 ;
          vtau = flag1 * tau1 + flag2 * tau2 + flag3 * tau3 ;
          vv = vnu + vtau .* randn(n,1) ;
          vv = exp(vv) ;
          
          vsimdat = vu .* vv ;

        end ;

        vsimdat = sort(vsimdat) ;
            %  sorted version

        vsimdat = vsimdat(vindkeep) ;


        %  Plot this curve
        %
        hold on ;
          plot(log10(vsimdat),log10(pgrid),'m-') ;
        hold off ;

      end ;    %  of isim loop



      %  replot most important parts
      %
      hold on ;
        plot(vxdat,vydat,'b-','LineWidth',2) ;
      hold off ;


    end ;    %  of ipart if-block




    %  overlay theoretical distributions
    %
    ngrid = 400 ;
    log10qgrid = linspace(left,right,ngrid) ;
    qgrid = 10.^log10qgrid ;
    if  ipart == 51  | ...
        ipart == 61  ;    %  Double-Pareto log-normal (fit by eye)

      cdf = normcdf((log(qgrid) - nu) / tau) - ...
              (beta / (alpha + beta)) .* ...
                (qgrid.^(-alpha)  .* ...
                    exp(alpha * nu + alpha^2 * tau^2 / 2) .* ...
                    normcdf((log(qgrid) - nu - alpha * tau^2) / tau)) + ...
              (alpha / (alpha + beta)) .* ...
                (qgrid.^beta  .* ...
                    exp(-beta * nu + beta^2 * tau^2 / 2) .* ...
                    (1 - normcdf((log(qgrid) - nu + beta * tau^2) / tau))) ;

    elseif  ipart == 52  | ...
            ipart == 62  ;    %  log-normal 2-mixture (fit by eye)

      cdf1 = normcdf((log(qgrid) - mu1) / sig1) ;
      cdf2 = normcdf((log(qgrid) - mu2) / sig2) ;
      cdf = w1 * cdf1 + w2 * cdf2 ;

    elseif  ipart == 53  | ...
            ipart == 63  ;    %  Double-Pareto log-normal 2-mixture (fit by eye)

      cdf1 = normcdf((log(qgrid) - nu1) / tau1) - ...
              (beta1 / (alpha1 + beta1)) .* ...
                (qgrid.^(-alpha1)  .* ...
                    exp(alpha1 * nu1 + alpha1^2 * tau1^2 / 2) .* ...
                    normcdf((log(qgrid) - nu1 - alpha1 * tau1^2) / tau1)) + ...
              (alpha1 / (alpha1 + beta1)) .* ...
                (qgrid.^beta1  .* ...
                    exp(-beta1 * nu1 + beta1^2 * tau1^2 / 2) .* ...
                    (1 - normcdf((log(qgrid) - nu1 + beta1 * tau1^2) / tau1))) ;
      cdf2 = normcdf((log(qgrid) - nu2) / tau2) - ...
              (beta2 / (alpha2 + beta2)) .* ...
                (qgrid.^(-alpha2)  .* ...
                    exp(alpha2 * nu2 + alpha2^2 * tau2^2 / 2) .* ...
                    normcdf((log(qgrid) - nu2 - alpha2 * tau2^2) / tau2)) + ...
              (alpha2 / (alpha2 + beta2)) .* ...
                (qgrid.^beta2  .* ...
                    exp(-beta2 * nu2 + beta2^2 * tau2^2 / 2) .* ...
                    (1 - normcdf((log(qgrid) - nu2 + beta2 * tau2^2) / tau2))) ;
      cdf = w1 * cdf1 + w2 * cdf2 ;

    elseif  ipart == 54  | ...
            ipart == 64  ;    %  log-normal 3-mixture (fit by eye)

      cdf1 = normcdf((log(qgrid) - mu1) / sig1) ;
      cdf2 = normcdf((log(qgrid) - mu2) / sig2) ;
      cdf3 = normcdf((log(qgrid) - mu3) / sig3) ;
      cdf = w1 * cdf1 + w2 * cdf2 + w3 * cdf3 ;

    elseif  ipart == 55  | ...
            ipart == 65  ;    %  Double-Pareto log-normal 3-mixture (fit by eye)

      cdf1 = normcdf((log(qgrid) - nu1) / tau1) - ...
              (beta1 / (alpha1 + beta1)) .* ...
                (qgrid.^(-alpha1)  .* ...
                    exp(alpha1 * nu1 + alpha1^2 * tau1^2 / 2) .* ...
                    normcdf((log(qgrid) - nu1 - alpha1 * tau1^2) / tau1)) + ...
              (alpha1 / (alpha1 + beta1)) .* ...
                (qgrid.^beta1  .* ...
                    exp(-beta1 * nu1 + beta1^2 * tau1^2 / 2) .* ...
                    (1 - normcdf((log(qgrid) - nu1 + beta1 * tau1^2) / tau1))) ;
      cdf2 = normcdf((log(qgrid) - nu2) / tau2) - ...
              (beta2 / (alpha2 + beta2)) .* ...
                (qgrid.^(-alpha2)  .* ...
                    exp(alpha2 * nu2 + alpha2^2 * tau2^2 / 2) .* ...
                    normcdf((log(qgrid) - nu2 - alpha2 * tau2^2) / tau2)) + ...
              (alpha2 / (alpha2 + beta2)) .* ...
                (qgrid.^beta2  .* ...
                    exp(-beta2 * nu2 + beta2^2 * tau2^2 / 2) .* ...
                    (1 - normcdf((log(qgrid) - nu2 + beta2 * tau2^2) / tau2))) ;
      cdf3 = normcdf((log(qgrid) - nu3) / tau3) - ...
              (beta3 / (alpha3 + beta3)) .* ...
                (qgrid.^(-alpha3)  .* ...
                    exp(alpha3 * nu3 + alpha3^2 * tau3^2 / 2) .* ...
                    normcdf((log(qgrid) - nu3 - alpha3 * tau3^2) / tau3)) + ...
              (alpha3 / (alpha3 + beta3)) .* ...
                (qgrid.^beta3  .* ...
                    exp(-beta3 * nu3 + beta3^2 * tau3^2 / 2) .* ...
                    (1 - normcdf((log(qgrid) - nu3 + beta3 * tau3^2) / tau3))) ;
      cdf = w1 * cdf1 + w2 * cdf2 + w3 * cdf3 ;
      
    end ;
    lccdf = log10(1 - cdf) ;
    hold on ;
      plot(log10qgrid,lccdf,'r-','LineWidth',2) ;
    hold off ;




      orient landscape ;
    print('-dpsc2',['UNC2001RS1CCDF' outdistr num2str(itimeb) '.ps']) ;



  end ;    %  of itimeb loop







elseif  ipart == 71  | ...
        ipart == 72  | ...
        ipart == 73  | ...
        ipart == 74  ;      %  Then do additions to part 42, which was:
                            %  then plot Log-Log Complementary CDF 
                            %  (full data version)


  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


  %  Load Quantiles
  %
  fid = fopen('DataUNC2001\responsesizequantilesRaw.txt','rt') ;
          %  'rt' is for "read only" and "text"
  fstr = '%f' ;
  [mq,cnt] = fscanf(fid,fstr) ;
  fclose('all') ;
  mq = reshape(mq,21,1000)' ;
          %  use this for 1000x21 matrix, since reshape works column wise


  %  loop through time blocks
  %
  for itimeb = vitimeb ;    %  loop through time blocks


    %  first get part that appeared earlier
    %  (i.e. for given quantiles)
    %
    vq = mq(:,virowraw(itimeb)) ;
        %  unpack vector of quantiles

    n = vnumresp(itimeb) ;
        %  unpack full sample size

    nq = length(vq) ;
    vccdf = flipud((1/n) + (0:(nq-1))' / nq) ;


    %  load large data values, and integrate
    %
    filename = filenamecell{itimeb} ;
    mtailresp = load(['DataUNC2001\' filename]) ;
    vtailsize = mtailresp(:,1) ;
    ntail = length(vtailsize) ;

    vq = [vq; sort(vtailsize)] ;
    vccdf = [vccdf; flipud((1:ntail)' / n)] ;
    
    [vq, vindsort] = sort(vq) ;
    vccdf = vccdf(vindsort) ;
    nq = length(vq) ;
    
    flagsame = vq(1:(nq-1)) == vq(2:nq) ;
        %  one where next entry is the same
    vq = vq(~flagsame) ;
    vccdf = vccdf(~flagsame) ;
        %  keep only entries that are different (always keeping last one)
    nq = length(vq) ;


    %  compute direct difference quotients
    %
    sbottom = -4.01 ;
    stop = 0.01 ;

    vcent = mean([log10(vq(2:nq)), log10(vq(1:(nq-1)))],2) ;
    vdenom = (log10(vq(2:nq)) - log10(vq(1:(nq-1)))) ;
    flag0 = vdenom <= 0 ;
    nflag0 = sum(flag0) ;
    if nflag0 > 0 ;
      vdenom(flag0) = ones(nflag0,1) ;
      vslope = (log10(vccdf(2:nq)) - log10(vccdf(1:(nq-1)))) ./ vdenom ;
      vslope(flag0) = sbottom * ones(nflag0,1) ;
    else ;
      vslope = (log10(vccdf(2:nq)) - log10(vccdf(1:(nq-1)))) ./ vdenom ;
    end ;

    left = 1 ;
    right = 9 ;
    bottom = -7 ;
    top = 0.01 ;


    %  Make gridded difference quotients
    %
    del1 = 0.3 ;
    del2 = 0.03 ;
        %  grid widths
    vgrid1 = (left:del1:right)' ;
    vgrid2 = (left:del2:right)' ;
        %  grid points covering full range
    minq = log10(min(vq)) ;
    maxq = log10(max(vq)) ;
        %  endpoints where have data beyond
    vgrid1 = vgrid1(minq < vgrid1) ;
    vgrid2 = vgrid2(minq < vgrid2) ;
    vgrid1 = vgrid1(vgrid1 < maxq) ;
    vgrid2 = vgrid2(vgrid2 < maxq) ;
        %  reduce to grid within the data
    vheightgrid1 = interp1(log10(vq),log10(vccdf),vgrid1) ;
    vheightgrid2 = interp1(log10(vq),log10(vccdf),vgrid2) ;
        %  interpolate to grid
    ngrid1 = length(vheightgrid1) ;
    ngrid2 = length(vheightgrid2) ;
    vcent1 = mean([vgrid1(2:ngrid1), vgrid1(1:(ngrid1-1))],2) ;
    vcent2 = mean([vgrid2(2:ngrid2), vgrid2(1:(ngrid2-1))],2) ;
        %  center points for plotting
    vslope1 = (vheightgrid1(2:ngrid1) - vheightgrid1(1:(ngrid1-1))) ./ del1 ;
    vslope2 = (vheightgrid2(2:ngrid2) - vheightgrid2(1:(ngrid2-1))) ./ del2 ;
        %  get difference quotients (slopes)



    %  Make graphics
    %
    subplot(2,2,1) ;
      plot(log10(vq),log10(vccdf),'b-','LineWidth',2) ;
        axis([left, right, bottom, top]) ;
        title(['log-log CCDF, HTTP R. S., n = ' num2str(n)],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('log_{10}(1 - prob.)','FontSize',18) ;
          tx = left + 0.1 * (right - left) ;
          ty = bottom + 0.1 * (top - bottom) ;
        text(tx,ty,[daystrcell{itimeb} ...
                    ', ' timestrcell{itimeb}],'FontSize',18) ;

        vprobplot = [0.5, 0.9, 0.999, 0.99999] ;
        vlq = [] ;
        for probplot = vprobplot ;
          ty = log10(1 - probplot) ;
          [temp,ind2min] = min(abs(vccdf - (1 - probplot))) ;
          tx = log10(vq(ind2min)) ;
          vlq = [vlq; tx] ;
          hold on ;
            plot([tx; tx],[bottom,top],'g--') ;
            tx = tx + 0.02 * (right - left) ;
            text(tx,ty,[num2str(probplot) ' quantile'], ...
                                  'FontSize',15,'Color','g') ;
          hold off ;
        end ;
        nlq = length(vlq) ;

    subplot(2,2,3) ;
      plot(vcent,vslope,'r-','LineWidth',2) ;
        axis([left,right, sbottom, stop]) ;
        title(['Corresponding Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([left; right],[-1; -1],'k--') ;
          plot([left; right],[-2; -2],'k--') ;
          plot([vlq'; vlq'],[sbottom * ones(1,nlq); stop * ones(1,nlq)],'g--') ;
        hold off ;

    subplot(2,2,2) ;
      plot(vcent1,vslope1,'m-','LineWidth',2) ;
        axis([left,right, sbottom, stop]) ;
        title(['\Delta = ' num2str(del1) ' Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([left; right],[-1; -1],'k--') ;
          plot([left; right],[-2; -2],'k--') ;
          plot([vlq'; vlq'],[sbottom * ones(1,nlq); stop * ones(1,nlq)],'g--') ;
        hold off ;


    %  add on simulated versions
    %
    if ipart == 71 ;    %  add simulated Paretos
      outstr = 'par' ;
      labelstr = 'Fit Pareto' ;
      alpha = 1.3923 ;
      sigma = 2226.2475 ;
    elseif ipart == 72 ;    %  add simulated double Pareto log normal
      outstr = 'dpln' ;
      labelstr = 'Fit DP-lN' ;
      alpha = 1.2 ;    %  "tail parameter" near origin
      beta = 1.2 ;    %  tail parameter
      nu = 7 ;    %  mean of log-normal
      tau = 1.0 ;    %  sd of log-normal
          %  these chosen by eye, for good fit in ipart = 10, Thurday morning
    elseif ipart == 73 ;    %  add simulated dPlN, 3 mixture
      outstr = '3dpln' ;
      labelstr = 'Fit DP-lN 3 mix' ;
      alpha1 = 2 ;    %  "tail parameter" near origin
      beta1 = 2 ;    %  tail parameter
      nu1 = 5.7 ;    %  mean of log-normal
      tau1 = 0.6 ;    %  sd of log-normal
      alpha2 = 10 ;    %  "tail parameter" near origin
      beta2 = 10 ;    %  tail parameter
      nu2 = 8.45 ;    %  mean of log-normal
      tau2 = 1.2 ;    %  sd of log-normal
      alpha3 = 1.3 ;    %  "tail parameter" near origin
      beta3 = 1.3 ;    %  tail parameter
      nu3 = 13.35 ;    %  mean of log-normal
      tau3 = 0.75 ;    %  sd of log-normal
      w1 = 0.6 ;
      w2 = 0.399 ;
      w3 = 1 - w1 - w2 ;
          %  these chosen by eye, for good fit in ipart = 10, Thurday morning
    elseif ipart == 74 ;    %  add simulated log gamma's
      outstr = 'lg' ;
      labelstr = 'Fit log Gamma' ;
      vqgrid = mq(:,virowraw(itimeb)) ;
          %  again unpack vector of quantiles (on grid)
      phat = gamfit(log(vqgrid + 2)) ;
    end ;    %    of ipart if-block

    nsim = 100 ;
    simseed = 79487520 ;
    rand('seed',simseed) ;
        %  set seed
    for isim = 1:nsim ;
      disp(['    working on isim = ' num2str(isim) ' of ' num2str(nsim)]) ;

      if ipart == 71 ;    %  add simulated Paretos
        vsimdat = sigma * ((1 - rand(n,1)).^(-1/alpha) - 1) ;
      elseif ipart == 72 ;    %  add simulated double Pareto log normal
        vunif = rand(n,1) ;
            %  Unif(0,1)
        flag = vunif < (alpha ./ (alpha + beta)) ;
            %  one with probability proportional to alpha
        vunif = rand(n,1) ;
            %  Unif(0,1)
        vu = flag .* vunif .^ (1 ./ beta) + ...
             (~flag) .* vunif .^ (-1 ./ alpha) ;
        vv = nu + tau .* randn(n,1) ;
        vv = exp(vv) ;
        vsimdat = vu .* vv ;
      elseif ipart == 73 ;    %  add simluated dPlN, 3 mixture
        vunif = rand(n,1) ;
            %  Unif(0,1)
        flag1 = vunif < w1 ;
            %  one with probability w1
        flag3 = (w1 + w2) < vunif ;
            %  one with probability w3 = 1 - (w1 + w2)
        flag2 = ~(flag1 | flag3) ;
            %  one when neither flag1 nor flag3 are
        valpha = flag1 * alpha1 + flag2 * alpha2 + flag3 * alpha3 ;
        vbeta = flag1 * beta1 + flag2 * beta2 + flag3 * beta3 ;
        vunif = rand(n,1) ;
            %  Unif(0,1)
        flag = vunif < (valpha ./ (valpha + vbeta)) ;
            %  one with probability proportional to alpha
        vunif = rand(n,1) ;
            %  Unif(0,1)
        vu = flag .* vunif .^ (1 ./ vbeta) + ...
             (~flag) .* vunif .^ (-1 ./ valpha) ;
        vnu = flag1 * nu1 + flag2 * nu2 + flag3 * nu3 ;
        vtau = flag1 * tau1 + flag2 * tau2 + flag3 * tau3 ;
        vv = vnu + vtau .* randn(n,1) ;
        vv = exp(vv) ;
        vsimdat = vu .* vv ;
      elseif ipart == 74 ;    %  add simulated log gamma's
        vsimdat = exp(gamrnd(phat(1),phat(2),n,1)) - 2 ;
      end ;    %    of ipart if-block

      vqsim = sort(vsimdat) ;
      vccdfsim = flipud((1:n)' / n) ;

      flagsame = vqsim(1:(n-1)) == vqsim(2:n) ;
          %  one where next entry is the same
      vqsim = vqsim(~flagsame) ;
      vccdfsim = vccdfsim(~flagsame) ;
          %  keep only entries that are different (always keeping last one)
      nqsim = length(vq) ;

      vgridsim = (left:del1:right)' ;
          %  grid points covering full range
      minq = log10(min(vqsim)) ;
      maxq = log10(max(vqsim)) ;
          %  endpoints where have data beyond
      vgridsim = vgridsim(minq < vgridsim) ;
      vgridsim = vgridsim(vgridsim < maxq) ;
          %  reduce to grid within the data
      vheightgridsim = interp1(log10(vqsim),log10(vccdfsim),vgridsim) ;
          %  interpolate to grid
      ngridsim = length(vheightgridsim) ;
      vcentsim = mean([vgridsim(2:ngridsim), vgridsim(1:(ngridsim-1))],2) ;
          %  center points for plotting
      vslopesim = (vheightgridsim(2:ngridsim) - vheightgridsim(1:(ngridsim-1))) ./ del1 ;
          %  get difference quotients (slopes)

      hold on ;
        plot(vcentsim,vslopesim,'b-') ;
      hold off ;


    end ;    %  of isim loop
    hold on ;
      plot(vcent1,vslope1,'m-','LineWidth',2) ;
        %  replot to bring back to top

        tx = left + 0.1 * (right - left) ;
        ty = sbottom + 0.1 * (stop - sbottom) ;
      text(tx,ty,[num2str(nsim) ' sim''d ' ...
                      labelstr ' copies'],'FontSize',18) ;
  hold off ;


    subplot(2,2,4) ;
      plot(vcent2,vslope2,'m-','LineWidth',2) ;
        axis([left,right, sbottom, stop]) ;
        title(['\Delta = ' num2str(del2) ' Slopes of log-log CCDF'],'FontSize',18) ;
        xlabel('log_{10}(quantile)','FontSize',18) ;
        ylabel('-\alpha','FontSize',18) ;

        hold on ;
          plot([left; right],[-1; -1],'k--') ;
          plot([left; right],[-2; -2],'k--') ;
          plot([vlq'; vlq'],[sbottom * ones(1,nlq); stop * ones(1,nlq)],'g--') ;
        hold off ;

      orient landscape ;
    print('-dpsc2',['UNC2001RS1allCCDFS' outstr num2str(itimeb) '.ps']) ;


  end ;    %  of itimeb loop




elseif  ipart == 81  | ...
        ipart == 82  | ...
        ipart == 83  | ...
        ipart == 84  ;    %  extremal index analysis,
                               %  based on the estimates (8.11) 
                               %  from Embrecths, Kueppelberg, Mikosch (1997)
                               %  page 429, section 8.1

  figure(1) ;
  clf ;


  %  Load summary data (gives filenames, etc.)
  %
  load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;
  
  
  %  loop through time blocks, read big data files, and make plots
  %
  for itimeb = vitimeb ;    %  loop through time blocks

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;

    mtailresp = load(['DataUNC2001\' filename]) ;

    mtailresp = mtailresp(:,[1 2 4]) ;
        %  keep only Response Sizes, start times and durations


    %  set up preliminaries
    %
    vstart = mtailresp(:,2) ;
        %  vector of start times
    vstart = vstart - min(vstart) ;
        %  set so that smallest is 0

    vsizes = mtailresp(:,1) ;
    vtimes = mtailresp(:,3) ;

    flagn0 = (vsizes > 0) & (vtimes > 0) ;
        %  ones where want to keep data
    vstart = vstart(flagn0) ;
    vsizes = vsizes(flagn0) ;
    vtimes = vtimes(flagn0) ;
    n = sum(flagn0) ;

    [vstart, vind] = sort(vstart) ;      
        %  sort according to start times
    vstart = linspace(1 / (2 * n), 1 - (1 / (2 * n)), n) ;
        %  normalize to be equally spaced on the interval [0,1] 


    if ipart == 81 ;    %  analyze response sizes, > 100k

      vdata = vsizes(vind) ;
          %  sizes in first column
          %  sorted by start times
      datastr = 'HTTP Resp. Sizes, > 100k' ;
      datoutstr = 'siz100k' ;

    elseif ipart == 82 ;    %  analyze response times, > 100k

      vdata = vtimes(vind) ;
          %  durations in third column
          %  sorted by start times
      datastr = 'HTTP Resp. Times, > 100k' ;
      datoutstr = 'tim100k' ;

    elseif ipart == 83 ;    %  analyze response rates, > 100k

      vdata = vsizes(vind) ./ vtimes(vind) ;
          %  durations in third column
          %  sorted by start times
      datastr = 'HTTP Resp. Rates, > 100k' ;
      datoutstr = 'rat100k' ;

    elseif ipart == 84 ;    %  analyze response inverse rates, > 100k

      vdata = vtimes(vind) ./ vsizes(vind) ;
          %  durations in third column
          %  sorted by start times
      datastr = 'HTTP Resp. Inverse Rates, > 100k' ;
      datoutstr = 'irat100k' ;

    end ;
    [vsortdata,vindsort] = sort(vdata) ;
        %  sorted in increasing order of data values
    [tempinvsort,vrank] = sort(vindsort) ;
        %  get ranks of data, by inverting sort indices


    vthreshprop = [0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001, ...
                     0.0005, 0.0002, 0.0001, 0.00005, 0.00002, 0.00001, ...
                     0.000005, 0.000002, 0.000001] ;
        %  vector of thresholds, as tail proportions
    flagthresh = round(n * vthreshprop) >= 5 ;
        %  one where can have at least 5 data points larger
    vthreshprop = vthreshprop(flagthresh) ;
        %  keep only these thresholds

    vnblocks = logspace(2,log10(n/100),30) ;
        %  vector of numbers of blocks


    %  Main data loop
    %
    vN = [] ;
    mtheta2 = [] ;
    for threshprop = vthreshprop ;

      N = round(n * threshprop) ;
          %  total number of points above threshold
      vN = [vN; N] ;
          %  
      flag =  (vrank > n - N)  ;
          %  ones at data values above threshold
      un = vsortdata(n - N) ;
          %  threshold value (value of last point before threshold)
      vstartthresh = vstart(flag) ;
          %  keep only start times where data is above threshold
%plot(vstart,vdata,'bo',vstartthresh,vdata(flag),'r+')
%pauseSM


      vtheta2 = [] ;
      for nblocks = vnblocks ;

        inblocks = round(nblocks) ;
            %  integer version of number of blocks

        disp(['    Working on nblocks = ' num2str(inblocks) ...
                      ' & threshprop = ' num2str(threshprop)]) ;

        binedges = linspace(0,1,inblocks+1) ;
            %  edges of blocks
        vcounts = histc(vstartthresh,binedges) ;
            %  count number of thresholded starting times within blocks
        vcounts = vcounts(1:inblocks) ;
            %  cut off last block (only counts # = last edge)
        K = sum(vcounts > 0) ;
            %  number of vcounts > 0
        theta2 = K / N ;
            %  from formula (8.11)

        vtheta2 = [vtheta2; theta2] ;


        disp(['        Results:  N = ' num2str(N) ...
                               ', K = ' num2str(K) ...
                               ', theta2 = ' num2str(theta2)]) ;

      end ;    %  of loop through nblocks


      mtheta2 = [mtheta2, vtheta2] ;


    end ;    %  of loop through threshprop


    %  Make main graphics result
    %
    subplot(2,2,1) ;    %  functions of log10(threshprop)
      plot(log10(vthreshprop'),mtheta2','-') ;
        title(['Extremal Index, ' datastr]) ;
        xlabel('log_{10}(Tail Threshold Proportion)') ;
        ylabel('{\theta}2') ;

    subplot(2,2,3) ;    %  functions of log10(number of blocks)
      plot(log10(vnblocks'),mtheta2,'-') ;
        title(['Extremal Index, ' daystrcell{itimeb} ', ' timestrcell{itimeb}]) ;
        xlabel('log_{10}(Number of Blocks)') ;
        ylabel('{\theta}2') ;

    subplot(2,2,2) ;
      surf(mtheta2)
        title(['Extremal Index, time block ' num2str(itimeb)]) ;

    subplot(2,2,4) ;
      plot(log10(vthreshprop'),log10(vN),'r-',...
                   'LineWidth',2) ;
        title(['# above threshold, N, n = ' num2str(n)]) ;
        xlabel('log_{10}(Tail Threshold Proportion)') ;
        ylabel('log_{10}(N)') ;

    pstr = ['UNC2001RS1ExtInd1T' num2str(itimeb) datoutstr] ;
    orient landscape ;
    print('-dpsc2',[pstr '.ps']) ;



  end ;



end ;    %  of ipart if-block






