PartiallyObservedMarkovProcesses.jl
The package is a Julia implementation of the pomp package for R.
Package Features
- Implementation of POMP models
- Simulation
- Particle filter
- Iterated filtering
- Trajectory matching
- Workhorses (low-level interface to basic model components)
- Helper functions
- Reproducibility tools
Function Documentation
Implementation of POMP models
Basic constructor
PartiallyObservedMarkovProcesses.pomp — Function
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. Ifdatais supplied as a DataFrame (or AbstractDataFrame),timesshould be a Symbol which denotes the time variable in the DfataFrame. NB: the types oftimesandt0must 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 att0. 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 (seeeuler,onestep,discrete_time, andvectorfield).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.
rprocess plugins
PartiallyObservedMarkovProcesses.euler — Function
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.
PartiallyObservedMarkovProcesses.discrete_time — Function
discrete_time(stepfun; dt = 1)The function stepfun should advance the state by one time unit. The magnitude of the time unit is dt.
PartiallyObservedMarkovProcesses.onestep — Function
onestep(stepfun)The function stepfun will be called once to advance the latent-state process over an interval of arbitrary duration.
PartiallyObservedMarkovProcesses.vectorfield — Function
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.
Simulation
PartiallyObservedMarkovProcesses.simulate — 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.
PartiallyObservedMarkovProcesses.simulate_array — Function
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.
Particle filter
PartiallyObservedMarkovProcesses.pfilter — Function
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.
Iterated filtering
PartiallyObservedMarkovProcesses.mif — Function
mif(object; Np = 1, Nmif = 1, perturbation_kernel, 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 perturbation kernel, cooling schedule, and number of iterations.
Arguments
object: ADataFrame,PompObject, or vector of data.Np: number of particles to use in the filtering.Nmif: number of MIF iterations to perform.perturbation_kernel: 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: seepfilter.params:NamedTupleof model parameters.rinit,rprocess,logdmeasure: necessary basic model components.kwargs...: other arguments are passed topomp.
Perturbation kernel
Cooling schedule
Trajectory matching
PartiallyObservedMarkovProcesses.traj_match_objfun — Function
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.
Workhorses
PartiallyObservedMarkovProcesses.rinit — Function
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 AbstractPompObjectparams: a NamedTuple of parameters or vector of NamedTuplest0: the time at whichrinitis to be simulated. This should be a single scalar.nsim: the number of simulations desired.
PartiallyObservedMarkovProcesses.rinit! — Function
rinit!(object, x0; t0=timezero(object), params = coef(object))rinit! is the in-place version of the rinit workhorse.
PartiallyObservedMarkovProcesses.rprocess — Function
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.
PartiallyObservedMarkovProcesses.rprocess! — Function
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.
PartiallyObservedMarkovProcesses.rmeasure — Function
rmeasure(object; x, times=times(object), params=coef(object))rmeasure is the workhorse for the simulator of the measurement distribution.
PartiallyObservedMarkovProcesses.logdmeasure — Function
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.
PartiallyObservedMarkovProcesses.logdmeasure! — Function
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.
PartiallyObservedMarkovProcesses.logdprior — Function
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.
PartiallyObservedMarkovProcesses.logdprior! — Function
logdprior!(object, ell; params=coef(object))logdprior! is the in-place version of the logdprior workhorse.
PartiallyObservedMarkovProcesses.rprior — Function
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.
Helper functions
PartiallyObservedMarkovProcesses.coef — Function
coef(object, names...)coef extracts the parameters stored in an AbstractPompObject or array thereof.
PartiallyObservedMarkovProcesses.obs — Function
obs(object)obs extracts the vector of observables from a PompObject.
PartiallyObservedMarkovProcesses.states — Function
states(object)states extracts the latent state trajectory of a PompObject.
PartiallyObservedMarkovProcesses.init_state — Function
init_state(object)init_state extracts the latent state at time t0.
PartiallyObservedMarkovProcesses.times — Function
times(object)times extracts the time vector from a PompObject.
PartiallyObservedMarkovProcesses.timezero — Function
timezero(object)timezero extracts the zero-time (t0) from a PompObject.
PartiallyObservedMarkovProcesses.melt — Function
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.
PartiallyObservedMarkovProcesses.logmeanexp — Function
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).
PartiallyObservedMarkovProcesses.traces — Function
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.
Reproducibility tools
PartiallyObservedMarkovProcesses.@freeze — Macro
@freeze [rng,] seed, codeFixes 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().
PartiallyObservedMarkovProcesses.@bake — Macro
@bake dish recipeA 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.