scaleonoff.m

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

function [rs_times, ircs_val] = scaleonoff(sc_type, maxtime, nproc, ...
   on_alpha, on_mu, off_alpha, off_mu, b_rescale_first, b_plot, b_verb)  

% SCALEONOFF Simulate a rescaled and centered sum of integrated
% on-off processes with heavy-tailed distribution of on-periods:
%
%      Y_m(t)=1/b (W(m,at)-EW(m,at)),
%
%      W(m,t)=int_0^t sum_{i=1}^m X_i(s) ds
%
% where X_i(t) are independent copies of an on-off process X(t),
% a=a(m) the time scale, b=b(m) the space scale. 
%
% On- and off-periods 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 1<alpha_on<2, and alpha_on<alpha_off. Hence, on-periods have
% finite expectation but infinte variance.
% 
% The process Y_m(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 on-distribution has heavy tails behaving as
% x^(-alpha) for 1<alpha<2, then the limit process of Y_m(t), when m
% grows to infinity, depends on the rate of growth of the time sequence
% a=a(m). Three scaling regimes exist and two known limit processes
% are fractional Brownian motion and stable Levy motion. The
% corresponding time and space scales are
% 
% 1. Fast (fBm): a(m) = m^(1/(alpha-1)-epsilon), epsilon>0
%                b(m) = (m*a(m)^(3-alpha))^0.5
%
% 2. Intermediate: a(m) = m^(1/(alpha-1))
%                  b(m) = a(m)
%
% 3. Slow (LM): a(m) = m^(1/(alpha-1)+epsilon), epsilon>0
%               b(m) = (m*a(m))^(1/alpha)
%
% Up to this date (Nov 2005), there are no publicated results on
% convergence in the intermediate scaling regime. You are kindly  
% requested to notify if you know any.
%
% 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. 
%
% Usage:
% [rs_times, ircs_val] = scaleonoff(sc_type, maxtime, nproc,
%   on_alpha, on_mu, off_alpha, off_mu, 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 
%        maxtime - maximal time for the process before rescaling
%        nproc - number of processes to superpose
%        on_alpha - tail parameter of the on-distribution:
%          normalized Pareto, 1<alpha<=2
%        on_mu - expected value of the on-distribution. 
%          Should be positive.
%        off_alpha - tail parameter of the off-distribution:
%          normalized Pareto, alpha_on<alpha_off
%        off_mu - expected value of the off-distribution. 
%          Should be positive.
%        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:
%        rs_times - a column vector of the break times of the
%          rescaled process(es) 
%        ircs_val - 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];

    maxtime = 100;
    nproc = 50; 
    
    on_alpha = 1.4;
    on_mu = 1;
    off_alpha = 2;
    off_mu = 2;

    b_rescale_first = 0;
    b_plot = 1;
    b_verb = 1;
  end
  
  
  % generate on-off processes as column vectors 
  [jmptimes, indproc, neff] = onoff(nproc, maxtime,  ...
              @simparetonrm, {on_alpha, 1/(on_mu*(on_alpha-1))}, ...
              @simparetonrm, {off_alpha, 1/(off_mu*(off_alpha-1))}, ...
	      b_verb);

  % add up the piecewise constant processes
  [s_times, s_val] = stairsum(jmptimes, indproc); 
  % get rid of accumulated zero times in the beginning
  zeri = find(s_times==0);
  s_times = s_times(zeri(end):end);
  s_val = s_val(zeri(end):end);
    
  if b_plot
    figure(1);
    stsumplot(jmptimes, indproc, s_times, s_val);
  end

  % center the superposition process
  indmu = on_mu/(on_mu+off_mu);
  cs_val = s_val-nproc*indmu;  

  % expected number of events   
  nevmu = nproc*maxtime*0.5*(1/on_mu+1/off_mu);

  if b_verb
    fprintf('##Expected number of events: %i\n', nevmu); 
  end  
  
  % set the time and space scales according to the scaling regime 
  alpha_min = min(on_alpha, off_alpha);   
  minev = NaN; % can be set explicitely, otherwise set in lracc
  acc = lracc(sc_type, alpha_min, nproc, neff, minev); % acceleration factor  
  [tm_scale, sp_scale] = lrscales(sc_type, alpha_min, nproc, ...
                                  acc, maxtime, b_verb);

    
  if b_rescale_first % first rescale, then integrate

    [rs_times, ircs_val, rcs_val] = strescint(s_times, cs_val, ...
                            tm_scale, sp_scale, b_verb);    
    
    if b_plot  % plot the rescaled processes
               % in the smallest common time window 
      figure(2); 
      clf;
      stairs(rs_times, rcs_val);
    end  
  else % first integrate, then rescale

    [rs_times, ircs_val, ics_val] = stintresc(s_times, cs_val, ...
                             tm_scale, sp_scale, b_verb);
    
    if b_plot  % plot the integrated system size
      figure(2);
      clf;
      hold on;
      stairs(s_times, cs_val);
      plot(s_times, ics_val, 'r');
      hold off;
    end    
  end   

  % plot parts of all processes in a largest common time 
  % window
  if b_plot  % plot the processes

    figure(3); 
    clf;
    plot(rs_times, ircs_val);
  end    
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%