中年男のぼやきというのは果てしないものだが、私の場合、仕事しててつくづく思うのは基礎学力の足りなさである。どこがどう足りないのかは面倒なので書きませんけど、それはもうね、辛いものです。
わざわざデータ・サイエンティスト(っていうんですか？)を目指そうという物好きな方々に、私は言いたい。あほかと。ばかかと。いい年こいて勉強を続けなきゃいけないってのがどういう地獄なのか、ほんとにわかってんのかと。
というわけで、時折発作的に「俺の知識には抜け漏れがある…！」「勉強しなきゃ…！」という焦燥感に駆られて、愛する布団からガバリと跳ね起き、小難しい本をメモをとりながら読んだりする次第である。
以下もそのようにしてとったメモである。
Gelman, A. et al.(2014) Bayesian Data Analysis, Third Edition. Chapter 10, Introduction to Bayesian computation.
この大著の10-13章は計算の話に割かれており、本章はそのイントロにあたる短めの章である。わたくし、苦手なんです、計算の話。
いわく。
ベイジアン計算は目標分布\(p(\theta|y)\)の計算と事後予測分布\(p(\tilde{y}|y)\)の計算を軸に回る。11-13章で計算アルゴリズムについて解説するが、そのまえに本章では積分を近似的に評価する方法について簡単に解説する。
以下では基準化されていない密度\(q(\theta|y)\)をなんらかの方法で簡単に計算でき、かつ \(q(\theta|y)/p(\theta|y)\)は\(y\)のみに依存する定数になると仮定する。たとえばベイズの定理の\(p(\theta)p(y|\theta)\)の部分がそれである。
断りない限り、事後密度は対数にして扱う。
10.1 数値積分
数値積分(求積法ともいう)とは、連続関数の積分を、有限個の点におけるその関数の値を求めることで評価する方法である。シミュレーションによる方法と決定論的方法がある。
任意の関数\(h(\theta)\)の事後期待値は $$ E(h(\theta|y)) = \int h(\theta) p(\theta|y) d \theta $$ と定義される。積分の次元数は\(\theta\)の次元数である。逆に、\(\theta\)の空間を通じた任意の積分を、適切に定義した\(h(\theta)\)の事後期待値とみることもできる。
仮に\(p(\theta|y)\)から\(\theta^s\)をドローできるとしよう。\( \frac{1}{S} \sum_s h(\theta^s) \)を求めれば積分を評価できる。その正確さは、おおまかには\( h(\theta^s) \) のSDでわかる。
シミュレーション法
上述のように、\(p(\theta)\)から無作為標本\(\theta^s\)を得て$$ E(h(\theta)|y) \approx \frac{1}{S} \sum_s h(\theta^s) $$ とする方法。[…]
決定論的方法
点\(\theta^s\)によって表現される空間のボリュームを重み \(w_s\)として、$$ E(h(\theta)|y) \approx \sum_s w_s h(\theta^s) p(\theta^s|y) $$ とする方法[原文では右辺の頭に\(\frac{1}{S}\)がついているけど、erattaによれば不要とのこと]。一般にシミュレーション法より分散が小さいが、高次元の場合は点の選び方が難しい。[…]
10.2 分布近似
分布近似とは、事後分布をなんらかの単純なパラメトリック分布で近似することである。[…]
一部の情報を無視した粗い推定
目標分布からのサンプリングをするまえに、目標分布の位置についてラフな推定値が手に入るとありがたい。その方法は問題によって異なるが、要するに、モデルとデータの一部を無視した単純な問題をつくる。たとえば階層モデルなら、ハイパーパラメータ \(\phi\)をおおざっぱに求めておいて、そこから条件付き事後分布 \(\gamma | \phi, y\)を求める、とか。[…]
10.3 直接シミュレーションと棄却サンプリング
[…]
格子点における計算による直接近似
もっとも単純な離散近似として、\(\theta\)の空間をカバーするように均等に打った点 \(\theta_1, \ldots, \theta_N\)について目標密度\(p(\theta|y)\)を求め、\(\theta_1, \ldots, \theta_N\)と\( \frac{p(\theta_i|y)}{\sum_j p(\theta_j|y)} \) によって連続的な\(p(\theta|y)\)を近似するという方法がある。この方法は高次元だと難しい。
予測分布からのシミュレーション
\(p(\theta|y)\)からのドローさえ手にはいれば、それぞれの\(\theta\)のドローについて\(p(\tilde{y}|\theta)\)からのドローを得れば、事後予測分布が手に入る。
棄却サンプリング
\(p(\theta|y)\)ないし\(q(\theta|y)\)からのドローを手に入れたいとしよう(以下では前者の場合について述べるが後者でも同様)。
棄却サンプリングとはこういう方法である。\(p(\theta|y) > 0\)であるすべての\(\theta\)について、以下の特徴を持つ正の関数 \(g(\theta)\)を定義する。
- \(g\)に比例した確率分布からドローできること。\(g(\theta)\)そのものの積分が1である必要はない。
- 重要度 \(\frac{p(\theta|y)}{g(\theta)}\)が既知の上界\(M\)を持つこと。
でもって、以下の手順でサンプリングする。
- \(g(\theta)\)に比例した確率密度から\(\theta\)をランダムドローして、
- 確率 \(\frac{p(\theta|y)}{M g(\theta)} \)で受理する。つまり、\(Mg(\theta)\)の確率密度曲線と\(p(\theta|y)\)の確率密度曲線を重ね描きした図上で(前者が上にある)、横軸\(\theta\)に縦線をひき、\(Mg(\theta)\)の高さに占める\(p(\theta|y)\)の高さを確率とみなして受理したりしなかったりするわけだ。
受理された\(\theta\)の分布は\(p(\theta|y)\)に従うわけである。この手法の効率は\(g\)の選び方で決まる(\(g \propto p\) だったら\(M\)を適切に選べば100%受理となる)。この手法の良いところは、受理率を監視できるという点である。
10.4 重点サンプリング
\(E(h(\theta)|y)\)に関心があるんだけど[← \(E(h(\theta|y))\)じゃないかなあ…]、\( p(\theta|y) \)からのランダムドローができない、という場合について考える。\(g(\theta)\)からのランダムドローならできるとしよう。$$ E(h(\theta|y)) = \frac{\int h(\theta) q(\theta|y) d \theta}{\int q(\theta|y) d \theta} = \frac{ \int \frac{h(\theta) q(\theta|y)}{g(\theta)} g(\theta) d \theta}{\int \frac{q(\theta|y)}{g(\theta)} g(\theta) d \theta} $$ と書けるから、\(g(\theta)\)から\(\theta^1, \ldots, \theta^S\)をドローして、$$ w(\theta^s) = \frac{q(\theta^s | y)}{g(\theta^s)} $$ として $$ \frac{ \frac{1}{S} \sum_s h(\theta^s) w(\theta^s)}{\frac{1}{S} \sum_s w(\theta^s)}$$ を求めれば推定できたことになる。[なるほど… 数式を手で書き写してみてはじめて腑に落ちた… \(w(\theta^s)\)ってのはつまり\(p(\theta|y)\)と\(g(\theta)\)の比で、これを重みにして\(h(\theta^s)\)を加重平均すりゃいいんじゃんってわけね]
この\(w(\theta^s)\)のことを重点比とか重点ウェイトという。
\(hq/g\)が一定になるような\(g(\theta)\)を選べれば正確な推定ができる。
重点サンプリング推定値の正確性と効率性
なので、ウェイトを $$ \tilde{w}(\theta^s) = \frac{w(\theta^s) S}{\sum_{s’} w(\theta^{s’})} $$ と基準化して、有効標本サイズ $$ S_{eff} = \frac{1}{\sum_s (\tilde{w}(\theta^s))^2 } $$ を調べるとよい。極端に高いウェイトがあると小さくなる。もっとも、ウェイトの分布がほんとに裾を引いているときは、この指標自体がノイジーになるのだけれど。
重点リサンプリング
[メモ省略。errataによれば、この方法よりもPareto smoothed importance samplingという方法がお勧めだそうな]
ベイジアン計算における重点サンプリングの使用
[メモ省略]
10.5 シミュレーションのドローはいくつ必要か？
[…]
多くの場合、\(S=100\)もあれば筋の通った事後分布の要約が手に入る。たとえば\(\theta\)がスカラーで、事後分布を\(N(\mu_\theta, \sigma_\theta^2)\)で近似できるとしよう。事後平均の推定値の正確性は近似的に\(s_\theta / \sqrt{S}\)である。モンテカルロ誤差(\(S\)が有限であることによる誤差)を考慮すれば\(s_\theta \sqrt{1+1/S}\)である[ここの理屈がわからない。きっと前の章を読んでいないからだろう]。\(S=100\)もあれば分母は1.005となるわけで、モンテカルロ誤差の影響はほとんどないことがわかる。
[… 事例による説明。省略]
10.6 計算環境
[略]
10.7 ベイジアン計算のデバッグ
偽データを使ったデバッグ
事後分布についての推論について確信を得るための通常の方法は、目指すモデルのいくつかのバージョンを当てはめて、推論が予期しない形で変わったりしないことを確かめることである。
ある特定のモデルの計算については、異なる初期値による複数のシミュレーションを並列で走らせ、収束をチェックする。
モデルが複雑だったり、予想しない推論が得られた場合には、偽データを使ったデバッグをするとよい。
- 真のパラメータベクトル\(\theta\)として、ある値を決める。厳密にいえばこの値は事前分布からのランダムドローであるべきだが、事前分布が無情報であればなにを選んでもよい。
- 階層モデルなら、上記をハイパーパラメータについて行い、他のパラメータをハイパーパラメータのもとでの事前分布からドローする。
- \(p(y|\theta)\)から偽データセット\(y^{fake}\)をドローする。
- \(p(\theta|y^{fake})\)から\(\theta\)についての事後推論を行う。
- 事後推論を真の\(\theta\)と比べる。
理屈から言えば上記を繰り返さないといけないのだが、デバッグとして１回でも役に立つ。
[…]
デバッグとしてのモデル・チェックと収束チェック
最後に、モデルのチェックと収束のチェックを行う。実務的には、モデルがデータに当てはまっていない場合、単純なプログラムのミスとかである場合も多い。
モデルが当てはまっているのに事後推論が変な場合、バグなのかモデルの根本的な問題なのか、区別しづらいこともある。モデルをもっと簡単にして様子をみるのがよいだろう。モデル構築とは多くの場合、単純なモデルと複雑なモデルの行き来の繰り返しである。
10.8 文献ノート
[略]
以上.
——–
この章はまだいいんですよ。いずれの話もわかりやすい。問題は次の章からだ。