disp('Running MATLAB script file VarHeavyTail2figG.m') ;

%    Makes Figure G for the Variable Heavy Tails paper

%    Copied from UNC2001RespSiz1.m
%    with parameters set to:
%        ipart = 63  -  Full Data Log-Log Complementary CDF
%                           plus Double-Pareto log-normal 3-mixture (fit by eye)
%                           plus lines showing variability
%        vitimeb = 10 ;    %  vector of indices of time blocks

%
%    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


icolor = 0 ;    %  0 - black and white only
                %  1 - original full color


itimeb = 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





figure(1) ;
clf ;
if icolor == 0 ;
  colorcell = {'k' 'k' [0.2 0.2 0.2]} ;
          %  1st - data
          %  2nd - model
          %  3rd - simulated versions
  modtypestr = '--' ;
  simltypestr = ':' ;
else ;
  colorcell = {'b' 'r' 'm'} ;
          %  1st - data
          %  2nd - model
          %  3rd - simulated versions
  modtypestr = '-' ;
  simltypestr = '-' ;
end ;


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


%  Load Quantiles
%
%fid = fopen('\Research\NetworkData\ResponseSize\DataUNC2001\responsesizequantilesRaw.txt','rt') ;
fid = fopen('..\ResponseSize\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


subplot(1,1,1) ;
  left = 1 ;
  right = 9 ;
  bottom = -7 ;
  top = 0.01 ;


disp(['       plotting for time block ' num2str(itimeb) ...
                       ', ' daystrcell{itimeb} ...
                      ', ' timestrcell{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(['\Research\NetworkData\ResponseSize\DataUNC2001\' filename]) ;
mtailresp = load(['..\ResponseSize\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) ;

vxdat = log10(vq) ;
vydat = log10(vccdf) ;

plot(vxdat,vydat,[colorcell{1} '-'],'LineWidth',2) ;
  axis([left, right, bottom, top]) ;

%  xlabel('log_{10}(quantile)','FontSize',12) ;
  xlabel('log_{10}(x)','FontSize',12) ;
%  ylabel('log_{10}(1 - prob.)','FontSize',12) ;
  ylabel('log_{10}(P[X > x])','FontSize',12) ;


hold on ;

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

%  vprobplot = [0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999, 0.99999] ;
  vprobplot = [0.5, 0.9, 0.99, 0.999, 0.9999, 0.99999] ;
    npp = length(vprobplot) ;
  vpph = [] ;
  vppth = [] ;
  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',12) ;
  end ;


  %  set parameters
  %
  diststr = 'D-Par log-N 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




  %  add simulated data
  %

  %  set up preliminaries
  %
  if icolor == 0 ;
    nsim = 40 ;
  else ;
    nsim = 100 ;
  end ;

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


    vsimdat = sort(vsimdat) ;
        %  sorted version

    vsimdat = vsimdat(vindkeep) ;


    %  Plot this curve
    %
    plot(log10(vsimdat),log10(pgrid),simltypestr,'Color',colorcell{3}) ;

  end ;    %  of isim loop



  %  replot most important parts
  %
  plot(vxdat,vydat,[colorcell{1} '-'],'LineWidth',2) ;



  %  overlay theoretical distributions
  %
  ngrid = 400 ;
  log10qgrid = linspace(left,right,ngrid) ;
  qgrid = 10.^log10qgrid ;
  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 ;
  
  lccdf = log10(1 - cdf) ;
  plot(log10qgrid,lccdf,[colorcell{2} modtypestr],'LineWidth',2) ;



hold off ;





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





if icolor == 0 ;
  pstr = 'VarHeavyTail2figGbw' ;
else ;
  pstr = 'VarHeavyTail2figG' ;
end ;
  orient portrait ;
  set(gcf,'PaperSize',[6.0, 4.5]) ; 
  set(gcf,'PaperPosition',[0.25, 0.25, 5.5, 4.0]) ; 
print('-dpsc', [pstr '.ps']) ;




