「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

【ベイズ統計学】入門 ベイズ統計学(07/X)

 Bayes統計学を今回は

を主に参照しつつ学んでいく。

5. Hamiltonian Monte Carlo法

 MCMCの実践的方法としてHamiltonian Monte Carlo法を扱った事後分布の評価方法を学ぶ。

5.1 Hamiltonian Monte Carlo法の必要性

 ランダムウォークMH法*1を成功させるためには、候補


\begin{aligned}
a=\theta^{(t)}+\varepsilon
\end{aligned}

において、\varepsilonの分散{\sigma_{\varepsilon}}^2を適切に指定する必要がある。
 母数が少なければ、事後分布を目視できるような状態であれば、\sigma_{\varepsilon}のの調整は比較的に容易である。しかし高次元となれば、そのコントロールはほぼ不可能になる。Bayesモデルの事後分布をオールマイティにシミュレーションを行うためには、高次元空間でも平均移動距離を長く取ることが出来、同時に受容率を高く保つことが出来るような方法が必須となる。このための方法がHMC法である。

5.2 Hamiltonianとは何か

 Hamiltonianは解析力学で導入された概念であり、まずは解析力学の世界を前提にこれを説明する。

 ある運動系において時間を表す変数として\tau、位置を表す変数として\theta(\tau)を導入する。このとき速度および加速度を


\begin{aligned}
v(\tau)&=\displaystyle{\frac{d\theta}{d\tau}},\\
a(\tau)&=\displaystyle{\frac{dv}{d\tau}}=\displaystyle{\frac{d^2\theta}{d\tau^2}}
\end{aligned}

で定義する。
 更に新たな概念として。動いている物体を止めるときに感じられる方向性を有した物理量として運動量p(\tau)を質量をmとして定義する。


\begin{aligned}
p(\tau)=mv(\tau)
\end{aligned}

5.2.1 運動エネルギー

 運動する物体はエネルギーを持つ。エネルギーは仕事をする能力であり、


\begin{aligned}
E(\tau)=F(\tau) \left(\theta(\tau)-\theta(0)\right)
\end{aligned}

で定義できる。
 物体が高所にあることで蓄えられているエネルギーがポテンシャルエネルギーである。
 ポテンシャルエネルギーU(\tau)は重力加速度をgとして


\begin{aligned}
U(\tau)=&F(\tau)\left(\theta(\tau)-\theta(0)\right)\\
=&mgh(\theta(\tau))
\end{aligned}

で表すことが出来る(h(\theta)は高さを表す変数である。)。
 ここから動いている物体が有する仕事をする能力を運動エネルギーという。運動エネルギーK(\tau)


\begin{aligned}
K(\tau)=&ma(\tau)\cdot\ \displaystyle{\frac{1}{2}a(\tau)\tau^2}\\
           =&\displaystyle{\frac{1}{2}mv(\tau)^{2}}\\
           =&\displaystyle{\frac{1}{2m}p(\tau)^2}
\end{aligned}

5.2.2 Hamiltonian

 運動エネルギーとポテンシャルエネルギーの総和を力学的エネルギーといい、


\begin{aligned}
H(\tau)=K(\tau)+U(\tau)
\end{aligned}

と書き、これをまたHamiltonianという。摩擦や空気抵抗が無い理想状態では力学的エネルギーは保存される。

5.2.3 Hamiltonianを用いた運動方程式

 力学的エネルギー保存の法則より、


\begin{aligned}
H(\tau)=K(\tau)+U(\tau)=K\in\mathbb{R}
\end{aligned}

が成り立つ。この両辺を\tauに関して微分してから移項すれば


\begin{aligned}
\displaystyle{\frac{d K(\tau)}{d\tau}}=-\displaystyle{\frac{d U(\tau)}{d\tau}}
\end{aligned}

が得られる。
 ここで運動エネルギーは運動量の関数、ポテンシャルエネルギーは位置の関数であるから、


\begin{aligned}
\                         &\displaystyle{\frac{d K(\tau)}{d\tau}}=-\displaystyle{\frac{d U(\tau)}{d\tau}}\\
\ \Leftrightarrow&\displaystyle{\frac{d p(\tau)}{d\tau}}\displaystyle{\frac{d K(\tau)}{d p(\tau)}}=-\displaystyle{\frac{d U(\tau)}{d\theta(\tau)}}\displaystyle{\frac{d \theta(\tau)}{d\tau}}\\
\ \Leftrightarrow&\displaystyle{\frac{p(\tau)}{m}}\displaystyle{\frac{d K(\tau)}{d p(\tau)}}=-\displaystyle{\frac{d U(\tau)}{d\theta(\tau)}}v(\tau)\\
\ \Leftrightarrow&\displaystyle{\frac{p(\tau)}{m}}\displaystyle{\frac{d K(\tau)}{d p(\tau)}}=-\displaystyle{\frac{d U(\tau)}{d\theta(\tau)}}\displaystyle{\frac{p(\tau)}{m}}
\end{aligned}

が成り立つ。
 したがって


\begin{aligned}
\displaystyle{\frac{d p(\tau)}{d\tau}}=&-\displaystyle{\frac{dU(\theta(\tau))}{d\theta(\tau)}},\\
\displaystyle{\frac{d\theta(\tau)}{d\tau}}=&\displaystyle{\frac{dK(p(\tau))}{dp(\tau)}}
\end{aligned}

が得られる。これをHamitonの運動方程式という。

5.2.4 リープフロッグ法

 Hamitonの運動方程式を具体的に解く。


\begin{aligned}
\                         &\displaystyle{\frac{d p(\tau)}{d\tau}}=-\displaystyle{\frac{dU(\theta(\tau))}{d\theta(\tau)}}\\
\Leftrightarrow  & F(\theta(\tau))=-\displaystyle{\frac{d }{d\theta(\tau)}}U(\theta(\tau))=-\displaystyle{\frac{d }{d\theta(\tau)}}mgh(\theta(\tau))
\end{aligned}

である。また


\begin{aligned}
\                         &\displaystyle{\frac{d\theta(\tau)}{d\tau}}=\displaystyle{\frac{dK(p(\tau))}{dp(\tau)}}\\
\Leftrightarrow  &v(\tau)= \displaystyle{\frac{d}{dp(\tau)}\frac{1}{2m}}p(\tau)^2=\displaystyle{\frac{1}{m}}p(\tau)
\end{aligned}

である。
 簡略化のためにm=g=1とすると


\begin{aligned}
\displaystyle{\frac{d p(\tau)}{d\tau}}&=-h^{\prime}(\theta(\tau)),\\
\displaystyle{\frac{d\theta(\tau)}{d\tau}}&=p(\tau)
\end{aligned}

が成り立つ。これは位置の微小な異動は運動量に等しいことを意味する。これを用いて経路を予測する方法がEuler法である。
 \varepsilon\in\mathbb{R}として


\begin{aligned}
p(\tau+1)&=p(\tau)-\varepsilon h^{\prime}(\theta(\tau)),\\
\theta(\tau+1)&=\theta(\tau)+\varepsilon p(\tau)
\end{aligned}

である。しかしHamiltonianの保存の精度は悪いことが知られているため、Hamiltonianをより精度よく保存するために、


\begin{aligned}
p\left(\tau+\displaystyle{\frac{1}{2}}\right)&=p(\tau)-\displaystyle{\frac{\varepsilon}{2}}h^{\prime}(\theta(\tau))\\
\theta(\tau+1)&=\theta(\tau)+\varepsilon p\left(\tau+\displaystyle{\frac{1}{2}}\right)\\
p(\tau+1)&=p\left(\tau+\displaystyle{\frac{1}{2}}\right)-\displaystyle{\frac{\varepsilon}{2}}h^{\prime}(\theta(\tau+1))
\end{aligned}

5.3 Hamiltonian Monte Carlo法のアルゴリズム

 ここまでで導入したHamiltonianを用いて、「遷移距離を大きくし、なおかつ候補の需要率を高め、事後分布に従う乱数を得る方法」を考える。
 事後分布f(\theta|\boldsymbol{x})と、それとは独立な標準正規分布確率密度関数f(p)との同時分布


\begin{aligned}
f(\theta,p|\boldsymbol{x})=f(\theta|\boldsymbol{x})f(p)
\end{aligned}

を考える。この方法では、この同時分布から乱数を発生させる。このときf(\theta|\boldsymbol{x})f(p)が互いに独立であるから、f(\theta,p|\boldsymbol{x})からの乱数\thetaf(\theta|\boldsymbol{x})からの乱数と等しい。
 標準正規分布確率密度関数f(p)カーネルを取り出して


\begin{aligned}
f(p)\propto \exp\left(-\displaystyle{\frac{1}{2}p^2}\right)
\end{aligned}

を考えておく。
 ここで正の値しか取らない関数f\gt0に対して


\begin{aligned}
f=\exp\left(\log{f}\right),\\
f=\log\left(\exp{(f)}\right)
\end{aligned}

が成り立つ。これを利用すると、


\begin{aligned}
f(\theta,p|\boldsymbol{x})&=\exp\left(\log{f(\theta,p|\boldsymbol{x})}\right)\\
                                         &=\exp\left\{\log{f(\theta|\boldsymbol{x})}+\log\left(f(p)\right)\right\}\\
                                         &\propto \exp\left[\log{f(\theta|\boldsymbol{x})}+\log\left\{\exp\left(-\displaystyle{\frac{1}{2}p^2}\right)\right\}\right]\\
                                         &=\exp\left(\log{f(\theta|\boldsymbol{x})}-\displaystyle{\frac{1}{2}p^2}\right)
\end{aligned}

ここでh(\theta)=-\log(f(\theta|\boldsymbol{x}))と仮におくと


\begin{aligned}
f(\theta,p|\boldsymbol{x})&\propto \exp\left(-h(\theta)-\displaystyle{\frac{1}{2}p^2}\right)=\exp\left(-H(\theta,p)\right)
\end{aligned}

となり、事後分布がHamiltonianと等価な式で書けることが導かれた。


Hamiltonian Monte Carlo法のアルゴリズム(母数が1つの場合)
(1) 初期値\theta^{(1)},\varepsilon,L,T,\ バーンイン期間を定め、t=1とする。
(2) 独立な標準正規乱数p^{(t)}を発生させる。
(3) リープフロッグ法で遷移させ、候補点\theta^{(a)},p^{(a)}を入手する。
(4) 確率\min(1,r)で受容(\theta^{(t+1)}=\theta^{(a)})し、さもなければ\theta^{(t+1)}=\theta^{(t)}とする。
(5)T=tならば終了し、T\lt tならばt=t+1として(2)に戻る。

 一般には、同じLに対して\varepsilonを大きくすると移動距離が長くなる代わりに受容率が下がる。逆に\varepsilonを小さくすると受容率が高くなる代わりに移動距離が短くなるため相関が大きくなる。Lを大きくする一方で\varepsilonを小さくすると、移動距離が長く受容率も上がるものの、1回の遷移のための計算コストが大きくなる。

5.4. 多次元の場合

 これまでは母数が1次元の場合を取り扱ってきた。今度は\boldsymbol{\theta}=(\theta_1,\cdots,\theta_d),d\gt1を推定する場合を考える。独立なd個の標準正規乱数\boldsymbol{p}=(p_1,\cdots,p_d)との同時事後分布は


\begin{aligned}
f(\boldsymbol{\theta},\boldsymbol{p}|\boldsymbol{x})=f(\boldsymbol{\theta}|\boldsymbol{x})f(\boldsymbol{p})=f(\boldsymbol{\theta}|\boldsymbol{x})\displaystyle{\prod_{i=1}^{d}f(p_i)}
\end{aligned}

である。独立なd次元標準正規乱数のカーネル


\begin{aligned}
f(\boldsymbol{p})=\displaystyle{\prod_{i=1}^{d}f(p_i)}\propto \exp\left(-\displaystyle{\frac{1}{2}\sum_{i=1}^{d}{p_i}^2}\right)
\end{aligned}

であるから、Hamiltonianは


\begin{aligned}
H(\boldsymbol{\theta},\boldsymbol{p})=h(\boldsymbol{\theta})+\displaystyle{\frac{1}{2}\sum_{i=1}^{d}{p_i}^2}
\end{aligned}

で与えられる。
 同時事後分布のカーネル


\begin{aligned}
f(\boldsymbol{\theta},\boldsymbol{p}|\boldsymbol{x})\propto \exp\left(-H(\boldsymbol{\theta},\boldsymbol{p})\right)
\end{aligned}

である。
 母数によるポテンシャル・エネルギーの微分を求めるべく、ベクトルの微分を確認しておくと、


\begin{aligned}
h^{\prime}(\boldsymbol{\theta})=\left(\displaystyle{\frac{\partial h(\boldsymbol{\theta})}{\partial \theta_1}},\cdots,\displaystyle{\frac{\partial h(\boldsymbol{\theta})}{\partial \theta_d}}\right)
\end{aligned}

である。
 補正係数は


\begin{aligned}
r=\exp(H(\boldsymbol{\theta}^{(t)},\boldsymbol{p}^{(t)})-\boldsymbol{\theta}^{(\alpha)},\boldsymbol{p}^{(\alpha)}))
\end{aligned}

とする。
 以上をまとめると以下のとおりとなる。


Hamiltonian Monte Carlo法のアルゴリズム(母数が複数の場合)
(1) 初期値\boldsymbol{\theta}^{(1)},\varepsilon,L,T,\ バーンイン期間を定め、t=1とする。
(2) 独立なd個の標準正規乱数\boldsymbol{p}^{(t)}を発生させる。
(3) リープフロッグ法で遷移させ、候補点\boldsymbol{\theta}^{(a)},\boldsymbol{p}^{(a)}を入手する。
(4) 確率\min(1,r)で受容(\boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(a)})し、さもなければ\boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(t)}とする。
(5)T=tならば終了し、T\lt tならばt=t+1として(2)に戻る。

プライバシーポリシー お問い合わせ