rencount.m

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

function [rntimes, rncount, nevents] = rencount(nproc, maxtime, distr1, ...
                              ren1_par, distr2, ren2_par, b_verb)
%
% RENCOUNT Simulate independent renewal counting processes
%
%    N(t) = max{n: S_n<t}, 
% 
% where 
%
%    S_n = \sum_{k=1}^{n} Y_k,
%
% and {Y_k, k>=2} are i.i.d. r.v. Y_1 may have a different
% distribution, e.g. the stationary one:
%  
%         G(t) = 1/mu_Y * int_0^t (1-F_Y(s)) ds.
% 
% Interrenewal distribution and distribution of the first renewal
% 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.
%
%
% [rntimes, rncount] = rencount(nproc, maxtime, distr1, ren1_par,
%                               distr2, ren2_par, b_verb) 
%
% Example:
% [rntimes, rncount] = rencount(10, 5, @simpareto, {0.6}, @simpareto, ...
%                              {1.6}, 1); 
% generates 10 stationary renewal counting processes on the interval [0, 5)
% with Pareto interrenewal distribution, alpha=1.6.
%
% Inputs:
%        nproc - number of processes to generate
%        maxtime - time window length
%        distr1 - distribution of the first renewal; a handle to
%          the external function generating random numbers from the
%          desired distribution, e.g.
%                @rand: uniform in (0,1)
%                @simexp: Exp(lambda)
%                @simpareto: Pareto(alpha) (see SIMPARETO)
%                @simparetonrm: Normalized Pareto(alpha, gamma) (see SIMPARETONRM)
%                @randn: standard normal
%        ren1_par - a cell array of parameters to the distribution
%          of the first renewal
%        distr2 - distribution of the remaining interrenewal times: a
%          handle to the external random number generator
%        ren2_par - a cell array of parameters to the distribution
%          of the remaining interrenewal times
%        b_verb - if 1 then print processing info
%
% Outputs:
%        rntimes - a matrix of jump times of the processes stored
%          columnwise
%        rncount - a matrix of values of the processes at the
%          jump points stored columnwise
%        nevents - the total number of events generated. It is not
%          the same as number of elements in rntimes, since some
%          vectors are padded with the last value.
%
% See also RENEWPP, RENREW.

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


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

    alpha = 1.4;
    renmu = 1;
    distr1 = @simparetonrm;
    ren1_par = {alpha-1, 1/(renmu*(alpha-1))};
    distr2 = @simparetonrm;
    ren2_par = {alpha, 1/(renmu*(alpha-1))};
    
    b_verb = 1;   
  end
  
  % generate renewal point processes as matrix columns
  [rntimes, ex_i] = renewpp(nproc, maxtime, distr1, ren1_par, ...
              distr2, ren2_par, b_verb);
  
  nren = size(rntimes, 1); % columns generated
  
  % generate the jumps of counting processes as matrix columns; 
  % do not add up as yet since we don't want values for times
  % greater than maxtime 
  rncount = [zeros(1, nproc); ones(nren-1, nproc)];
  % set the counts of the exceeding times to zero
  rncount(ex_i) = 0;
  % add up the jumps
  rncount = cumsum(rncount);    
  
  % plot the counting process
%  stairs(rntimes, rncount);
  
  nevents = nren*nproc-length(ex_i);
  
  if b_verb
    fprintf('##Events generated: %i\n', nren*nproc-length(ex_i)); 
  end