nhppp
nhppp is a package for simulating events from one dimensional nonhomogeneous Poisson point processes (NHPPPs). Its functions are based on three algorithms that provably sample from a target NHPPP: the time-transformation of a homogeneous Poisson process (of intensity one) via the inverse of the integrated intensity function; the generation of a Poisson number of order statistics from a fixed density function; and the thinning of a majorizing NHPPP via an acceptance-rejection scheme.
Installation
You can install the release version of nhppp from CRAN with:
install.packages("nhppp")You can install the development version of nhppp from GitHub with:
# install.packages("devtools")
devtools::install_github("bladder-ca/nhppp-fast")Example
These examples use the generic function draw(), which is a wrapper for
the packages specific functions.
Consider the time varying intensity function $\lambda(t) = e^{(0.2t)} (1 + \sin t)$, which is a sinusoidal intensity function with an exponential amplitude. To draw samples over the interval $(0, 6\pi]$ execute
l <- function(t) (1 + sin(t)) * exp(0.2 * t)
nhppp::draw(lambda = l, lambda_maj = l(6 * pi), range_t = c(0, 6 * pi)) |>
head(n = 20)
#> [1] 1.197587 1.238620 1.497499 1.713629 1.761914 2.256739 2.537528 3.622938
#> [9] 5.822574 6.064265 6.645696 6.651551 6.684603 6.875765 6.891348 7.130680
#> [17] 7.446557 7.453139 7.545474 7.557381where lambda_maj is a majorizer constant.
When available, the integrated intensity function $\Lambda(t) = \int_0^t \lambda(s) \ ds$ and its inverse $\Lambda^{-1}(z)$ result in faster simulation times. For this example, $\Lambda(t) = \frac{e^{0.2t}(0.2 \sin t - \cos t)+1}{1.04} + \frac{e^{0.2t} - 1}{0.2}$; $\Lambda^{-1}(z)$ is constructed numerically upfront (or can be calculated numerically by the function, at a computational cost).
L <- function(t) {
exp(0.2 * t) * (0.2 * sin(t) - cos(t)) / 1.04 +
exp(0.2 * t) / 0.2 - 4.038462
}
Li <- stats::approxfun(x = L(seq(0, 6 * pi, 10^-3)), y = seq(0, 6 * pi, 10^-3), rule = 2)
nhppp::draw(Lambda = L, Lambda_inv = Li, range_t = c(0, 6 * pi)) |>
head(n = 20)
#> [1] 0.01152846 0.23558627 0.32924742 0.49921843 0.63509297 1.36677413
#> [7] 2.38941548 3.19511655 3.28049866 4.62140995 5.96916564 6.37504015
#> [13] 6.68283108 6.76577784 7.12919141 7.29249262 7.38665270 7.92953383
#> [19] 7.94791744 7.96591106Function naming conventions
- All functions whose name start with
ppporztpppsample from constant or piecewise constant intensity functions, as described below:
Functions whose names start with
ppp_[sequential|orderstats]sample event times in an interval with constant intensity functions with the sequential and order statistics algorithms.Function
ztppp()samples one or more event times in an interval with constant intensity, i.e., from a zero-truncated Poisson process.Functions
ppp_n()andppp_next_n()samplenevents in an interval and the nextnevent times after a timet0.
- All functions whose name starts with
draworvdrawsample from NHPPPs.
Functions with names starting with
draw_ztsample at least one event in the interval, i.e., from a zero-truncated NHPPP.Functions with names starting with
[draw|draw_zt]_intensity[_majorizer]expect an intensity argument. The third part ([_majorizer]) denotes what, if any, majorizer function is used.Functions with names starting with
[draw|draw_zt]_cumulative_intensity[_algorithm]expect a cumulative (integrated) intensity argument. The third part ([_algorithm]) denotes the algorithm used, if more than one algorithms are pertinent.Functions with names starting with
vdraware vectorized.