Metropolis Hastings Sampling

Metropolis Hastings sampling is a widely used Markov Chain Monte Carlo (MCMC) algorithm for generating a sequence of samples from a target probability distribution where direct sampling is difficult. The samples can be used to estimate properties of the target distribution, like its mean, variance, and higher moments. The samples can also be used to approximate the target distribution, for example via a histogram. The algorithm works by constructing a Markov chain with a stationary distribution equal to the target distribution.

Steps of the algorithm are as follows:

  1. Choose an initial state for the Markov chain, usually from the target distribution.

  2. At each iteration, propose a new state by sampling from a proposal distribution.

  3. Calculate the acceptance probability.

  4. Accept or reject the proposed state based on the acceptance probability.

  5. Repeat steps 2-4 for a large number of iterations to obtain a sequence of samples.

  6. Apply “burn-in” and “thining”.

Final samples from the Markov chain can then be used to estimate the target distribution or its properties.

MetropolisHastings Class

The MetropolisHastings class is imported by:

from psimpy.sampler.metropolis_hastings import MetropolisHastings

Methods

class MetropolisHastings(ndim, init_state, f_sample, target=None, ln_target=None, bounds=None, f_density=None, symmetric=True, nburn=0, nthin=1, seed=None, args_target=None, kwgs_target=None, args_f_sample=None, kwgs_f_sample=None, args_f_density=None, kwgs_f_density=None)[source]

Metropolis Hastings sampling.

Parameters
  • ndim (int) – Dimension of parameter x.

  • init_state (numpy array) – Contains ndim numeric values representing the starting point of the Metropolis hastings chain. Shape (ndim,).

  • f_sample (Callable) – A function which proposes a new state x' given a current state x. Call with f_sample(x, *args_f_sample, **kwgs_f_sample). Return the proposed state x' as a numpy.ndarray of shape (ndim,). If \(ndim=1\), f_sample may also return a scalar value.

  • target (Callable, optional) – Target probability density function, up to a normalizing constant. Call with target(x, *args_target, **kwgs_target). Return the probability density value at x. If target is not given, ln_target must be provided.

  • ln_target (Callable, optional) – Natural logarithm of target probability density function. Call with ln_target(x, *args_target, **kwgs_target). Return the value of natural logarithm of target probability density value at x. If ln_target is not given, target must be provided.

  • bounds (numpy array, optional) – Upper and lower boundaries of each parameter. Shape (ndim, 2). bounds[:, 0] corresponds to lower boundaries of each parameter and bounds[:, 1] to upper boundaries of each parameter.

  • f_density (Callable, optional) – Call with f_density(x1, x2, *args_f_density, **kwgs_f_density). Return the probability density at state x1 given state x2. It must be provided if symmetric is set to False.

  • symmetric (bool, optional) – Whether f_density is symmetric or asymmetric.

  • nburn (int, optional) – Number of burnin. Default no burnin.

  • nthin (int, optional) – Number of thining to reduce dependency. Default no thining.

  • seed (int, optional) – Seed to initialize the pseudo-random number generator. Default None.

  • args_target (list, optional) – Positional arguments for target or ln_target.

  • kwgs_target (dict, optional) – Keyword arguments for target or ln_target.

  • args_f_sample (list, optional) – Positional arguments for f_sample.

  • kwgs_f_sample (dict, optional) – Keyword arguments for f_sample.

  • args_f_density (list, optional) – Positional arguments for f_density.

  • kwgs_f_density (dict, optional) – Keyword arguments for f_density.

sample(nsamples)[source]

Draw samples.

Parameters

nsamples (int) – Number of samples to be drawn.

Return type

tuple[ndarray, ndarray]

Returns

  • mh_samples (numpy array) – Samples drawn from target. Shape (nsamples, ndim).

  • mh_accept (numpy array) – Shape of (nsamples,). Each element indicates whether the corresponding sample is the proposed new state (value 1) or the old state (value 0). np.mean(mh_accept) thus gives the overall acceptance ratio.