randomNpp¶

Generates pseudorandom numbers from a nonhomogeneous Poisson process.

Synopsis¶

randomNpp (tbegin, tend, ftheta(), thetaMin, thetaMax, neub, ne)

Required Arguments¶

float tbegin (Input): Lower endpoint of the time interval of the process. tbegin must be nonnegative. Usually, tbegin = 0.
float tend (Input): Upper endpoint of the time interval of the process. tend must be greater than tbegin.
float ftheta(float t) (Input): 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 tbegin to tend and must be nonnegative in that interval.
float thetaMin (Input): Minimum value of the rate function ftheta() in the interval (tbegin, tend). If the actual minimum is unknown, set thetaMin = 0.0.
float thetaMax (Input): Maximum value of the rate function ftheta() in the interval (tbegin, tend). If the actual maximum is unknown, set thetaMax to a known upper bound of the maximum. The efficiency of randomNpp is less the greater thetaMax exceeds the true maximum.
int neub (Input): Upper bound on the number of events to be generated. In order to be reasonably sure that the full process through time tend is generated, calculate neub as neub = X + 10.0 × sqrt(X), where X = thetaMax × (tend – tbegin).
int ne (Output): Number of events actually generated. If ne is less that neub, the time tend is reached before neub events are realized.

Return Value¶

An array of length neub containing the times to events in the first ne elements.

Description¶

Function randomNpp simulates a one-dimensional nonhomogeneous Poisson process with rate function ftheta in a fixed interval (tbegin, tend].

Let \(\lambda(t)\) be the rate function and \(t_0\) = tbegin and \(t_1\) = tend. Function randomNpp uses a method of thinning a nonhomogeneous Poisson process \(\{ N*(t),t\geq t_0\ \}\) with rate function \(\lambda(t)\geq\lambda *(t)\) in \(\left(t_0,t_1\right]\), where the number of events, N*, in the interval \(\left(t_0,t_1\right]\) has a Poisson distribution with parameter

\[\mu_0 = \int_{t_0}^{t_1} \lambda(t) dt\]

The function

\[\mathit{\Lambda}(t) = \int_0^{t'} \lambda(t) dt\]

is called the integrated rate function. In randomNpp, \(\lambda*(t)\) is taken to be a constant λ*(= thetaMax) so that at time \(t_i\), the time of the next event \(t_{i+1}\) is obtained by generating and cumulating exponential random numbers

\[E_{1,i}^*, E_{2,i}^*, \ldots\]

with parameter λ*, until for the first time

\[u_{j,i} \leq \left(t_i + E_{1,i}^* + \ldots + E_{j,i}^*\right) / \lambda^*\]

where the \(u_{j,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, tend, 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 thetaMax is actually greater than the maximum of \(\lambda(t)\) in \((t_0,t_1]\), the function will work, but less efficiently. Also, if \(\lambda(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, tend 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 tbegin to the time of the last event (the sum of the elements in R) and calling randomNpp again. Cox and Lewis (1966) discuss modeling applications of nonhomogeneous Poisson processes.

Example¶

In this example, randomNpp 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 \exp(0.001427t)\]

for \(0<t\leq 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 \exp(0.02854)=0.652561\).

from __future__ import print_function
from numpy import *
from pyimsl.stat.randomNpp import randomNpp
from pyimsl.stat.randomSeedSet import randomSeedSet


def ftheta(t):
    return 0.6342 * exp(0.001427 * t)


neub = 5
tmax = .652561
tmin = .6342
tbeg = 0.
tend = 20.
ne = []
randomSeedSet(123457)
r = randomNpp(tbeg, tend, ftheta, tmin, tmax, neub, ne)
print("Inter-event times for the first ", ne[0], "events in the process:")
for i in range(0, ne[0]):
    print("%15.6f\t" % r[i])

Output¶

Inter-event times for the first  5 events in the process:
       0.052660	
       0.407978	
       0.258399	
       0.019766	
       0.167641

Fatal Errors¶

IMSLS_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”.