Implements the Hamiltonian Monte Carlo (HMC) algorithm with the No-U-Turn Sampler (NUTS) and Dual Averaging.
More...
|
| std::default_random_engine & | rng_ |
| | Reference to the random number generator.
|
| |
| Eigen::VectorXd | currentPosition_ |
| | The current position (parameter state).
|
| |
| Eigen::VectorXd | currentMomentum_ |
| | The current momentum auxiliary variable.
|
| |
| double | currentScore_ |
| | The current log-posterior score value.
|
| |
| size_t | dim_ |
| | Dimension of the parameter space.
|
| |
| Eigen::VectorXd | mean_ |
| | Running mean of the parameter samples.
|
| |
| Eigen::MatrixXd | cov_ |
| | Running covariance of the parameter samples.
|
| |
| double | epsilon_ |
| | Step size for the leapfrog integrator.
|
| |
| size_t | nSteps_ = 0 |
| | Number of steps taken in the Markov chain.
|
| |
| double | sumAcceptanceProb_ = 0.0 |
| | Sum of acceptance probabilities for step size adaptation.
|
| |
| double | mu_ |
| | Target value for log step size.
|
| |
| double | gamma_ = 0.05 |
| | Adaptation shrinkage parameter.
|
| |
| double | t0_ = 10.0 |
| | Adaptation stabilization parameter.
|
| |
| double | kappa_ = 0.75 |
| | Adaptation exponential decay parameter.
|
| |
| double | H_bar_ = 0.0 |
| | Running average of difference between target and accept probability.
|
| |
| double | log_epsilon_bar_ = 0.0 |
| | Log of the adapted step size.
|
| |
| double | target_accept_ = 0.80 |
| | Target acceptance probability.
|
| |
| std::normal_distribution< double > | distN_ {0.0, 1.0} |
| | Normal generator helper for momentum initialization.
|
| |
| std::uniform_real_distribution< double > | distU_ {0.0, 1.0} |
| | Uniform generator helper for transition accept checks.
|
| |
Implements the Hamiltonian Monte Carlo (HMC) algorithm with the No-U-Turn Sampler (NUTS) and Dual Averaging.
Mathematical Foundations
Hamiltonian Monte Carlo (HMC) maps the target probability density \(p(\mathbf{q})\) to the potential energy of a physical system: \(U(\mathbf{q}) = -\log p(\mathbf{q})\). We introduce auxiliary momentum variables \(\mathbf{p} \sim \mathcal{N}(\mathbf{0}, \mathbf{M})\) with kinetic energy \(K(\mathbf{p}) = \frac{1}{2} \mathbf{p}^T \mathbf{M}^{-1} \mathbf{p}\).
The total energy is given by the Hamiltonian:
\[ H(\mathbf{q}, \mathbf{p}) = U(\mathbf{q}) + K(\mathbf{p}) \]
The joint system evolves according to Hamilton's equations:
\[ \frac{d\mathbf{q}}{dt} = \frac{\partial H}{\partial \mathbf{p}} = \mathbf{M}^{-1} \mathbf{p}, \quad \frac{d\mathbf{p}}{dt} = -\frac{\partial H}{\partial \mathbf{q}} = -\nabla U(\mathbf{q}) \]
Implementation Algorithms
- Leapfrog Integrator: To integrate the system numerically with step size \(\epsilon\), we use the symplectic Leapfrog scheme:
\[ \mathbf{p}\left(t + \frac{\epsilon}{2}\right) = \mathbf{p}(t) - \frac{\epsilon}{2} \nabla U(\mathbf{q}(t)) \]
\[ \mathbf{q}(t + \epsilon) = \mathbf{q}(t) + \epsilon \mathbf{M}^{-1} \mathbf{p}\left(t + \frac{\epsilon}{2}\right) \]
\[ \mathbf{p}(t + \epsilon) = \mathbf{p}\left(t + \frac{\epsilon}{2}\right) - \frac{\epsilon}{2} \nabla U(\mathbf{q}(t + \epsilon)) \]
- No-U-Turn Sampler (NUTS): Avoids manual tuning of the number of integration steps \(L\). It recursively builds a binary tree of leapfrog steps forward and backward in time. The tree stops growing when the trajectory starts to turn back on itself:
\[ (\mathbf{q}_{\text{plus}} - \mathbf{q}_{\text{minus}})^T \mathbf{p}_{\text{minus}} < 0 \quad \text{or} \quad (\mathbf{q}_{\text{plus}} - \mathbf{q}_{\text{minus}})^T \mathbf{p}_{\text{plus}} < 0 \]
- Dual Averaging Step Size Adaptation: Uses Nesterov's dual averaging scheme to adaptively tune the step size \(\epsilon\) to match a target acceptance probability \(\delta\) (default 0.8):
\[ H_t = \delta - \alpha_t, \quad \bar{H}_t = \left(1 - \frac{1}{t + t_0}\right) \bar{H}_{t-1} + \frac{1}{t + t_0} H_t \]
\[ \log \epsilon_t = \mu - \frac{\sqrt{t}}{\gamma} \bar{H}_t, \quad \log \bar{\epsilon}_t = t^{-\kappa} \log \epsilon_t + (1 - t^{-\kappa}) \log \bar{\epsilon}_{t-1} \]
Constraints & Invariants
- Symplectic Flow: The Leapfrog scheme must preserve volume in phase space and hold total energy approximately constant.
- U-turn Condition: Tree growth is terminated immediately if a U-turn is detected to ensure detailed balance is maintained.
1.9.8