Define some bespoke functions for the different simulation approaches
below. The trick is to define functions so that they work in a
vectorized form. For an example, take a look at the intensity function
below. The other functions have the same bahavior.
Intensity function \(\lambda()\)
Define a vectorized form of the intensity function.
l <- function(t, alpha = pop$alpha, beta = pop$beta, ...) exp(alpha + beta * t)
Arguments alpha and beta can be scalars or
vectors (or column matrices). Argument t can be a scalar, a
vector of \(K\) ages (times), or a
\(K \times s\) matrix. See some
examples of the behavior of the vectorized function. The ellipses
(...) allow l() to ignore extra arguments
without breaking the execution of the script. While we won’t do that on
purpose, some sampling functions in the nhppp package that
take function arguments need the ability to pass optional extraneous
arguments.
Scalar arguments, evaluated at t = 45
l(t = 45, alpha = -4, beta = 0.03)
#> [1] 0.07065121
Scalar arguments, evaluated at t = 45:50 (vector) or as
a column matrix The function returns results in the format of the
t argument.
l(t = 45:50, alpha = -4, beta = 0.03)
#> [1] 0.07065121 0.07280286 0.07502004 0.07730474 0.07965902 0.08208500
l(t = matrix(45:50, ncol = 1), alpha = -4, beta = 0.03)
#> [,1]
#> [1,] 0.07065121
#> [2,] 0.07280286
#> [3,] 0.07502004
#> [4,] 0.07730474
#> [5,] 0.07965902
#> [6,] 0.08208500
Vector arguments, using the first 5 people and evaluating all people
at t = 45
l(t = 45, alpha = pop$alpha[1:5], beta = pop$beta[1:5])
#> [1] 0.04943966 0.06380303 0.05296605 0.03496243 0.06470944
Matrix arguments are convenient: the rows are people and the columns
are the ages (times) – and they can differ across persons. For the first
5 people, we evaluate a different age for each person. The results are
returned in a matrix format.
l(t = matrix(c(45, 50, 45, 47.4, 30), ncol = 1), alpha = pop$alpha[1:5], beta = pop$beta[1:5])
#> [,1]
#> [1,] 0.04943966
#> [2,] 0.07494128
#> [3,] 0.05296605
#> [4,] 0.03691928
#> [5,] 0.04142702
For the first 5 people, we evaluate three different ages for each
person. We arrange a t_mat matrix of ages (times) with
people in the rows and ages (times) for each person in the columns. The
results are returned in a matrix format.
t_mat <- matrix(c(
45, 50, 45, 47.4, 30,
45.1, 50.1, 45.5, 47.8, 38,
48, 52.7, 60.1, 70.1, 99.9
), ncol = 3, byrow = TRUE)
t_mat
#> [,1] [,2] [,3]
#> [1,] 45.0 50.0 45.0
#> [2,] 47.4 30.0 45.1
#> [3,] 50.1 45.5 47.8
#> [4,] 38.0 48.0 52.7
#> [5,] 60.1 70.1 99.9
l(t = t_mat, alpha = pop$alpha[1:5], beta = pop$beta[1:5])
#> [,1] [,2] [,3]
#> [1,] 0.04943966 0.05719509 0.04943966
#> [2,] 0.06892608 0.03937330 0.06400868
#> [3,] 0.06145926 0.05374399 0.05747231
#> [4,] 0.02982757 0.03742537 0.04163741
#> [5,] 0.10137780 0.13647889 0.33101750
The defaults for the alpha and beta
arguments of the l() functions are the respective columns
of the whole population. We can then evaluate the intensity function
\(\lambda()\) for the whole population
passing only the t argument.
t_mat <- matrix(rep(c(45, 50, 55), each = K), ncol = 3, byrow = TRUE)
l(t = t_mat) |> head()
#> [,1] [,2] [,3]
#> [1,] 0.04943966 0.04943966 0.04943966
#> [2,] 0.06380303 0.06380303 0.06380303
#> [3,] 0.05296605 0.05296605 0.05296605
#> [4,] 0.03496243 0.03496243 0.03496243
#> [5,] 0.06470944 0.06470944 0.06470944
#> [6,] 0.08214702 0.08214702 0.08214702
Cumulative intensity function \(\Lambda(t)\)
The cumulative intensity function is \(\Lambda(a, t) = \int_a^t \lambda(s) \
\textrm{d}s\). It does not matter what we choose for the lower
limit of integration in the definition of \(\Lambda\). The lower limit cancels out in
the mathematics of the sampling algorithms. Thus we are free to use any
lower limit (say 0) for any antiderivative of \(\lambda()\) that is convenient (for any
integration constant).
With a slight abuse of notation, define 0 as the lower integration
limit and write \({\Lambda(t) := \Lambda(0, t)
= \int_0^t \lambda(s) \ \textrm{d}s}\). For the log-linear
intensity in this example, \(\Lambda(t) =
\frac{1}{\beta} (e^{\alpha + \beta t} - e^\alpha)\).
The vectorized version of \(\Lambda\) for our example is
L <- function(t, alpha = pop$alpha, beta = pop$beta, ...) {
(exp(alpha + beta * t) - exp(alpha)) / beta
}
Inverse cumulative intensity function \(\Lambda^{-1}(z)\)
By its construction, \(\Lambda\) is
a strictly positive monotone function in \(t\), and thus invertible. The inverse of
the cumulative intensity function \(\Lambda^{-1}\) is defined as the function
that recovers the \(t\) when you pass
it the \(\Lambda(t)\). The definition
is that it satisfies \({\Lambda^{-1} \big (
\Lambda(t) \big ) = t}\). In our example, \(\Lambda^{-1}(z) = \big(\log(\beta z + e^\alpha) -
\alpha\big)/\beta\), which is easily derived from the formula of
\(\Lambda(t)\).
The vectorized implementation of \(\Lambda^{-1}\) for our example is
Li <- function(z, alpha = pop$alpha, beta = pop$beta, ...) {
(log(beta * z + exp(alpha)) - alpha) / beta
}
