pymc.smc.sample_smc#

pymc.smc.sample_smc(draws=2000, kernel=<class 'pymc.smc.kernels.IMH'>, *, start=None, model=None, random_seed=None, chains=None, cores=None, compute_convergence_checks=True, return_inferencedata=True, idata_kwargs=None, progressbar=True, **kernel_kwargs)[source]#

Sequential Monte Carlo based sampling.

Parameters
drawsint, default 2000

The number of samples to draw from the posterior (i.e. last stage). And also the number of independent chains. Defaults to 2000.

kernelSMC Kernel, optional

SMC kernel used. Defaults to pymc.smc.smc.IMH (Independent Metropolis Hastings)

startdict or array of dict, optional

Starting point in parameter space. It should be a list of dict with length chains. When None (default) the starting point is sampled from the prior distribution.

modelModel (optional if in with context).
random_seedint, array_like of int, RandomState or Generator, optional

Random seed(s) used by the sampling steps. If a list, tuple or array of ints is passed, each entry will be used to seed each chain. A ValueError will be raised if the length does not match the number of chains.

chainsint, optional

The number of chains to sample. Running independent chains is important for some convergence statistics. If None (default), then set to either cores or 2, whichever is larger.

coresint, default None

The number of chains to run in parallel. If None, set to the number of CPUs in the system.

compute_convergence_checksbool, default True

Whether to compute sampler statistics like R hat and effective_n. Defaults to True.

return_inferencedatabool, default True

Whether to return the trace as an InferenceData (True) object or a MultiTrace (False). Defaults to True.

idata_kwargsdict, optional

Keyword arguments for pymc.to_inference_data().

progressbarbool, optional, default True

Whether or not to display a progress bar in the command line.

**kernel_kwargsdict, optional

Keyword arguments passed to the SMC_kernel. The default IMH kernel takes the following keywords:

thresholdfloat, default 0.5

Determines the change of beta from stage to stage, i.e. indirectly the number of stages, the higher the value of threshold the higher the number of stages. Defaults to 0.5. It should be between 0 and 1.

correlation_thresholdfloat, default 0.01

The lower the value the higher the number of MCMC steps computed automatically. Defaults to 0.01. It should be between 0 and 1.

Keyword arguments for other kernels should be checked in the respective docstrings.

Notes

SMC works by moving through successive stages. At each stage the inverse temperature \(\beta\) is increased a little bit (starting from 0 up to 1). When \(\beta\) = 0 we have the prior distribution and when \(\beta = 1\) we have the posterior distribution. So in more general terms, we are always computing samples from a tempered posterior that we can write as:

\[p(\theta \mid y)_{\beta} = p(y \mid \theta)^{\beta} p(\theta)\]

A summary of the algorithm is:

  1. Initialize \(\beta\) at zero and stage at zero.

  2. Generate N samples \(S_{\beta}\) from the prior (because when :math beta = 0 the tempered posterior is the prior).

  3. Increase \(\beta\) in order to make the effective sample size equal some predefined value (we use \(Nt\), where \(t\) is 0.5 by default).

  4. Compute a set of N importance weights W. The weights are computed as the ratio of the likelihoods of a sample at stage i+1 and stage i.

  5. Obtain \(S_{w}\) by re-sampling according to W.

  6. Use W to compute the mean and covariance for the proposal distribution, a MvNormal.

  7. Run N independent MCMC chains, starting each one from a different sample in \(S_{w}\). For the IMH kernel, the mean of the proposal distribution is the mean of the previous posterior stage and not the current point in parameter space.

  8. The N chains are run until the autocorrelation with the samples from the previous stage stops decreasing given a certain threshold.

  9. Repeat from step 3 until \(\beta \ge 1\).

  10. The final result is a collection of N samples from the posterior.

References

Minson2013

Minson, S. E., Simons, M., and Beck, J. L. (2013). “Bayesian inversion for finite fault earthquake source models I- Theory and algorithm.” Geophysical Journal International, 2013, 194(3), pp.1701-1726. link

Ching2007

Ching, J., and Chen, Y. (2007). “Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging.” J. Eng. Mech., 2007, 133(7), pp. 816-832. doi:10.1061/(ASCE)0733-9399(2007)133:7(816). link