Generates pseudorandom numbers from a nonhomogeneous Poisson process.
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),
w:here
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.
Generic: CALL RNNPP (TIMBEG, TIMEND, FTHETA, THEMIN, THEMAX, NEUB, NE, R)
Specific: The specific interface names are S_RNNPP and D_RNNPP.
Single: CALL RNNPP (TIMBEG, TIMEND, FTHETA, THEMIN, THEMAX, NEUB, NE, R)
Double: The double precision name is DRNNPP.
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 came 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.
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.
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
Inter-event times for the first 5 events
in the process: 0.0527 0.4080 0.2584 0.0198 0.1676
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
Only 12 events occurred before the time
limit expired. The last event occurred at
time = 18.809
PHONE: 713.784.3131 FAX:713.781.9260 |