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
!
REAL FUNCTION FTHETA (T)
REAL T
!
REAL EXP
INTRINSIC EXP
!
FTHETA = 0.6342*EXP(0.001427*T)
RETURN
END
Output
Only 12 events occurred before the time
limit expired. The last event occurred at
time = 18.809
