scalemginfty.m

用matalab开发的随机数生成程序包

function [r_times, irc_sys] = scalemginfty(sc_type, lambda, ...
            alpha, servmu, maxtime, b_rescale_first, b_plot, b_verb)

% SCALEMGINFTY Simulate the rescaled and centered cumulative
% workload process in a stationary M/G/infinity queue:
%
%      Y(lambda,t)=1/b (W(lambda,at)-EW(lambda,at)),
%
%      W(lambda,t)=int_0^t H(lambda, s) ds,
%
% where H(lambda,t) is the system size process, W(lambda,t) is
% cumulative workload process, lambda arrival intensity,
% a=a(lambda) the time scale, b=b(lambda) the space scale. 
%
% Service times follow normalized Pareto distribution with
% 
% pdf     f(x) = alpha*gamma/(1+gamma*x)^(1+alpha), x>0,
% cdf     F(x) = 1-(1+gamma*x)^(-alpha),
%
% for alpha>1, which implies finite expectation. If 1<alpha<2,
% which is the interesting case, the variance is infinite.
%
% The model is also known as infinite source Poisson model. It
% describes traffic sources which arrive according to a 
% a homogeneous stationary Poisson process and during iid
% random times produce traffic at a constant rate.
% 
% The process Y(lambda,t) is a random piecewise linear function. It
% is represented by two vectors of the break points and the values
% at these points respectively. 
% 
% It is known that if the service-time distribution has heavy tails 
% behaving as x^(-alpha) for 1<alpha<2, then depending on the rate
% of growth of the time sequence a=a(lambda), as lambda grows to
% infinity, the process Y(lambda,t) converges to one of the three
% the possible limits. The limit processes are fractional Brownian
% motion, stable Levy motion and the third process, which can be
% expressed as an integral with respect to a Poisson random
% measure. The corresponding time and space scales are
% 
% 1. Fast (fBm): a(lambda) = lambda^(1/(alpha-1)-epsilon), epsilon>0
%                b(lambda) = (lambda*a(lambda)^(3-alpha))^0.5
%
% 2. Intermediate: a(lambda) = lambda^(1/(alpha-1))
%                  b(lambda) = a(lambda)
%
% 3. Slow (LM): a(lambda) = lambda^(1/(alpha-1)+epsilon), epsilon>0
%               b(lambda) = (lambda*a(lambda))^(1/alpha)
%
% For details see
%
% T.Mikosch, S.Resnick, H.Rootz閚 and A.Stegeman (2002) Is network
% traffic approximated by stable L関y motion or fractional Brownian
% motion? Ann. Appl. Probab. 12(1), 23--68. 
% 
% I.Kaj, M.S.Taqqu (2004) Convergence to fractional Brownian motion 
% and to the Telecom process: the integral representation
% approach. Department of Mathematics, Uppsala University,
% Technical report U.U.D.M. 2004:16.
%
% Usage:
%   [r_times, irc_sys] =  scalemginfty(sc_type, lambda, alpha,
%               servmu, maxtime, b_rescale_first, b_plot, b_verb) 
%
% Inputs:
%        sc_type - type of scaling: 1-fBm, 2-Intermediate, 3-Levy
%          motion, can be a vector containing combinations of the
%          above 
%        lambda - arrival intensity
%        alpha - tail parameter of the service-time distribution:
%          normalized Pareto, 1<alpha<=2 (see above)
%        servmu - expected value of the service-time distribution. 
%          Should be positive.
%        maxtime - maximal time for the process before rescaling
%        b_rescale_first - if 1 then first rescale, then
%          integrate. If 0 then the other way around.
%        b_plot - if 1 then plot the processes
%        b_verb - if 1 then print processing info
%
% Outputs:
%        r_times - a column vector of the break times of the
%          rescaled process(es) 
%        irc_sys - a column vector of values of the rescaled
%          process(es) at the break points
%
% See also SIMSTMGINFTY, SCALEREN, SCALEONOFF, LRSCALES.

% Authors: R.Gaigalas, I.Kaj
% v2.1 17-Dec-05


  if (nargin<1) % default parameter values
    sc_type = [1 2 3];
    lambda = 50;
    alpha = 1.2;    
    servmu = 1;
    maxtime = 5000;
    b_rescale_first = 1;
    b_plot = 1;
    b_verb = 1;
  end
  

  % generate the system size process in a stationary M/G/infinity
  % queue 
  [jmptimes, syssize] = simstmginfty(maxtime, lambda, ...
              @simparetonrm, {alpha, 1/(servmu*(alpha-1))}, b_verb);  

  if b_plot
    figure(1);
    clf;
    stairs(jmptimes, syssize);
  end

  % subtract the expected value
  c_sys = syssize-lambda*servmu;

  % expected number of events   
  nevmu = lambda*(1+1/servmu)*maxtime;   
  
  % set the time and space scales according to the scaling regime 
  minev = NaN; % can be set explicitely, otherwise set in lracc
  acc = lracc(sc_type, alpha, lambda, length(c_sys), minev); % acceleration factor  
  [tm_scale, sp_scale] = lrscales(sc_type, alpha, lambda, ...
                                  acc, maxtime, b_verb);
  
  if b_rescale_first % first rescale, then integrate

    [r_times, irc_sys, rc_sys] = strescint(jmptimes, c_sys, ...
                            tm_scale, sp_scale, b_verb);
    
    if b_plot  % plot the centered and rescaled system size
               % in the smallest common time window 
      figure(2); 
      clf;
      stairs(r_times, rc_sys);
    end  
  else % first integrate, then rescale
    
    [r_times, irc_sys, ic_sys] = stintresc(jmptimes, c_sys, ...
                             tm_scale, sp_scale, b_verb);

    % plot the centered and integrated system size process
    if b_plot  
      figure(2);
      clf;
      hold on;
      stairs(jmptimes, c_sys);
      plot(jmptimes, ic_sys, 'r');      
      hold off;      
    end
  end   

  % plot parts of the rescaled versions in a largest common 
  % time  window
  if b_plot

    figure(3); 
    clf;
    plot(r_times, irc_sys);
  end  

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%