Exported functions and types

AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, dynamics_density, initial_density; p = SciMLBase.NullParameters(), threads = false, kwargs...)

This type represents a standard particle filter but affords extra flexibility compared to the ParticleFilter type, e.g., non-additive noise in the dynamics and measurement functions.

See the docs for more information: https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/#AdvancedParticleFilter-1

Arguments:

• N: Number of particles
• dynamics: A discrete-time dynamics function (x, u, p, t, noise=false) -> x⁺. It's important that the noise argument defaults to false.
• measurement: A measurement function (x, u, p, t, noise=false) -> y. It's important that the noise argument defaults to false.
• measurement_likelihood: A function (x, u, y, p, t)->logl to evaluate the log-likelihood of a measurement.
• dynamics_density: This field is not used by the advanced filter and can be set to nothing.
• initial_density: The distribution of the initial state.
• threads: use threads to propagate particles in parallel. Only activate this if your dynamics is thread-safe. SeeToDee.SimpleColloc is not thread-safe by default due to the use of internal caches, but SeeToDee.Rk4 is.

AuxiliaryParticleFilter(args...; kwargs...)

Takes exactly the same arguments as ParticleFilter, or an instance of ParticleFilter.

DAEUnscentedKalmanFilter(ukf; g, get_x_z, build_xz, xz0, threads=false)

An Unscented Kalman filter for differential-algebraic systems (DAE).

Ref: "Nonlinear State Estimation of Differential Algebraic Systems", Mandela, Rengaswamy, Narasimhan

Arguments

• ukf is a regular UnscentedKalmanFilter that contains dynamics(xz, u, p, t) that propagates the combined state xz(k) to xz(k+1) and a measurement function with signature (xz, u, p, t)
• g(x, z, u, p, t) is a function that should fulfill g(x, z, u, p, t) = 0
• get_x_z(xz) -> x, z is a function that decomposes xz into x and z
• build_xz(x, z) is the inverse of get_x_z
• xz0 the initial full state.
• threads: If true, evaluates dynamics on sigma points in parallel. This typically requires the dynamics to be non-allocating (use StaticArrays) to improve performance.

Assumptions

• The DAE dynamics is index 1 and can be written on the form

\[\begin{aligned} ẋ &= f(x, z, u, p, t) \quad &\text{Differential equations}\ 0 &= g(x, z, u, p, t) \quad &\text{Algebraic equations}\ y &= h(x, z, u, p, t) \quad &\text{Measurements} \begin{aligned}\]

the measurements may be functions of both differential state variables x and algebraic variables z. Please note, the actual dynamics and measurement functions stored in the internal ukf should have signatures (xz, u, p, t), i.e., they take the combined state (descriptor) containing both x and z in a single vector as dictated by the function build_xz. It is only the function g that is assumed to actually have the signature g(x,z,u,p,t).

ExtendedKalmanFilter(kf, dynamics, measurement)
ExtendedKalmanFilter(dynamics, measurement, R1,R2,d0=MvNormal(Matrix(R1)); nu::Int, p = SciMLBase.NullParameters(), α = 1.0, check = true)

A nonlinear state estimator propagating uncertainty using linearization.

The constructor to the extended Kalman filter takes dynamics and measurement functions, and either covariance matrices, or a KalmanFilter. If the former constructor is used, the number of inputs to the system dynamics, nu, must be explicitly provided with a keyword argument.

The filter will internally linearize the dynamics using ForwardDiff.

The dynamics and measurement function are on the following form

x(t+1) = dynamics(x, u, p, t) + w
y      = measurement(x, u, p, t) + e

where w ~ N(0, R1), e ~ N(0, R2) and x(0) ~ d0

See also UnscentedKalmanFilter which is typically more accurate than ExtendedKalmanFilter. See KalmanFilter for detailed instructions on how to set up a Kalman filter kf.

KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1); p = SciMLBase.NullParameters(), α=1)

The matrices A,B,C,D define the dynamics

x' = Ax + Bu + w
y  = Cx + Du + e

where w ~ N(0, R1), e ~ N(0, R2) and x(0) ~ d0

The matrices can be time varying such that, e.g., A[:, :, t] contains the $A$ matrix at time index t. They can also be given as functions on the form

Afun(x, u, p, t) -> A

For maximum performance, provide statically sized matrices from StaticArrays.jl

α is an optional "forgetting factor", if this is set to a value > 1, such as 1.01-1.2, the filter will, in addition to the covariance inflation due to $R_1$, exhibit "exponential forgetting" similar to a Recursive Least-Squares (RLS) estimator. It is thus possible to get a RLS-like algorithm by setting $R_1=0, R_2 = 1/α$ and $α > 1$ ($α$ is the inverse of the traditional RLS parameter $α = 1/λ$). The exact form of the covariance update is

\[R(t+1|t) = α AR(t)A^T + R_1\]

ParticleFilter(N::Integer, dynamics, measurement, dynamics_density, measurement_density, initial_density; threads = false, p = SciMLBase.NullParameters(), kwargs...)

See the docs for more information: https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/#Particle-filter-1

Arguments:

• N: Number of particles
• dynamics: A discrete-time dynamics function (x, u, p, t) -> x⁺
• measurement: A measurement function (x, u, p, t) -> y
• dynamics_density: A probability-density function for additive noise in the dynamics. Use AdvancedParticleFilter for non-additive noise.
• measurement_density: A probability-density function for additive measurement noise. Use AdvancedParticleFilter for non-additive noise.
• initial_density: Distribution of the initial state.

SqKalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1); p = SciMLBase.NullParameters(), α=1)

A standard Kalman filter on square-root form. This filter may have better numerical performance when the covariance matrices are ill-conditioned.

The matrices A,B,C,D define the dynamics

x' = Ax + Bu + w
y  = Cx + Du + e

where w ~ N(0, R1), e ~ N(0, R2) and x(0) ~ d0

The matrices can be time varying such that, e.g., A[:, :, t] contains the $A$ matrix at time index t. They can also be given as functions on the form

Afun(x, u, p, t) -> A

The internal fields storing covariance matrices are for this filter storing the upper-triangular Cholesky factor.

α is an optional "forgetting factor", if this is set to a value > 1, such as 1.01-1.2, the filter will, in addition to the covariance inflation due to $R_1$, exhibit "exponential forgetting" similar to a Recursive Least-Squares (RLS) estimator. It is thus possible to get a RLS-like algorithm by setting $R_1=0, R_2 = 1/α$ and $α > 1$ ($α$ is the inverse of the traditional RLS parameter $α = 1/λ$). The form of the covariance update is

\[R(t+1|t) = α AR(t)A^T + R_1\]

Ref: "A Square-Root Kalman Filter Using Only QR Decompositions", Kevin Tracy https://arxiv.org/abs/2208.06452

TupleProduct(v::NTuple{N,UnivariateDistribution})

Create a product distribution where the individual distributions are stored in a tuple. Supports mixed/hybrid Continuous and Discrete distributions

UnscentedKalmanFilter(dynamics, measurement, R1, R2, d0=MvNormal(Matrix(R1)); p = SciMLBase.NullParameters(), ny, nu)

A nonlinear state estimator propagating uncertainty using the unscented transform.

The dynamics and measurement function are on the following form

x' = dynamics(x, u, p, t) + w
y  = measurement(x, u, p, t) + e

where w ~ N(0, R1), e ~ N(0, R2) and x(0) ~ d0

The matrices R1, R2 can be time varying such that, e.g., R1[:, :, t] contains the $R_1$ matrix at time index t. They can also be given as functions on the form

Rfun(x, u, p, t) -> R

For maximum performance, provide statically sized matrices from StaticArrays.jl

ny, nu indicate the number of outputs and inputs.

commandplot(pf, u, y, p=parameters(pf); kwargs...)

Produce a helpful plot. For customization options (kwargs...), see ?pplot. After each time step, a command from the user is requested.

• q: quit
• s n: step n steps

ll, e = correct!(f, u, y, p = parameters(f), t = index(f))

Update state/covariance/weights based on measurement y, returns loglikelihood and prediction error (the error is always 0 for particle filters).

correct!(kf::SqKalmanFilter, u, y, p = parameters(kf), t::Real = index(kf); R2 = get_mat(kf.R2, kf.x, u, p, t))

For the square-root Kalman filter, a custom provided R2 must be the upper triangular Cholesky factor of the covariance matrix of the measurement noise.

(; ll, e, S, Sᵪ, K) = correct!(kf::AbstractKalmanFilter, u, y, p = parameters(kf), t::Integer = index(kf), R2)

The correct step for a Kalman filter returns not only the log likelihood ll and the prediction error e, but also the covariance of the output S, its Cholesky factor Sᵪ and the Kalman gain K.

If R2 stored in kf is a function R2(x, u, p, t), this function is evaluated at the state before the correction is performed. The measurement noise covariance matrix R2 stored in the filter object can optionally be overridden by passing the argument R2, in this case R2 must be a matrix.

debugplot(pf, u, y, p=parameters(pf); runall=false, kwargs...)

Produce a helpful plot. For customization options (kwargs...), see ?pplot.

• runall=false: if true, runs all time steps befor displaying (faster), if false, displays the plot after each time step.

The generated plot becomes quite heavy. Initially, try limiting your input to 100 time steps to verify that it doesn't crash.

densityplot(x,[w])

Plot (weighted) particles densities

sol = forward_trajectory(pf, u::AbstractVector, y::AbstractVector, p=parameters(pf))

Run the particle filter for a sequence of inputs and measurements. Return a solution with x,w,we,ll = particles, weights, expweights and loglikelihood

If MonteCarloMeasurements.jl is loaded, you may transform the output particles to Matrix{MonteCarloMeasurements.Particles} using Particles(x,we). Internally, the particles are then resampled such that they all have unit weight. This is conventient for making use of the plotting facilities of MonteCarloMeasurements.jl.

sol can be plotted

plot(sol::ParticleFilteringSolution; nbinsy=30, xreal=nothing, dim=nothing)

sol = forward_trajectory(kf::AbstractKalmanFilter, u::Vector, y::Vector, p=parameters(kf))

Run a Kalman filter forward

Returns a KalmanFilteringSolution: with the following

• x: predictions
• xt: filtered estimates
• R: predicted covariance matrices
• Rt: filter covariances
• ll: loglik

sol can be plotted

plot(sol::KalmanFilteringSolution; plotx = true, plotxt=true, plotu=true, ploty=true)

ll(θ) = log_likelihood_fun(filter_from_parameters(θ::Vector)::Function, priors::Vector{Distribution}, u, y, p)

returns function θ -> p(y|θ)p(θ)

ll = loglik(filter, u, y, p=parameters(filter))

Calculate loglikelihood for entire sequences u,y

ll = logsumexp!(w, we [, maxw])

Normalizes the weight vector w and returns the weighted log-likelihood

https://arxiv.org/pdf/1412.8695.pdf eq 3.8 for p(y) https://discourse.julialang.org/t/fast-logsumexp/22827/7?u=baggepinnen for stable logsumexp

x,ll = mean_trajectory(pf, u::Vector{Vector}, y::Vector{Vector}, p=parameters(pf))

This method resets the particle filter to the initial state distribution upon start

mean_trajectory(sol::ParticleFilteringSolution)
mean_trajectory(x::AbstractMatrix, we::AbstractMatrix)

Compute the weighted mean along the trajectory of a particle-filter solution. Returns a matrix of size T × nx. If x and we are supplied, the weights are expected to be in the original space (not log space).

metropolis(ll::Function(θ), R::Int, θ₀::Vector, draw::Function(θ) = naive_sampler(θ₀))

Performs MCMC sampling using the marginal Metropolis (-Hastings) algorithm draw = θ -> θ' samples a new parameter vector given an old parameter vector. The distribution must be symmetric, e.g., a Gaussian. R is the number of iterations. See log_likelihood_fun

Example:

filter_from_parameters(θ) = ParticleFilter(N, dynamics, measurement, MvNormal(n,exp(θ[1])), MvNormal(p,exp(θ[2])), d0)
priors = [Normal(0,0.1),Normal(0,0.1)]
ll     = log_likelihood_fun(filter_from_parameters,priors,u,y,1)
θ₀ = log.([1.,1.]) # Initial point
draw = θ -> θ .+ rand(MvNormal(0.1ones(2))) # Function that proposes new parameters (has to be symmetric)
burnin = 200 # If using threaded call, provide number of burnin iterations
# @time theta, lls = metropolis(ll, 2000, θ₀, draw) # Run single threaded
# thetam = reduce(hcat, theta)'
@time thetalls = LowLevelParticleFilters.metropolis_threaded(burnin, ll, 5000, θ₀, draw) # run on all threads, will provide (2000-burnin)*nthreads() samples
histogram(exp.(thetalls[:,1:2]), layout=3)
plot!(thetalls[:,3], subplot=3) # if threaded call, log likelihoods are in the last column

reset!(kf::AbstractKalmanFilter; x0)

Reset the initial distribution of the state. Optionally, a new mean vector x0 can be provided.

Reset the filter to initial state and covariance/distribution

reset!(kf::SqKalmanFilter; x0)

Reset the initial distribution of the state. Optionally, a new mean vector x0 can be provided.

x,u,y = simulate(f::AbstractFilter, T::Int, du::Distribution, p=parameters(f), [N]; dynamics_noise=true, measurement_noise=true)
x,u,y = simulate(f::AbstractFilter, u, p=parameters(f); dynamics_noise=true, measurement_noise=true)

Simulate dynamical system forward in time T steps, or for the duration of u, returns state sequence, inputs and measurements du is a distribution of random inputs.

A simulation can be considered a draw from the prior distribution over the evolution of the system implied by the selected noise models. Such a simulation is useful in order to evaluate whether or not the noise models are reasonable.

If MonteCarloMeasurements.jl is loaded, the argument N::Int can be supplied, in which case N simulations are done and the result is returned in the form of Vector{MonteCarloMeasurements.Particles}.

xT,RT,ll = smooth(kf::KalmanFilter, u::Vector, y::Vector, p=parameters(kf))
xT,RT,ll = smooth(kf::ExtendedKalmanFilter, u::Vector, y::Vector, p=parameters(kf))

Returns smoothed estimates of state x and covariance R given all input output data u,y

xb,ll = smooth(pf, M, u, y, p=parameters(pf))
xb,ll = smooth(pf, xf, wf, wef, ll, M, u, y, p=parameters(pf))

Perform particle smoothing using forward-filtering, backward simulation. Return smoothed particles and loglikelihood. See also smoothed_trajs, smoothed_mean, smoothed_cov

smoothed_cov(xb)

Helper function to calculate the covariance of smoothed particle trajectories

smoothed_mean(xb)

Helper function to calculate the mean of smoothed particle trajectories

smoothed_trajs(xb)

Helper function to get particle trajectories as a 3-dimensions array (N,M,T) instead of matrix of vectors.

ll, e = update!(f::AbstractFilter, u, y, p = parameters(f), t = index(f))

Perform one step of predict! and correct!, returns loglikelihood and prediction error

weighted_cov(x,we)

Similar to weighted_mean, but returns covariances

x̂ = weighted_mean(x,we)

Calculated weighted mean of particle trajectories. we are expweights.

x̂ = weighted_mean(pf)
x̂ = weighted_mean(s::PFstate)

predict!(f, u, p = parameters(f), t = index(f))

Move filter state forward in time using dynamics equation and input vector u.

predict!(kf::SqKalmanFilter, u, p = parameters(kf), t::Real = index(kf); R1 = get_mat(kf.R1, kf.x, u, p, t))

For the square-root Kalman filter, a custom provided R1 must be the upper triangular Cholesky factor of the covariance matrix of the process noise.

predict!(kf::AbstractKalmanFilter, u, p = parameters(kf), t::Integer = index(kf); R1)

Perform the prediction step (updating the state estimate to $x(t+1|t)$). If R1 stored in kf is a function R1(x, u, p, t), this function is evaluated at the state before the prediciton is performed. The dynamics noise covariance matrix R1 stored in kf can optionally be overridden by passing the argument R1, in this case R1 must be a matrix.

prediction_errors!(res, f::AbstractFilter, u, y, p = parameters(f), λ = 1)

Calculate the prediction errors and store the result in res. Similar to sse, this funciton is useful for sum-of-squares optimization. In contrast to sse, this function returns the residuals themselves rather than their sum of squares. This is useful for Gauss-Newton style optimizers, such as LeastSquaresOptim.LevenbergMarquardt.

Arguments:

• res: A vector of length ny*length(y). Note, for each datapoint in u and u, there are ny outputs, and thus ny residuals.
• f: Any filter
• λ: A weighting factor to minimize dot(e, λ, e). A commonly used metric is λ = Diagonal(1 ./ (mag.^2)), where mag is a vector of the "typical magnitude" of each output. Internally, the square root of W = sqrt(λ) is calculated so that the residuals stored in res are W*e.

See example in Solving using Gauss-Newton optimization.