Hamiltonian Monte Carlo法

記事の内容


前に一度とりあげましたが, HMCをもう一度. 説明も割と被ってます. 一応動機を述べておくと, Bayes深層学習や確率的勾配Langevin動力学法のための準備のつもりです.

基本的なアイデア

概要

Hamiltonian Monte Carlo法は, 主にBayes推論で用いられるサンプリング近似法MCMCの一種です. Hybrid Monte Carloとも呼ばれ, HMCと省略されます. MCMCの中でも特に効率的な方法です. HMCの元のアイデアは物理から借りてきたものです. 直感的にはわかりやすい反面, 数理的に解析するには幾何的なアプローチをとるため, 難解です. HMCの拡張が提案される際には, Riemann幾何が登場します. また, Hamilton系の文脈ではシンプレクティック幾何が絡んできます. さらに, GelmanとHoffmanが提案したNUTSも正直訳がわかりません. 今回はあくまで今後の準備のための記事なので, 難しい話は無しです. そもそも説明できませんし...

本当はうまくいかない例を記事にしたかったのですが, なかなか大変だったのでまた機会があれば紹介します. 特に階層モデルに対するHMCは全然うまくいきません. 今, 試行錯誤しているところです. とりあえず今回は綺麗で簡単な例だけです.

物理的イメージ

HMCの基本的なアイデアを述べます. そもそもMCMCでやりたいことは, ある分布(以下, 目的の分布)からのサンプルを発生させることです. そのサンプルを使って期待値などを近似します. 目的の分布からのサンプルを得るためには, 分布の中でも確率の大きい部分のサンプルが十分得られれば, それなりの精度で近似できると期待できます. 典型的には, 密度関数は次のような形状をしています(現実にはもっと複雑な形状かもしれません).

この密度関数の上の部分ほど確率が大きいです. したがって, 上の部分周辺をうろちょろしたい訳です. 突然ですが, ここでこの密度関数を上下ひっくり返してみます.

ひっくり返したので, 下の窪んだ部分からのサンプルが欲しいですね. したがって, 次のように玉を転がすイメージでサンプルできれば理想的です. 窪んだ部分の周辺を満遍なく転がせば, 分布は再現できると期待します. このアイデアをもとに, HMCのアルゴリズムを構成していきます.

アルゴリズムの大まかな流れを説明します. 目的の分布が次のように書けるとします.

\begin{equation} p(\boldsymbol{x}) = \frac{1}{Z}\tilde{p}(\boldsymbol{x}) \propto \mathrm{exp}\left( -U(\boldsymbol{x} )\right) \end{equation}

\(U\)は上の図のように窪んだ形をしているとします. この\(U\)を, ポテンシャルエネルギーと呼びます. 今, 前回のサンプルが手元にあるとします. ここから出発して, 先ほどの玉転がしを物理的にシミュレートします. ここでいう物理的とは, 微分方程式でシミュレートすることを指します. 当然このシミュレートには微分方程式の数値解法を用います. 適当な反復回数と適当な刻み幅で軌道を計算します. 軌道を計算して行き着いた先が, 新しいサンプルの候補です. この候補を棄却するか, 採択するかはMetropolis-Hastings法と同様に吟味します. 要するに, ポテンシャルエネルギーが小さくなる方へ下っていき, エネルギーが小さくなれば採択します.

アルゴリズムの構成

この節では, 玉転がしのシミュレート方法を具体的に考えます.

Hamilton系の導入

先ほどと同じ記号を用います. 関数\(U\)は, 目的の分布に関する情報を集約したものです. ここで, 新しい変数\(\boldsymbol{p}\)を導入して, 次のような関数を定めます. もちろん欲しいサンプルは\(\boldsymbol{x}\)だけです. \(\boldsymbol{p}\)は便宜上導入して, あとで無視すれば大丈夫です.

\begin{equation} H(\boldsymbol{x}, \boldsymbol{p}) = U(\boldsymbol{x}) + \frac{1}{2}\|\boldsymbol{p}\|^2 \end{equation}

これは, 物理ではHamiltonianとして知られています.

【(可分な)Hamilton系】

次のHamiltonian \begin{equation} H(\boldsymbol{x},\boldsymbol{p}) = U(\boldsymbol{x}) + \frac{1}{2}\| \boldsymbol{p}\|^2 \end{equation} に対応するHamilton系は, 次式で表される.

\begin{align} \frac{d\boldsymbol{x}}{d\tau} &= \boldsymbol{p}\\ \frac{d\boldsymbol{p}}{d\tau} &= -\nabla U(\boldsymbol{x}) \end{align}

ということで, 上で定義した関数\(H\)をHamiltonianとするHamilton系を反復ごとにシミュレートすることにします. サンプルしたい\(\boldsymbol{x}\)は位置に相当し, \(\boldsymbol{p}\)は運動量に相当します. Hamiltonianは系の全エネルギーを表し, 時間に依存せず, 少なくともHMCの文脈では 保存されます. したがって, Hamilton系をシミュレートしても到達できる位置は限定されます. そこで, 反復ごとに運動量\(\boldsymbol{p}\)を変えて確率的な揺らぎを加えます. 普通は, 運動量\(\boldsymbol{p}\)は標準正規分布からとってきます. 要するに, 同じような転がり方だけでは満遍なくサンプルできないので, 時々玉を別の方向に弾いてやるってことです. あとは, Hamiltonianが定める微分方程式を如何に解くか, という問題が残っています. 前回のサンプルと, 正規分布からサンプルした運動量を初期値として, 以下の微分方程式を解きます.

\begin{equation} H(\boldsymbol{x},\boldsymbol{p}) = U(\boldsymbol{x}) + \frac{1}{2}\| \boldsymbol{p}\|^2 \end{equation}

これを解くのに, leap-frog法を用います.

leap-frog法による離散化

微分方程式を数値的に解くためには, 離散化する必要があります. 離散化の仕方によって数値解法が決まります. HMCで用いるleap-frog法は, 次のように離散化します.

【leap-frog法(Störmer-Verlet法)の1ステップ】

\begin{align} \boldsymbol{p}_{n+\frac{1}{2}} &= \boldsymbol{p}_n - \frac{h}{2}\nabla U(\boldsymbol{x}_n) \\ \boldsymbol{x}_{n+1} &= \boldsymbol{x}_n + h \boldsymbol{p}_{n+\frac{1}{2}} \\ \boldsymbol{p}_{n+1} &= \boldsymbol{p}_{n+\frac{1}{2}} - \frac{h}{2} \nabla U(\boldsymbol{x}_{n+1}) \end{align}

leap-frog法は, 半ステップごとに更新します. Störmer-Verlet法とも言います. このスキームは, 天文学の分野でStörmerが, 分子動力学の分野でVerletがそれぞれ独立に発見しました. しかし, どうやらNewtonが遥か前に発見していたようです... さらに, Feynmanも発見していたとか.

微分方程式の数値解法

leap-frog法に関して補足しておきます. なぜleap-frog法が必要なのでしょうか. まずは数値解法を整理しておきましょう. 次の1階微分方程式を考えます.

\begin{equation} \frac{d}{dt}\boldsymbol{x} = f(\boldsymbol{x}),\quad \boldsymbol{x}:\mathbb{R}\to \mathbb{R}^d,\quad f:\mathbb{R}^d\to \mathbb{R}^d \end{equation}

このような微分方程式の数値解法として有名なものを2つほど導入します. 以下, \(h\)を離散化の幅, \(t_n=nh\)とし, \(\boldsymbol{x}_n\)を\(\boldsymbol{x}(t_n)\)の近似値とします. 1つ目は, 陽的Euler法です.

【スキーム1: 陽的Euler法】

\begin{equation} \boldsymbol{x}_{n+1} = \boldsymbol{x}_n + hf(\boldsymbol{x}_n) \end{equation}

2つ目は4段4次Runge-Kutta法です.

【スキーム2: 4段4次Runge-Kutta法】

\begin{equation} \boldsymbol{x}_{n+1} = \boldsymbol{x}_n + \frac{h}{6}(\boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4) \end{equation}

ただし, \begin{align} \boldsymbol{k}_1 &= f(\boldsymbol{x}_n) \\ \boldsymbol{k}_2 &= f\left(\boldsymbol{x}_n + \frac{h}{2}\boldsymbol{k}_1\right) \\ \boldsymbol{k}_3 &= f\left(\boldsymbol{x}_n + \frac{h}{2}\boldsymbol{k}_2\right) \\ \boldsymbol{k}_4 &= f(\boldsymbol{x}_n + h\boldsymbol{k}_3) \end{align} とする.

微分方程式の数値解法の適用範囲を広げるため, 階数を1つあげます. つまり, 2階の微分方程式も考えることにします.

\begin{equation} \frac{d^2}{dt^2} \boldsymbol{x} = f(\boldsymbol{x}) \end{equation}

ここで, 次のような変数を導入します.

\begin{equation} \boldsymbol{p} = \frac{d}{dt}\boldsymbol{x} \end{equation}

このとき, 上の2階常微分方程式は以下のように書けます.

\begin{equation} \frac{d}{dt}\begin{bmatrix} \boldsymbol{x}\\ \boldsymbol{p} \end{bmatrix} = \begin{bmatrix} \boldsymbol{p} \\ f(\boldsymbol{x}) \end{bmatrix} \end{equation}

このように, 2階の微分方程式であっても, 1階に帰着します. したがって, 先ほどの数値解法, 陽的Euler法, 4段4次Runge-Kutta法も適用できます. もう1つだけ数値解法を導入します. 3つ目の数値解法が, leap-frog法です.

【スキーム3: leap-frog法(Störmer-Verlet法)】

\begin{align} \boldsymbol{p}_{n+\frac{1}{2}} &= \boldsymbol{p}_n + \frac{h}{2}f(\boldsymbol{x}_n) \\ \boldsymbol{x}_{n+1} &= \boldsymbol{x}_n + h \boldsymbol{p}_{n+\frac{1}{2}} \\ \boldsymbol{p}_{n+1} &= \boldsymbol{p}_{n+\frac{1}{2}} + \frac{h}{2} f(\boldsymbol{x}_{n+1}) \end{align}

数値解法を3つほど導入したところで, これらを比較します. 例題として, Kepler問題を考えます.

【Kepler問題】

次の微分方程式系を考える.

\begin{equation} \ddot{x_1} = -\frac{x_1}{(x_1^2 + x_2^2)^{\frac{3}{2}}}, \quad \ddot{x_2} = -\frac{x_2}{(x_1^2 + x_2^2)^{\frac{3}{2}}} \\ \end{equation}

この系の解の様子をそれぞれの解法について, プロットして確かめる. さらに, この系はHamilton系である. そこで, Hamiltonianの値の変化も確かめる.

上つきのドットは時間微分を表します. Kepler問題も1階の系に帰着します. 次のようにおきます.

\begin{equation} \boldsymbol{x} = \begin{bmatrix} x_1\\ x_2\end{bmatrix},\quad \boldsymbol{p} = \begin{bmatrix} \dot{x_1} \\ \dot{x}_2 \end{bmatrix} \end{equation}

よって, 次のように書けます.

\begin{equation} \frac{d}{dt}\boldsymbol{x} = \boldsymbol{p}, \quad \frac{d}{dt}\boldsymbol{p} = -\frac{\boldsymbol{x}}{\|\boldsymbol{x}\|^3} \end{equation}

以上より, \begin{equation} \frac{d}{dt}\begin{bmatrix} \boldsymbol{x}\\ \boldsymbol{p} \end{bmatrix} = \begin{bmatrix} \boldsymbol{p} \\ f(\boldsymbol{x}) \end{bmatrix}, \quad f(\boldsymbol{x}) = -\frac{\boldsymbol{x}}{\| \boldsymbol{x} \|^3} \end{equation} と書けます. 上式に様々な解法を適用します. また, Hamiltonianは, \begin{equation} H(\boldsymbol{x},\boldsymbol{p}) = - \frac{1}{\|\boldsymbol{x}\|} + \frac{1}{2}\| \boldsymbol{p} \|^2 \end{equation} で与えられます. 下図は, ステップ数150で, ステップ幅0.05とした場合の軌道のアニメーションです. 初期値は, 以下のように設定しました.

\begin{equation} \boldsymbol{x}_0 = \begin{bmatrix} 0.4 \\ 0.0 \end{bmatrix},\quad \boldsymbol{p}_0 = \begin{bmatrix} 0.0 \\ 2.0 \end{bmatrix} \end{equation}
【コード3の実行結果】

Euler法のE軌道はかなり外れています. その他の解法は一定の軌道を描いています. 以下に, Hamiltonianの時間変化を示します. 10000サンプル分計算しました.

【コード4の実行結果】

長時間軌道計算すると, Runge-Kutta法の誤差は大きくなります. 一方, leap-frog法(Störmer-Verlet法)は誤差が一定です. 特に, 長時間計算した際に, 顕著な差が現れます. この解法は, エネルギー保存というHamilton系の性質を保っています.

アルゴリズム

以下にHMCのアルゴリズムを示します. ここで, 目的の分布の密度関数は以下のように書けるとします.

\begin{equation} p(\boldsymbol{x}) = \frac{1}{Z}\tilde{p}(\boldsymbol{x}) \end{equation}

また, ポテンシャルエネルギーを以下のように定義します.

\begin{equation} U(\boldsymbol{x}) = -\log \tilde{p}(\boldsymbol{x}) \end{equation}

さらに, Hamiltonianを以下のように定義します.

\begin{equation} H(\boldsymbol{x}, \boldsymbol{p}) = U(\boldsymbol{x}) + \frac{1}{2}\| \boldsymbol{p}\|^2 \end{equation}
【Hamiltonian Monte Carlo法】

Initialize \(\boldsymbol{x}^{(0)}\) and set step size \(h\) and number of iteration of \(T\) of leap-frog.
for \(s=1,2, \cdots\)
Sample \(\boldsymbol{p}\) from \(\mathrm{N}(0,I)\).
Compute \((\boldsymbol{x}^{\prime},\boldsymbol{p}^{\prime}) = \mathrm{leapfrog}(\boldsymbol{x}^{(s-1)}, \boldsymbol{p}, h, T)\).
Compute \(\Delta H = H(\boldsymbol{x}^{\prime}, \boldsymbol{p}^{\prime})-H(\boldsymbol{x}^{(s-1)}, \boldsymbol{p})\).
Set \(\alpha = \min\left\{ 1,e^{-\Delta H}\right\} \).
Sample \(u\) from uniform distribution.
if \(u \leq \alpha\)
Set \(\boldsymbol{x}^{(s)}=\boldsymbol{x}^{\prime}\).
else
Set \(\boldsymbol{x}^{(s)}=\boldsymbol{x}^{(s-1)}\).
end if
end for

数値実験

2つほど実験結果を載せます.

例1: 2次元正規分布からのサンプル

これは以前にも用いた例です. 次の正規分布からサンプルします.

\begin{equation} \mathrm{N}(\boldsymbol{0}, \Sigma),\quad \Sigma = \begin{bmatrix} 1 & 0.99 \\ 0.99 & 1\end{bmatrix} \end{equation}

密度関数は次式で書けます.

\begin{equation} p(\boldsymbol{x}) \propto \mathrm{exp}\left( -\frac{1}{2}\boldsymbol{x}^{\mathrm{T}}\Sigma^{-1}\boldsymbol{x}\right) \end{equation}

ポテンシャルエネルギーは以下の2次形式です. ちなみにプログラムの方では, ポテンシャルエネルギーやその勾配を求めなくてもサンプルできるようにしています.

\begin{equation} U(\boldsymbol{x}) = \frac{1}{2}\boldsymbol{x}^{\mathrm{T}}\Sigma^{-1}\boldsymbol{x} \end{equation}

これくらい斜めに歪んでいると, Gibbs samplerではうまくいきません. 原点を初期値, サンプル数を5000, leap-frogの総ステップ数を100, 刻み幅を0.05としました. また, 最初の500サンプル分はburn-in期間として除去しました. 以下にサンプルの結果を示します.

【コード6の実行結果】

例2: 混合正規分布からのサンプル

次に, 混合正規分布からサンプルします. 平均が5角形を描くように正規分布を5つ混合しました.

\begin{equation} p(\boldsymbol{x}) = \frac{1}{5}\sum_{i=1}^{5} \mathrm{N}\left( \boldsymbol{\mu}_i, \frac{1}{2}I\right), \quad \boldsymbol{\mu}_i = \begin{bmatrix} 2\cos\left( i\theta+\frac{\pi}{10}\right) \\ 2\sin\left(i\theta+\frac{\pi}{10} \right)\end{bmatrix},\quad \theta = \frac{2\pi}{5} \end{equation}
【コード7の実行結果】

コード

【Juliaコード1; インポート】
using LinearAlgebra
using ForwardDiff
using Random 
using Statistics
using Distributions
using Plots
pyplot()
【Juliaコード2; 微分方程式の数値解法の実装】
#one step of explicit Euler method
function myEuler(qvec, pvec, h, f)
    q_new = qvec + h * pvec
    p_new = pvec + h * f(qvec)
    return q_new, p_new
end

#one step of RK method
function myRK(qvec, pvec, h, f)
    #内部段の計算
    kq1 = pvec
    kp1 = f(qvec)
    kq2 = pvec + h * kp1/2 #kq₂ = v_old + h*kv₁/2;
    kp2 = f(qvec + h * kq1/2) #kv₂ = f(q_old + h*kq₁/2);
    kq3 = pvec + h * kp2/2 #kq₃ = v_old + h*kv₂/2;
    kp3 = f(qvec + h * kq2/2) #kv₃ = f(q_old + h*kq₂/2);
    kq4 = pvec + h * kp3 #kq₄ = v_old + h*kv₃;
    kp4 = f(qvec + h * kq3) #kv₄ = f(q_old + h*kq₃);
    
    #pとqの更新
    q_new = qvec + h*(kq1+2*kq2+2*kq3+kq4)/6
    p_new = pvec + h*(kp1+2*kp2+2*kp3+kp4)/6
    return q_new, p_new
end

#one step of Störmer-Verlet method
function myStörmerVerlet(qvec, pvec, h, f)
    p_mid = pvec + h * f(qvec)/2;
    q_new = qvec + h * p_mid;
    p_new = p_mid + h * f(q_new)/2;
    return q_new, p_new
end

#gradient of potentiol energy
∇Urev(U, x) = -ForwardDiff.gradient(U, x)

#simulate Hamilton dynamics
function simHamilton(T, h, f, qvec_init, pvec_init, ODEsolver)
    #initialize
    d = length(qvec_init)
    qvecs = zeros(d, T)
    pvecs = zeros(d, T)
    qvecs[:,1] = qvec_init
    pvecs[:,1] = pvec_init
    
    #solve ODE iteratively
    for t in 2:T
        qvecs[:,t], pvecs[:,t] = ODEsolver(qvecs[:,t-1], pvecs[:,t-1], h, f)
    end
    return qvecs, pvecs
end
【Juliaコード3; Kepler問題】
#potentiol energy
U(q) = -1/norm(q)

#kinetic energy
K(p) = norm(p)^2/2

#calculate Hamiltonian
H(q, p, U, K) = U(q) + K(p)

#set the parameters and initialize
T = 150
h = 0.05
qvec_init = [0.4, 0.0]
pvec_init = [0.0, 2.0]

#simulate Hamilton dynamics using explicit Euler method
qvecs_E, pvecs_E = simHamilton(T, h, x->∇Urev(U,x), qvec_init, pvec_init, myEuler)

#simulate Hamilton dynamics using RK method
qvecs_RK, pvecs_RK = simHamilton(T, h, x->∇Urev(U,x), qvec_init, pvec_init, myRK)

#simulate Hamilton dynamics using Srörmer-Verlet method
qvecs_SV, pvecs_SV = simHamilton(T, h, x->∇Urev(U,x), qvec_init, pvec_init, myStörmerVerlet)

#animation
anim_KP = @animate for t in 1:T
    plot(title="Kepler problem", xlim=[-2,1], ylim=[-1,1.5], xlabel="q₁", ylabel="q₂")
    
    #Euler
    plot!(qvecs_E[1,1:t-1], qvecs_E[2,1:t-1], label="Euler")
    plot!([qvecs_E[1,t]], [qvecs_E[2,t]], st=:scatter, color=:red, markersize=5, label=false)
    
    #Runge-Kutta
    plot!(qvecs_RK[1,1:t-1], qvecs_RK[2,1:t-1], label="RK")
    plot!([qvecs_RK[1,t]], [qvecs_RK[2,t]], st=:scatter, color=:red, markersize=5, label=false)
    
    #Störmer-Verlet
    plot!(qvecs_SV[1,1:t-1], qvecs_SV[2,1:t-1], label="SV")
    plot!([qvecs_SV[1,t]], [qvecs_SV[2,t]], st=:scatter, color=:red, markersize=5, label=false)
end
#gif(anim_KP, "figs-HMC/anim_KP.gif")
【Juliaコード4; Hamiltonianの誤差】
#set the parameters and initialize
T = 10000
h = 0.05
qvec_init = [0.4, 0.0]
pvec_init = [0.0, 2.0]

#simulate Hamilton dynamics using explicit Euler method
qvecs_E, pvecs_E = simHamilton(T, h, x->∇Urev(U,x), qvec_init, pvec_init, myEuler)

#simulate Hamilton dynamics using RK method
qvecs_RK, pvecs_RK = simHamilton(T, h, x->∇Urev(U,x), qvec_init, pvec_init, myRK)

#simulate Hamilton dynamics using Srörmer-Verlet method
qvecs_SV, pvecs_SV = simHamilton(T, h, x->∇Urev(U,x), qvec_init, pvec_init, myStörmerVerlet)

#calculate Hamiltonian
exactH = H(qvec_init, pvec_init, U, K)
Herror_E = zeros(T) 
Herror_RK = zeros(T)
Herror_SV = zeros(T)
for t in 1:T
    Herror_E[t] = abs(H(qvecs_E[:,t], pvecs_E[:,t], U, K)-exactH)
    Herror_RK[t] = abs(H(qvecs_RK[:,t], pvecs_RK[:,t], U, K)-exactH)
    Herror_SV[t] = abs(H(qvecs_SV[:,t], pvecs_SV[:,t], U, K)-exactH)
end

#plot Hamiltonian
fig1 = plot(title="Error of Hamiltonian", xlabel="iter", ylabel="error of H(q,p)", lw=0.1, xscale=:log10, yscale=:log10, ylim=[1e-6,1e0])
plot!(1:T, Herror_E, label="Euler")
plot!(1:T, Herror_RK, label="RK")
plot!(1:T, Herror_SV, label="SV")
#savefig(fig1, "figs-HMC/fig1.png")
【Juliaコード5; HMCの実装】
#update the position
function update(T, h, f, qvec, pvec)
    qvec_new = qvec
    pvec_new = pvec
    for t in 1:T
        qvec_new, pvec_new = myStörmerVerlet(qvec_new, pvec_new, h, f)
    end
    return qvec_new, pvec_new
end

#MH acceptance and rejection
function accept_or_reject(xvec, xvec_old, pvec, pvec_old, H)
    ΔH = H(xvec, pvec)-H(xvec_old, pvec_old)
    α = min(1.0, exp(-ΔH))
    u = rand()
    if u≤α
        return xvec, pvec
    else
        return xvec_old, pvec_old
    end
end

#Hamiltonian Monte Carlo
function myHMC(x₀, n_samps, n_burnin, pdf_func, T, h)
    #initialization
    d = length(x₀)
    xsamps = zeros(d, n_samps)
    xsamps[:,1] = x₀
    xvec = zeros(d)
    pvec = zeros(d)
    
    #Hamiltonian
    U(xvec) = -log(pdf_func(xvec))
    ∇Uneg(xvec) = -ForwardDiff.gradient(U, xvec)
    H(xvec, pvec) = U(xvec) + norm(pvec)^2/2
    
    #sample
    xvec_old = x₀
    pvec_old = randn(d)
    for i in 2:n_samps
        pvec = randn(d)
        xvec, pvec = update(T, h, ∇Uneg, xvec, pvec)
        xvec, pvec = accept_or_reject(xvec, xvec_old, pvec, pvec_old, H)
        xsamps[:,i] = xvec
        xvec_old = xvec
        pvec_old = pvec
    end
    return xsamps[:,n_burnin:end]
end

#HMC for hierarchical model
function HMC_hm(λ₀, β₀, n_samps, n_burnin, T, h, Y, N, α, μβ)
    #initialization
    λsamps = zeros(N,n_samps)
    βsamps = zeros(n_samps)
    λsamps[:,1] = λ₀
    βsamps[1] = β₀
    λvec = zeros(N)
    pvec = zeros(N)
    
    #Hamiltonian
    pdf_func(xvec) = exp(dot(Y+(α+1)*ones(N), log.(xvec)))
    U(xvec) = -log(pdf_func(xvec))
    ∇Uneg(xvec) = -ForwardDiff.gradient(U, xvec)
    H(xvec, pvec) = U(xvec) + norm(pvec)^2/2
    
    #sample
    λvec_old = λ₀
    pvec_old = randn(N)
    for i in 2:n_samps
        #sample using leap frog
        pvec = randn(N)
        λvec, pvec = update(T, h, ∇Uneg, λvec, pvec)
        λvec, pvec = accept_or_reject(λvec, λvec_old, pvec, pvec_old, H)
        λsamps[:,i] = λvec
        λvec_old = λvec
        pvec_old = pvec
        
        #saple using Gibbs sampler
        βsamps[i] = rand(Gamma(α*N+1, 1/(N+μβ)))
    end
    return λsamps[:,n_burnin:end], βsamps[n_burnin:end]
end
【Juliaコード6; 正規分布からのサンプル】
#multivariate Normal distribution
μ = [0.0, 0.0]
Σ = [1.0 0.99; 0.99 1.0]
mv_normal = MvNormal(μ, Σ)
pdf_func(x) = pdf(mv_normal, x)

#sample from multivariate normal distribution
x₀ = zeros(2)
n_samps = 5000
n_burnin = div(n_samps,10)
T = 100
h = 0.05
xsamps = myHMC(x₀, n_samps, n_burnin, pdf_func, T, h)

#visualize
x1s = -3:0.1:3
x2s = -3:0.1:3
fig2 = plot(x1s, x2s, (x,y)->pdf_func([x,y]), st=:contour, xlim=[-3,3], ylim=[-3,3], xlabel="x₁", ylabel="x₂", title="Normal sample")
plot!(xsamps[1,:], xsamps[2,:], st=:scatter, label=false, alpha=0.3, markerstrokewidth=0, color=:green)
#savefig(fig2, "figs-HMC/fig2.png")
【Juliaコード7; 混合正規分布からのサンプル】
#Guassian mixture 
θ = 2*π/5
μ_func(i, θ) = 2*[cos(π/10+i*θ), sin(π/10+i*θ)]
Σ = I(2)/2
μ₀ = μ_func(0, θ)
μ₁ = μ_func(1, θ)
μ₂ = μ_func(2, θ)
μ₃ = μ_func(3, θ)
μ₄ = μ_func(4, θ)
mixed_normal = MixtureModel([MvNormal(μ₀, Σ), MvNormal(μ₁, Σ), MvNormal(μ₂, Σ), MvNormal(μ₃, Σ), MvNormal(μ₄, Σ)])
pdf_func(x) = pdf(mixed_normal, x)

#sample from Gussian mixture distribution
x₀ = zeros(2)
n_samps = 20000
n_burnin = div(n_samps,10)
T = 100
h = 0.05
xsamps = myHMC(x₀, n_samps, n_burnin, pdf_func, T, h)

#visualize
x1s = -3:0.1:3
x2s = -3:0.1:3
fig3 = plot(x1s, x2s, (x,y)->pdf_func([x,y]), st=:contour, xlim=[-3,3], ylim=[-3,3], xlabel="x₁", ylabel="x₂", title="Gaussian Mixture")
plot!(xsamps[1,:], xsamps[2,:], st=:scatter, label=false, alpha=0.3, markerstrokewidth=0, color=:green)
#savefig(fig3, "figs-HMC/fig3.png")
参考文献

      [1]R.Neal, Bayesian learning for Neural Network, 1995, Springer