renrew.m

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

function [rntimes, pval] = renrew(nproc, maxtime, ren1_dist, ...
         ren1_par, ren2_dist, ren2_par, b_verb, rew_dist, rew_par)
%
% RENREW Generate N independent renewal reward processes
%
%     Z(t)= \sum_{k=1}^{N(t)} X_k,
%
% where N(t) is the counting process of a renewal process {S_n}, 
%
%    N(t) = max{n: S_n<t},     S_n = \sum_{k=1}^{n} Y_k,
%
% {Y_k, k>=2} are i.i.d. r.v. while 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.
%
% Distribution of rewards {X_k} is arbitrary. It is passed into
% the function as a MATLAB function handle to an external random
% number generator together with a cell-array 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, pval] = renrew(alpha, mu, maxtime, nproc, b_verb,
%		             [rew_dist, rew_par])
%
% Inputs:
%        nproc - number of processes to superpose
%        maxtime - time window length
%        ren1_dist - distribution of the first renewal; a handle to
%          the external function generating random numbers from the
%          distribution:
%                @rand: uniform in (0,1)
%                @simexp: Exp(lambda)
%                @simparetonrm: Pareto(alpha, gamma) (see SIMPARETONRM)
%                @randn: standard normal
%        ren1_par - a cell array of parameters to the distribution
%          of the first renewal
%        ren2_dist - distribution of the remaining renewals: a
%          handle to the external random nuber generator
%        ren2_par - a cell array of parameters to the distribution
%          of the remaining renewals
%        b_plot - if 1 then plot the processes
%        b_verb - if 1 then print processing info
%        rew_dist - optional, distribution of the rewards; a handle
%          to the external random number generator. The default
%          value @ones gives the renewal counting process. 
%        rew_par - optional, additional parameters for the
%          distribution of rewards 
%
% Outputs:
%        rntimes - a matrix of jump times of the processes stored
%           columnwise
%        pval - a matrix of values of the processes at the
%             jump points stored columnwise
%
% See also RENCOUNT, RENEWPP.

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


 if (nargin<1) % default parameter values
   alpha = 1.4;
   mu = 1;
   maxtime = 5;
   nproc = 2;
   b_verb = 1;

   rew_dist = @simparetonrm;
   rew_par = {1.6, 2};
  end

  if (nargin==5) % if the rewards distribution is not given
   rew_dist = @ones; % the default is deterministic
   rew_par = {};
  end
  
  % generate stationary renewal point processes as matrix columns
  [rntimes, ex_i] = renewpp(nproc, maxtime, ren1_dist, ren1_par, ...
              ren2_dist, ren2_par, b_verb);

  if (b_verb==1)
    fprintf('##Generating reward processes\n');
    fprintf('  Rewards, distribution: %s\n', func2str(rew_dist));
    fprintf('    distr. params: %f\n', rew_par{:});
    fprintf('\n');
   end

  nren = size(rntimes, 1); % columns generated   
   
  % parameters for the random number generator
  rnd_par = {nren-1 nproc rew_par{:}};
  % generate rewards as matrix columns; 
  % do not add up as yet since we don't want values for times
  % greater than maxtime 
  pval = [zeros(1, nproc); feval(rew_dist, rnd_par{:})];
	       
  % set the rewards of the exceeding times to zero
  pval(ex_i) = 0;
  
  % sum up the rewards
  pval = cumsum(pval);  

  if b_verb
    fprintf('##Events generated: %i\n', ...
            nren*nproc-length(ex_i)); 
  end