pomp provides C entry points to a number of facilities for model specification. Below, we describe these, providing the prototypical function calls and descriptions of their arguments.
The final section describes the prototypes for the basic model components. Users wishing to write libraries to hold basic model components must furnish functions of these prototypes that perform the basic model component computations.
double rbetabinom(double size, double prob, double theta);
double dbetabinom(double x, double size, double prob,
double theta, int give_log);
\(X\) is said to be Beta-binomially distributed with size \(n\), mean probability \(p\), and dispersion parameter \(\theta\) if \[P \sim \mathrm{Beta}\left(\theta\,p,\theta\,(1-p)\right)\] and \[X|P \sim \mathrm{Binomial}\left(n,P\right).\] If \(X\sim\mathrm{BetaBinomial}\left(n,p,\theta\right)\), then \[\mathbb{E}\left[X\right]=n\,p\qquad\text{and}\qquad\mathrm{Var}\left[X\right]=n\,p\,(1-p)\,\frac{\theta+n}{\theta+1}.\]
The R
C API provides a simulator
for the multinomial distribution, rmultinom
. See the
“Rmath.h” header file for information on this facility.
pomp provides an evaluator for the probability mass
function (dmultinom
) for this distribution.
Input:
m
is a positive integer, is the dimension of the random
variable.prob
is a pointer to an m
-vector of
probabilities.x
is a pointer to an m
-vector containing
the data.give_log
is an integer:
give_log=1
if the log probability is desired.give_log=0
if probability is desired.The return value is the probability or log probability (as requested).
The Euler multinomial approximation of a continuous-time, stochastic compartmental model is as follows. Suppose a compartment has occupancy \(N_t\) at time \(t\) and that there are \(K\) ways of exiting the compartment, with per capita rates (hazards) \(\mu_1,\dots,\mu_K\), respectively.
Diagram: A single compartment within a compartmental model. Here, there are \(K=2\) paths out of the compartment.
To make the Euler multinomial approximation, we approximate the total exit rate as constant over a small interval \([t,t+\Delta{t})\). Let the random variable \(\Delta{n_k}\), \(k=1,\dots,K\), be the number that exit by path \(k\) in this time interval and \(\Delta{n_0}\) be the number that remain. Under this assumption, the vector of numbers of exits, \((\Delta{n_{0}},\Delta{n_{1}},\dots,\Delta{n_{K}})\) is multinomially distributed with size \(N_t\) and probabilities \((p_k)_{k=0}^K\), where \[p_0 = \exp\left(-\sum\!\mu_i\,\Delta{t}\right),\] and \[p_k = \frac{\mu_k}{\sum\!\mu_i}\,\left(1-p_0\right),\qquad k=1,\dots,K.\] By way of shorthand, we say that \(\Delta{n}=(\Delta{n_k})_{k=1}^K\) is Euler-multinomially distributed with size \(N_t\), rates \(\mu=(\mu_k)_{k=1}^K\), and time-step \(\Delta{t}\) and we write \[\Delta{n} \sim \mathrm{Eulermultinom}\left(N_t,\mu,\Delta{t}\right).\]
The pomp C API provides three functions that relate to the Euler-multinomial distribution. Their descriptions follow.
The reulermultinom
function draws a random sample from
this distribution. Using the notation above, one has to pack the \(K\) rates \(\mu_1,\dots,\mu_K\) into contiguous memory
locations and retrieve the results in (a different set of) contiguous
memory locations. For example, if rate
is a pointer
to \(K\) contiguous memory locations
holding the rates and dn
is a pointer to \(K\) contiguous memory locations ready to
hold the results, then
will result in a random sample from the Euler multinomial
distribution (with timestep dt) being stored in dn[0]
, …,
dn[K-1]
. In the foregoing, we’ve assumed that the
quantities \(N_t\) and \(K\) are stored in the integer variables
N
and K
, respectively, and that the double
precision variable dt
holds the timestep.
Input:
m
, a positive integer, is number of potential
transitions (“deaths”).size
, a positive integer, is the number of individuals
at risk.rate
is a pointer to the vector of transition (“death”)
rates.dt
, a positive real number, is the duration of time
interval.trans
is a pointer to the vector that will hold the
random deviate.Output:
On return, trans[0]
, …, trans[m-1]
will be
the numbers of individuals making each of the respective
transitions.
See ?reulermultinom
for more on the Euler-multinomial distributions.
NB: reulermultinom
does not call
GetRNGstate()
or PutRNGstate()
internally.
This must be done by the calling program. But note that when
reulermultinom
is called inside a pomp
rprocess, there is no need to call either GetRNGState()
or
PutRNGState()
; this is handled by
pomp.
If \(\Delta{n} \sim
\mathrm{Eulermultinom}\left(N_t,\mu,\Delta{t}\right)\), then the
probability it takes a specific value can be computed using the C
function deulermultinom
. Its prototype is:
double deulermultinom(int m, double size, const double *rate,
double dt, double *trans, int give_log);
Input:
m
, a positive integer, is the number of potential
transitions (“deaths”).size
, a positive integer, is the number of individuals
at risk.rate
is a pointer to vector of transition (“death”)
rates.dt
, a positive real number, is the duration of time
interval.trans
is pointer to vector containing the data, which
are numbers of individuals making the respective transitions.give_log
is an integer:
give_log=0
requests that the probability be
returned.give_log=1
requests that the log probability to be
returned.Output:
The value returned is the probability or log probability (as requested).
See also ?deulermultinom
.
If \(\Delta{n} \sim
\mathrm{Eulermultinom}\left(N_t,\mu,\Delta{t}\right)\), then the
expectation of its \(i\)-th component
is \[\mathbb{E}\left[\Delta{n}_i\right]=p_k
N_t,\] where \(p_k\) is as defined above. The C function
eeulermultinom
computes this. Its prototype is:
Input:
The parameters m
, size
, rate
,
and dt
have the same meaning as above.
Output:
After a call to eeulermultinom
, trans
points to an array of double
s holding the expected
values of the Euler-multinomial random variables.
Corresponding to the R function
rgammawn
, this C function draws a single random increment
of a gamma white-noise process. This will have expectation equal to
dt
and variance sigma^2*dt
.
In particular, when dW = rgammawn(sigma,dt);
is
executed, mu*dW/dt
is a candidate for a random rate process
within an Euler-multinomial context, i.e., mu*dW
will have
expectation mu*dt
and variance
mu*sigma^2*dt
.
void bspline_basis_eval(double x, double *knots, int degree,
int nbasis, double *y);
void bspline_basis_eval_deriv(double x, double *knots, int degree,
int nbasis, int deriv, double *y);
void periodic_bspline_basis_eval(double x, double period, int degree,
int nbasis, double *y);
void periodic_bspline_basis_eval_deriv(double x, double period,
int degree, int nbasis, int deriv, double *y);
These functions work with ordinary and periodic B-spline basis
functions. For ordinary splines, one passes the spline knots in
knots
, together with the degree (degree
)
desired. Note that knots
must point to an array of length
at least nbasis + degree + 1
. The first form evaluates the
nbasis
B-spline basis functions at x
, the
values being returned in the array pointed to by y
. The
second evaluates the order-deriv
derivative of each basis
function.
For periodic splines, the period and smoothness of the functions are
given by period
and degree
, respectively. The
fundamental domain is always assumed to be \([0,T]\), where \(T=\)period
. The first form
evaluates nbasis
periodic B-spline basis functions at
x
, the values being returned in the array pointed to by
y
. The second evaluates the order-deriv
derivative of each basis function.
The logit transformation is defined by \[x = \log\frac{p}{1-p}.\] Its inverse is therefore \[p = \frac{1}{1+e^{-x}}.\]
void to_log_barycentric(double *xt, const double *x, int n);
void from_log_barycentric(double *xt, const double *x, int n);
The log-barycentric transformation takes a vector \(X\in\mathbb{R}^n_+\) to a vector \(Y\in\mathbb{R}^n\), where \[Y_i = \log\frac{X_i}{\sum_j\!X_j}.\] For every \(c>0\), this transformation maps the simplex \(\{X\in\mathbb{R}^n_+\;\vert\;\sum_i\!X_i = c\}\) bijectively onto \(\mathbb{R}^n\).
The pseudo-inverse transformation takes \(\mathbb{R}^n\) to the unit simplex \(S=\{X\in\mathbb{R}^n_+\;\vert\;\sum_i\!X_i=1\}\). Specifically, \[X_i = \frac{e^{Y_i}}{\sum_j\!e^{Y_j}}.\]
Note that if \(T:\mathbb{R}^n_+\to\mathbb{R}^n\) is the log-barycentric transformation so defined, \(U\) is the pseudo-inverse, and \(Id\) denotes the identity map, then \(T\circ U=Id:\mathbb{R}^n\to\mathbb{R}^n\) but \(U\circ T\ne Id\). However, if \(T\) is restricted to the unit simplex S, then \(U\circ{T\vert_{S}}=Id:S\to S\).
Input:
x
is a pointer to vector of parameters to be tranformed
either to or from log barycentric coordinates.n
is the length of this vector.xt
is a pointer to the vector that will hold the
results.On return, xt[0]
, …, xt[n-1]
will contain
the transformed coordinates.
The return value is the dot (inner) product of the
n
-vectors x
and y
. By definition,
\[x{\cdot}y = \sum_i\!x_i\,y_i.\]
This function computes \(r\) such that if \[N \sim \mathrm{Geometric}\left(p=1-e^{-r\,\Delta{t}}\right)\] and \[T \sim \mathrm{Exponential}\left(\mathrm{rate}=R\right),\] then \(\mathbb{E}\left[N\,\Delta{t}\right] = \mathbb{E}\left[T\right]\). That is, \(r = \log{(1+R\,\Delta{t})}/{\Delta{t}}\) is the rate for an Euler process that gives the same expected waiting time as the exponential process it approximates. In particular \(r \to R\) as \(\Delta{t} \to 0\).
const int *get_userdata_int(const char *name);
const double *get_userdata_double(const char *name);
const SEXP get_userdata(const char *name);
The first function returns a pointer to the integer element of the
userdata
with the given name
. The second
retrieves a pointer to a double-precision element of
userdata
by name. The third form returns a general S
expression (SEXP
). See ?userdata
for more information.
pomp provides a facility whereby model codes can be compiled into a dynamically linked library for use in pomp objects. Specifically, basic model components are coded as C functions with the following prototypes.
NB: These functions should not be used within C snippets!
Each of the following functions is supplied one or more of the
stateindex
, parindex
, covindex
,
obsindex
, vmatindex
arguments. Each of these
is an integer vector: the integers within are indices giving the
positions of specific model variables, according to the user’s
specification, the latter being given by means of the
statenames
, paramnames
,
covarnames
, and obsnames
arguments. See ?pomp
for more explanation. Thus, for example, within the body of a function
of prototype pomp_rinit
(see below),
refer to the first state variable, the fourth state variable, the third parameter, and the second covariate, respectively.
void pomp_rinit (double *x, const double *p, double t0,
const int *stateindex, const int *parindex, const int *covindex,
const double *covars);
Description:
p
is a pointer to parameter vector.t0
is the zero time.stateindex
, parindex
,
covindex
: see Indices, above.covars
is a pointer to a vector containing the
(possibly interpolated) values of the covariates at the current
time.x
is a vector that will, on return, contain a draw from
the initial-state distribution.NB: There is no need to call
GetRNGstate()
or PutRNGstate()
in the body of
the user-defined function. The RNG is initialized before any call to
this function, and the RNG state is written afterward. Inclusion of
these calls in the user-defined function may result in significant
slowdown.
void pomp_dinit (double *loglik, const double *x, const double *p, double t0,
const int *stateindex, const int *parindex, const int *covindex,
const double *covars);
Description:
loglik
is a pointer to the scalar that will, on return,
contain the log probability density.x
is the state vector at time t
.p
is a pointer to parameter vector.t0
is the zero time.stateindex
, parindex
,
covindex
: see Indices, above.covars
is a pointer to a vector containing the
(possibly interpolated) values of the covariates at the current
time.step.fun
as used by euler
and
onestep
void pomp_onestep_sim (double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
const double *covars, double t, double dt);
Description:
p
is the parameter vector.stateindex
, parindex
,
covindex
: see Indices, above.covars
is the vector of covariates.t
is the time at the beginning of the step.dt
is the step size (duration of the interval).x
is the vector, that will, on return, contain a draw
from the state process at time t
+dt
.NB: There is no need to call
GetRNGstate()
or PutRNGstate()
in the body of
the user-defined function. The RNG is initialized before any call to
this function, and the RNG state is written afterward. Inclusion of
these calls in the user-defined function may result in significant
slowdown.
rate.fun
as used by gillespie
double pomp_ssa_rate_fn (int event, double t, const double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
const double *covars);
Description:
event
is an integer specifying the number of the
reaction whose rate is desired (the first is event is 1, not 0).t
is the current time.x
is the vector of state variables.p
is the vector of parameters.stateindex
, parindex
,
covindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.The function returns the rate of the requested reaction.
void pomp_dprocess (double *loglik,
const double *x1, const double *x2, double t1, double t2, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
const double *covars);
Description:
t1
, t2
are the times at the beginning and
end of the interval, respectively.x1
, x2
are the state vectors at time
t1
and t2
, respectively.p
is the parameter vector.stateindex
, parindex
,
covindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.loglik
is a pointer to the scalar that will, on return,
contain the log probability density.void pomp_skeleton (double *f, const double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
const double *covars, double t);
Description:
t
is the time.x
is the state vector at time t
.p
is the parameter parameter vector.stateindex
, parindex
,
covindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.f
is a vector, of the same length as x
,
that will, on return, contain the value of the map or vectorfield.void pomp_rmeasure (double *y, const double *x, const double *p,
const int *obsindex, const int *stateindex, const int *parindex,
const int *covindex, const double *covars, double t);
Description:
t
is the time at the beginning of the Euler step.x
is the state vector at time t
.p
is the parameter vector.obsindex
, stateindex
,
parindex
, covindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.y
is a vector that will, on return, contain the
simulated observations.NB: There is no need to call
GetRNGstate()
or PutRNGstate()
in the body of
the user-defined function. The RNG is initialized before any call to
this function, and the RNG state is written afterward. Inclusion of
these calls in the user-defined function may result in significant
slowdown.
void pomp_dmeasure (double *lik, const double *y, const double *x, const double *p,
int give_log, const int *obsindex, const int *stateindex, const int *parindex,
const int *covindex, const double *covars, double t);
Description:
y
is the vector of observables at time
t
.x
is the state vector at time t
.p
is the parameter vector.give_log
is an integer:
give_log=1
if log probability is desired;give_log=0
if probability is desired.obsindex
, stateindex
,
parindex
, covindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.lik
is a pointer to a scalar that will, on return,
contain the requested likelihood or log likelihood.void pomp_emeasure (double *e, const double *x, const double *p,
const int *obsindex, const int *stateindex, const int *parindex, const int *covindex,
const double *covars, double t);
Description:
t
is the time.x
is the state vector at time t
.p
is the parameter parameter vector.stateindex
, parindex
,
covindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.e
is a vector, of the same length as x
,
that will, on return, contain the expected value of the observed
variables.void pomp_vmeasure (double *v, const double *x, const double *p,
const int *vmatindex, const int *stateindex, const int *parindex, const int *covindex,
const double *covars, double t);
Description:
t
is the time.x
is the state vector at time t
.p
is the parameter parameter vector.stateindex
, parindex
,
covindex
, vmatindex
: see Indices, above.covars
is the vector of (possibly interpolated)
covariates at time t
.v
points to a square matrix that will, on return,
contain the covariance matrix of the observed variables. In particular,
if \(X\) is the latent state and the
nobs
observables are \(Y_i\), then
v[vmatindex[i+nobs*j]]
contains \(\mathrm{Cov}[Y_i,Y_j\;\vert\;X]\).Description:
p
is the parameter vector.parindex
: see Indices, above.On return, p
will contain a new random draw from the
prior distribution.
NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function. The RNG is initialized before any call to this function, and the RNG state is written afterward. Inclusion of these calls in the user-defined function may result in significant slowdown.
Description:
p
is the parameter vector.give_log
is an integer:
give_log=1
if log probability is desired;give_log=0
if probability is desired.parindex
: see Indices, above.lik
is a pointer to a scalar that will, on return,
contain the requested probability density or log probability
density.Description:
p
is the parameter vector.parindex
: see Indices, above.pt
is the vector wherein the results will be
returned.