9.1 例題：種子数のポアソン回帰（個体差なし）

例題：架空植物20個体のサイズ$x_i$，種子数$y_i$．
目的：平均種子数が体サイズにどう依存しているか調べる．
個体差なし$\to$ポアソン分布で表現できる．
第3章の例題のようにすれば切片と傾きの最尤推定値が得られる．
このGLMをベイズ統計モデル化する方法を考える．

9.2 GLMのベイズモデル化

ベイズモデル化したGLMでも，中核部分はポアソン回帰のGLM．
個体$i$の種子数$y_i$．平均$\lambda_i$のポアソン分布$p(y_i|\lambda_i)$にしたがう．
線形予測子：$\lambda_i = \exp(\beta_1+\beta_2 x_i)$
この$\beta_1$と$\beta_2$の事後分布を推定する．

尤度関数
$L(\beta_1, \beta_2)
= \displaystyle \prod_i p(y_i | \lambda_i)
= \displaystyle \prod_i p(y_i | \beta_1, \beta_2, x_i) $

パラメータ{$\beta_1$, $\beta_2$}がある値を取っているときに$Y$が得られる確率は$p(Y|\beta_1, \beta_2) = L(\beta_1, \beta_2)$．

ベイズモデルでは
（事後分布）$\propto$ （尤度）$\times$（事前分布）
なので，
$p(\beta_1, \beta_2 | Y) \propto p(Y|\beta_1, \beta_2) p(\beta_1) p(\beta_2)$

右辺の$p(\beta_1)$と$p(\beta_2)$は事前分布．これらを適切に指定すればベイズモデル化したGLMとなる．

9.3 無情報事前分布

無情報事前分布：$[-\infty, \infty]$で任意の値を取ることのできる事前分布．
以下の二つがよく使われる¹．

9.4 ベイス統計モデルの事後分布の推定

$\beta$の事前分布を無情報事前分布とする．
例題データにもとづいて$\beta_1$と$\beta_2$の事後分布をMCMCサンプリングで推定する．

WinBUGSというソフトウェアを使う．

9.4.1 ベイズ統計モデルのコーディング

ベイズモデル化したGLMをBUGSコードで表現．

forの中身の1行目：
個体$i$の種子数$y_i$が平均$\lambda_i$のポアソン分布にしたがう．
$\sim$は「左辺は右辺の確率分布にしたがう」という確率論的な関係を表す．

2行目：
個体$i$の平均種子数$\lambda_i$の指定には対数リンク関数を使う．
右辺のような線形予測子を使う．

forの下，beta1とbeta2で平均ゼロで分散の大きい正規分布を指定している．

9.4.2 事後分布推定の準備

Rを使うと便利．

9.4.3 どれだけ長くMCMCサンプリングすればいいのか？

MCMCサンプリングの際にパラメータを設定する必要がある．

n.iter：MCMCサンプリングのステップ数
n.burnin：事後分布推定用に使わない初期区間の長さ
n.chains：MCMCサンプリングの反復回数

MCMCサンプリングを複数回実行することにより，データや統計モデル，MCMCサンプリングなどの妥当性がわかる．

図9.4に示すように，異なるMCMCサンプリングの間の乖離の大小を調べることを収束診断と呼ぶ．
収束診断には$\hat R$指数が使われる．3つ以上のMCMCサンプリングを比較する．

\hat R
= \sqrt{\frac{\hat{\rm var}^+}{W}}

$W$：MCMCサンプリングごとの分散の平均
$\hat{\rm var}^+$：周辺事後分布の分散で，次式で定義される．

\hat{\rm var}^+
= \frac{n-1}{n} W + \frac{1}{W} B

ただし：$B$はMCMCサンプリング間の分散．詳細は参考文献をみないとわからない．

$\hat R$指数の推定値が$1.0$に近いと，$B/W$が小さいことに対応して，MCMCサンプリング間より各MCMCサンプリングでのばらつきの方が大きいことに対応するので，収束しているとみなす．
経験的には，$\hat R > 1.1$のときはMCMCサンプリング間のばらつきが大きく，定常分布・事後分布は推定できないと判断する．

収束しない原因：
・n.iterやn.burninが小さすぎる．
・不適切な統計モデリング
・コーディングの間違い
・データの間違い
・パラメータの初期値が不適切

9.5 MCMCサンプルから事後分布を推定

9.4.2の結果を図示．

図9.5 (A)・(C)：$\beta_1$と$\beta_2$のMCMCサンプリング結果．それぞれ3本あって$\hat R \simeq 1$．
図9.5 (B)・(D)：確率密度関数$p(\beta_1|Y)$と$p(\beta_2|Y)$．周辺事後分布．

MCMCサンプリングにもとづいて，たとえば$\beta_1$の周辺事後分布について調べたい場合は，単に$\beta_2$を無視して$\beta_1$のサンプリング値だけを使って標本統計量などを評価する．

9.5.1 事後分布の統計量

MCMCサンプリング結果の統計量（平均，標準偏差，分位点，収束の良さを表す指数）．

95%信用区間：あるパラメータの値が95%の事後確率で含まれる範囲．
事後分布の95%信用区間を推定したい場合，上下の両端2.5%に十分な個数のほぼ独立なMCMCサンプルが入るようにサンプリング回数を設定する必要がある．

9.6 複数パラメータのMCMCサンプリング

この章の例題の事後分布$p(\beta_1, \beta_2 | Y)$は二つのパラメータ$\beta_1$と$\beta_2$の同時確率分布．
この場合，次のように$\beta_1$と$\beta_2$の値を交互に更新する方法がある²．

9.6.1 ギブスサンプリング：この章の例題の場合

メトロポリス法：新しい値の候補をあげて，それに更新するかどうかを決める方式
ギブスサンプリング：「新しい値の確率分布」を作り，その確率分布のランダムサンプルを新しい値とする

「新しい値の確率分布」：多変量確率分布からひとつの変量を除いて他の変量すべてを定数とする一変量確率分布（全条件付き分布，FCD）

たとえば$\beta_1 = 1.5$，$\beta_2 = 0$とおく．

$\beta_1$のサンプリングから始める．$\beta_2=0$とおいているので，事前分布$p(\beta_2=0)$は定数となり，$\beta_1$のFCDは

p(\beta_1 | Y , \beta_2=0) 
\propto \prod_i \frac{\lambda_i \exp (-\lambda_i)}{y_i!} p (\beta_1)

ここで，$\lambda_i = \exp (\beta_1 + 0)$．
このFCDにしたがう乱数をひとつ発生させて，$\beta_1^{\rm new} = 2.052$が得られた．

次に$\beta_2$．$\beta_1=2.052$として，$\beta_2$のFCDは

p(\beta_2 | Y , \beta_1=2.052) 
\propto \prod_i \frac{\lambda_i \exp (-\lambda_i)}{y_i!} p (\beta_2)

ここで，$\lambda_i = \exp (2.052 + \beta_2 x_i)$．
このFCDにしたがう乱数をひとつ発生させて，$\beta_2^{\rm new} = -0.016$が得られた．

以上により第1MCMCステップは終了．この後も同様にして$\beta_1$と$\beta_2$を交互にサンプリングしていく．十分な数が得られるまで繰り返す．

9.6.2 WinBUGSの挙動はどうなっている？

パラメータ$\beta_\ast$の事前分布$p(\beta_\ast)$が，$\beta_\ast$がある値をとるときにデータ$Y$が得られる確率$p(Y|\beta_\ast)$の共役事前分布であるとき，その事後分布は事前分布と同じ種類の確率分布になる．

事前分布が共役でない場合は，さまざまな数値的な方法によりMCMCサンプリングを行う．

9.7 この章のまとめと参考文献

ベイズ統計モデルでも第4章で説明したようなモデル選択ができればよいが，その方法についてはまだ決着が付いていない．

ToDo: 具体的に問題に対して実装して図9.5や図9.6のように結果を図示してみる．

$p(\beta)=$定数となる事前分布（improper prior）も使われるとある．
・$10^{-9} < \beta < 10^9$の範囲を取る一様分布
・平均ゼロで標準偏差がとても大きい正規分布
ここでは後者を採用． ↩
多変量確率分布を利用して，同時に複数の更新をする方法もあるとか．
・いったん$\beta_2$を定数だとみなして$\beta_1$を変更．
・変更された$\beta_1$を定数とみなして$\beta_2$を変更．
このような部分的なサンプリングが許されるのは，MCMCアルゴリズムによって定められる定常分布を変えないため． ↩