覚え書き: BDA3 Chap.12 計算効率の良いマルコフ連鎖シミュレーション

　引き続き、BDA3のメモ。
　Gelman, A. et al.(2014) Bayesian Data Analysis, Third Edition. Chapter 12. Computationally efficient Markov chain simulation.

　計算に関わる解説10-13章のうち中盤戦。いよいよHMCが登場する。
　実をいうと、これまで何度かHMCについて勉強しようと試みたことがあったのだが、そのたびに挫折していたのである。とにかくですね、物理学の用語が出てくると心が折れるのである。高校時代のトラウマかしらん。
　幸い本章は、HMCについても物理学ぬきで説明して下さる模様である。はい深呼吸して…

12.1 効率的なGibbsサンプラー
変換と再パラメータ化
　Gibbsサンプラーは、要素が独立であるようにパラメータ化されているときに効率が良くなる。再パラメータ化のもっとも簡単な方法として、パラメータを線形変換することが挙げられるが、事後分布が近似的に正規分布に従わない場合は特別な手法が必要になる。
　メトロポリス法のジャンプについても同じことがいえる。[…]

補足変数
　Gibbsサンプラーの計算は、補足変数を追加することでシンプルになったり収束が加速したりする。
　たとえば、t分布は正規分布の混合として表現できる。いま$t_\nu(\mu, \sigma^2)$ から独立に得た$n$個のデータ点があるとして、$\mu, \sigma^2$について推論するとしよう。話を簡単にするため$\nu$は既知とする。$\mu, \log \sigma$に一様事前分布を与える。
　それぞれのデータ点の尤度は以下と等価になる。$V_i$が補足変数。$$ y_i \sim N(\mu, V_i) $$ $$ V_i \sim \mathrm{Inv-} \chi^2 (\mu, \sigma^2) $$ 同時分布$p(\mu, \sigma^2, V | y)$を使って推論し、$\mu, \sigma$のことだけ考えれば良い。
　では$\mu, \sigma^2$をどうサンプリングすれば良いかというと…

それぞれの$V_i$の条件付き事後分布を得る。$y$で条件付けると$V_i$の事前分布は逆カイ二乗分布なので、事後分布も逆カイ二乗分布となり、$$ V_i | \mu, \sigma^2, \nu, y \sim \mathrm{Inv-}\chi^2 \left( \nu+1, \frac{\nu \sigma^2 + (y_i – \mu)^2}{\nu + 1} \right) $$
$\mu$の条件付き事後分布を得る。$$ \mu | \sigma^2, V, \nu, y \sim N \left( \frac{\sum_i (1/V_i) y_i}{\sum_i (1/V_i)} , \frac{1}{\sum_i (1/V_i)} \right) $$
$\sigma^2$の条件付き事後分布を得る。$$ p(\sigma^2 | \mu, V, \nu, y) \propto \sigma^{-2} \prod_i \sigma^\nu \exp(-\frac{\nu \sigma^2}{2V_i}) $$ [この辺になると式の成り立ちがさっぱりわかない。省略するけど、結局右辺はガンマ分布になるのだそうだ]

パラメータ拡張
　パラメータ間で事後分布が依存していてGibbsサンプラーの収束が遅いとき、パラメータを追加してもっと広い空間をランダムウォークさせてやるとなぜか早くなることがある。
　上の例だと、$\sigma$のドローがゼロに近いとき、$V_i$のドローもゼロに近くなり、また$\sigma$のドローもゼロに近くなり… でなかなか抜け出せなくなる。そこで、$\tau$と$V_i$の間の依存性を破り脱出させるためだけに、それ自身には特に意味の無い新パラメータ$\alpha > 0$を付け加えよう。$$ y_i \sim N(\mu, \alpha^2 U_i) $$ $$ U_i \sim \mathrm{Inv-}\chi^2(\nu, \tau^2)$$ $\alpha$には対数スケール上で一様事前分布を与える。

$U_i$を更新する: $$ U_i | \alpha, \mu, \tau^2, \nu, y \sim \mathrm{Inv-}\chi^2 \left( \nu+1, \frac{\nu \tau^2 + (\frac{y_i – \mu}{\alpha})^2}{\nu+1} \right) $$
$\mu$を更新する。[上の例の式と同じなのでメモ省略]
$\tau^2$を更新する。[上の例でいう$\sigma^2$と同じなので省略。結局ガンマ分布になる]
$\alpha^2$を更新する。$$ \alpha^2 | \mu, \tau^2, U, \nu, y \sim \mathrm{Inv-}\chi^2 \left(n, \frac{1}{n} \sum_i \frac{(y_i – \mu)^2}{U_i} \right) $$

$\alpha^2, U, \tau$は識別できない。しかし、$\mu, \sigma = \alpha \tau, V_i = \alpha^2 U_i$の収束を監視する限りにおいて、モデル全体は識別できる。[なるほど…]

12.2 効率的なメトロポリス・アルゴリズム
　メトロポリス・アルゴリズムのジャンピング・ルールには無数のバージョンがあるが、単純なルールについてみると大きくふたつにわかれる。

本質的に、パラメータ空間をランダムウォークする奴。正規カーネルを使い、その平均はパラメータの現在の値とし、その分散は工夫する、という場合が多い。
目標分布を近似する提案分布を作る奴。ここでいう目標分布はGibbsサンプラーにおける一部についての条件付き分布だったり、同時事後分布だったりする。一度に変えるパラメータがひとつだけだということには自然な利点はなく、単に計算の節約に過ぎない。

　ジャンピング・ルールを効率よくする一般的な方法を考えるのは難しいが、ランダム・ウォークについては次のことがわかっている。$d$個のパラメータの事後分布$\theta = (\theta_1, \ldots, \theta_d)$が、なんらかの適切な変換をすれば、既知の分散行列$\Sigma$のMVNに従う、としよう。正規ジャンピングカーネルによるメトロポリス・アルゴリズム $$ J(\theta^* | \theta^{t-1}) = N(\theta^* | \theta^{t-1}, c^2 \Sigma)$$ を使うとする。このときもっとも効率的なのは、$c \approx 2.4 \sqrt{d}$とすることだ。このとき効率性はおよそ$0.3/d$となる。ちなみにGibbsサンプラーの効率は、パラメータの事後分布が互いに独立なら$1/d$である。$d$回のドローでまっさらな$\theta$が得られるわけだから。

　メトロポリス・アルゴリズムのもうひとつの特徴として、ジャンプの受容率に注目することもできる。
　MVNランダムウォークでジャンピング・カーネルが目標分布と同じ形状であるとき、最適なルールは、受容率がある次元でおよそ0.44、高次元(およそ$d > 5$)ではおよそ0.23となる。ここから以下の適応的アルゴリズムを考えることができる。Gibbsサンプラーなりメトロポリスなりでシミュレーションを始めて、しばらく続けたのち、以下のメトロポリス・ジャンピング・ルールを使う。

ジャンピング分布の共分散が、シミュレーションから推定される事後共分散行列と比例するように調整し、
受容率が高すぎたらジャンピング分布のスケールを増やし、受容率が低すぎたら減らす。目標は、一次元なら0.44, 多次元なら0.23に近くなるようにすることである。

[えええ… 最適なルールでは受容率がこれこれになるからといって、受容率をそれに近づけたら最適ルールに近づくといえるのだろうか？いえるってことか…]

適応的アルゴリズム
　反復シミュレーション・アルゴリズムを走らせながらチューニングする際は、変な分布に収束するのを避けるように気をつけなければならない。
　たとえば、分布の平たいエリアではより早く動き、事後密度が素早く変わったら遅くする、という適応的工夫をしたとしよう。これは目標分布の探索方法としてはいいけれど、得られたドローの分布は目標分布と一致しない。[あー、そりゃそうだな]
　安全のため、適応的アルゴリズムは適応的フェイズと固定的フェイズに分けておき、最後の推論には後者のドローのみを使う、というやりかたにすることが多い。

12.3 Gibbs・メトロポリスのさらなる拡張

スライス・サンプリング
　$d$次元目標分布$p(\theta | y)$からの$\theta$のランダムサンプルは、その分布の下にあるエリアからのランダムサンプルと等価である。実は$(\theta, u)$の$d+1$次元分布からのサンプリングであり、それは$u \in [0, p(\theta|y)]$のときに一様分布$p(\theta, u|y) \propto 1$であり、そうでなかったら0である、と考えても良い。
　この一様分布上でのシミュレーション・アルゴリズムのことをスライス・サンプリングという。実装は複雑だが、応用範囲は広い。特に、Gibbsサンプリング構造における一次元条件つき分布のサンプリングにおいて有用である。

次元の異なる空間の間を移動するためのリバーシブル・ジャンプ・サンプリング
　trans-dimensionalなMCMCをやりたくなることも多い。つまり、パラメータ空間の次元がイテレーションごとに変わるようなMCMCである。たとえばモデル平均では、多数の可能なモデルたち(たとえば予測子セットが異なる回帰モデルたち)の間を動くマルコフ連鎖をつくりたくなる。また有限混合モデルでは、混合する素数を変動させてやりたい。
　こういうメトロポリス・アルゴリズムは、リバーシブル・ジャンプ・サンプリングという方法でつくれる。

　マルコフ連鎖が多数の候補モデル$M_k (k=1, \ldots, K)$の間を移動するという場合について考える。パラメータの次元数を$d_k$、パラメータ・ベクトルを$\theta_k$とする。
　補足変数$u$を考え、$k$から$k^*$に移動したときのジャンピング分布を$J(u|k, k^*, \theta_k)$とする。さて、異なる次元を1対1にマッチさせる決定論的な関数を定義する。$$ (\theta_{k^*}, u^*) = g_{k,k^*} (\theta_k, u) $$ $$ d_k + dim(u) = d_{k^*} + dim(u^*) $$ [$dim(u)$というのは$u$の次元数じゃなくて、$u$によって決まる次元数増減関数のことだろう]
　たとえば…
　[以下、事例による説明が続くが、いまあんまし関心無いのでスキップする]

simulated temperingと並列tempering
　[焼きなまし法のこと? まあいいや、いつか必要になったら読もう。スキップ]

粒子フィルタリング、ウェイティング、遺伝的アルゴリズム
　[それぞれについて1段落づつ紹介。スキップ]

12.4 ハミルトニアン・モンテカルロ
　ハミルトニアン・モンテカルロ(HMC)は、メトロポリス・アルゴリズムの局所ランダムウォークを抑制するために、物理学からあるアイデアを借りてきた方法である。目標空間の要素$\theta_j$ごとに、「モメンタム」変数$\phi_j$を追加する。$\theta$と$\phi$を一緒に更新するが、その際$\theta$のジャンピング分布は主に$\phi$で決まる。MCMCと決定論的シミュレーションを結合しているので、ハイブリッド・モンテカルロともいう。
　HMCでは、事後密度$p(\theta|y)$はモメンタムの独立な分布$p(\phi)$によってaugmentされ、同時分布$ p(\theta, \phi | y) = p(\phi) p(\theta |y) $を定義する。つまり、ベクトル$\phi$は補足変数で、パラメータ空間のより速く移動できるように追加しただけである。
　HMCでは事後密度の他に、対数事後密度の勾配が必要である。実務的には勾配は分析的に求める必要がある。数値微分していると計算効率が悪くなるからだ。たいていのモデルで、勾配のベクトルは容易に求まる。

モメンタム分布
　通常、$\phi$には平均0, 共分散行列$M$のMVNを与える。物理学のアナロジーで、$M$を質量行列という。ふつう$M$は対角行列にする(つまり$\phi_j \sim N(0, M_{jj})$)。

HMCイテレーションの3ステップ
　HMCの各イテレーションは次の3ステップからなる。

$\phi$を$N(0, M)$からドローする。
以下を$L$回繰り返す。この1回を「リープフロッグ・ステップ」という。モメンタム$\phi$の更新を2回に分けているのでこう呼ばれている。
1. $\theta$の対数事後密度の勾配を使って$\phi$を更新する。$$ \phi \leftarrow \phi + \frac{1}{2}\epsilon \frac{d \log p(\theta|y)}{d\theta} $$
2. $\phi$を使って$\theta$の位置を更新する。$$ \theta \leftarrow \theta + \epsilon M^{-1} \phi$$ 質量行列$M$が対角なら、結局これは$\theta$の更新をスケーリングしていることになる。なのにわざわざ$\epsilon$を入れている理由は、$M$を固定にして$\epsilon$をチューニングしたほうが便利だからである。
3. ふたたび$\phi$を更新する。$$ \phi \leftarrow \phi + \frac{1}{2}\epsilon \frac{d \log p(\theta|y)}{d\theta} $$
リープフロッグの前の値を$\theta^{t-1}, \phi^{t-1}$, 後の値を$\theta^*, \phi^*$として、以下を求める。$$ r = \frac{p(\theta^* | y) p(\phi^*)}{p(\theta^{t-1}|y)p(\phi^{t-1})} $$
確率 $\mathrm{min}(r, 1)$で$\theta^*$を採用しこれを$\theta^t$とする。不採用の場合は$\theta^{t-1}$を引き続き$\theta^t$とする。[ああそうか。散々リープフロッグしたけど不採用ってこともあるわけね]

棄却されたパラメータと事後密度がゼロのエリア
　HKCはすべて正の目標密度のために設計されている。もしアルゴリズムが事後密度ゼロのところにいっちゃったら、諦めて前の$\theta$に戻るわけである。
　なお、そうではなく、モメンタムの符号を反転させるという手もある。
　[…]

チューニング・パラメータの設定
　HMCには、チューニング・パラメータが3つある。(1)$\phi$の確率分布, (2)スケーリング・ファクター$\epsilon$, (3)イテレーションあたりのリープフロッグ数$L$。
　メトロポリス法一般と同じく、これらは事前に決めておくかランダムに決めるべきである。適応的にやるんならウォームアップでのみ変えるべし。
　どう設定するか？我々は、まずは$\epsilon$は目標分布のスケールの粗い推定値に合わせ、質量行列は単位行列にし、$\epsilon L$が1になるように$L$を決める。

　理論的には、HMCがもっとも効率的になるのは受容率がおよそ65%のときである。受容率が低すぎる場合、リープフロッグが野心的過ぎるわけだから、$\epsilon$を下げ、そのぶん$L$を上げるべし。高すぎる場合はその逆。[へー、なるほどね]

実行中にチューニング・パラメータを変動させる
　[…] $\epsilon$と$L$を実行中にランダムに変えるという手もある。[…]

局所適応的HMC
　難しいHMC問題では、アルゴリズムが事後分布の全体を動けるように、チューニング・パラメータを変動させるのが望ましい。HMCの拡張の試みを紹介しよう。

no-U-turnサンプラー。ステップ数$L$をイテレーションごとに適応的に決める。具体的にいうと、モメンタム変数$\phi$と、イテレーションの開始以来の$\theta$の移動距離とのドット積が負になるまで続ける)。このルールによって、軌跡は行けるところまで伸びることになる。とはいえ、こういうルールを導入するとシミュレーションが目標分布に収束しなくなるかもしれないので、軌跡をいったり戻ったりして、ジャンプのバランスを保つ。また質量行列$M$とステップサイズ$\epsilon$もウォームアップ中は適応的にチューニングする。
リーマン適応。質量行列を、それぞれのステップの対数事後密度の局所曲率に合わせる。

　いずれの拡張も、HMCの問題の全てを解決してくれるわけではない。no-U-turnサンプラーは裾がすごく短い分布や
長い分布では裾に移動しにくいという問題を依然として抱えているし、リーマン適応は高次元では非現実的である。

HMCとGibbsサンプリングの結合
GibbsサンプラーとHMCを併用する方法がふたつある。

変数をブロックにわける。たとえば$J$群の階層モデルがあって、群$j$のパラメータベクトルが$\eta^j$[原文では$\eta^{(j)}$だが略記する], ハイパーパラメータが$\phi$だとしよう。事後分布は次のように分解できる: $$ p(\theta | y) \propto p(\phi) \prod_j p(\eta^j | \phi) p(y^j | \eta^j) $$ ここで、すべてのパラメータのベクトル$\theta = (\eta^1, \ldots, \eta^J, \phi)$についてHMCをやるんじゃなくて、$\eta^1, \ldots, \eta^J, \phi$についてそれぞれ更新したほうが速いだろう。[それぞれについてHMCってこと？]
離散変数の更新。ハミルトニアン・ダイナミクスは連続分布についてのみ定義されているので、混合要素の潜在インジケータのように、離散的空間上で定義されるパラメータについては、Gibbsとかメトロポリスとかスライス・サンプリングとかで更新したい。こういうときは、空間を離散パラメータと連続パラメータに分けて、連続のほうはHMCで更新すればよい。

12.5 単純な階層モデルにおけるハミルトニアン・ダイナミクス
　5.5節の8 schools実験のデータを例に説明しよう。[8つの学校の生徒さんを対象になにかの実験をやったという例だと思う。ここでは個体は生徒じゃなくて学校で、変数として効果の推定値$y_j$とそのSE $\sigma_j$を持っている]
　学校$j$での効果を$\alpha_j$, ハイパーパラメータを$\mu, \tau$とする。10個からなるベクトルを$\theta$とする。

まずは対数事後密度の勾配を求める。それぞれの$j$について $$ \frac{d \log p(\theta|y)}{d \alpha_j} = – \frac{\alpha_j -y_j}{\sigma^2_j} – \frac{\alpha_j – \mu}{\tau^2} $$ また $$ \frac{d \log p(\theta|y)}{d \mu} = – \sum_j \frac{\mu-\alpha_j}{\tau^2} $$ $$ \frac{d \log p(\theta|y)}{d \tau} = – \frac{J}{\tau} + \sum_j \frac{(\mu -\alpha_j)^2}{\tau^3} $$ [ぶひー。こんなん自力で導出できる気がしない]
質量行列を決める。データを眺めた後、それぞれのパラメータにおおまかにスケール15を設定する。つまり$M$はすべての対角に15がはいった対角行列とする。
初期値を決める。4チェーン走らせることにする。各チェーンの10個のパラメータの初期値を$N(0, 15^2)$からドローする。
$\epsilon, L$のチューニング。ここでは、$\epsilon_0 L_0 = 1, L_0 = 10$として、$\epsilon$は$(0, 2\epsilon_0)$の一様分布から、$L$は$[1, 2L_0]$の離散一様分布から、それぞれステップのたびに独立にドローすることにしよう。とりあえず20イテレーションやってみて、プログラムがうまく走ることを確認した。こんどは100イテレーションやってみたら、収束せず、受容率が低かった。$\epsilon_0 = 0.05, L_0 = 20$に再設定したらこんどはうまくいった。そこでこれらのチューニングパラメータで10000イテレーション回した。

log tauへの変換
　パラメータに制約がかかっているとき、アルゴリズムは境界の外に出てしまい、イテレーションの無駄遣いになる。一番簡単な方法は、たとえば$\tau > 0$ならそれを$\log \tau$に変換しちゃうというやり方である。以下の手順となる。

いまや$\theta = (\alpha_1, \ldots, \alpha_8, \mu, \log \tau)$である。
事後分布$p(\theta | y)$にヤコビアン$\tau$を掛け、計算に用いる対数事後密度に$\log \tau$を加える。[ううう… わからない… 学力不足だ…]
対数事後密度の勾配を変える。まず、上記で追加した新しい項を説明する必要がある。次に、勾配の最後の要素の導関数は$\log \tau$に関してではなくて$\tau$に関してになっているので、ヤコビアン$\tau$を掛ける必要がある。$$ \frac{d \log p(\theta | y)}{d \log \tau} = (J-1) + \sum_j \frac{(\mu – \alpha_j)^2}{\tau^2} $$
質量行列を変える。$\log \tau$の質量は1にする。
初期値を変える。$\log \tau$については$N(0, 1)$からドローする。

12.6 Stan: 計算環境の構築
　というわけで、我々は与えられたベイジアン・モデルにHMCを自動的に適応するプログラムStanを開発しました。
　以下、Stanの概要。使い方はAppendix Cをみてね。

データとモデルの入力。モデルの各行では、データとパラメータの対数確率密度を定義する。Stanは勝手に自動微分してくれる。[…] Stanをコールする際にチェーン数などを指定する。初期値は指定してもいいけどStanが勝手に生成してくれる。
ウォームアップにおけるパラメータ・チューニングの設定。個別の事例についてHMCをチューニングするのは難しいが、no-U-turnサンプラーのおかげで、ステップ数$L$を決める必要はない。密度行列$M$, ステップサイズ$\epsilon$はウォームアップ中に、Stanが確率的最適化によって調整してくれる。もっともいつもうまくいくわけではないけれど。[…]
No-U-turnサンプラーによるHMC。複雑な制約があるときは再パラメータ化も重要。Stanは受容確率を監視しているので、収束しないときに中でなにが起きているのかを知るのに役立つだろう。
推論と事後処理。ウォームアップ期のイテレーションを捨てて、$\hat{R}$とか$n_{eff}$とかを求めましょう。

12.7 文献ノート
[略]
————
　やれやれ… 数学の能力がからきしないもので、本章は悔し涙に暮れながら写経した箇所が多い。特に最後の、パラメータを対数変換したらどうなるかという話とか。
　まあいいや、HMCがおおまかになにをやっているのか理解できたような気がするので、よしとしましょう。

読書日記

読んだ本を淡々と記録します

覚え書き: BDA3 Chap.12 計算効率の良いマルコフ連鎖シミュレーション