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

%    Makes Figure I for the Variable Heavy Tails paper

%    Copied from NZIX2000RespSiz1.m
%    with parameters set to:
%        ipart = 5  -  log-log CCDF, with vis. fit mix of 3 dP-LN
               %            (as in VarHeavyTail1figG.m)
%        vitimeb = 1 ;    %  vector of indices of time blocks

%
%    For making analysis of Response Data,
%    based on traces from New Zealand, Internet Exchange
%    gathered in July 2000,
%    and filtered by Felix Hernandez Campos.
%
%    This works with the full data sets, 
%    in the directory:
%        \Research\NetworkData\ResponseSize\DataNZIX2000
%
%    For further details, see the file:    ReadMe.txt


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


vitimeb = 1 ;    %  vector of indices of time blocks
                 %  1 - July 6, 2000, Thursday, 8:00AM - 12:00noon
                 %  2 - July 6, 2000, Thursday, 1:00PM - 5:00PM


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 ;



%  Set up basics
%
%indatstr = '\Research\NetworkData\ResponseSize\DataNZIX2000\20000706_800_1200_rsps.txt' ;
indatstr = '..\ResponseSize\DataNZ\20000706_800_1200_rsps.txt' ;
timestr = 'Thurs. Morn.' ;
%    vpardpln = [1.6, 1.6 6.8 1.7] ;
    %  needs to be less steep
%    vpardpln = [1.3, 1.6 6.8 1.2] ;
    %  reduce the s. d.
%    vpardpln = [1.3, 1.3 6.8 1.2] ;
    %  move th mean right
vpardpln = [1.3, 1.3 6.9 1.2] ;
    %  vector of parameters for double pareto log normal
    %  chosen by eye
    %    alpha, beta, nu, tau
%    vpar3mdpln = [2 2 6.0 0.6 10 10 8.84 1.05 1.6 1.6 12.4 1 0.62 0.3775] ;
    %  lower the 2nd prob
%    vpar3mdpln = [2 2 6.0 0.6 10 10 8.84 1.05 1.6 1.6 12.4 1 0.67 0.3275] ;
    %  lower the 3rd prob
%    vpar3mdpln = [2 2 6.0 0.6 10 10 8.84 1.05 1.6 1.6 12.4 1 0.67 0.32935] ;
    %  move 3rd mean right
%    vpar3mdpln = [2 2 6.0 0.6 10 10 8.84 1.05 1.6 1.6 13.6 1 0.67 0.32935] ;
    %  move 2nd mean left
%    vpar3mdpln = [2 2 6.0 0.6 10 10 8.7 1.05 1.6 1.6 13.6 1 0.67 0.32935] ;
    %  move first mean left
vpar3mdpln = [2 2 5.7 0.6 10 10 8.7 1.05 1.6 1.6 13.6 1 0.67 0.32935] ;
        %  vector of parameters for 3-mixture of double pareto log normal
    %  chosen by eye
    %    alpha1, beta1, nu1, tau1, alpha2, beta2, nu2, tau2, 
    %    alpha3, beta3, nu3, tau3, w1, w2, 
    %    (w3 is automatically chosen as 1 - w1 - w2)



%  Load big data set
%
mdata = load(indatstr) ;

%size(mdata)
%pauseSM
    %  used these line to check load


vsize = mdata(:,1) ;
    %  keep response sizes
n = length(vsize) ;
mdata = [] ;
    %  to save space




%  Remove 0 duration responses
%
flag0 = (vsize == 0) ;
n0 = sum(flag0) ;
if n0 > 0 ;
  vsize = vsize(~flag0) ;
  disp(['  Removed ' num2str(n0) ' Zero Duration responses, out of ' num2str(n)]) ;
  n = n - n0 ;
  disp(['  Have ' num2str(n) ' left']) ;
end ;




%  Make main graphic
%
vq = sort(vsize) ;
vccdf = flipud((1:n)' / n) ;
nq = length(vq) ;


subplot(1,1,1) ;
  left = 0.01 ;
  right = 8 ;
  bottom = -6 ;
  top = 0.01 ;


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


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

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

%  titstr = ['log-log CCDF, HTTP R. S., n = ' ...
%                     num2str(n) ' & DP-lN 3 mix'] ;
%  tith = title(titstr,'FontSize',18) ;

%  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] ;
  vprobplot = [0.5, 0.9, 0.99, 0.999, 0.9999] ;
    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
  alpha1 = vpar3mdpln(1) ;    %  "tail parameter" near origin
  beta1 = vpar3mdpln(2) ;    %  tail parameter
  nu1 = vpar3mdpln(3) ;    %  mean of log-normal
  tau1 = vpar3mdpln(4) ;    %  sd of log-normal
  alpha2 = vpar3mdpln(5) ;    %  "tail parameter" near origin
  beta2 = vpar3mdpln(6) ;    %  tail parameter
  nu2 = vpar3mdpln(7) ;    %  mean of log-normal
  tau2 = vpar3mdpln(8) ;    %  sd of log-normal
  alpha3 = vpar3mdpln(9) ;    %  "tail parameter" near origin
  beta3 = vpar3mdpln(10) ;    %  tail parameter
  nu3 = vpar3mdpln(11) ;    %  mean of log-normal
  tau3 = vpar3mdpln(12) ;    %  sd of log-normal
  w1 = vpar3mdpln(13) ;
  w2 = vpar3mdpln(14) ;
  w3 = 1 - w1 - w2 ;




  %  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',['NZIX00RS1allCCDFfullT' num2str(itimeb) ...
%                                'p' num2str(ipart) '.ps']) ;


if icolor == 0 ;
  pstr = 'VarHeavyTail2figIbw' ;
else ;
  pstr = 'VarHeavyTail2figI' ;
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']) ;





