MCMC

Dimension-independent samplers, tracers and diagnostics.

hippymfem.mcmc.chain

Running an MCMC chain.

class hippymfem.mcmc.chain.SampleStruct(kernel)

Bases: object

One chain state: the parameter, the state, and whatever the kernel needs.

assign(other)

class hippymfem.mcmc.chain.MCMC(kernel)

Bases: object

Drive a transition kernel, with burn-in and a tracer.

Parameters are number_of_samples, burn_in, print_progress (how many progress lines to print) and print_level.

run(m0, qoi=None, tracer=None)

Run the chain from m0; returns the number of accepted proposals.

consume_random()

Advance the random stream as a full run would, without running it.

hippymfem.mcmc.kernels

MCMC transition kernels for the posterior over the parameter field.

All of these are dimension-robust: their acceptance rates do not degrade as the mesh is refined, because the proposals are built from the prior rather than from a fixed step in the discrete parameter space.

pCN	preconditioned Crank-Nicolson: propose \(m' = \bar m + s\,\xi + \sqrt{1-s^2}(m-\bar m)\) with \(\xi \sim N(0, R^{-1})\), which leaves the prior invariant, so the acceptance ratio involves only the misfit.
gpCN	the same, but with the Laplace approximation in place of the prior, so the proposal already knows where the posterior mass is; far more efficient when the data is informative.
MALA	Metropolis-adjusted Langevin: a gradient-informed proposal.
IS	importance sampling from the Laplace approximation (independent proposals, weights returned rather than accepted/rejected).

References

• S.L. Cotter, G.O. Roberts, A.M. Stuart, D. White, MCMC methods for functions, Statistical Science 28(3), 2013 (pCN).

• F.J. Pinski, G. Simpson, A.M. Stuart, H. Weber, Algorithms for Kullback-Leibler approximation of probability measures in infinite dimensions, SIAM J. Sci. Comput. 37(6), 2015, Alg. 5.2 (gpCN).

class hippymfem.mcmc.kernels.pCNKernel(model)

Bases: _KernelBase

Preconditioned Crank-Nicolson with the prior as reference measure.

name()

proposal(current)

sample(current, proposed)

One Metropolis-Hastings step; returns 1 if accepted, 0 otherwise.

consume_random()

Advance the stream as one step would, without taking the step.

class hippymfem.mcmc.kernels.gpCNKernel(model, nu)

Bases: _KernelBase

pCN with the Laplace approximation nu as reference measure.

name()

delta(sample)

The log-density difference that the gpCN ratio needs.

\(\Phi(m) + \tfrac12\|m\|_R^2 - \tfrac12\|m-\bar m_\nu\|^2_{H_\nu}\): the true negative log-posterior minus the reference measure’s, so the ratio corrects the proposal’s approximation of the posterior.

proposal(current)

sample(current, proposed)

consume_random()

Advance the stream as one step would, without taking the step.

class hippymfem.mcmc.kernels.MALAKernel(model)

Bases: _KernelBase

Metropolis-adjusted Langevin, preconditioned by the prior covariance.

The proposal is \(m' = m + \tfrac{\tau}{2}\,(-R^{-1}g) + \sqrt{\tau}\,\xi\), with \(g\) the gradient of the negative log-posterior and \(\xi\sim N(0,R^{-1})\); the Metropolis correction makes the chain exact despite the discretization of the Langevin dynamics.

name()

derivativeInfo()

Highest derivative the kernel needs (0 = cost only, 1 = gradient).

init_sample(sample)

Forward and adjoint solves, the cost, and the preconditioned gradient.

proposal(current)

sample(current, proposed)

consume_random()

Advance the stream as one step would, without taking the step.

class hippymfem.mcmc.kernels.ISKernel(model, nu)

Bases: _KernelBase

Importance sampling from the Laplace approximation.

Proposals are independent draws from nu, and each carries the weight \(\exp(-(\delta(m) - \delta(\bar m)))\), with \(\bar m\) the mean of nu; log_weight holds \(-\delta(m)\) itself. The factor \(\exp(\delta(\bar m))\) common to all draws is left out of weight because \(\delta\) is about the misfit at the MAP point, half the number of observations, and \(\exp(-\delta)\) is exactly zero beyond about 745: every weight of a realistic problem would be 0 and a self-normalized estimate 0/0. Self-normalized estimates do not depend on the factor. \(\delta(\bar m)\) costs one forward solve, at the first draw. The chain is not Markov: every “step” is accepted and the weight is what corrects the approximation, so the usual acceptance rate is meaningless and always reported as 1.

delta_reference()

\(\delta\) at the mean of nu, where its quadratic term vanishes.

name()

delta(sample)

proposal(current)

sample(current, proposed)

consume_random()

Advance the stream as one step would, without taking the step.

hippymfem.mcmc.tracers

Recording what a chain visits.

class hippymfem.mcmc.tracers.NullTracer

Bases: object

Records nothing.

append(current, q)

class hippymfem.mcmc.tracers.QoiTracer(n)

Bases: object

Records the quantity of interest at every sample.

append(current, q)

trim()

class hippymfem.mcmc.tracers.FullTracer(n, space=None, basename=None, every=1)

Bases: object

Records the quantity of interest and, optionally, the parameter fields.

Writing every field of a long chain is expensive, so every subsamples.

append(current, q)

trim()

hippymfem.mcmc.diagnostics

Chain diagnostics.

hippymfem.mcmc.diagnostics.integratedAutocorrelationTime(samples, max_lag=None)

Integrated autocorrelation time, by the initial-positive-sequence rule.

Returns (iact, lags, autocorrelation). The sum is truncated at the first negative autocorrelation, which is the standard way to keep the estimator from accumulating noise at large lags.

hippymfem.mcmc.diagnostics.effective_sample_size(samples, max_lag=None)

Number of effectively independent samples in the chain.

hippymfem.mcmc.diagnostics.chain_summary(samples, max_lag=None)

Mean, standard error, IACT and effective sample size.