HamiltonianSampler class

Hamiltonian Monte Carlo (HMC) sampler

Description

A Hamiltonian Monte Carlo (HMC) sampler is a gradient-based Markov Chain Monte Carlo sampler that you can use to generate samples from a probability density P(x). HMC sampling requires specification of log P(x) and its gradient.

The parameter vector x must be unconstrained, meaning that every element of x can be any real number. To sample constrained parameters, transform these parameters into unconstrained variables before using the HMC sampler.

After creating a sampler, you can compute MAP (maximum-a-posteriori) point estimates, tune the sampler, draw samples, and check convergence diagnostics using the methods of this class. For an example of this workflow, see Bayesian Linear Regression Using Hamiltonian Monte Carlo.

Construction

hmc = hmcSampler(logpdf,startpoint) creates a Hamiltonian Monte Carlo (HMC) sampler, returned as a HamiltonianSampler object. logpdf is a function handle that evaluates the logarithm of the probability density of the equilibrium distribution and its gradient. The column vector startpoint is the initial point from which to start HMC sampling.

hmc = hmcSampler(___,Name,Value) specifies additional options using one or more name-value pair arguments. Specify name-value pair arguments after all other input arguments.

Input Arguments

logpdf — Logarithm of target density and its gradient

function handle

Logarithm of target density and its gradient, specified as a function handle.

logpdf must return two output arguments: [lpdf,glpdf] = logpdf(X). Here, lpdf is the base-e log probability density (up to an additive constant), glpdf is the gradient of the log density, and the point X is a column vector with the same number of elements as startpoint.

The input argument X to logpdf must be unconstrained, meaning that every element of X can be any real number. Transform any constrained sampling parameters into unconstrained variables before using the HMC sampler.

If the 'UseNumericalGradient' value is set to true, then logpdf does not need to return the gradient as the second output. Using a numerical gradient can be easier since logpdf does not need to compute the gradient, but it can make sampling slower.

Data Types: function_handle

startpoint — Initial point to start sampling from

numeric column vector

Initial point to start sampling from, specified as a numeric column vector.

Data Types: single | double

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'VariableNames',{'Intercept','Beta'},'MassVectorTuningMethod','hessian' specifies sampling variable names and the mass vector tuning method to be 'hessian'.

'StepSize' — Step size of Hamiltonian dynamics

0.1 (default) | positive scalar

Step size of Hamiltonian dynamics, specified as the comma-separated pair consisting of 'StepSize' and a positive scalar.

To propose a new state for the Markov chain, the HMC sampler integrates the Hamiltonian dynamics using leapfrog integration. This argument controls the step size of that leapfrog integration.

You can automatically tune the step size using tuneSampler.

Example: 'StepSize',0.2

'NumSteps' — Number of steps of Hamiltonian dynamics

50 (default) | positive integer

Number of steps of Hamiltonian dynamics, specified as the comma-separated pair consisting of 'NumSteps' and a positive integer.

To propose a new state for the Markov chain, the HMC sampler integrates the Hamiltonian dynamics using leapfrog integration. This argument controls the number of steps of that leapfrog integration.

You can automatically tune the number of steps using tuneSampler.

Example: 'NumSteps',20

'MassVector' — Mass vector of momentum variables

ones(size(startpoint,1),1) (default) | numeric column vector

Mass vector of momentum variables, specified as the comma-separated pair consisting of 'MassVector' and a numeric column vector with positive values and the same length as startpoint.

The “masses” of the momentum variables associated with the variables of interest control the Hamiltonian dynamics in each Markov chain proposal.

You can automatically tune the mass vector using tuneSampler.

Example: 'MassVector',rand(3,1)

'JitterMethod' — Method for jittering step size and number of steps

'jitter-both' (default) | 'jitter-numsteps' | 'none'

Method for jittering the step size and number of steps, specified as the comma-separated pair consisting of 'JitterMethod' and one of the following:

Value	Description
'jitter-both'	Randomly jitter the step size and number of steps for each leapfrog trajectory.
'jitter-numsteps'	Jitter only the number of steps of each leapfrog trajectory.
'none'	Perform no jittering.

With jittering, the sampler randomly selects the step size or the number of steps of each leapfrog trajectory as values smaller than the 'StepSize' and 'NumSteps' values. Use jittering to improve the stability of the leapfrog integration of the Hamiltonian dynamics.

Example: 'JitterMethod','jitter-both'

'StepSizeTuningMethod' — Method for tuning sampler step size

'dual-averaging' (default) | 'none'

Method for tuning the sampler step size, specified as the comma-separated pair consisting of 'StepSizeTuningMethod' and 'dual-averaging' or 'none'.

If the 'StepSizeTuningMethod' value is set to 'dual-averaging', then tuneSampler tunes the leapfrog step size of the HMC sampler to achieve a certain acceptance ratio for a fixed value of the simulation length. The simulation length equals the step size multiplied by the number of steps. To set the target acceptance ratio, use the 'TargetAcceptanceRatio' name-value pair argument of the tuneSampler method.

Example: 'StepSizeTuningMethod','none'

'MassVectorTuningMethod' — Method for tuning sampler mass vector

'iterative-sampling' (default) | 'hessian' | 'none'

Method for tuning the sampler mass vector, specified as the comma-separated pair consisting of 'MassVectorTuningMethod' and one of the following values

Value	Description
'iterative-sampling'	Tune the MassVector via successive approximations by drawing samples using a sequence of mass vector estimates.
'hessian'	Set the MassVector equal to the negative diagonal Hessian of the logpdf at the startpoint.
'none'	Perform no tuning of the MassVector.

To perform the tuning, use the tuneSampler method.

Example: 'MassVectorTuningMethod','hessian'

'CheckGradient' — Flag for checking analytical gradient

true (or 1) (default) | false (or 0)

Flag for checking the analytical gradient, specified as the comma-separated pair consisting of 'CheckGradient' and either true (or 1) or false (or 0).

If 'CheckGradient' is true, then the sampler calculates the numerical gradient at the startpoint and compares it to the analytical gradient returned by logpdf.

Example: 'CheckGradient',true

'VariableNames' — Sampling variable names

{'x1','x2',...} (default) | string array | cell array of character vectors

Sampling variable names, specified as the comma-separated pair consisting of 'VariableNames' and a string array or cell array of character vectors. Elements of the array must be unique. The length of the array must be the same as the length of startpoint.

Supply a 'VariableNames' value to label the components of the vector you want to sample using the HMC sampler.

Example: 'VariableNames',{'Intercept','Beta'}

'UseNumericalGradient' — Flag for using numerical gradient

false (or 0) (default) | true (or 1)

Flag for using numerical gradient, specified as the comma-separated pair consisting of 'UseNumericalGradient' and either true (or 1) or false (or 0).

If you set the 'UseNumericalGradient' value to true, then the HMC sampler numerically estimates the gradient from the log density returned by logpdf. In this case, the logpdf function does not need to return the gradient of the log density as the second output. Using a numerical gradient makes HMC sampling slower.

Example: 'UseNumericalGradient',true

Properties

StepSize — Step size of Hamiltonian dynamics

0.1 (default) | positive scalar

Step size of Hamiltonian dynamics, specified as a positive scalar.

To propose a new state for the Markov chain, the HMC sampler integrates the Hamiltonian dynamics using leapfrog integration. The value of this property controls the step size of that leapfrog integration.

NumSteps — Number of steps of Hamiltonian dynamics

50 (default) | positive integer

Number of steps of Hamiltonian dynamics, specified as a positive integer.

To propose a new state for the Markov chain, the HMC sampler integrates the Hamiltonian dynamics using leapfrog integration. The value of this property controls the number of steps of that leapfrog integration.

MassVector — Mass vector of momentum variables

ones(size(startpoint,1),1) (default) | numeric column vector

Mass vector of momentum variables, specified as a numeric column vector with positive values and the same length as startpoint.

The “masses” of the momentum variables associated with the variables of interest control the Hamiltonian dynamics in each Markov chain proposal.

JitterMethod — Method for jittering step size and number of steps

'jitter-both' (default) | 'jitter-numsteps' | 'none'

Method for jittering the step size and the number of steps, specified as one of the following values.

Value	Description
'jitter-both'	Randomly jitter the step size and number of steps of each leapfrog trajectory.
'jitter-numsteps'	Jitter only the number of steps of each leapfrog trajectory.
'none'	Perform no jittering.

With jittering, the sampler randomly selects the step size or the number of steps of each leapfrog trajectory as values smaller than the 'StepSize' and 'NumSteps' values. Use jittering to improve the stability of the leapfrog integration of the Hamiltonian dynamics.

StepSizeTuningMethod — Method for tuning sampler step size

'dual-averaging' (default) | 'none'

Method for tuning the sampler step size, specified as 'dual-averaging' or 'none'.

If StepSizeTuningMethod equals 'dual-averaging', then tuneSampler tunes the leapfrog step size of the HMC sampler to achieve a certain acceptance ratio for a fixed value of the simulation length. The simulation length equals the step size multiplied by the number of steps. To set the target acceptance ratio, use the 'TargetAcceptanceRatio' name-value pair argument of the tuneSampler method.

MassVectorTuningMethod — Method for tuning sampler mass vector

'iterative-sampling' (default) | 'hessian' | 'none'

Method for tuning the sampler mass vector, specified as one of the following values.

Value	Description
'iterative-sampling'	Tune the MassVector via successive approximations by drawing samples using a sequence of mass vector estimates.
'hessian'	Set the MassVector equal to the negative diagonal Hessian of the logpdf at the startpoint.
'none'	Perform no tuning of the MassVector.

To perform the tuning, use the tuneSampler method.

LogPDF — Logarithm of target density and its gradient

function handle

Logarithm of target density and its gradient, specified as a function handle.

LogPDF returns two output arguments: [lpdf,glpdf] = LogPDF(X). Here, lpdf is the base-e log probability density (up to an additive constant) and glpdf is the gradient of the log density at the point X. The input argument X must be a column vector with the same number of elements as the StartPoint property.

If you set the 'UseNumericalGradient' value to true when creating the sampler, then LogPDF returns the numerical gradient in glpdf.

StartPoint — Initial point to start sampling from

numeric column vector

Initial point to start sampling from, specified as a numeric column vector.

VariableNames — Sampling variable names

{'x1','x2',...} (default) | cell array of unique character vectors

Sampling variable names, specified as a cell array of unique character vectors.

Methods

diagnostics	Markov Chain Monte Carlo diagnostics
drawSamples	Generate Markov chain using Hamiltonian Monte Carlo (HMC)
estimateMAP	Estimate maximum of log probability density
tuneSampler	Tune Hamiltonian Monte Carlo (HMC) sampler

Examples

Create Hamiltonian Monte Carlo Sampler

Create a Hamiltonian Monte Carlo (HMC) sampler to sample from a normal distribution.

First, save a function normalDistGrad on the MATLAB® path that returns the multivariate normal log probability density and its gradient (normalDistGrad is defined at the end of this example). Then, call the function with arguments to define the logpdf input argument to the hmcSampler function.

means = [1;-3];
standevs = [1;2];
logpdf = @(theta)normalDistGrad(theta,means,standevs);

Choose a starting point for the HMC sampler.

startpoint = randn(2,1);

Create the HMC sampler and display its properties.

smp = hmcSampler(logpdf,startpoint);

smp

smp = 
  HamiltonianSampler with properties:

                  StepSize: 0.1000
                  NumSteps: 50
                MassVector: [2x1 double]
              JitterMethod: 'jitter-both'
      StepSizeTuningMethod: 'dual-averaging'
    MassVectorTuningMethod: 'iterative-sampling'
                    LogPDF: @(theta)normalDistGrad(theta,means,standevs)
             VariableNames: {2x1 cell}
                StartPoint: [2x1 double]

The normalDistGrad function returns the logarithm of the multivariate normal probability density with means in Mu and standard deviations in Sigma, specified as scalars or columns vectors the same length as startpoint. The second output argument is the corresponding gradient.

function [lpdf,glpdf] = normalDistGrad(X,Mu,Sigma)
Z = (X - Mu)./Sigma;
lpdf = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2));
glpdf = -Z./Sigma;
end