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

　計算に関わる解説10-13章のうち序盤戦、11章である。MHアルゴリズムまで出てくるが、HMCはまだ出てこない。いやー、計算の話って苦手なもので、どこまでついていけるものか、冷や汗が出ますね。
　原文はマルコフ連鎖シミュレーションという表現を好んで使っているんだけど、字数が多いので、このメモではMCMCと書く。またメトロポリス-ヘイスティングスはMHと書く。

　いわく。

　MCMCは、近似的な分布から\(\theta\)をドローし、目標事後分布\( p(\theta|y)\)をより良く近似するようにそれらを修正する一般的な方法である。サンプリングは系列的に行われ、あるドローが従う分布はその前のドローのみに依存する(マルコフ連鎖からのドロー)。しかしこの手法の成功の鍵はマルコフ性ではなくて、近似分布が徐々に改善し目標分布に収束するという点である。[ああそうか… マルコフ性自体はことの本質ではないということか]
　[図による説明…]
　\(t\)番目のドロー\(\theta^t\)のドロー元の分布を遷移分布 \(T_t(\theta^t | \theta^{t-1})\)とする。マルコフ連鎖が一意な定常分布に収束するように遷移分布を構築したいわけである。
　本章では基本的なMCMC手法であるGibbsサンプラーとMHアルゴリズムを紹介する。

11.1 Gibbsサンプラー
　パラメータ・ベクトルが\(d\)個のサブベクトルに分かれる、すなわち \(\theta = (\theta_1, \ldots, \theta_d)\)だとしよう。
　Gibbサンプラーでは、イテレーション\(t\)は\(d\)ステップからなる。まず、\(\theta\)からサブベクトルを取り出す順番を決める。次に、各サブベクトル\(\theta_j\)を、他のサブベクトルの下での条件つき分布 \( p(\theta_j | \theta^{t-1}_{-j}, y)\)からドローして更新する。

　たとえば、二変量正規分布に従う母集団から観察\((y_1, y_2)\)を得たとしよう。母平均\(\theta = (\theta_1, \theta_2)\)は未知, 共分散行列\( \left(\begin{array}{cc} 1 & \rho \\ \rho & 1 \end{array} \right) \)は既知としよう。\(\theta\)の事前分布を一様分布として、事後分布は $$ \left( \begin{array}{c} \theta_1 \\ \theta_2 \end{array} \right) = N \left( \left( \begin{array}{c} y_1 \\ y_2 \end{array} \right) , \left( \begin{array}{cc} 1 & \rho \\ \rho & 1 \end{array} \right) \right) $$ MVNの性質より、条件付き事後分布は $$ \theta_1 | \theta_2, y \sim N(y_1 + \rho(\theta_2 – y_2), 1-\rho^2) $$ $$ \theta_2 | \theta_1, y \sim N(y_2 + \rho(\theta_1 – y_1), 1-\rho^2) $$ である。Gibbsサンプラーはこれらの2つの分布からサンプリングする。[つまり、\((\theta_1, \theta_2)\)の2次元空間上を上下左右にカクカクと動いていくわけね]

11.2 メトロポリス・アルゴリズムとMHアルゴリズム
　GibbsサンプラーはMHアルゴリズムの特殊ケースとみることができる。

メトロポリス・アルゴリズム
　まず、初期分布\(p_0(\theta) > 0\)をどうにかして用意する。10章で紹介したようなラフな推定値の周りをちらばるような分布でもいいし、なんでもいい。

1. 初期分布\(p_0(\theta)\)から\(\theta^0\)をドローする。
2. \(t = 1, 2, \ldots\)について
   1. 提案分布 \(J_t(\theta^* | \theta^{t-1})\)から提案 \(\theta^*\)をドローする。メトロポリス・アルゴリズムの場合、提案分布は対称的、つまり\(J_t(\theta_a | \theta_b) = J(\theta_b | \theta_a)\)である。
      　たとえば、目標分布が二変量正規分布\(p(\theta|y) = N(\theta|\mathbf{0}, \mathbf{I})\)だとして、提案分布はたとえば\(J_t(\theta^* | \theta^{t-1}) = N(\theta^* | \theta^{t-1}, 0.2^2 \mathbf{I}) \)とする。
   2. 以下の密度比を求める。$$ r = \frac{p(\theta^* | y)}{p(\theta^{t-1} | y)} $$
   3. \(\mathrm{min}(r, 1)\)の確率で\(\theta^*\)を、そうでなければ\(\theta^{t-1}\)を、\(\theta^t\)とする。[たとえば\(p(\theta^* | y)\)が\(p(\theta^{t-1} | y)\)の半分だったら、確率0.5で\(\theta^*\)に乗り換え、もしも\(p(\theta^{t-1} | y)\)より大きかったら問答無用で乗り換えるわけね]

遷移分布\(T_t(\theta^t|\theta^{t-1})\)は\(\theta^{t-1}\)の点密度と提案分布の混合となる。
　[\((\theta_1, \theta_2)\)の2次元空間上をうねうねと動いていくことになるわけね]

最適化との関係
　[略]

メトロポリス・アルゴリズムはなぜうまくいくのか？
　\(\theta^1, \theta^2, \ldots\)が目標分布に収束するという証明は2段階に分かれる。(1)シミュレートされた系列が一意な定常分布に従うマルコフ連鎖であるという証明。(2)その定常分布が目標分布と等しいという証明。
　(1)は、マルコフ連鎖が既約irreducibleであり、非周期的aperiodicであり、一時的transientでないならば、成立する。非周期性と非一次性は、瑣末な例外を除き、任意の正則な分布上のランダム・ウォークにおいて成立する。既約性は、ランダム・ウォークが任意の状態から他の任意の状態へと到達する確率がすべて正である限り成立する。
　(2)について。時点\(t-1\)に目標分布\(p(\theta|y)\)から\(\theta^{t-1}\)をドローしたとしよう。さて、\(p(\theta|y)\)からドローできる任意の点\(\theta_a, \theta_b\)があって\(p(\theta_a|y) \leq p(\theta_b |y)\)だとしよう。\(\theta_b\)から\(\theta_a\)への無条件な遷移確率は$$ p(\theta^t = \theta_a, \theta^{t-1} = \theta_b) = p(\theta_b|y) J_t(\theta_a | \theta_b) \left( \frac{p(\theta_a|y)}{p(\theta_b|y)} \right) = p(\theta_a | y) J_t(\theta_a | \theta_b) $$ である。\(J_t(\cdot|\cdot)\)が対称的なのだから、これは\(\theta_a\)から\(\theta_b\)への遷移確率と同じである。\(\theta^t\)と\(\theta^{t-1}\)の同時分布も対称、つまり周辺分布は同一、つまり\(p(\theta|y)\)は\(\theta\)のマルコフ遷移の定常分布である。
　[ええええ？ わからない。ここで示さないといけないのは、マルコフ連鎖の定常分布\(p(\theta^{t-1}) = p(\theta^t) = \ldots \)が目標分布\(p(\theta|y)\)と等しいということでしょう？ なぜこれがその証明になっているの？ 私がなにかを根本的にわかってないんだろうけど…]

MHアルゴリズム
　MHアルゴリズムはメトロポリス・アルゴリズムを以下の2点について一般化したものである。
　その1。提案分布\(J_t(\theta^* | \theta^{t-1})\)はもはや対称でなくてもよい。
　その2。非対称性を修正するために、密度比を以下のように書き換える。$$ r = \frac{ \frac{p(\theta^*|y)}{J_t(\theta^* | \theta^{t-1})} }{ \frac{p(\theta^{t-1}|y)}{J_t(\theta^{t-1}|\theta^*)} } $$
　収束はメトロポリス・アルゴリズムより早くなる。

提案分布とシミュレーションの効率性の関係
　理想的な提案分布とは、すべての\(\theta\)について\(J(\theta^* | \theta) = p(\theta^* | y)\)であることである。もっともそれができないからこうしているわけですけど。
　提案分布の特性として以下が望まれる。

• 任意の\(\theta\)について\(J(\theta^* | \theta)\)からのドローが容易であること。
• 密度比\(r\)が簡単に計算できること。
• パラメータ空間上を適切な距離だけジャンプできること。
• ジャンプが棄却されすぎないこと。

11.3 Gibbsとメトロポリスをbuilding blocksとして使う
　Gibbsサンプラーは条件付き共役モデル(条件付き事後分布から直接ドローできるモデル)での第一候補である。メトロポリス・アルゴリズムはそうでないときに優れている。
　条件付き事後分布のうちいくつかが直接ドローできて他のはできない、というようなときには組み合わせて使うこともできる。
　条件付きサンプリング・アルゴリズムの一般的な問題のひとつとして、目標分布においてパラメータが高い相関を持っているときに遅くなりうるという点がある。これは再パラメータ化(12章で述べる)とかもっと進んだアルゴリズムとかで解決できる。

GibbsサンプラーはMHアルゴリズムの特殊ケースとして解釈できる
　[略]

近似を伴うGibbsサンプラー
　[略]

11.4 推論と収束評価
　反復シミュレーションであれ、それ以外のベイジアン・シミュレーションであれ、推論の方法は同じで、\(p(\theta|y)\)のドローの集まりを使って事後密度を要約したり、関心ある量を求めたり、\(\tilde{y}\)の事後予測シミュレーションをしたりするわけである。以下では反復シミュレーションに特有な注意点を述べよう。

　難題はふたつある。(1)反復が十分に長くないとき、シミュレーションは目標分布を表さない。また初期の反復は目標分布でなく初期の近似を表す。(2)系列内相関。それ自体には問題ないのだが(つまり収束しているんならドロー順序は無視して良いのだが)、そのぶんシミュレーションは非効率になる。
　対策は3つある。(a)初期値をパラメータ空間上で散らばらせて複数の系列を走らせ、収束を監視する。(b)関心あるすべての量について、系列間の変動と系列内の変動がだいたい同じになるかどうかを監視する。(c)シミュレーションの効率が悪すぎるときはアルゴリズムを変える。

　初期値の影響を減らすため、系列の前半の半分を捨てることが多い。捨てる部分をウォームアップという(バーンインという言い方もあるが我々はあまり好きでない)。半分というのは保守的な選択で、場面によってはもっと少なくても良い。

　近似的な収束に到達したら、ドローのうち毎\(k\)個目だけを抜き出すという方法もある。これをthinという。容量的に助かる。1000個くらい取り出せれば十分である。

　異なる初期値からの複数の系列をつくるのがお勧め。そうすれば系列間のちらばりと系列内のちらばりを比較できる。

　関心ある全てのパラメータ、ならびに関心あるスカラー量を別々に監視しよう。事後密度の対数を監視するのもよい(メトロポリス・アルゴリズムの場合はもう計算してあるはずである)。
　それぞれの系列が定常になり、かつ系列間でmixedであることが求められる。

　収束監視の方法はいろいろある。シンプルな方法は以下である。まずウォームアップを捨てる。次に、それぞれの系列を前半と後半に分ける(5系列なら10本の半系列ができる)。で…
　上記でつくった半系列\(j(=1,\ldots,m)\)のイテレーション\(i(=1,\ldots,n)\)の、関心あるスカラーestimandの値を\(\psi_{ij}\)としよう。系列間分散\(B\)と系列内分散\(W\)を求める。$$ \bar{\psi}_{.j} = \frac{1}{n} \sum_i \psi_{ij} $$ $$ \bar{\psi}_{..} = \frac{1}{m} \sum_j \bar{\psi}_{.j} $$ $$ s^2_j = \frac{1}{n-1} \sum_i (\psi_{ij} – \bar{\psi}_{.j})^2 $$ $$ B = \frac{n}{m-1} \sum_j (\bar{\psi}_{.j} – \bar{\psi}_{..})^2 $$ $$ W = \frac{1}{m} \sum_j s^2_j $$ でもって、\(\psi\)の周辺事後分散\( \mathrm{var}(\psi|y)\)を推定する。$$ \hat{\mathrm{var}}^+(\psi|y) = \frac{n-1}{n} W + \frac{1}{n}B $$ これは初期値の分布が適切にoverdispersedだと仮定しているという点で過大評価になっているのだが、定常性のもとでは不偏である。クラスタ抽出における古典的な分散推定値みたいなものである。
　いっぽう、\(W\)は系列内分散の推定値として過小である。なぜなら目標分布全体をカバーできてはいないからである(\(n\)が無限大に使づくと\(W\)の期待値は\( \mathrm{var}(\psi|y)\)に近づく)。
　反復シミュレーションの収束を監視するために、現在の\(\psi\)の分布が、\(n\)が無限大に使づくとしたらどこまで減りうるか (potential scale reduction) を推定する。すなわち$$ \hat{R} = \sqrt{ \frac{\hat{\mathrm{var}}^+(\psi|y)}{W} } $$ これは\(n\)が無限大に近づくにつれて1へと減少する。
　[…]

11.5 有効ドロー数
　[考えてみると、11.4-11.5は前にメモを取って読んでいた。道理で既視感があると思った。省略する]

11.6 事例: 階層正規モデル
　動物24頭を4条件に割り当て(条件あたり4～8頭)、なにかを測ったとしよう。
　群\(j(=1,\ldots,J)\)の個体\(i(=1,\ldots,n_j)\)の値を\(y_{ij}\)とする。総個体数を\(n\)とする。
　\(y_{ij}\)は\(N(\theta_j, \sigma^2)\)にIIDに従うと仮定する。\(\theta_j\)は\(N(\mu, \tau^2)\)に従うと仮定する。
　\((\mu, \log \sigma, \tau)\)について一様事前分布を仮定する(ただし\(\sigma > 0, \tau > 0\))。これは\(p(\mu, \log \sigma, \log \tau) \propto \tau\)と等価である。\(\log \tau\)について一様事前分布でないのは、そうすると事後分布が非正則になるからである。5章をみよ。[5章についてのメモをこのエントリの末尾に添える]
　すべてのパラメータの同時事後密度は下式となる。 $$ p(\theta, \mu, \log \sigma, \log \tau | y) \propto \tau \prod_j N(\theta_j | \mu, \tau^2) \prod_j \prod_i N(y_{ij} | \theta_j, \sigma^2)$$

　\(\theta_j\)の初期値として、各群の\(y_{ij}\)からランダムに10個の点を取ろう。\(\theta_j\)の平均を\(\mu\)の初期値とする。\(\sigma, \tau\)の初期値はいらない。
　では、Gibbsサンプリングをはじめよう。

1. それぞれの\(\theta_j\)についての条件つき事後分布からドローする。$$ \theta_j | \mu, \sigma, \tau, y \sim N(\hat{\theta}_j, V_{\theta_j} ) $$ ただし $$ \hat{\theta}_j = \frac{ \frac{1}{\tau^2} \mu + \frac{n_j}{\sigma^2} \bar{y}_{.j}}{\frac{1}{\tau^2} + \frac{n_j}{\sigma^2} } $$ $$ V_{\theta_j} = \frac{1}{ \frac{1}{\tau^2} + \frac{n_j}{\sigma^2} } $$ [おちつけ、そんなに複雑な式ではない]
2. \(\mu\)の条件付き事後分布からドローする。$$ \mu | \theta, \sigma, \tau, y \sim N(\hat{\mu}, \frac{\tau^2}{J} ) $$ ただし $$ \hat{\mu} = \frac{1}{J} \sum_j \theta_j $$
3. \(\sigma^2\)の条件つき事後分布からドローする。$$ \sigma^2 | \theta, \mu, \tau, y \sim \mathrm{Inv-}\chi^2(n, \hat{\sigma}^2) $$ ただし $$ \hat{\sigma}^2 = \frac{1}{n} \sum_j \sum_i (y_{ij} – \theta_j)^2 $$
4. \(\tau^2\)の条件付き事後分布からドローする。$$ \tau^2 | \theta, \mu, \sigma, y \sim \mathrm{Inv-}\chi^2(J-1, \hat{\tau}^2) $$ ただし $$ \hat{\tau}^2 = \frac{1}{J-1} \sum_j (\theta_j – \mu)^2$$ \(p(\tau) \propto 1 \)なので自由度は\(J\)ではなく\(J-1\)になる。[5章で出てきた話ね]

　ここでは\(p(\theta, \mu, \sigma, \tau | y)\)についてGibbsサンプラーでシミュレーションしたが、\(p(\mu, \log \sigma, \log \tau | y)\)についてメトロポリスアルゴリズムでシミュレーションしたほうが効率的である。
　メトロポリスアルゴリズムで10個の並列な系列を走らせたところ、500イテレーションで、全パラメータが\(\hat{R} < 1.1\)に到達した。

11.7 文献ノート
[略]

以上.
———————————–
補足. \(\tau\)の一様事前分布
　11.6節の、\(\log \tau\)ではなく\(\tau\)について一様事前分布を想定するというくだり、あいにく頭が悪いもので、なにをいうとるのか、という感じである。5章p.129に戻ってメモをとってみよう。私には難しい話なのでほぼ逐語訳。

　階層モデルの分散パラメータについて、まずは一様事前分布を与える場合を考えてみよう。ただし、分布を定義する尺度については明確でなければならない。分散パラメータのモデリングについてはいろいろな方法が提案されている。自然に見える方法としては、\(\log \tau\)について一様分布を与えるという方法がある(パラメータは正でなければならないので対数をとる)。しかしこの方法では、事後分布は非正則になってしまう。もう一つの方法として、事前分布をあるcompact set上に定義する(たとえば\(A\)をなんらかの大きな値として\([-A, A]\)について定義する)という方法がある。しかしこの方法では、事後分布が\(-A\)に強く依存してしまう。
　問題が生じる理由は、\(\theta, \mu\)を積分した後の周辺尤度\(p(y|\tau)\)が、\(\tau \rightarrow 0\)とともに有限の非ゼロ値に近づくという点にある[←そうならないと困る、という主旨かなあ?]。従って、仮に\(\log \tau\)の事前分布が一様だとすると、the posterior will have infinite mass integrating to the limit \(\log \tau \rightarrow -\infty\). 言い換えると、階層モデルにおいて群レベル分散が0である可能性をデータからは排除できないので、事後分布はこのエリアにan infinite massを持っていてはいけないのである。[ううう。\(\tau = 0\)ってことは十分にあるうるのに、\(p(\log \tau)\)について事前分布を与えると\(\tau = 0\)のときに\(\log \tau\)が-Inf に飛んじゃうから困る、しかし\(\sigma = 0\)はありえないから\(\log \sigma\)について事前分布を与えてもOK、ってこと？]

　もうひとつの選択肢は、\(\tau\)そのものに一様事前分布を与えることである。\(\tau = 0\)付近で has a finite integralであり、よって上記の問題を回避できる。我々は応用場面ではふつうこの無情報分布を使う。しかし、正の値に向けてわずかに合意できないmiscalibrationを持つ。\(\tau \rightarrow \infty\)のエリアではinfinite prior massを持つからである。もし群の数が1とか2とかだったら、これは実際に非正則な事後密度をもたらす。要するに\(\tau = \infty\)となりプーリングが起きない。ある意味でこれは筋の通った挙動である。たった1～2群しかないデータについてどれだけプーリングすべきかを決めることはデータからは困難だと思われるからだ。しかしながら、ベイジアンの観点からいって、問題についてデータが何も言ってくれないときに前もって決定してしまうのは良いやり方ではない。さらに、4～5群のとき、事後分布の右の裾が重いせいで\(\tau\)の過大推定が生じ、個々の\(\theta_j\)の推定において最適とは言えないやり方でのプーリングが起きてしまうという心配がある。
　これらの非正則事前分布密度は、弱情報的な条件つき共役事前分布の極限[?]と解釈することができる。\(\log \tau\)についての一様事前分布は \(p(\tau) \propto \tau^{-1}\)ないし\(p(\tau^2) \propto \tau^{-2}\)と等価であり、自由度0の逆カイ二乗密度の形式を持ち、正則な逆ガンマ事前分布の極限とみることができる。

　\(\tau\)についての一様密度は\(p(\tau^2) \propto \tau^{-1}\)と等価であり、自由度-1の逆カイ二乗密度の形式を持つ。この密度を正則な逆カイ二乗密度の極限と見ることは容易ではないが(そうであるなら自由度は正でなければならない)、\(\tau\)についての半t分布族の、スケールが\(\infty\)であり\(\nu\)が任意の値であるような極限と解釈することができる。