scaleren.m

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

function [crs_times, crs_val] = scaleren(sc_type, alpha, ren_mu, ...
          maxtime, nproc, b_rescale_first, b_plot, b_verb)  

% SCALEREN Simulate a rescaled and centered sum of independent
% stationary renewal counting processes with heavy-tailed
% interrenewal distribution:
%
%     Z_m(t) = 1/b * \sum_{i=1}^{m} (N_i(a*t) - a*t/ren_mu),
%
% where a=a(m) is the time scale, b=b(m) the space scale, N_i(t)
% are independent copies of
%
%     N(t) = max{n: S_n<t},      S_n = \sum_{k=1}^{n} Y_k,
%
% {Y_k, k>=2} are iid normalized Pareto random variables 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.
%
% Y_1 has the stationary distribution:
%  
%         G(t) = 1/mu_Y * int_0^t (1-F_Y(s)) ds,
%
% which is also normalized Pareto with alpha1=alpha-1 and
% gamma1=gamma. 
%
% The process Z_m(t) is a random piecewise linear function with
% jumps (a "saw"). It is represented by two vectors of the jump or
% break points and the values at these points respectively.
% 
% It is known that if the interrenewal 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(m), as m grows to infinity,
% the process Z_m(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(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)
%
% For details see
%
% R.Gaigalas, I.Kaj (2003) Convergence of scaled renewal processes
% and a packet arrival model. Bernoulli  9(4), 671--703.
%
% Usage:
% [crs_times, crs_val] = scaleren(sc_type, alpha, ren_mu, 
%            nproc, 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 
%        alpha - tail parameter of the interrenewal distribution:
%          normalized Pareto, 1<alpha<=2 (see above)
%        ren_mu - expected value of the interrenewal
%          distribution. Should be positive
%        nproc - number of processes to superpose
%        maxtime - maximal time for the process before rescaling
%        b_rescale_first - if 1 then first rescale, then center. If
%          0 then the other way around.
%        b_plot - if 1 then plot the processes
%        b_verb - if 1 then print processing info
%
% Outputs:
%        crs_times - a column vector of the jump times of the
%          rescaled process(es) 
%        crs_val - a column vector of values of the rescaled
%          superposition process(es) at the jump times
%
% See also SCALEMGINFTY, SCALEONOFF, LRSCALES.

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


  if (nargin<1) % default parameter values
    sc_type = [1 2 3];

    alpha = 1.4;
    ren_mu = 1; 

    maxtime = 100;
    nproc = 100; 
    b_rescale_first = 1;
    b_plot = 1;
    b_verb = 1;
  end

  % interrenewal distributions: normalized Pareto
  distr1 = @simparetonrm;
  ren1_par = {alpha-1, 1/(ren_mu*(alpha-1))};
  distr2 = @simparetonrm;
  ren2_par = {alpha, 1/(ren_mu*(alpha-1))};  

  % generate renewal counting processes as column vectors   
  [a_times, a_val, nevents] = rencount(nproc, maxtime, distr1, ...
                          ren1_par, distr2, ren2_par, b_verb);
  
  % sum up the piecewise constant processes
  [s_times, s_val] = stairsum(a_times, a_val); 

  if b_plot
    figure(1);
    stsumplot(a_times, a_val, s_times, s_val);
  end

  % expected number of events 
  nevmu = nproc*maxtime/ren_mu;
  
  if b_verb
    fprintf('##Expected number of events: %i\n', nevmu); 
  end
    
  % 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, nproc, nevents, minev); % acceleration factor
  [tm_scale, sp_scale] = lrscales(sc_type, alpha, nproc, ...
                                  acc, maxtime, b_verb);
  
  if b_rescale_first % first rescale, then center

    [crs_times, crs_val, rs_times, rs_val] = ren_rescale_center(...
    s_times, s_val, tm_scale, sp_scale, nproc, ren_mu, b_verb); 

    if b_plot  % plot the rescaled versions of the sum
               % in the largest common time window 
      figure(2); 
      clf;
      stairs(rs_times, rs_val);
    end  
  else % first center, then rescale
    
    [crs_times, crs_val, cs_times, cs_val] = ren_center_rescale(...
     s_times, s_val, tm_scale, sp_scale, nproc, ren_mu, b_verb);
    
    % plot the centered renewal counting processes together with
    % the centered sum
    if b_plot  
      figure(2);
      plot_centered(a_times, a_val, cs_times, cs_val, ren_mu);
    end
  end   

  % plot parts of the rescaled and centered versions in the
  % largest common time window
  if b_plot
    figure(3); 
    clf;
    plot(crs_times, crs_val);
  end  

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [crs_times, crs_val, rs_times, rs_val] = ren_rescale_center(...
s_times, s_val, tm_scale, sp_scale, nproc, ren_mu, b_verb)
%
% first rescale, then center
%

  % rescale time
  % N(at) jumps at S_n/a
  rs_times = s_times*(1./tm_scale)';
  % rescale space 
  rs_val = s_val*(1./sp_scale)';
  
  % cut the piecewise constant processes at the smallest common time
  cut_time = min(rs_times(end, 1:end));
  [rs_times rs_val] = staircut(rs_times, rs_val, cut_time);

  if b_verb % print info
    fprintf(1, '##Largest common time window: %f\n', cut_time);
  end  

  % center the rescaled sum
  % convert the stair processes into piecewise linear processes
  [crs_times crs_val] = stairs(rs_times, rs_val);
  crs_times = crs_times(1:end-1, :); % delete the copy of the last point
  crs_val = crs_val(1:end-1, :);  
  % subtract the expected value
  % OBS! crs_times*diag(tm_scale) should be used instead of s_times
  % since crs_times is truncated and doubled
  crs_val = crs_val - nproc/ren_mu* ...
              crs_times*diag(tm_scale./sp_scale);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rcs_times, rcs_val, cs_times, cs_val] = ren_center_rescale(...
s_times, s_val, tm_scale, sp_scale, nproc, ren_mu, b_verb)
%
% first center, then rescale
%

  % convert the stair processes into piecewise linear processes
  [cs_times cs_val] = stairs(s_times, s_val);
  % subtract the expected value
  cs_val = cs_val-nproc/ren_mu*cs_times; 
  
  % rescale time
  % N(at) jumps at S_n/a
  rcs_times = cs_times*(1./tm_scale)';
  % rescale space
  rcs_val = cs_val*(1./sp_scale)';

  % cut the piecewise linear processes at the smallest common time
  cut_time = min(rcs_times(end, 1:end));
  [rcs_times rcs_val] = linjpcut(rcs_times, rcs_val, cut_time);

  
  if b_verb % print info
    fprintf(1, '##Largest common time window: %f\n', cut_time);
  end  
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = plot_centered(a_times, a_val, cs_times, cs_val, ...
                            ren_mu) 
%
% Subtract the expectation from the renewal-rewards processes and plot
% together with their sum
%

  % convert the stair processes into piecewise linear processes
  [ca_times, ca_val] = stairs(a_times, a_val);
  % subtract the expected value  
  ca_val = ca_val-ca_times/ren_mu;

      
  clf;
  plot(ca_times, ca_val, 'b', cs_times, cs_val, 'r');   

%  hold on;
  
%  plot(ca_times, ca_val, 'b'); 
%  plot(cs_times, cs_val, 'r'); 
  
%  hold off;

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