はじめに

その2では、rjagsでMCMCするところまでを紹介しました。今回はMCMCサンプリングの結果をプロットし、推定結果が妥当かどうか確かめます。

推定値を取り出す

MCMCサンプリングの結果から推定値のデータを取り出す方法には、以下のようにします。

# 座標のサンプリング結果の中央値を取り出す
# grep("x", varnames(res)) <= xと名前が付く行名を抽出
res.xy <- summary(res[, grep("x", varnames(res))])$quantiles[, "50%"]
res.x <- res.xy[1:100]
res.y <- res.xy[101:200]
# 移動モードのサンプリング結果の中央値を取り出す
res.bmode <- summary(res[, grep("bmode", varnames(res))])$quantiles[, "50%"]

サンプリング結果の中央値をプロットしたものが図1です。左から真の位置、観測値、推定値となっており、プロットした結果からはうまく推定できたように見えます。ただ、テストデータの移動モードと推定結果があまり一致していないようなのが気になります...

推定値が妥当かどうか確かめる

MCMCの収束判断の際に使用するのRhat値は1.1以下であればモデルが収束したと言えます*1。
以下のコードを実行することで、Rhat値を参照することができます。

gelman.diag(res)

今回の推定結果に対して、Rhat値は全て1.1以下となっており、収束しているようです。

おわりに

rjagsでMCMCして動物の移動モードと真の位置を推定してみました。
しかし、全てRhat値の値は1.1以下となっているにも関わらず、実際のデータの移動モードと推定値はあまり一致せず... 原因がサンプルデータなのか、モデルの定義が正しくないのか、それ以外なのかよく分かりませんでした。
文献を見てみても*2 *3、収束に関して言及されておらず、論文になっているデータに対して適用するとどのような結果だったのか気になる次第です*4。なので、実際のデータで試してみて、どうなるか確かめてみたいなと思いました。

*1:http://takehiko-i-hayashi.hatenablog.com/entry/20100512/1273643029

*2:https://www.int-res.com/articles/meps2007/337/m337p255.pdf

*3:https://www.jstage.jst.go.jp/article/mammalianscience/54/1/54_19/_pdf

*4:前者の論文ではサンプリングの平均値が1.25の時は移動、1.75以上は索餌としていました。1.25から1.75までの間は不明と判定していました。