Reference

AbstractPompObject

An AbstractPompObject must, at a minimum, have a method for pomp() defined. This should return a proper PompObject.

MifdPompObject

Created by a call to mif, this struct holds the results of the iterated filtering computation. In particular, one can re-run a mif computation, optionally with modifications to the model, the model parameters, or algorithm settings, by calling mif on a MifdPompObject.

PfilterdPompObject

Created by a call to pfilter, this struct holds the results of the particle-filter computation. In particular, one can re-run a pfilter computation, optionally with modifications to the model, the model parameters, or algorithm settings, by calling pfilter on a PfilterdPompObject.

PompObject

This is the basic data structure for the package. It contains data and model components.

coef(object, names...)

coef extracts the parameters stored in an AbstractPompObject or array thereof.

cond_logLik(object::PfilterdPompObject)

Returns the time series of conditional log likelihoods from the particle filter computation stored in object.

discrete_time(stepfun; dt = 1)

The function stepfun should advance the state by one time unit. The magnitude of the time unit is dt.

eff_sample_size(object::PfilterdPompObject)

Returns the time series of effective sample sizes from the particle filter computation stored in object.

euler(stepfun; dt)

The function stepfun should advance the state by an arbitrary time-increment. The time-increment will be chosen so that equal-sized steps of duration at most dt are taken over any desired interval.

geometric_cooling(frac)

Returns a geometric cooling schedule under which the perturbations are at a fraction frac of their original magnitude after 50 iterations.

init_state(object)

init_state extracts the latent state at time t0.

logLik(object::MifdPompObject)

Returns the esimated log likelihood obtained using a pfilter computation after the final mif iteration.

logLik(object::PfilterdPompObject)

Returns the estimated log likelihood resulting from the particle filter computation stored in object.

logdmeasure!(object, ell; times=times(object), y=obs(object),
             x=states(object), params=coef(object))

logdmeasure! is the in-place version of the logdmeasure workhorse. If no logdmeasure component has been specified, this returns 0 for all inputs.

logdmeasure(object; times=times(object), y=obs(object),
            x=states(object), params=coef(object))

logdmeasure is the workhorse for the evaluator of the log measurement density.

logdprior!(object, ell; params=coef(object))

logdprior! is the in-place version of the logdprior workhorse.

logdprior(object; params=coef(object))

logdprior is the workhorse for the evaluator of the log prior density. If no prior is specified, the logdprior returns 0 for all inputs.

logmeanexp(x; se = false, ess = false)

Compute the log-mean-exp of x. Optionally, return a jack-knife estimate of the standard error (se = true) and/or the effective sample size (ess = true).

melt(x::AbstractPompObject; id...)

Convert an AbstractPompObject to a data frame, with columns for time, observables, and latent states, if present.

melt(x::PfilterdPompObject; id...)

Convert a PfilterdPompObject to a data frame, with columns for time, observables, and latent states, if present, as well as effective sample size and conditional log likelihood.

melt(x::AbstractArray, names::Symbol...)

Convert an array into a data frame, using the variables named in names as identifiers. If x is an array of NamedTuples, the tuples will be arranged into columns of the resulting data frame.

mif(object; Np = 1, Nmif = 1, perturbations, cooling_schedule,
    trigger, target, params, rinit, rprocess, logdmeasure, kwargs...)

Iterated filtering. In addition to the components needed for a pfilter (i.e., Np, trigger, target), one must specify a perturbations function, cooling schedule, and number of iterations. After performing the requested Nmif iterations, mif runs a particle filter using the estimated parameters.

Arguments

• object: A DataFrame, PompObject, or vector of data.
• Np: number of particles to use in the filtering.
• Nmif: number of MIF iterations to perform.
• perturbations: a function that returns perturbed versions of some or all of the model parameters. See below for details.
• cooling_schedule: a function that specifies the MIF cooling schedule. See below for details.
• trigger, target: see pfilter.
• params: NamedTuple of model parameters.
• rinit, rprocess, logdmeasure: necessary basic model components.
• kwargs...: other arguments are passed to pomp.

Return value

mif returns a MifdPompObject. This contains the PfilterdPompObject containing the results of the final particle filter computation. It also records the algorithmic parameters (i.e., Np, Nmif, perturbations and cooling_schedule functions, trigger and target).

Perturbation kernel

Iterated filtering works by applying a random perturbation to some or all of the model parameters during a particle-filter computation. The perturbations argument specifies the nature of the perturbations to be applied. Specifically, it should be function with signature f(scale, lag; params...) where scale is the relative scale of the perturbations (a fraction between 0 and 1) and lag is an integer indicating the observation number (ranging from 1 to n if there are n observations); at the initial time, lag = 0. The named arguments params should be the model parameters to be perturbed. In particular, when the function is called, these arguments will contain the values of the corresponding model parameters.

The function should return a NamedTuple containing the perturbed parameters. Thus, for example, if one is attempting to estimate parameters α and β, while leaving parameters γ and δ fixed, one might furnish a function such as the following as the perturbations argument to mif:

    p(scale, lag; α, β, _...) = begin
        α=rand(LogNormal(log(α),scale*0.02))
        β=rand(LogNormal(log(β),scale*0.02))
        (;α,β)
    end

Note that this function allows for, but ignores, additional arguments (_...).

Initial value parameters

For certain types of parameters, one does not wish to apply the perturbations at every lag. For example, parameters that specify the initial conditions of the latent state process should have perturbations applied only at lag 0. One can use the lag argument of the perturbations function to accommodate these cases. For example, suppose the perameters α, β, γ, and δ mentioned above are regular parameters, but that x₀ is a parameter that fixes the value of the latent state at the zero-time. Then the following perturbations function might be appropriate

    p(scale, lag; α, β, x₀, _...) = begin
        α=rand(LogNormal(log(α),scale*0.02))
        β=rand(LogNormal(log(β),scale*0.02))
        x₀ = (lag == 0) ? rand(LogNormal(log(x₀),scale*0.05)) : x₀
        (;α,β,x₀)
    end

Note that the perturbations are only applied to x₀ at lag 0, i.e., at the zero-time.

Cooling schedule

The cooling schedule is specified by a function which, when furnished a non-negative integer n, returns the fractional reduction of perturbation scale (relative to that determined by the perturbation kernel) that are to be applied in the n-th mif iteration. The function geometric_cooling produces such a function, but the user is free to specify alternative cooling schedules.

obs(object)

obs extracts the vector of observables from a PompObject.

onestep(stepfun)

The function stepfun will be called once to advance the latent-state process over an interval of arbitrary duration.

paramsymbs(object)

Attempt to determine the names of parameters. Return the result as a Vector{Symbol}. If object is an AbstractPompObject without init_state or states, this function may mistakenly include state variables.

pfilter(object; kwargs...)

Calling pfilter on the result of a mif computation runs a particle filter. By default, the parameters Np, trigger, and target used in the mif computation are re-used; one can optionally modify these.

pfilter(object; Np, trigger, target, kwargs...)

Running pfilter on a PfilterdPompObject re-runs the particle filter. Additional arguments adjust the number of particles (Np), trigger and/or target specifications, model parameters, or basic model components.

pfilter(object; Np = 1, params, rinit, rprocess, logmeasure,
        trigger, target, kwargs...)

pfilter runs a sequential Monte Carlo computation, also known as a particle filter. At least the rinit, rprocess, and logdmeasure basic components are needed. Resampling is triggered whenever the effective sample size falls below trigger*Np. The resampling is performed so that the weights are renormalized to the power target, i.e., if target = β, w is a particle weight, and W is the corresponding renormalized weight, then W ∝ wᵝ. One must have 0 ≤ trigger ≤ 1 and 0 ≤ target < 1. trigger = missing is a synonym for trigger = 1 and target = missing is equivalent to target = 0.

As in other POMP.jl functions, kwargs... can be used to modify or unset additional fields in the AbstractPompObject object.

pomp(object::AbstractPompObject; params=missing, accumvars=missing,
     rinit=missing, rprocess=missing, rmeasure=missing,
     logdmeasure=missing)

This form returns a modified version of object. Individual basic components can be modified or removed. The default is to leave them unchanged.

pomp(
    data;
    t0, times, timevar,
    params,
    accumvars,
    rinit, rprocess,
    rmeasure, logdmeasure,
    rprior, logdprior,
    userdata
    )

pomp is the constructor for the PompObject class.

Arguments

• data: observations. The default constructor takes a vector of NamedTuples as data. One can also supply an AbstractDataFrame.
• t0: zero time, t₀.
• times: observation times. If data is supplied as a DataFrame (or AbstractDataFrame), times should be a Symbol which denotes the time variable in the DataFrame. NB: the types of times and t0 must match, or an error will be generated.
• timevar: optional symbol. Name of the time variable.
• params: parameters. A NamedTuple or vector of NamedTuples.
• accumvars: a NamedTuple of state variables to be reset (usually to zero) immediately before each simulation stage.
• rinit: simulator of the latent-state distribution at t0. This component should be a function that takes parameters and, optionally, t0, the initial time. It should return a NamedTuple of state variables.
• rprocess: simulator of the latent-state process. This component should be a plugin (see euler, onestep, discrete_time, and vectorfield).
• rmeasure: simulator of the measurement process. This component should be a function that takes states, parameters, and, optionally, t, the current time. It should return a NamedTuple of observable variables.
• logdmeasure: log pdf of the measurement process. This component should be a function that takes data, states, parameters, and, optionally, t, the current time. It should return a scalar.
• rprior: simulator of the prior distribution on parameters. This component should be a function that takes parameters and returns a NamedTuple of parameters. It should return a NamedTuple of parameters.
• logdprior: log pdf of the prior distribution on parameters. This component should be a function that takes parameters. It should return a scalar.
• userdata: an optional NamedTuple containing elements that will be furnished to each of the basic model components.

rand!(trans, d::EulerMultinomial)

Fill the vector trans with a random draw from the Euler-multinomial distribution defined by d.

rinit!(object, x0; t0=timezero(object), params = coef(object))

rinit! is the in-place version of the rinit workhorse.

rinit(object; t0=timezero(object), params=coef(object), nsim=1)

rinit is the workhorse for the simulator of the initial-state distribution.

Arguments

• object: the AbstractPompObject
• params: a NamedTuple of parameters or vector of NamedTuples
• t0: the time at which rinit is to be simulated. This should be a single scalar.
• nsim: the number of simulations desired.

rmeasure(object; x, times=times(object), params=coef(object))

rmeasure is the workhorse for the simulator of the measurement distribution.

rprior(object; params=coef(object), nsim = 1)

rprior is the workhorse for the simulator of the prior distribution. If nsim > 1, then a matrix is returned.

rprocess!(object, x; x0 = init_state(object),
          t0 = timezero(object), times=times(object),
          params = coef(object))

rprocess! is the in-place version of the rprocess workhorse.

rprocess(object; x0, t0 = timezero(object), times=times(object),
         params = coef(object))

rprocess is the workhorse for the simulator of the process.

If there is no user-supplied rprocess component, the dynamics are trivial.

simulate(object; nsim = 1, params, rinit, rprocess, rmeasure, args...)

Simulate the POMP. Returns an array of PompObjects. At least the rinit, rprocess, and rmeasure basic components, are needed.

simulate_array(object; nsim = 1, params, rinit, rprocess, rmeasure, args...)

Simulate the POMP. At least the rinit, rprocess, and rmeasure basic components, are needed. Return an array containing the simulated sample paths.

states(object)

states extracts the latent state trajectory of a PompObject.

times(object)

times extracts the time vector from a PompObject.

timezero(object)

timezero extracts the zero-time (t0) from a PompObject.

traces(object)

Returns a DataFrame containing the traces of an iterated filtering computation. Specifically, these are the iteration-by-iteration trajectory of the model parameters and the estimated mif likelihood.

traj_match_objfun(object, estimvars; rinit, rprocess,
    logdmeasure, logdprior, bigvalue = Inf, whitelist,
    kwargs...)

Returns a function that compares a model trajectory with the data, returning minus the sum of the log likelihood and log prior density. Trajectory matching consists of the minimization of such an objective function. At least the rinit, rprocess, and logdmeasure basic components are needed. kwargs... can be used to modify or unset these fields. rprocess must be a VectorfieldPlugin (i.e., constructed via a call to vectorfield).

The returned function takes as input a vector or tuple of length equal to that of estimvars. This vector is associated, element-for-element, with the symbols in estimvars. By default, estimvars = keys(coef(object)).

Note that failure of the integrator (signaled via an unsuccessful retcode) will result in an objective-function value of bigvalue and a warning message. Depending on the nature of the objective function, and on the integration algorithm, additional exceptions may arise. In some cases, one may wish to trap these and return an objective-function value of bigvalue. The whitelist argument specifies an exception type (or Union of exception types) that should be handled in this way.

val_array(x)

Stop makin' a fool outta me! Why don'cha come on over...?

vectorfield(vf, integration_alg; integrator_args...)

The function vf should return a vector containing the components of the vectorfield in the same order in which they are passed to vf.

The integration_alg argument specifies the integration algorithm to be used. See OrdinaryDiffEq.jl for more on the choices. The keyword arguments in integrate_args... are passed to the OrdinaryDiffEq.solve() function.

@bake dish recipe

A facility for caching results of computations in a file and retrieving them when needed. @bake checks to see whether the code used to produce the results matches that of the call, and recomputes as needed.

@bake first parses, then computes a hash of the expression recipe. It then checks to see if the file at path dish exists. If it does, and if the hash stored in this file matches the computed hash, the results saved in dish are loaded and returned. If either dish does not exist, or if the digests do not match, it evaluates recipe and stores the result (and the digest) in dish.

@bake uses a serialized representation to store the results.

@freeze [rng,] seed, code

Fixes the state of the pseudorandom number generator stream rng to seed for the duration of the evaluation of code, then restores the state of rng to its original value. By default rng = Random.default_rng().

brownian_motion(;times, t0 = 0, x0, σ, τ)

returns a PompObject for a multivariate Brownian motion process with observations at times. The zero-time is t0 and the starting state at that time is x0. The intensity of the Brownian motion is given by the matrix σ. The measurement error is multivariate normal with variance τ^⊤ τ.

drmca(
     r = 1, K = 1e4, A = 1e3,
     b = 1e-3, c = 1, m = 0.8,
     σ = 0.01,
     N0 = 3000, P0 = 4, t0 = 0.0,
     times = range(start=0,stop=500,step=0.2),
     integrator = AutoTsit5(Rosenbrock23())
   )

Parameters

• r: intrinsic growth rate of prey
• K: carrying capacity for prey
• A: half-saturation prey density
• c: predator foraging rate
• b: predator yield (predators born per prey item killed)
• m: predator death rate
• σ: measurement noise magnitude
• N0, P0: initial densities
• t0: zero-time
• times: vector of observation times
• integrator: integration algorithm.

Observables

• n: prey density
• p: predator density

State variables

• $X = \log(N)$
• $Y = \log(P)$

Details

drmca returns a PompObject containing simulated data from a deterministic Rosenzweig-MacArthur model implemented as a vectorfield. Specifically, the model is the classical Rosenzweig-MacArthur model

\[\begin{aligned} \frac{dN}{dt} &= r N \left(1-\frac{N}{K}\right) - \frac{c N P}{1+N/A} \\ \frac{dP}{dt} &= \frac{b c N P}{1+N/A} - m P \\ \end{aligned}\]

There is also a stochastic version of this model available as an example (see rmca).

In this system, the predator is inviable unless $R = \frac{bcA}{m} > 1$. Even if the predator is viable, the environment is too impoverished to support predators unless $R>1+\frac{A}{K}$. If the environment is rich enough, and if moreover $R>\frac{1+\frac{A}{K}}{1-\frac{A}{K}}$, then the nontrivial equilibrium of the system is unstable. For the default parameters, we have $R = 1.25$ and $\frac{A}{K} = 0.1$, so the latter condition holds.

gompertz()

gompertz is a PompObject containing Parus major data and a simple Gompertz population model. The population model has a single scalar state variable, $X_t$, which obeys

\[X_t = X_{t-1}^S\,K^{1-S}\,\varepsilon_t,\]

where $S = e^{-r\delta{t}}$ and $\varepsilon_t \sim \mathrm{LogNormal}(0,\sigma_p)$. The time-step is one unit: $\delta{t}=1$. The data are assumed to be drawn from a log-normal distribution. In particular,

\[\mathrm{pop}_t \sim \mathrm{LogNormal}(\log{X_t},\sigma_m).\]

Parameters

• r: the growth rate
• K: the equilibrium population density
• X0: the initial population density
• σₚ: process noise s.d.
• σₘ: measurement noise s.d.

rmca(
     r = 1, K = 1e4, A = 1e3,
     b = 1e-3, c = 1, m = 0.8,
     V = 100, σ = 0.01,
     N0 = 3000, P0 = 4, t0 = 0.0,
     δt = 0.01,
     times=range(start=0,stop=500,step=0.2)
   )

Parameters

• r: intrinsic growth rate of prey
• K: carrying capacity for prey
• A: half-saturation prey density
• c: predator foraging rate
• b: predator yield (predators born per prey item killed)
• m: predator death rate
• V: system size
• σ: measurement noise magnitude
• N0, P0: initial densities
• t0: zero-time
• δt: Euler stepsize
• times: vector of observation times

Observables

• n: prey density
• p: predator density

State variables

• $X = \log(N)$
• $Y = \log(P)$

Details

rmca returns a PompObject containing simulated data from a Rosenzweig-MacArthur model implemented as an Itô diffusion. Specifically, if $N$ and $P$ are prey and predator densities, respectively, then $dN = dG - dC - dS$ and $dP = b dS - dM$, where

\[\begin{aligned} dG &= r N dt + \sqrt{\frac{1}{V} r N} dW_1 \\ dC &= \frac{r N^2}{K} dt + \sqrt{\frac{1}{V} \frac{r N^2}{K}} dW_2 \\ dS &= \frac{c N P}{1+N/A} dt + \sqrt{\frac{1}{V} \frac{c N P}{1+N/A}} dW_3 \\ dM &= m P dt + \sqrt{\frac{1}{V} m P} dW_4 \\ \end{aligned}\]

Here, the $dW_i$ are increments of independent standard Wiener processes. Thus, the process noise scales demographically. Specifically, the system size, $V$, converts the densities $N$, $P$ into numbers. It controls the relative magnitude (coefficient of variation) of the demographic process noise. Moreover, $V$ determines a lower threshold on the population sizes, such that if ever $N V < 1$ or $P V < 1$, the population is taken to be extinct. Otherwise, it plays no role in the dynamics. The measurement error is assumed to scale environmentally:

\[\begin{aligned} n &\sim \mathrm{LogNormal}(\log{N},\sigma) \\ p &\sim \mathrm{LogNormal}(\log{P},\sigma) \\ \end{aligned}\]

Note that, in the limit $V\to\infty$, the Itô diffusion becomes the classical Rosenzweig-MacArthur ordinary differential equation model (see drmca):

\[\begin{aligned} \frac{dN}{dt} &= r N \left(1-\frac{N}{K}\right) - \frac{c N P}{1+N/A} \\ \frac{dP}{dt} &= \frac{b c N P}{1+N/A} - m P \\ \end{aligned}\]

sir(
    β = 0.5, γ = 0.25, N = 10000,
    ρ = 0.3, k = 10,
    S0 = 0.9, I0 = 0.01, R0 = 0.1,
    δt = 0.1, t0 = 0.0,
    times = range(start=1.0,stop=90,step=1.0)
   )

sir returns a PompObject containing simulated SIR data.

Parameters

• β: transmission rate
• γ: recovery rate
• N: population size
• ρ: reporting rate
• k: overdispersion coefficient (negative binomial size parameter)
• S0, I0, R0: relative proportions of susceptible, infected, recovered (respectively) in the population at t=t0.
• δt: Euler stepsize
• t0: zero-time
• times: vector of observation times

parus_data

Parus major (Great Tit) census (all individuals), Wytham Wood, Oxfordshire. Global Population Dynamics Database dataset #10163. (NERC Centre for Population Biology, Imperial College (2010) The Global Population Dynamics Database Version 2. http://www.sw.ic.ac.uk/cpb/cpb/gpdd.html).

Original source: McCleery, R. & Perrins, C. (1991) Effects of predation on the numbers of Great Tits, Parus major. In: Bird Population Studies, edited by Perrins, C.M., Lebreton, J.-D. & Hirons, G.J.M. Oxford. Univ. Press. pp. 129–147.