Models and conditioning

A state-space model describes an unobserved process and the measurements we make of it. For example, a sensor might record noisy measurements of a signal that changes over time. Inference uses those measurements to estimate the signal and, when needed, learn the parameters that govern its evolution.

GeneralisedFilters uses the model interface defined by SSMProblems, with three components:

A prior describes the initial state at time zero.
Dynamics describe how the state changes between consecutive times.
An observation process describes a measurement given the current state.

You combine these components in a StateSpaceModel and choose an inference algorithm separately. This lets you reuse a model with different algorithms, or change one component without rewriting the rest of the model.

A linear-Gaussian model

Consider a scalar signal $z_t$ observed with Gaussian noise:

\[z_0 \sim \mathcal{N}(0, 1), \qquad z_t = 0.9 z_{t-1} + w_t, \qquad y_t = z_t + v_t,\]

where $w_t \sim \mathcal{N}(0, 0.1)$ and $v_t \sim \mathcal{N}(0, 0.2)$ are independent. The second arguments here are variances. The corresponding model uses three Gaussian components. Their vector and matrix arguments also support multivariate states.

using GeneralisedFilters, StaticArrays
model = StateSpaceModel(
    GaussianPrior(SA[0.0], SMatrix{1,1}(1.0)),
    LinearGaussianDynamics(SMatrix{1,1}(0.9), SA[0.0], SMatrix{1,1}(0.1)),
    LinearGaussianObservation(SMatrix{1,1}(1.0), SA[0.0], SMatrix{1,1}(0.2)),
)
ys = [SA[0.2], SA[-0.1], SA[0.3]]
state, ll = GeneralisedFilters.filter(model, KF(), ys)
@assert ll ≈ marginal_loglikelihood(model, KF(), ys)
state

GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.15171270461579556], [0.0966783382594036;;])

KF() selects the Kalman filter. It returns the Gaussian distribution of the final state given all observations, along with the log marginal likelihood of the observations. Use marginal_loglikelihood when you only need that scalar likelihood.

This example uses StaticArrays for its small, fixed-size state. Ordinary Julia vectors and matrices work too; see the static arrays example for a performance comparison. Observations start at time one: the filter first propagates the initial state, then conditions on the first observation.

Distribution functions and analytical components

For particle filtering, it is often convenient to describe each conditional distribution with a function. The following model has the same scalar dynamics as the example above:

using Distributions
scalar_model = StateSpaceModel(
    DistributionPrior(Normal(0.0, 1.0)),
    DistributionDynamics((t, z_prev) -> Normal(0.9z_prev, sqrt(0.1))),
    DistributionObservation((t, z) -> Normal(z, sqrt(0.2))),
)

StateSpaceModel{DistributionPrior{Distributions.Normal{Float64}}, DistributionDynamics{Main.var"#3#4"}, DistributionObservation{Main.var"#5#6"}}(DistributionPrior{Distributions.Normal{Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0)), DistributionDynamics{Main.var"#3#4"}(Main.var"#3#4"()), DistributionObservation{Main.var"#5#6"}(Main.var"#5#6"()))

These adapters belong to SSMProblems and are also exported by GeneralisedFilters. DistributionDynamics calls the supplied function with the time and previous state. SSMProblems uses the returned distribution to sample a state or evaluate its log density. DistributionObservation works the same way with the current state. DistributionPrior wraps the initial distribution directly.

The two model definitions expose different information to an algorithm. A function returning a distribution can describe almost any transition, but it does not tell a Kalman filter whether the transition is linear. LinearGaussianDynamics(A, b, Q) exposes that structure: its conditional mean is A * z_prev + b and its covariance is Q. Use these analytical components when you want an algorithm to exploit their structure. There is no need to define a new Julia type for either form.

Time-varying models

A distribution function already receives the time index, so it can look up coefficients or external inputs directly:

coefficients = [0.9, 0.8, 0.7]
time_varying_scalar = StateSpaceModel(
    DistributionPrior(Normal()),
    DistributionDynamics((t, z_prev) -> Normal(coefficients[t] * z_prev, sqrt(0.1))),
    DistributionObservation((t, z) -> Normal(z, sqrt(0.2))),
)

StateSpaceModel{DistributionPrior{Distributions.Normal{Float64}}, DistributionDynamics{Main.var"#9#10"}, DistributionObservation{Main.var"#11#12"}}(DistributionPrior{Distributions.Normal{Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0)), DistributionDynamics{Main.var"#9#10"}(Main.var"#9#10"()), DistributionObservation{Main.var"#11#12"}(Main.var"#11#12"()))

To use the same time dependence with a Kalman filter, return a linear-Gaussian component for each time. TimeVaryingDynamics receives a named tuple containing t:

time_varying_gaussian = StateSpaceModel(
    GaussianPrior(SA[0.0], SMatrix{1,1}(1.0)),
    TimeVaryingDynamics(ctx -> LinearGaussianDynamics(
        SMatrix{1,1}(coefficients[ctx.t]), SA[0.0], SMatrix{1,1}(0.1),
    )),
    LinearGaussianObservation(SMatrix{1,1}(1.0), SA[0.0], SMatrix{1,1}(0.2)),
)
GeneralisedFilters.filter(time_varying_gaussian, KF(), ys)

(GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.13405109323683329], [0.08558267356951614;;]), -2.0358455519811454)

TimeVaryingObservation accepts the corresponding function for the observation process. The package evaluates these functions as it visits each time step. This is called resolving a component. A dynamics function must return a dynamics component, rather than a distribution. For Kalman filtering, that component must expose linear-Gaussian parameters.

Rao-Blackwellised models

Some models become linear and Gaussian only after conditioning on another latent process. For example, suppose a signal's process variance changes over time. Let $x_t$ describe its log variance and $z_t$ the signal:

\[\begin{aligned} x_0 &\sim \mathcal{N}(0, 1), & x_t \mid x_{t-1} &\sim \mathcal{N}(0.9x_{t-1}, q), \\ z_0 &\sim \mathcal{N}(0, 1), & z_t \mid z_{t-1}, x_t &\sim \mathcal{N}(0.8z_{t-1}, \exp(x_t)), \\ && y_t \mid z_t &\sim \mathcal{N}(z_t, 0.2). \end{aligned}\]

The joint model is not linear-Gaussian. However, if we fix the log-variance trajectory $x_{0:T}$, a Kalman filter can integrate out the signal $z_{0:T}$ exactly. A Rao-Blackwellised particle filter uses this structure: it samples trajectories of $x$ and keeps a Gaussian filtering distribution for $z$ within each particle. This reduces the part of the state space that must be represented by samples.

The interface calls $x$ the outer state and $z$ the inner state. The five-component model constructor takes the outer prior and dynamics, followed by the inner prior, dynamics, and observation process. Inner components can depend on the outer state through functions that return the appropriate model component.

using Distributions
function build(q)
    inner_dyn((; t, x_prev, x_new)) = LinearGaussianDynamics(
        SMatrix{1,1}(0.8), SA[0.0], SMatrix{1,1}(exp(x_new)),
    )
    return StateSpaceModel(
        DistributionPrior(Normal()),
        DistributionDynamics((t, x) -> Normal(0.9x, sqrt(q))),
        GaussianPrior(SA[0.0], SMatrix{1,1}(1.0)),
        inner_dyn,
        LinearGaussianObservation(SMatrix{1,1}(1.0), SA[0.0], SMatrix{1,1}(0.2)),
    )
end
hierarchical = build(0.1)

StateSpaceModel{HierarchicalPrior{DistributionPrior{Distributions.Normal{Float64}}, GaussianPrior{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, HierarchicalDynamics{DistributionDynamics{Main.var"#build##1#build##2"{Float64}}, Main.var"#inner_dyn#build##0"}, HierarchicalObservation{LinearGaussianObservation{StaticArraysCore.SMatrix{1, 1, Float64, 1}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}}(HierarchicalPrior{DistributionPrior{Distributions.Normal{Float64}}, GaussianPrior{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(DistributionPrior{Distributions.Normal{Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0)), GaussianPrior{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.0], [1.0;;])), HierarchicalDynamics{DistributionDynamics{Main.var"#build##1#build##2"{Float64}}, Main.var"#inner_dyn#build##0"}(DistributionDynamics{Main.var"#build##1#build##2"{Float64}}(Main.var"#build##1#build##2"{Float64}(0.1)), Main.var"#inner_dyn#build##0"()), HierarchicalObservation{LinearGaussianObservation{StaticArraysCore.SMatrix{1, 1, Float64, 1}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(LinearGaussianObservation{StaticArraysCore.SMatrix{1, 1, Float64, 1}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([1.0;;], [0.0], [0.2;;])))

Combine a bootstrap particle filter for the outer states with a Kalman filter for the inner states to filter this model:

using Random
rng = Xoshiro(42)
particles, loglikelihood_estimate = GeneralisedFilters.filter(
    rng, hierarchical, RBPF(BF(100), KF()), ys,
)
loglikelihood_estimate

-3.193559666799898

Here q is a model parameter, and inner_dyn uses the current outer state to set the signal's process variance. A builder such as build is an ordinary Julia function. It can compute quantities shared across time once and capture data or controls from its enclosing scope.

Each inner component receives a named tuple with the conditioning values available to it:

Component	Input to its function	Returned component
Initial prior	`(; x0)`	A prior for the inner initial state
Dynamics	`(; t, x_prev, x_new)`	Inner dynamics for this transition
Observations	`(; t, x)`	Inner observation process at this time

A component that does not depend on these values can be supplied directly, as the inner prior and observation process are above. If dynamics need a longer outer history, include that history in the outer state so each transition has the information it needs.

Time-varying hierarchical dynamics

Both the outer dynamics and the inner dynamics can depend on time. The outer distribution function receives t as above. The inner function receives t together with the adjacent outer states, so it does not need an additional TimeVaryingDynamics wrapper:

time_varying_hierarchical = StateSpaceModel(
    DistributionPrior(Normal()),
    DistributionDynamics((t, x_prev) -> Normal(coefficients[t] * x_prev, sqrt(0.1))),
    GaussianPrior(SA[0.0], SMatrix{1,1}(1.0)),
    ctx -> LinearGaussianDynamics(
        SMatrix{1,1}(coefficients[ctx.t]), SA[0.0], SMatrix{1,1}(exp(ctx.x_new)),
    ),
    LinearGaussianObservation(SMatrix{1,1}(1.0), SA[0.0], SMatrix{1,1}(0.2)),
)
GeneralisedFilters.filter(Xoshiro(42), time_varying_hierarchical, RBPF(BF(100), KF()), ys)

(GeneralisedFilters.ParticleDistribution{Float64, GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}, Vector{GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}}}(GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}[GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.012153289518600278, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.24719239226471512], [0.16909451166015363;;])), -4.690679854648092, 1), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.035150435361452396, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.24905250326640144], [0.16953610816250364;;])), -4.695362058402319, 2), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.15897550714725048, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2531193700464192], [0.172475698706786;;])), -4.743421031451653, 3), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.40847028304362965, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2612842347120196], [0.1776356190707522;;])), -4.840807412333722, 4), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.42793100581608207, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.22742776892897576], [0.15705924894474566;;])), -4.542445969379019, 6), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.45786085184780373, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.26430685087909983], [0.17844652759170634;;])), -4.855317035095881, 7), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.668686924301909, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.21588871933302162], [0.14920901596714464;;])), -4.466816959128179, 7), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.6252504491859481, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.26704662672735024], [0.1814704451419298;;])), -4.9306500779624685, 9), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.8977580350971313, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.20665829637191085], [0.14048027395815244;;])), -4.3910686657017886, 10), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.0105315868247041, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2465365396035749], [0.169116479713326;;])), -4.69219588346262, 11)  …  GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.5701365900692221, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2236387648726554], [0.151911728375343;;])), -4.485140409046853, 91), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.39005307113583676, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.22932920209074853], [0.15817184600104392;;])), -4.553947210541805, 92), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.7724837676527981, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2715712844545032], [0.18367852164212675;;])), -4.989548651845211, 93), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.4089168651722875, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2602890294162649], [0.17772266314315613;;])), -4.844639985584576, 94), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.6939415278175882, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.26908404986730655], [0.18253651947453878;;])), -4.958340706783532, 95), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-1.0692118520884137, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.19642984229854316], [0.13426634330938603;;])), -4.349818277815521, 96), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.6419512056232445, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.21789069773918343], [0.14999828294836268;;])), -4.472419466971339, 97), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(-0.03337210572568717, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2465744326341681], [0.1678197364435457;;])), -4.669702889422215, 98), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.962532936695221, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2750801749512006], [0.18626586580072904;;])), -5.073696759632735, 99), GeneralisedFilters.Particle{RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}, Float64, Int64}(RBState{Float64, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}}(0.263842324727391, GaussianState{StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SMatrix{1, 1, Float64, 1}}([0.2556115148106051], [0.1748460805387547;;])), -4.787689834317775, 100)], 0.0), -3.306998762456484)

The particle filter resolves the inner dynamics separately for each particle. When the outer trajectory is fixed for a parameter update, the conditional likelihood resolves the same function using the states from that trajectory. The model therefore specifies the same transition for both calculations.

Conditioning on a trajectory

condition_inner turns a hierarchical model and a fixed outer trajectory into an ordinary state-space model. You can then apply the same Kalman likelihood calculation used above.

outer = ReferenceTrajectory(0.1, [0.2, -0.3, 0.4])
inner = condition_inner(hierarchical, outer)
@assert inner_loglikelihood(KF(), hierarchical, outer, ys) ≈
    marginal_loglikelihood(inner, KF(), ys)
trajectory_logdensity(hierarchical, KF(), outer, ys)

-7.19323212253257

ReferenceTrajectory stores the initial state separately and uses time indices 0:T. You may also pass an ordinary vector [x0, x1, ..., xT]. Either representation must contain one more state than there are observations.

The two density evaluations serve different purposes:

inner_loglikelihood evaluates $\log p(y_{1:T} \mid x_{0:T})$, integrating out the inner states.
trajectory_logdensity also adds the outer initial and transition densities, giving $\log p(x_{0:T}, y_{1:T})$. It does not include a prior on model parameters or a parameter-transformation Jacobian.

For parameter inference, differentiate an objective such as θ -> trajectory_logdensity(build(exp(θ[1])), KF(), outer, ys). The outer trajectory stays fixed while parameter changes affect the model. See Particle Gibbs and Turing for a sampler that alternates parameter and trajectory updates.

The conditional model borrows its trajectory. Keep that trajectory unchanged while using the model, and rebuild the conditional model when the trajectory or parameters change. The prior is resolved at construction, while dynamics and observations are resolved as the filter visits each time step.

Custom process and model types

Functions are convenient for individual models. For reusable components or specialised methods, you can instead define subtypes of SSMProblems' StatePrior, LatentDynamics and ObservationProcess. Implement SSMProblems.distribution to obtain the default sampling and density methods, or implement simulate and logdensity directly when no suitable distribution object exists. Algorithms require the methods and analytical structure they use, so a component that works with a bootstrap filter need not support Kalman filtering or backward simulation.

StateSpaceModel is the standard container. Ordinary model operations also accept custom subtypes of AbstractStateSpaceModel through the SSMProblems.prior, SSMProblems.dyn and SSMProblems.obs accessors. Implement these accessors to return the three process components. The abstract model interface does not depend on AbstractMCMC.

Use the five-component StateSpaceModel constructor for the hierarchical models described above. A custom container can also return HierarchicalPrior, HierarchicalDynamics and HierarchicalObservation through its accessors. Conditioning and Rao-Blackwellised algorithms use these component objects. Additional fields on the container do not enter the likelihood unless its components use them.

Covariance storage and differentiation

Gaussian model components retain the arrays you supply, including Diagonal, Symmetric, and real Hermitian covariances. During Kalman filtering, the state distribution uses a full covariance because prediction and conditioning can introduce correlations. SVector means use SMatrix covariances. Dynamic means use ordinary vectors and matrices. This conversion respects the selected triangle of symmetric wrappers and promotes mean and covariance scalar types together, including ForwardDiff dual numbers.

Automatic differentiation follows the parameterisation of structured model covariances. Mooncake also provides a specialised Kalman derivative rule for plain static arrays. Activity annotations are not required. If you use with_activity to mark inactive inputs manually, check those declarations against finite differences.

Covariances must be symmetric, and innovation covariances must be positive definite for the Kalman likelihood. For numerical stability and square-root filtering, see Numerical stability.

The smoother allocates its history after the first update, so later states must retain that storage and scalar type. The Kalman marginal_loglikelihood requires at least one observation and starts its total with the first likelihood increment. It does not force Float32 calculations to accumulate in Float64.

Particle weights and numeric types

Particle weights take their numeric type from the evaluated log-density contributions, which need not have the same type as the latent state. For example, an integer-valued state can have a floating-point likelihood, and ForwardDiff parameters can produce dual weights. Initialization and the first complete filtering step may change these types. Subsequent steps must retain the numeric log-weight type established by that first update. Use consistent scalar types in your model's density calculations across time.

If you implement custom particle-level updates, combine an existing weight and a new contribution with add_logweight. Initial particles and bootstrap proposal corrections can contain an internal marker for an exact zero contribution. This operation handles the marker without choosing a floating-point type. Model density methods should return ordinary real scalars, including a numeric zero when appropriate.

CSMC history preserves the numeric weight type established at the first update. Passing weights of a different type to that history raises an error rather than silently converting them. For differentiation through particle weights, disable resampling. This interface does not define a gradient estimator for discrete ancestor selection.