MCMCサンプリング練習

概要

・前回の最尤推定の練習からの続き
・二項分布のパラメータ最尤推定
・参考文献(後述)の二項分布のパラメータ「ふらふら最尤推定」
・同じく参考文献の二項分布のパラメータをMCMCサンプリング(メトロポリス法)

まずは普通に最尤推定

ということで前回は正規分布のパラメータ推定をoptimメソッドを用いて行いましたが,
正規分布はパラメータが2つなので,もっと単純にということで二項分布のパラメータ推定

今回は最適化アルゴリズムとして,準ニュートン法を使用しています.
実行すると

推定結果は0.44ということで,真値(0.45)が推定されてそうですというのが最尤推定

ふらふら最尤推定

参考文献P.176からの抜粋ですが,ふらふら最尤推定とは

1. パラメータ\(q\)の初期値を選ぶ
2. \(q\)を増やすか減らすかをランダムに決める(新しく選んだ\(q\)の値を\(q^新\)としましょう)
3. \(q^新\)において尤度が大きくなる(あてはまりが良くなる)なら\(q\)の値を\(q^新\)に変更する

この手順を\(q\)の値が変化しなくなるまで実行するものとなります.

Rでの再帰の練習もかねてfor文は使わずに書いています(後々面倒なことになるのですが…)
また,最尤推定の時に使った対数尤度を計算する関数Lも使いまわしています.
これを実行すると,

131229_1

というような感じに,\(q\)が徐々に収束していきます.
最初の条件分岐のところでランダムにどちらに進むか(尤度が改善される方向か,されない方向か)を決めるので,
ステップを進むにつれて段々と\(q\)の値が更新されていきます.

MCMCサンプリング

上記のふらふら最尤推定のアルゴリズムを少し改変したのがメトロポリス法によるMCMCサンプリングとなります(らしい)
上述したアルゴリズムに4番目が追加されて以下のようになります.

1. パラメータ\(q\)の初期値を選ぶ
2. \(q\)を増やすか減らすかをランダムに決める(新しく選んだ\(q\)の値を\(q^新\)としましょう)
3. \(q^新\)において尤度が大きくなる(あてはまりが良くなる)なら\(q\)の値を\(q^新\)に変更する
4. \(q^新\)で尤度が小さくなる(あてはまりが悪くなる)場合であっても,確率\(r\)で\(q\)の値を変更する

ここで確率\(r\)は尤度比
$$r = \frac{L\left( q^新 \right)}{L\left( q \right)}$$
で与えられます.

これをRで実装すると

となります.尤度が改善されない方向に進もうとした場合,尤度比(likelihood ratio)を計算して,
尤度比の確率で,尤度が改善されない方向にも\(q\)を更新していきます.

131229_2

初期値0.05から始めた場合も何やら0.45周辺でサンプリングされてるっぽいですね!!
試しに初期値を色々動かした場合にも同じようにサンプリングされるかというと,

131229_3

ということで大体100ステップ以降くらいからはqの定常分布になってるっぽいですね!何かすごい!

余談

で,

調子に乗って10万ステップくらいで計算しようとしたら

どうやら再帰で計算しているため,関数のネストが深くなりすぎて処理できなくなっている様子.
教えてジェネラル!

参考

Amazon.co.jp: データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学): 久保 拓弥: 本
http://www.amazon.co.jp/dp/400006973X