metrop

Purpose

Markov Chain Monte Carlo sampling with Metropolis algorithm.

Synopsis


samples = metrop(f, x, options)
samples = metrop(f, x, options, [], P1, P2, ...)
[samples, energies, diagn] = metrop(f, x, options)
s = metrop('state')
metrop('state', s)

Description

samples = metrop(f, x, options) uses the Metropolis algorithm to sample from the distribution p ~ exp(-f), where f is the first argument to metrop. The Markov chain starts at the point x and each candidate state is picked from a Gaussian proposal distribution and accepted or rejected according to the Metropolis criterion.

samples = metrop(f, x, options, [], p1, p2, ...) allows additional arguments to be passed to f(). The fourth argument is ignored, but is included for compatibility with hmc and the optimisers.

[samples, energies, diagn] = metrop(f, x, options) also returns a log of the energy values (i.e. negative log probabilities) for the samples in energies and diagn, a structure containing diagnostic information (position and acceptance threshold) for each step of the chain in diagn.pos and diagn.acc respectively. All candidate states (including rejected ones) are stored in diagn.pos.

s = metrop('state') returns a state structure that contains the state of the two random number generators rand and randn. These are contained in fields randstate, randnstate.

metrop('state', s) resets the state to s. If s is an integer, then it is passed to rand and randn. If s is a structure returned by metrop('state') then it resets the generator to exactly the same state.

The optional parameters in the options vector have the following interpretations.

options(1) is set to 1 to display the energy values and rejection threshold at each step of the Markov chain. If the value is 2, then the position vectors at each step are also displayed.

options(14) is the number of samples retained from the Markov chain; default 100.

options(15) is the number of samples omitted from the start of the chain; default 0.

options(18) is the variance of the proposal distribution; default 1.

Examples

The following code fragment samples from the posterior distribution of weights for a neural network.


w = mlppak(net);
[samples, energies] = metrop('neterr', w, options, 'netgrad', net, x, t);

Algorithm