# hmc

## Purpose

Hybrid Monte Carlo sampling.

## Synopsis

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

## Description

`samples = hmc(f, x, options, gradf)` uses a hybrid Monte Carlo algorithm to sample from the distribution `p ~ exp(-f)`, where `f` is the first argument to `hmc`. The Markov chain starts at the point `x`, and the function `gradf` is the gradient of the `energy' function `f`.

`hmc(f, x, options, gradf, p1, p2, ...)` allows additional arguments to be passed to `f()` and `gradf()`.

`[samples, energies, diagn] = hmc(f, x, options, gradf)` 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, momentum and acceptance threshold) for each step of the chain in `diagn.pos`, `diagn.mom` and `diagn.acc` respectively. All candidate states (including rejected ones) are stored in `diagn.pos`.

`[samples, energies, diagn] = hmc(f, x, options, gradf)` also returns the `energies` (i.e. negative log probabilities) corresponding to the samples. The `diagn` structure contains three fields:

`pos` the position vectors of the dynamic process.

`mom` the momentum vectors of the dynamic process.

`acc` the acceptance thresholds.

`s = hmc('state')` returns a state structure that contains the state of the two random number generators `rand` and `randn` and the momentum of the dynamic process. These are contained in fields `randstate`, `randnstate` and `mom` respectively. The momentum state is only used for a persistent momentum update.

`hmc('state', s)` resets the state to `s`. If `s` is an integer, then it is passed to `rand` and `randn` and the momentum variable is randomised. If `s` is a structure returned by `hmc('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(5)` is set to 1 if momentum persistence is used; default 0, for complete replacement of momentum variables.

`options(7)` defines the trajectory length (i.e. the number of leap-frog steps at each iteration). Minimum value 1.

`options(9)` is set to 1 to check the user defined gradient function.

`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(17)` defines the momentum used when a persistent update of (leap-frog) momentum is used. This is bounded to the interval [0, 1).

`options(18)` is the step size used in leap-frogs; default 1/trajectory length.

## Examples

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

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

## Algorithm

The algroithm follows the procedure outlined in Radford Neal's technical report CRG-TR-93-1 from the University of Toronto. The stochastic update of momenta samples from a zero mean unit covariance gaussian.

## See Also

`metrop`

Pages: Index

Copyright (c) Ian T Nabney (1996-9)