RNNPP

Generates pseudorandom numbers from a nonhomogeneous Poisson process.

Required Arguments

TIMBEG — Lower endpoint of the time interval of the process. (Input)
TIMBEG must be nonnegative. Usually, TIMBEG = 0.

TIMEND — Upper endpoint of the time interval of the process. (Input)
TIMEND must be greater than TIMBEG.

FTHETA — User‑supplied FUNCTION to provide the value of the rate of the process as a function of time. This function must be defined over the interval from TIMBEG to TIMEND and must be nonnegative in that interval. The form is FTHETA(TIME),
where

TIME — Time at which the rate function is evaluated. (Input)

FTHETA — Value of the rate function. (Output)

FTHETA must be declared EXTERNAL in the calling program.

THEMIN — Minimum value of the rate function FTHETA in the interval (TIMBEG, TIMEND). (Input)
If the actual minimum is unknown, set THEMIN = 0.0.

THEMAX — Maximum value of the rate function FTHETA in the interval (TIMBEG, TIMEND). (Input)
If the actual maximum is unknown, set THEMAX to a known upper bound of the maximum. The efficiency of RNNPP is less the greater THEMAX exceeds the true maximum.

NEUB — Upper bound on the number of events to be generated. (Input)
In order to be reasonably sure that the full process through time TIMEND is generated, calculate NEUB as NEUB = X + 10.0 * SQRT(X), where X = THEMAX * (TIMEND ‑ TIMBEG). The only penalty in setting NEUB too large is that the output vector must be dimensioned of length NEUB.

NE — Number of events actually generated. (Output)
If NE is less that NEUB, the time TIMEND is reached before NEUB events are realized.

R — Vector of length NE containing the times to events. (Output)
R must be dimensioned to be of length NEUB.

FORTRAN 90 Interface

Generic: CALL RNNPP (TIMBEG, TIMEND, FTHETA, THEMIN, THEMAX, NEUB, NE, R)

Specific: The specific interface names are S_RNNPP and D_RNNPP.

FORTRAN 77 Interface

Single: CALL RNNPP (TIMBEG, TIMEND, FTHETA, THEMIN, THEMAX, NEUB, NE, R)

Double: The double precision name is DRNNPP.

Description

Routine RNNPP simulates a one‑dimensional nonhomogeneous Poisson process with rate function THETA in a fixed interval (TIMBEG, TIMEND].

Let λ(t) be the rate function and t0 = TIMBEG and t1 = TIMEND. Routine RNNPP uses a method of thinning a nonhomogeneous Poisson process {N*(t), t ≥ t0} with rate function λ*(t) ≥ λ(t) in (t0, t1], where the number of events, N*, in the interval (t0, t1] has a Poisson distribution with parameter

The function

is called the integrated rate function. In RNNPP, λ*(t) is taken to be a constant λ*(= THEMAX) so that at time ti, the time of the next event ti+1 is obtained by generating and cumulating exponential random numbers

with parameter λ*, until for the first time

where the uj,i are independent uniform random numbers between 0 and 1. This process is continued until the specified number of events, NEUB, is realized or until the time, TIMEND, is exceeded. This method is due to Lewis and Shedler (1979), who also review other methods. The most straightforward (and most efficient) method is by inverting the integrated rate function, but often this is not possible.

If THEMAX is actually greater than the maximum of λ(t) in (t0, t1], the routine will work, but less efficiently. Also, if λ(t) varies greatly within the interval, the efficiency is reduced. In that case, it may be desirable to divide the time interval into subintervals within which the rate function is less variable. This is possible because the process is without memory.

If no time horizon arises naturally, TIMEND must be set large enough to allow for the required number of events to be realized. Care must be taken, however, that FTHETA is defined over the entire interval.

After simulating a given number of events, the next event can be generated by setting TIMBEG to the time of the last event (the sum of the elements in R) and calling RNNPP again. Cox and Lewis (1966) discuss modeling applications of nonhomogeneous Poisson processes.

Comments

The routine RNSET can be used to initialize the seed of the random number generator. The routine RNOPT can be used to select the form of the generator.

Examples

Example 1

In this example, RNNPP is used to generate the first five events in the time 0 to 20 (if that many events are realized) in a nonhomogeneous process with rate function

λ(t) = 0.6342 e0.001427t

for 0 < t ≤ 20.

Since this is a monotonically increasing function of t, the minimum is at t = 0 and is 0.6342, and the maximum is at t = 20 and is 0.6342 e0.02854 = 0.652561.

USE RNNPP_INT

USE UMACH_INT

USE RNSET_INT

IMPLICIT NONE

INTEGER NEUB

PARAMETER (NEUB=5)

INTEGER I, ISEED, NE, NOUT

REAL FTHETA, R(NEUB), THEMAX, THEMIN, TIMBEG, TIMEND

EXTERNAL FTHETA

CALL UMACH (2, NOUT)

TIMBEG = 0.0

TIMEND = 20.0

THEMIN = 0.6342

THEMAX = 0.652561

ISEED = 123457

CALL RNSET (ISEED)

CALL RNNPP (TIMBEG, TIMEND, FTHETA, THEMIN, THEMAX, NEUB, NE, R)

WRITE (NOUT,99999) NE, (R(I),I=1,NE)

99999 FORMAT (' Inter-event times for the first ', I1, ' events', /, &

' in the process: ', 5F7.4)

END

REAL FUNCTION FTHETA (T)

REAL T

REAL EXP

INTRINSIC EXP

FTHETA = 0.6342*EXP(0.001427*T)

RETURN

END

Output

Inter-event times for the first 5 events

in the process: 0.0527 0.4080 0.2584 0.0198 0.1676

Example 2

As it turns out in the simulation above, the first five events are realized before time equals 20. If it is desired to continue the simulation to time equals 20, setting NEUB to 49 (that is,

would likely ensure that the time is reached. In the following example, we see that there are twelve events realized by time equals 20.

USE UMACH_INT

USE RNSET_INT

USE RNNPP_INT

USE SSUM_INT

IMPLICIT NONE

INTEGER NEUB

PARAMETER (NEUB=49)

INTEGER ISEED, NE, NOUT

REAL FTHETA, R(NEUB), T, THEMAX, THEMIN, TIMBEG, TIMEND

EXTERNAL FTHETA

CALL UMACH (2, NOUT)

TIMBEG = 0.0

TIMEND = 20.0

THEMIN = 0.6342

THEMAX = 0.652561

ISEED = 123457

CALL RNSET (ISEED)

CALL RNNPP (TIMBEG, TIMEND, FTHETA, THEMIN, THEMAX, NEUB, NE, R)

T = TIMBEG + SSUM(NE,R,1)

IF (NE .LT. NEUB) THEN

WRITE (NOUT,99998) NE, T

99998 FORMAT (' Only ', I2, ' events occurred before the time', &

/, ' limit expired. The last event occurred at', /, &

' time = ', F6.3)

ELSE

WRITE (NOUT,99999) NE, T

99999 FORMAT (' Possibly more than ', I2, ' events would have', &

/, ' occurred before the time limit expired.', /, &

' The last event occurred at time = ', F6.3)

END IF

END