GitHub

Particle sampling algorithms

Particle sampling algorithms are a class of Monte Carlo methods that approximate a target distribution through a set of "particles" or points in the parameter space which are weighted or evolved according to some update rule. The stationary distribution of the particles after some number of updates then typically corresponds to the target. In the setting of Bayesian inference, the target distribution is the posterior.

The simplest form of particle method is importance sampling which reweights samples from a so-called proposal distribution. If the proposal is taken to be the prior and the weighting function is simply the likelihood, the resulting weighted samples will approximately follow the posterior distribution. This is implemented in SimulationBasedInference as the "ensemble importance sampling" (EnIS) algorithm:

SimulationBasedInference.EnIS — Type
EnIS

Basic ensemble importance sampling (EnIS) inference algorithm. Also sometimes referred to as the "particle batch smoother" (PBS) or "generalized likelihood uncertainty estimation" (GLUE) depending on the context.

source

where "ensemble" here refers to the collection of particles used for the Monte Carlo approximation.

Naive importance sampling algorithms tend to be very inefficient when applied to high-dimensional or highly nonlinear systems.

For fully Gaussian problems (i.e., Gaussian prior and likelihood), Kalman-type methods are a much more efficient alternative.

SimulationBasedInference currently provides two ensemble Kalman-type algorithms; the first is the Ensmeble Kalman Sampler (EKS; Garbuno-Inigo et al. 2020):

SimulationBasedInference.EKS — Type
EKS <: EnsembleInferenceAlgorithm

Represents a proxy for the Ensemble Kalman Sampler (Garbuno-Inigo et al. 2020) implementation provided by EnsembleKalmanProcesses.

source

and the second is the Ensemble Smoother with Multiple Data Assimilation (ES-MDA; Emerick and Reynolds 2013):

SimulationBasedInference.ESMDA — Type
ESMDA <: EnsembleInferenceAlgorithm

Implementation of the ensemble-smother multiple data assimilation algorithm of Emerick et al. 2013.

Emerick, Alexandre A., and Albert C. Reynolds. "Ensemble smoother with multiple data assimilation." Computers & Geosciences 55 (2013): 3-15.

source

Ensemble interface methods

All EnsembleInferenceAlgorithms implement the following method interface:

SimulationBasedInference.get_ensemble — Method
get_ensemble(state::EnsembleState)

Retrieves the current ensemble matrix from the given EnsembleState. Must be implemented for each ensemble algorithm state type.

source
SimulationBasedInference.isiterative — Method
isiterative(alg::EnsembleInferenceAlgorithm)

Returns true if the given ensemble inference algorithm is iterative, false otherwise. Default implementation returns false (non-iterative).

source
SimulationBasedInference.hasconverged — Method
hasconverged(alg::EnsembleInferenceAlgorithm, state::EnsembleState)

Should return true when alg at the current state has converged, false otherwise. Must be implemented for each ensemble algorithm state type.

source
SimulationBasedInference.initialstate — Method
initialstate(
    alg::EnsembleInferenceAlgorithm,
    ens::AbstractMatrix,
    obs::AbstractVector,
    obscov::AbstractMatrix;
    rng::AbstractRNG=Random.GLOBAL_RNG
)

Constructs the initial ensemble state for the given algorithm and observations.

source
SimulationBasedInference.ensemblestep! — Method
ensemblestep!(solver::EnsembleSolver{algType}) where {algType}

Executes a single ensemble step (forward solve + update) for the given algorithm type. Must be implemented by all ensemble algorithm implementations.

source
SimulationBasedInference.finalize! — Method
finalize!(solver::EnsembleSolver)

Finalizes the solver state after iteration has completed. Default implementation runs ensemble_solve on the current ensemble state and pushes the results to sol.outputs.

source

EnsembleState is an abstract type which must be subtyped by each ensemble algorithm implementation.

EnsembleSolver acts as a container for all of the common components shared by all ensemble algorithms.

SimulationBasedInference.EnsembleSolver — Type
EnsembleSolver{algType,probType,ensalgType,stateType<:EnsembleState,kwargTypes}

Generic implementation of an iterative solver for any ensemble-based algorithm. Uses the SciML EnsembleProblem interface to automatically parallelize forward runs over the ensemble.

source

Utility methods

The following utility methods are also provided

SimulationBasedInference.get_transformed_ensemble — Method
get_transformed_ensemble(sol::SimulatorInferenceSolution{<:EnsembleInferenceAlgorithm}, iter::Int=-1)

Fetches the transformed ensemble from the given solution object. For iterative algorithms, the optinal argument iter may be provided, which then retrieves the ensemble at the given iteration.

source

Internally, the algorithms use ensemble_solve to construct and solve an EnsembleProblem over the parameter ensemble.

SimulationBasedInference.ensemble_solve — Method
ensemble_solve(
    ens::AbstractMatrix,
    initial_prob::SciMLBase.AbstractSciMLProblem,
    ensalg::SciMLBase.BasicEnsembleAlgorithm,
    dealg::Union{Nothing,SciMLBase.AbstractSciMLAlgorithm},
    param_map;
    iter::Integer=1,
    prob_func,
    output_func,
    pred_func,
    solve_kwargs...
)

Performs a single step/iteration for the given ensemble and returns a named tuple (; pred, sol) where sol are the full ensemble forward solutions and pred is the prediction matrix produced by pred_func.

source