renewpp.m

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

function [rntimes, ex_i] = renewpp(nproc, maxtime, distr1, ren1_par, ...
                                   distr2, ren2_par, b_verb)

% RENEWPP Generate a matrix of N independent renewal point
% processes {S_n} stored columnwise, i.e. 
%
%    S_n = \sum_{k=1}^{n} Y_k,
%
% {Y_k, k>=2} are i.i.d. r.v., while Y_1 may follow a different
% distribution e.g. the stationary distribution. Distributions are
% provided as handles to the external functions and parameters are
% stored as cell-arrays.
%
% The processes are truncated at maximal time and may contain
% different number of points. To obtain equal lenghts, each vector
% is padded with the maximal time.
%
% [rntimes, ex_i] = renewpp(nproc, maxtime, distr1, ren1_par,
%                           distr2, ren2_par, b_verb)
%
% Inputs:
%        nproc - number of processes to generate
%        maxtime - time window length
%        distr1 - distribution of the first renewal: a pointer 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
%        ren1_par - a cell array of distribution parameters
%        distr2 - distribution of the remaining renewals: a pointer to
%           the external function generating random numbers 
%        ren2_par - a cell array of distribution parameters
%        b_verb - if 1 then print processing info
%
% Outputs:
%        rntimes - a matrix containing renewal times of the
%          processes columnwise 
%        ex_i - indices of the discarded renewal times exceeding
%          the time window
%
% See also RENCOUNT, RENREW.


% Authors: R.Gaigalas, I.Kaj  
% v1.0 13-Nov-05

  
  if (b_verb==1)
    fprintf('##Generating renewal point processes\n');
    fprintf('  %i processes\n', nproc);
    fprintf('  max time: %f\n', maxtime);
    fprintf('  the first renewal, distribution: %s\n', func2str(distr1));
    fprintf('    distr. parameters: %f\n', ren1_par{:});
    fprintf('  the remaining renewals, distribution: %s\n', func2str(distr2));
    fprintf('    distr. parameters: %f\n', ren2_par{:});    
  end

  % parameters for the random number generators
  rnd_params1 = {1 nproc ren1_par{:}};
  rnd_params2 = {1 nproc ren2_par{:}};

  % the first renewal
  rntimes = [zeros(1, nproc); feval(distr1, rnd_params1{:})];

  % generate renewal processes as matrix columns
  i = 2; 
  while (min(rntimes(i, :))<=maxtime)
    rntimes = [rntimes; rntimes(i, :)+ feval(distr2, ...
               rnd_params2{:})];
    i = i+1;
  end

%  n_points = i; % length of the trajectories  

  % cut off renewal times greater than maxtime
  ex_i = find(rntimes>maxtime);
  rntimes(ex_i) = maxtime;
  
  
  if (b_verb==1)
    fprintf(1, 'Maximal number of jumps: %i\n\n', i); 
  end