onoff.m

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

function [jmptimes, indproc, neff] = onoff(nproc, maxtime, on_dist, on_par, ...
                                     off_dist, off_par, b_verb)

% ONOFF generate N independent stationary on-off processes. An
% on-off process on the real line alternates between values 1 and 0 
% depending if the process is in the on-state or off-state. 
% On-periods and off-periods are iid positive random variables. 
%
%   X(t)=B 1_[0, U_00)(t) + \sum_{n=0}{\infty} 1_[T_n, T_n+U_{n+1})(t)
%
% where 
%       T_0 = B(U_00+V_0)+(1-B)V_00,
%
%       T_n = \sum_{k=1}^{n} (U_n+V_n), n>=1
%
% {U_n, n>=1} are positive iid on-periods, {V_n, n>=1} are positive iid
% off-periods. B is a Bernouli rv such that P(B=1)=EU/(EU+EV), U_00 has
% the stationary on-period distribution, V_00 has the stationary
% off-period distribution. Hence, if B=1 then the process starts with a
% stationary on-period followed by the usual off-period. If B=0 then the
% process starts with a stationary off-period. 
%
% On- and off-distributions are passed into the function as MATLAB
% function handles to external random number generators together with
% two cell-arrays containing distribution parameters.

% The output, two matrices, contain the jump points and the values of
% processes at these points stored columnwise. The processes are
% truncated at the maximal time and may contain different number of
% jumps. To obtain equal lenghts, each process is padded with the
% maximal time and the last value not exceeding the time window.
%
% [jmptimes, indproc] = onoff(nproc, maxtime, on_dist, on_par, 
%                             off_dist, off_par, b_verb)
% Example:
%   [jmptimes, indproc] = onoff(10, 2, @simpareto, {1.4}, ...
%                               @rand, {}, 1);
% generates 10 stationary on-off processes on the interval [0, 2)
% with Pareto(1.6) on-distribution and standard normal
% off-distribution. 
%
% Inputs:
%        nproc - number of processes to generate
%        maxtime - time window length
%        on_dist - distribution of the on-periods: a handle to the
%           external function generating random numbers from the desired
%           distribution, e.g. 
%                @rand: uniform in (0,1)
%                @simexp: Exp(lambda)
%                @simparetonrm: Pareto(alpha, gamma) (see SIMPARETONRM)
%                @randn: standard normal
%        on_par - a cell array of distribution parameters
%        off_dist - distribution of the off-periods: a handle to the
%           external random number generator
%        off_par - a cell array of distribution parameters
%        b_verb - if 1 then print processing info
%
% Outputs:
%        jmptimes - a matrix containing jump times of the indicator
%          processes columnwise
%        indproc - values of the indicator process at the jump
%          points
%        neff - the effective number of periods generated. It is not
%          the same as number of elements in jmptimes, since some
%          vectors are padded with the last value.
%

% Authors: R.Gaigalas, I.Kaj 
% v1.0 12-Oct-05

  % default parameter values
  if (nargin==0)
    nproc = 10;
    maxtime = 5;

    onalpha = 1.4;
    onmu = 1;
    on_dist = @simparetonrm;
    on_par = {onalpha, 1/(onmu*(onalpha-1))};

    offalpha = 2;
    offmu = 1;
    off_dist = @simparetonrm;
    off_par = {offalpha, 1/(offmu*(offalpha-1))};    
    
    b_verb = 1;
  end

  if (b_verb==1)
    fprintf('##Generating on-off processes\n');
    fprintf('  %i processes\n', nproc);
    fprintf('  max time: %f\n', maxtime);
    fprintf('  the on-periods distribution: %s\n', func2str(on_dist));
    fprintf('    distr. parameters: %f\n', on_par{:});
    fprintf('  the off-periods distribution: %s\n', func2str(off_dist));
    fprintf('    distr. parameters: %f\n', off_par{:});    
  end

  % extract the expected values from the parameters
  on_mu = distrmu(on_dist, on_par);
  off_mu = distrmu(off_dist, off_par);
  
  % generate the first on- and off-periods according to the stationary
  % distribution 
  % a sample from Bernoulli: do we start with an on-period?
  oni = find(rand(nproc, 1)<on_mu/(on_mu+off_mu));
  bernrv = zeros(1, nproc);
  bernrv(oni) = 1;
  non = length(oni);

  % get handles to the stationary distributions
  [staton_dist, staton_par] = distrstat(on_dist, on_par);
  [statoff_dist, statoff_par] = distrstat(off_dist, off_par);  

  % parameters for the random number generators
  rndston_par = {1 non staton_par{:}};
  rndoff_par = {1 non off_par{:}};
  rndstoff_par = {1 nproc-non statoff_par{:}};  

  % initialize variables
  onfirst = zeros(1, nproc);
  offfirst = zeros(1, nproc); 
  % some processes start with a stationary on-period followed by 
  % a normal off-period 
  onfirst(oni) = feval(staton_dist, rndston_par{:});
  offfirst(oni) = feval(off_dist, rndoff_par{:});
  % the remaining processes start with a stationary off-period  
  offfirst(find(bernrv==0)) = feval(statoff_dist, rndstoff_par{:});  

  % jump times for the indicator process
  jmptimes = [zeros(1, nproc); onfirst; onfirst+offfirst];
  
  % generate the main on- and off-periods
  % parameters for the random number generators
  rndon_par = {1 nproc on_par{:}};
  rndoff_par = {1 nproc off_par{:}};

  i = 3; 
  while (min(jmptimes(i, :))<=maxtime)
  
    if (rem(i, 2) == 1) % on-period    
      jmptimes = [jmptimes; jmptimes(i, :)+ feval(on_dist, ...
                  rndon_par{:})];
    else % off-period
      jmptimes = [jmptimes; jmptimes(i, :)+ feval(off_dist, ...
                  rndoff_par{:})];
    end
    i = i+1;
  end

  % number of periods generated
  nperiods = i;
  
  % cut off times greater than maxtime
  ex_i = find(jmptimes>maxtime);
  jmptimes(ex_i) = maxtime;
  
  % create the indicator process
  % start with the increment process - avoid assignments at 
  % the exceeding jump times 
  indproc = zeros(nperiods, nproc);
  indproc(1:2:nperiods, :) = 1;
  indproc(2:2:nperiods, :) = -1;
  % set the increments of the exceeding times to zero
  indproc(ex_i) = 0;
  indproc = cumsum(indproc);
  % the starting period is different
  indproc(1, :) = bernrv;

  % the "real" number of periods generated
  neff = nperiods*nproc-length(ex_i);
  
  if (b_verb==1)
    fprintf('##Periods in the longest path: %i\n', nperiods); 
    fprintf('##Overall periods generated: %i\n', nperiods*nproc-length(ex_i)); 
  end