Mcmc.Proposal.Hamiltonian

Creation date: Mon Jul 5 12:59:42 2021.

The Hamiltonian Monte Carlo (HMC) proposal.

For references, see:

• [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.), Handbook of Markov Chain Monte Carlo (2011), CRC press.
• [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data analysis (2014), CRC Press.
• [3] Review by Betancourt and notes: Betancourt, M., A conceptual introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).

NOTE on implementation:

• The implementation assumes the existence of the Gradient. Like so, the user can use automatic or manual differentiation, depending on the problem at hand.
• The Hamiltonian proposal acts on a vector of storable Values. Functions converting the state to and from this vector have to be provided. See HSettings.
• The desired acceptance rate is 0.65, although the dimension of the proposal is high.
• The speed of this proposal can change drastically when tuned because the leapfrog trajectory length is changed.
• The Hamiltonian proposal is agnostic of the actual prior and likelihood functions, and so, points with zero posterior probability cannot be detected. This affects models with constrained parameters. See Gelman p. 303. This problem can be ameliorated by providing a Validate function so that the proposal can gracefully fail as soon as the state becomes invalid.

Documentation

type Values = Vector Double

The Hamiltonian proposal acts on a vector of floating point values.

type Gradient a = a -> a

Gradient of the log posterior function.

The gradient has to be provided for the complete state. The reason is that the gradient may change if parameters untouched by the Hamiltonian proposal are altered by other proposals.

type Validate a = a -> Bool

Function validating the state.

Useful if parameters are constrained.

Also the validity of the state may depend on parameters untouched by the Hamiltonian proposal.

type Masses = Herm Double

Parameter mass matrix.

The masses roughly describe how reluctant the particles move through the state space. If a parameter has higher mass, the momentum in this direction will be changed less by the provided gradient, than when the same parameter has lower mass. Off-diagonal entries describe the covariance structure. If two parameters are negatively correlated, their generated initial momenta are likely to have opposite signs.

The proposal is more efficient if masses are assigned according to the inverse (co)-variance structure of the posterior function. That is, parameters changing on larger scales should have lower masses than parameters changing on lower scales. In particular, the optimal entries of the diagonal of the mass matrix are the inverted variances of the parameters distributed according to the posterior function.

Of course, the scales of the parameters of the posterior function are usually unknown. Often, it is sufficient to

• set the diagonal entries of the mass matrix to identical values roughly scaled with the inverted estimated average variance of the posterior function; or even to
• set all diagonal entries of the mass matrix to 1.0, and all other entries to 0.0, and trust the tuning algorithm (see HTune) to find the correct values.

type LeapfrogTrajectoryLength = Int

Mean leapfrog trajectory length \(L\).

Number of leapfrog steps per proposal.

To avoid problems with ergodicity, the actual number of leapfrog steps is sampled per proposal from a discrete uniform distribution over the interval \([\text{floor}(0.8L),\text{ceiling}(1.2L)]\).

For a discussion of ergodicity and reasons why randomization is important, see [1] p. 15; also mentioned in [2] p. 304.

Usually set to 10, but larger values may be desirable.

NOTE: To avoid errors, the left bound has an additional hard minimum of 1, and the right bound is required to be larger equal than the left bound.

NOTE: Call error if value is less than 1.

type LeapfrogScalingFactor = Double

Mean of leapfrog scaling factor \(\epsilon\).

Determines the size of each leapfrog step.

To avoid problems with ergodicity, the actual leapfrog scaling factor is sampled per proposal from a continuous uniform distribution over the interval \((0.8\epsilon,1.2\epsilon]\).

For a discussion of ergodicity and reasons why randomization is important, see [1] p. 15; also mentioned in [2] p. 304.

Usually set such that \( L \epsilon = 1.0 \), but smaller values may be required if acceptance rates are low.

NOTE: Call error if value is zero or negative.

data HTuneLeapfrog

Tune leapfrog parameters?

Constructors

HNoTuneLeapfrog

HTuneLeapfrog

We expect that the larger the leapfrog scaling factor the lower the acceptance ratio. Consequently, if the acceptance rate is too low, the leapfrog scaling factor is decreased and vice versa. Further, the leapfrog trajectory length is scaled such that the product of the leapfrog scaling factor and leapfrog trajectory length stays roughly constant.

Instances

Eq HTuneLeapfrog

Defined in Mcmc.Proposal.Hamiltonian

Methods

(==) :: HTuneLeapfrog -> HTuneLeapfrog -> Bool

(/=) :: HTuneLeapfrog -> HTuneLeapfrog -> Bool

Show HTuneLeapfrog

Defined in Mcmc.Proposal.Hamiltonian

Methods

showsPrec :: Int -> HTuneLeapfrog -> ShowS

show :: HTuneLeapfrog -> String

showList :: [HTuneLeapfrog] -> ShowS

data HTuneMasses

Tune masses?

The masses are tuned according to the (co)variances of the parameters obtained from the posterior distribution over the last auto tuning interval.

Constructors

HNoTuneMasses

HTuneDiagonalMassesOnly

Diagonal only: The variances of the parameters are calculated and the masses are amended using the old masses and the inverted variances. If, for a specific coordinate, the sample size is 60 or lower, or if the calculated variance is out of predefined bounds [1e-6, 1e6], the mass of the affected position is not changed.

HTuneAllMasses

All masses: The covariance matrix of the parameters is estimated and the inverted matrix (sometimes called precision matrix) is used as mass matrix. This procedure is error prone, but models with high correlations between parameters it is necessary to tune off-diagonal entries. The full mass matrix is only tuned if more than 200 samples are available. For these reasons, when tuning all masses it is recommended to use tuning settings such as

BurnInWithCustomAutoTuning ([10, 20 .. 200] ++ replicate 5 500)

Instances

Eq HTuneMasses

Defined in Mcmc.Proposal.Hamiltonian

Methods

(==) :: HTuneMasses -> HTuneMasses -> Bool

(/=) :: HTuneMasses -> HTuneMasses -> Bool

Show HTuneMasses

Defined in Mcmc.Proposal.Hamiltonian

Methods

showsPrec :: Int -> HTuneMasses -> ShowS

show :: HTuneMasses -> String

showList :: [HTuneMasses] -> ShowS

data HTune

Tuning settings.

Constructors

HTune HTuneLeapfrog HTuneMasses

Instances

Eq HTune

Defined in Mcmc.Proposal.Hamiltonian

Methods

(==) :: HTune -> HTune -> Bool

(/=) :: HTune -> HTune -> Bool

Show HTune

Defined in Mcmc.Proposal.Hamiltonian

Methods

showsPrec :: Int -> HTune -> ShowS

show :: HTune -> String

showList :: [HTune] -> ShowS

data HSettings a

Specifications of the Hamilton Monte Carlo proposal.

Constructors

HSettings

Fields

• hToVector :: a -> Values

  Extract values to be manipulated by the Hamiltonian proposal from the state.

• hFromVectorWith :: a -> Values -> a

  Put those values back into the state.

• hGradient :: Gradient a
• hMaybeValidate :: Maybe (Validate a)
• hMasses :: Masses
• hLeapfrogTrajectoryLength :: LeapfrogTrajectoryLength
• hLeapfrogScalingFactor :: LeapfrogScalingFactor
• hTune :: HTune

hamiltonian

Arguments

:: Eq a

=> a

The sample state is used for error checks and to calculate the dimension of the proposal.

-> HSettings a

-> PName

-> PWeight

-> Proposal a

Hamiltonian Monte Carlo proposal.