PartiallyObservedMarkovProcesses.jl

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.

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.

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.

discrete_time(stepfun; dt = 1)

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

onestep(stepfun)

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

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.

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.

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.

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; 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.

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.

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.

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.

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

rinit! is the in-place version of the rinit 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.

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.

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

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

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.

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.

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.

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

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

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.

coef(object, names...)

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

obs(object)

obs extracts the vector of observables from a PompObject.

states(object)

states extracts the latent state trajectory of a PompObject.

init_state(object)

init_state extracts the latent state at time t0.

times(object)

times extracts the time vector from a PompObject.

timezero(object)

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

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.

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.

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).

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.

@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().

@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.