lampe.inference.mcmc

Monte Carlo Markov chain (MCMC) components.

Classes

MetropolisHastings

Creates a batched Metropolis-Hastings sampler.

Descriptions

class lampe.inference.mcmc.MetropolisHastings(x_0, f=None, log_f=None, sigma=1.0)

Creates a batched Metropolis-Hastings sampler.

Metropolis-Hastings is a Markov chain Monte Carlo (MCMC) sampling algorithm used to sample from intractable distributions \(p(x)\) whose density is proportional to a tractable function \(f(x)\), with \(x \in \mathcal{X}\). The algorithm consists in repeating the following routine for \(t = 1\) to \(T\), where \(x_0\) is the initial sample and \(q(x' | x)\) is a pre-defined transition distribution.

\[\begin{split}1. ~ & x' \sim q(x' | x_{t-1}) \\ 2. ~ & \alpha \gets \frac{f(x')}{f(x_{t-1})} \frac{q(x_{t-1} | x')}{q(x' | x_{t-1})} \\ 3. ~ & u \sim \mathcal{U}(0, 1) \\ 4. ~ & x_t \gets \begin{cases} x' & \text{if } u \leq \alpha \\ x_{t-1} & \text{otherwise} \end{cases}\end{split}\]

Asymptotically, i.e. when \(T \to \infty\), the distribution of samples \(x_t\) is guaranteed to converge towards \(p(x)\). In this implementation, a Gaussian transition \(q(x' | x) = \mathcal{N}(x'; x, \Sigma)\) is used, which can be modified by sub-classing MetropolisHastings.

Wikipedia

https://wikipedia.org/wiki/Metropolis-Hastings_algorithm

Parameters:
  • x_0 (Tensor) – A batch of initial points \(x_0\), with shape \((*, L)\).

  • f (Callable[[Tensor], Tensor]) – A function \(f(x)\) proportional to a density function \(p(x)\).

  • log_f (Callable[[Tensor], Tensor]) – The logarithm \(\log f(x)\) of a function proportional to \(p(x)\).

  • sigma (float | Tensor) – The standard deviation of the Gaussian transition. Either a scalar or a vector.

Example

>>> x_0 = torch.randn(128, 7)
>>> log_f = lambda x: -(x**2).sum(dim=-1) / 2
>>> sampler = MetropolisHastings(x_0, log_f=log_f, sigma=0.5)
>>> samples = [x for x in sampler(256, burn=128, step=4)]
>>> samples = torch.stack(samples)
>>> samples.shape
torch.Size([32, 128, 7])