【MCMC法】ハミルトニアンモンテカルロをわかりやすく解説｜ベイズ統計学
ハミルトニアンモンテカルロ（Hamiltonian Monte Carlo
ベイズ統計学において、ハミルトニアンモンテカルロ（Hamiltonian Monte Carlo, HMC）を利用するアプローチは、主に複雑な事後確率分布からの効率的なサンプリングに用いられます。
このプロセスを理解するためには、ベイズの定理とハミルトニアンの動力学の基本概念を組み合わせる必要があり、結構むずいと思います。
まず、国語的にMHとの差分を語ると以下のようになります
ひとまずMH法との差を理解
メトロポリスヘイスティングスは「提案分布 」から候補点をサンプリングし、それを受容・棄却する確率を計算して次の点を選択。受容されたら次のサンプルはその候補点となり、受容されなければ現状の点が再びサンプルとなった。一方ハミルトニアンモンテカルロは、物理学の力学系に着想を得たもので、補助変数 (運動量) を導入し、ハミルトニアン力学系に従ってサンプルを「飛ばす」ことで、より効率よくサンプリングを行う手法
まあいずれにせよ、両者とも最終的にはターゲットとする分布に従うマルコフ連鎖を生成しています。
事前知識
ベイズ推定
ベイズの定理は以下のような形でした。
$${P(\theta|D)=\frac{P(D|\theta)P(\theta)}{P(D)}}$$
- 事後確率 (${P(θ|D)}$): データを観測した後に、パラメータ${\theta}$について得られる新しい知識を表します。
- 尤度 (${P(D|θ)}$): あるパラメータ${\theta}$のもとで、観測されたデータDが生じる確率を表します。つまり、データがどれだけそのパラメータ${\theta}$を支持しているかを示します。
- 事前確率 (${P(θ)}$): データを観測する前に、パラメータ${\theta}$について持っていた信念を表します。
- 周辺尤度 (${P(D)}$): すべての可能なパラメータ${\theta}$を考慮した上で、データDが生じる確率を表します。これは、通常、直接計算が難しい場合が多いです。別名エビデンスですね。
ベイズの定理は、尤度と事前確率を使って事後確率を計算するための公式ということですね。
ベイズ推定の手順
- 事前分布の設定: パラメータθに関する事前知識に基づいて、事前分布(P(θ))を仮定します。
- 尤度の計算: 尤度関数(P(D|θ))を計算します。
- 事後分布の計算: ベイズの定理を用いて、事後分布(P(θ|D))を計算します。
- 事後分布からのサンプリング: 事後分布からサンプルを生成し、パラメータθの推定や予測を行います。
で、上の手順４にある事後分布が複雑な形をしている場合、直接計算することが難しいことがあります。
このような場合に、マルコフ連鎖モンテカルロ法(MCMC)を用いて、事後分布からサンプルを生成する、という話でしたね。
MCMC法は、マルコフ連鎖を用いて、事後分布の高い確率密度を持つ領域を探索する手法です。
MCMCに関する基礎的な解説は、WEB版青の統計学の以下のコンテンツがおすすめです。
MCMC法
さて、MCMC法（ハミルトニアンモンテカルロもその一つ）についても少し補足をしておきます。
解析的に解けると思いますがコイン投げの例で説明します。
コインを${N}$回投げた結果、表が${K}$回出たとします。このとき、コインの裏が出る確率${\theta}$をベイズ推定で求めてみましょう。
事前分布の設定
θの事前分布として、ベータ分布(${Beta(\alpha, \beta)}$)を仮定します。
ベータ分布は、確率の範囲(0から1)をとる変数に対してよく用いられる分布です。αとβは、事前知識に基づいて設定するハイパーパラメータです。
尤度の計算
ベルヌーイ分布を用いて、尤度を計算します。
$${P(D∣θ)=θ^K(1−θ)^{N−K}}$$
事後分布の導出
ベイズの定理より、事後分布は次のようになります。
$${P(θ∣D)∝P(D∣θ)P(θ)=θ^K(1−θ){N−K}×θ^{α−1}(1−θ)^{β−1} }$$
これは、(Beta(α+K, β+N-K))に比例します。
つまり、事後分布もベータ分布になります。
事後分布が解析的に求まったため、厳密なサンプリングは不要ですが、MCMCを用いてサンプリングする手順を示します。
ここでは、最も簡単なメトロポリスヘイスティングス法を使います。
メトロポリスヘイスティングス法
初期値の設定
${θ}$の初期値(${θ^{(0)}}$)を、事前分布からランダムにサンプリングします。
提案分布からサンプリングしないの？と疑問に思った方は鋭いです。初期値は、事前分布からランダムにサンプリングするのが一般的です。事前知識に基づいた最も合理的な出発点と考えることができるからですね。
新しいθの提案
現在の${θ^{(t)}}$から、少しだけランダムに変化させた新しい${θ^{(t+1)}}$を提案します。
例えば、正規分布などを用いて提案できます。
受容確率の計算
メトロポリス・ヘイスティングス法の受容確率を計算します。
$${\alpha = \min \left( 1, \frac{P(\theta^{t+1} \mid D) \, q(\theta^t \mid \theta^{t+1})}{P(\theta^t \mid D) \, q(\theta^{t+1} \mid \theta^t)} \right)}$$
${q(\theta^{(t+1)}|\theta^{(t)})}$は提案分布です。
状態の更新
一様乱数uを生成し、$${u \leq α}$$であれば、${θ^{(t+1)}}$を受け入れ、${θ^{(t+1)} = θ^{(t)}}$とします。
で、上記のステップを繰り返し行います。
ハミルトニアン力学の話
さて、HMCに話を戻します。
ハミルトニアンモンテカルロ（Hamiltonian Monte Carlo, HMC）に関してきっちり理解を得るためには、まず物理学におけるハミルトニアンの概念についての基礎知識が必要です。
ハミルトニアン（系の全エネルギー）は、物体のエネルギーを記述するための関数であり、力学的エネルギーの総和、すなわち運動エネルギーと位置エネルギー（ポテンシャルエネルギー）の和として表されます。
ハミルトニアン \(H\) は、運動量 \(p\) と位置 \(q\) を使って次のように表されます
$$H(p,q)=T(p)+V(q)$$
ここで、\(T(p)\)は運動エネルギー（例えば\(\frac{p^2}{2m}\)）、\(V(q)\)は位置エネルギーを示します。
物体の運動を追跡するには、ハミルトンの正準方程式を解く必要があります。
これらは、時間に関して運動量と位置の導関数を与えるもので、以下の形で表されます。
$$\frac{dq}{dt}=\frac{∂H}{∂p}$$
$$\frac{dp}{dt}=-\frac{∂H}{∂q}$$
ハミルトニアンモンテカルロでは、このハミルトンの正準方程式を用いて、高次元の確率分布に沿った効率的な探索を実現します。
もう少し詳細に説明すると、この正準方程式は系の状態が時間と共にどのように変化するかを記述しているので、正準方程式を解くことで、物理的に意味のある方法（エネルギー保存の原則に従う）でパラメータ空間を探索できます。
これにより、効率的かつ効果的にターゲット分布からのサンプリングが可能になります。
HMCのアプローチ
MCMCと同様にやりたいことは、事後確率分布 \(P(\theta|D)\)からサンプリングです。
しかし、直接サンプリングするのは難しいため、ハミルトニアン力学を用いてこの問題を解決します。
先ほどのハミルトニアン方程式は、ベイズ統計の文脈だと以下のように解釈します。
$$H(\theta|p)=T(p)+V(\theta)=T(p)-log(P(\theta|D))$$
ベイズ文脈では、ポテンシャルエネルギーは通常、負の対数事後確率（尤度と事前確率の負の対数の合計）として定義されます。
また、\(T(p)\)は運動エネルギーを表し、通常、運動量の二次形式として定義されます。
ここで\(p\)は運動量のベクトルです。
次に正準方程式についてです。
パラメータ \(\theta\)と運動量\(p\) の時間発展は、以下の正準方程式によって記述されます：
$$\frac{d \theta}{dt}=\frac{∂H}{∂p}$$
$$\frac{dp}{dt}=-\frac{∂H}{∂ \theta}$$
このあとは、数値積分を行います。
上の方程式は通常、リープフロッグ法と呼ばれる数値積分法を用いて近似的に解かれます。
これにより、新しいパラメータと運動量の候補が提案されます。
leap flog法
リープフロッグ法は、数値解析のアルゴリズムで、特に微分方程式の解を近似的に求める手法として重要です。
リープフロッグ法は、オイラー法よりも高精度な方法であり、特に振動系や保存系のシミュレーションに適しています。
この方法は、半ステップずつ進むことで、オイラー法の数値不安定性を改善しています。
リープフロッグ法は次のような手順で進行します。
まず、位置 \(q\) と運動量 \(p\)を初期化し、次に以下のステップを繰り返します。
1：半ステップの時間 \(\frac{h}{2}\)で運動量 \(p\)を更新する。
２：完全なステップの時間 \(h\)で位置 \(q\) を更新する。
３：さらに半ステップの時間 \(\frac{h}{2}\)で運動量 \(p\)を再び更新する。
これにより、次の時間ステップにおける \(q\)と\(p\) の値を得ます。数学的には、以下の式で表されます。
$$p_{n+\frac{1}{2}}=p_n+\frac{1}{2}×f_q(q_n,p_n)$$
$$q_{n+1}=q_{n}+h×f_p(q_{n},p_{n+1})$$
$$p_{n+1}=p_{n+\frac{1}{2}}+\frac{h}{2}×f_q(q_{n+1},p_{n+1})$$
正準方程式を解くことで得られた新しい状態\(p,q\)が得られた後、それを受容するかどうかというプロセスに移ります。
受容確率は以下のようになります。
$$min(1,exp(-H(\theta’,p’)+H(\theta,p)))$$
新しい状態を受容するかどうかは、メトロポリス基準（一定の確率で受容・拒否）によって決定されます。
この確率は、新しい状態が古い状態よりも低いハミルトニアン（すなわち低いエネルギー）を持つ場合に高くなります。
そして上のプロセスを繰り返し、多数のサンプルを生成します。
これにより、事後分布の表現が得られます。
→というのも詳細釣り合いの原則があるからですね。
メトロポリス受容ステップは、マルコフ連鎖が詳細釣り合いの原則を満たすことを保証します。
これにより、長期的にはターゲット分布に従ったサンプリングが保証されます。
この方法では、ランダムウォークメトロポリスヘイスティングスのように単純なランダムサンプリングに依存せず、複雑な確率分布、特に高次元での探索が効率的に行えます。
HMCの基本的なアイディアは、各反復でランダムに運動量を初期化し、それを使用してハミルトンの正準方程式を数値的に解くことです。
この過程で生成された軌道は、新しいサンプル候補を提供します。
この新しいサンプルは、特定の条件（詳細釣り合い条件など）に基づいて受け入れられるか、拒否されます。
このプロセスは、サンプリングを目指す分布に従って運動量を「誘導」し、より効率的な探索を可能にします。
ただHMCは計算コストが高く、数値積分のためのステップサイズやステップ数などの調整が必要な場合があります。
また、運動量の初期化や数値積分のエラーによっては、提案されたサンプルが拒否される可能性があります。
しかし、これらのデメリットを上回る効率と精度を提供するため、ベイズ統計学において広く利用されています。
