MCMC
The samplers here are dimension independent: their acceptance rate does not degrade as the mesh is refined, because the proposals respect the prior.
kernel = hm.pCNKernel(model)
kernel.parameters["s"] = 0.1
chain = hm.MCMC(kernel)
chain.parameters["number_of_samples"] = 5000
chain.parameters["burn_in"] = 500
tracer = hm.QoiTracer(chain.parameters["number_of_samples"])
n_accept = chain.run(m0, qoi=qoi, tracer=tracer)
preconditioned Crank-Nicolson; proposes along the prior |
|
pCN around the Laplace approximation, with much better mixing when it is a good fit |
|
Metropolis-adjusted Langevin, using the gradient |
|
importance sampling from the Laplace approximation |
Reporting a chain honestly
An acceptance rate is not a convergence diagnostic. Report the integrated autocorrelation time and the effective sample size:
iact, lags, acorr = hm.integratedAutocorrelationTime(tracer.data)
ess = len(tracer.data) / iact
A chain with 1500 samples and an IACT of 70 has an effective sample size of about
21, and a mean computed from it is not a result. Tune the step size to an acceptance
rate in the 20-40% range first, then run long enough that the effective sample size
justifies the precision being quoted. integratedAutocorrelationTime()
is checked against an AR(1) process whose IACT is known analytically.
Cost
Every proposal requires a forward solve, and for a nonlinear problem that means reassembling the Jacobian. This is where the element kernels dominate the runtime and where Running on a GPU pays off most directly.