renuni.m

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

function [rntimes, rncount] = renuni(nproc, maxtime)

% RENUNI Generate a matrix of N independent renewal counting
% processes 
%
%    N(t) = max{n: S_n<t},
%
% where 
%
%    S_n = \sum_{k=1}^{n} Y_k,
%
% and {Y_k, k>=1} are i.i.d. uniformly on (0,1) distributed random
% variables. 
%
% 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] = renuni(nproc, maxtime)
%
% Inputs:
%        nproc - number of processes to generate
%        maxtime - time window length
%
% Outputs:
%        rntimes - a matrix containing renewal times of the
%          processes columnwise 
%        rncount - a matrix containing values of the counting
%          processes at the renewal times columnwise 
%
% See also RENCOUNT, RENREW.

% Author: R.Gaigalas 
% v1.0 12-Nov-05

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

  % add zero to the renewal times for nicer plots 
  rntimes = zeros(1, nproc);

  % generate renewal processes as matrix columns
  i = 1; 
  while (min(rntimes(i, :))<=maxtime)
    rntimes = [rntimes; rntimes(i, :)+rand(1, nproc)];
    i = i+1;
  end

  % cut off renewal times greater than maxtime
  ex_i = find(rntimes>maxtime);
  rntimes(ex_i) = maxtime;  
  
  % 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(size(rntimes, 1)-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);