GitHub

MarkovChainMonteCarlo

Top-level class and methods

CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCWrapper — Type
struct MCMCWrapper{VV1<:(AbstractVector), VV2<:(AbstractVector), VV3<:(AbstractVector)}

Top-level class holding all configuration information needed for MCMC sampling: the prior, emulated likelihood and sampling algorithm (DensityModel and Sampler, respectively, in AbstractMCMC's terminology).

Fields

  • prior::EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution: ParameterDistribution object describing the prior distribution on parameter values.

  • encoded_prior::EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution: ParameterDistribution object describing the encoded prior distribution on the encoded parameter values.

  • observations::AbstractVector: [outputdim x Nsamples] matrix, of given observation data.

  • encoded_observations::AbstractVector: Vector of observations describing the data samples to actually used during MCMC sampling (that have been transformed into a space consistent with emulator outputs).

  • log_posterior_map::AbstractMCMC.AbstractModel: AdvancedMH.DensityModel object, used to evaluate the posterior density being sampled from.

  • mh_proposal_sampler::AbstractMCMC.AbstractSampler: Object describing a MCMC sampling algorithm and its settings.

  • sample_kwargs::NamedTuple: NamedTuple of other arguments to be passed to AbstractMCMC.sample().

  • encoder_schedule::AbstractVector: Vector of encoders dictating how to encode/decode data.

source
CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCWrapper — Method
MCMCWrapper(
    mcmc_alg::CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCProtocol,
    observation::Union{AbstractMatrix, AbstractVector},
    prior::EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution,
    em_or_fmw::Union{CalibrateEmulateSample.Emulators.Emulator, CalibrateEmulateSample.Emulators.ForwardMapWrapper};
    init_params,
    burnin,
    kwargs...
) -> CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCWrapper

Constructor for MCMCWrapper that will perform an MCMC sampling in the encoded space given by em_or_fmw (Emulator or ForwardMapWrapper). It creates and wraps an instance of EmulatorPosteriorModel, for sampling from the Emulator (predicting means and covariances), and MetropolisHastingsSampler, for generating the MC proposals.

  • mcmc_alg: MCMCProtocol describing the MCMC sampling algorithm to use. Currently implemented algorithms are:

    • RWMHSampling: Metropolis-Hastings sampling from a vanilla random walk with fixed stepsize.
    • pCNMHSampling: Metropolis-Hastings sampling using the preconditioned Crank-Nicholson algorithm, which has a well-behaved small-stepsize limit.
    • BarkerSampling: Metropolis-Hastings sampling using the Barker proposal, which has a robustness to choosing step-size parameters.

  • observation: Vector (for one sample) or matrix with columns as samples from the observation. Can, e.g., be picked from an Observation struct using get_obs_sample.

  • prior: ParameterDistribution object containing the parameters' prior distributions.

  • em_or_fmw: Emulator or ForwardMapWrapper to sample from.

  • init_params: Starting parameter values for MCMC sampling. (defined in unconstrained parameter coordinates)

  • burnin: Initial number of MCMC steps to discard from output (pre-convergence).

source
StatsBase.sample — Function

sample([rng,] mcmc::MCMCWrapper, args...; kwargs...)

Extends the sample methods of AbstractMCMC (which extends StatsBase) to sample from the emulated posterior, using the MCMC sampling algorithm and Emulator configured in MCMCWrapper. Returns a MCMCChains.Chains object containing the samples.

Supported methods are:

  • sample([rng, ]mcmc, N; kwargs...)

    Return a MCMCChains.Chains object containing N samples from the emulated posterior.

  • sample([rng, ]mcmc, isdone; kwargs...)

    Sample from the model with the Markov chain Monte Carlo sampler until a convergence criterion isdone returns true, and return the samples. The function isdone has the signature

        isdone(rng, model, sampler, samples, state, iteration; kwargs...)

    where state and iteration are the current state and iteration of the sampler, respectively. It should return true when sampling should end, and false otherwise.

  • sample([rng, ]mcmc, parallel_type, N, nchains; kwargs...)

    Sample nchains Monte Carlo Markov chains in parallel according to parallel_type, which may be MCMCThreads() or MCMCDistributed() for thread and parallel sampling, respectively.

source
CalibrateEmulateSample.MarkovChainMonteCarlo.get_posterior — Function
get_posterior(
    mcmc::CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCWrapper,
    chain::MCMCChains.Chains;
    noise_injector_threshold,
    noise_injector_scaling
) -> EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution

Returns a ParameterDistribution object corresponding to the empirical distribution of the samples in chain.

keyword args

  • noise_injector_threshold[=0.001]: If the encoded space is lossy, and the lost variability due to encoding exceeds a threshold noise_injector_threshold, then in place of decoding posterior samples, additional noise consistent with the prior is injected into the null space of the encoder. See decode_and_add_noise() for more detail.
  • noise_injector_scaling[=1.0]: Scales the injected noise; though 1.0 is the only "consistent" value, reduction may be necessary if noise injection causes posterior samples to be unstable in simulations.
Note

This method does not currently support combining samples from multiple Chains.

source
CalibrateEmulateSample.MarkovChainMonteCarlo.optimize_stepsize — Function
optimize_stepsize(
    rng::Random.AbstractRNG,
    mcmc::CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCWrapper;
    init_stepsize,
    N,
    max_iter,
    target_acc,
    sample_kwargs...
) -> Float64

Uses a heuristic to return a stepsize for the mh_proposal_sampler element of MCMCWrapper which yields fast convergence of the Markov chain.

The criterion used is that Metropolis-Hastings proposals should be accepted between 15% and 35% of the time.

source

See AbstractMCMC sampling API for background on our use of Turing.jl's AbstractMCMC API for MCMC sampling.

Sampler algorithms

CalibrateEmulateSample.MarkovChainMonteCarlo.RWMHSampling — Type
struct RWMHSampling{T<:CalibrateEmulateSample.MarkovChainMonteCarlo.AutodiffProtocol} <: CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCProtocol

MCMCProtocol which uses Metropolis-Hastings sampling that generates proposals for new parameters as as vanilla random walk, based on the covariance of prior.

source
CalibrateEmulateSample.MarkovChainMonteCarlo.pCNMHSampling — Type
struct pCNMHSampling{T<:CalibrateEmulateSample.MarkovChainMonteCarlo.AutodiffProtocol} <: CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCProtocol

MCMCProtocol which uses Metropolis-Hastings sampling that generates proposals for new parameters according to the preconditioned Crank-Nicholson (pCN) algorithm, which is usable for MCMC in the stepsize → 0 limit, unlike the vanilla random walk. Steps are based on the covariance of prior.

source
CalibrateEmulateSample.MarkovChainMonteCarlo.MetropolisHastingsSampler — Function
MetropolisHastingsSampler(
    _::CalibrateEmulateSample.MarkovChainMonteCarlo.RWMHSampling{T<:CalibrateEmulateSample.MarkovChainMonteCarlo.AutodiffProtocol},
    encoded_prior::EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution,
    encoder_schedule::AbstractVector
) -> CalibrateEmulateSample.MarkovChainMonteCarlo.RWMetropolisHastings{AdvancedMH.RandomWalkProposal{false, _A}} where _A

Constructor for all Sampler objects, with one method for each supported MCMC algorithm.

Warning

Both currently implemented Samplers use a Gaussian approximation to the prior: specifically, new Metropolis-Hastings proposals are a scaled combination of the old parameter values and a random shift distributed as a Gaussian with the same covariance as the prior.

This suffices for our current use case, in which our priors are Gaussian (possibly in a transformed space) and we assume enough fidelity in the Emulator that inference isn't prior-dominated.

source

Emulated posterior (Model)

CalibrateEmulateSample.MarkovChainMonteCarlo.EmulatorPosteriorModel — Function
EmulatorPosteriorModel(
    encoded_prior::EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution,
    em_or_fmw::Union{CalibrateEmulateSample.Emulators.Emulator, CalibrateEmulateSample.Emulators.ForwardMapWrapper},
    obs_vec::AbstractVector
) -> AdvancedMH.DensityModel{F} where F<:(CalibrateEmulateSample.MarkovChainMonteCarlo.var"#13#14"{EnsembleKalmanProcesses.ParameterDistributions.ParameterDistribution{PDType, CType, ST}, _A, <:AbstractVector{T}} where {PDType<:EnsembleKalmanProcesses.ParameterDistributions.ParameterDistributionType, CType<:EnsembleKalmanProcesses.ParameterDistributions.ConstraintType, ST<:AbstractString, _A, T})

Factory which constructs AdvancedMH.DensityModel objects given an (Encoded) prior on the model parameters (encoded_prior) and an Emulator of the log-likelihood of the data given parameters. Together these yield the log posterior density we're attempting to sample from with the MCMC, which is the role of the DensityModel class in the AbstractMCMC interface.

source

Internals - MCMC State

CalibrateEmulateSample.MarkovChainMonteCarlo.MCMCState — Type
struct MCMCState{T, L<:Real} <: AdvancedMH.AbstractTransition

Extends the AdvancedMH.Transition (which encodes the current state of the MC during sampling) with a boolean flag to record whether this state is new (arising from accepting a Metropolis-Hastings proposal) or old (from rejecting a proposal).

Fields

  • params::Any: Sampled value of the parameters at the current state of the MCMC chain.

  • log_density::Real: Log probability of params, as computed by the model using the prior.

  • accepted::Bool: Whether this state resulted from accepting a new MH proposal.

source