JAGSで正規分布のパラメータの推定

投稿日:2016年9月3日 更新日:

作業療法や信念対立解明アプローチは、難しいデータを扱うことが多いので、ベイズ推定が役立つだろうと考えています。

来週の日本作業療法学会では、ぼくたちはベイズ推定を用いた研究を2つ発表します。

簡単に表すと、ベイズ推定は「事後分布∝尤度×事前分布」と表現できます。

尤度は、手元にあるデータがもっとも得られやすいモデルです。

事前分布は、データを得る前の分布です。

事後分布は、尤度と事前分布に比例する確率分布です。

ベイズ推定は、最尤法のように1つの真値を推定するのではなく、この事後分布(パラメータの確率分布)を推定します。

ベイズ推定は、多変量正規性からの逸脱という問題を解消し、サンプルサイズが小さくてもおおよそ妥当な推定ができ、不適解になりにくく、複雑なモデリングができるというプラグマティックな利点があります。

作業療法や信念対立解明アプローチでは、複雑なモデルで難しいデータで何とかしないといけないことがままあるので、ベイズ推定が広まればよいなぁと思っています。

これを実行するために、ぼくたちの研究室ではMplusを活用していましたが、もっと自由にやるにはStanが良さそうだと思い、コードの写経を行い続ける強化月間を過ごしたわけです。

インプットしたらアウトプットしないと忘れるので、Stanで実データを解析した英語論文を書きつつあります。

いろいろ調べるとJAGSもかなり良さそうだと思いはじめて地道に独習しつつありまする。

作業療法の世界でも広まればよいなぁと思いつつ、ぼく自身の勉強をかねて、めっちゃ簡単なモデルの備忘録を残しておきます。

なお、JAGSコードに変なところがあるかもなので、気づいたら遠慮なくご指摘ください。

よろしくお願いいたします。

さて、慣れないツールの独習は、めっちゃくちゃ簡単なことからはじめるのが鉄板です。

なので、今回のお題は「正規分布のパラメータの推定」というごくごく簡単なものにします。

以下のソフトを使用するため、初めての人は適宜ダウンロードしてください(全部無料です)。

参考にしたサイトは以下です。

充実した内容でとても勉強になります(感謝!)。

まず、サンプルサイズは100、平均は10、標準偏差は5というシミュレーションデータを生成します。

平均と標準偏差を確認すると、以下のようになります。

おおよそ平均は10、標準偏差は5のデータが作れました。

ヒストグラムは以下のような感じ。

NewImage
JAGSはStanと同様に、データをリストで渡す必要があるので、以下のコードを実行しておきます。


次に、JAGSコードの作成です。

JAGSコードは model{} の中にコードを書きます。

さて今回のJAGSコードは以下になります。

muは平均、tauは精度(分散の逆数)、sigmaは標準偏差で事前分布を指定しています。

for文は、正規分布の確率モデルを宣言しています。

またこれは繰り返しの設定であり、for(i in 1:N)と記述すると、iという変数にN回代入する処理が繰り返されるという意味になります。

今回はNは100なので、平均がmu、精度がtauになる正規分布の確率モデルを推定する処理を、100回繰り返すように指定していることになります。

このコードは、txtファイルに書き、ここではnormal1.txtと名づけて保存します。

以下のコードも同じく、お題に対応したものなので、結果を見比べるためにnormal2.txtと名づけて保存しておきましょう。

なお例では、sigmaの事前分布は一様分布ですが、半コーシー分布にするときはt分布(dt)を使うようです。

まずnormal1.txtのモデルの半コーシー分布版。これはnormal3.txtで保存します。

次にnormal2.txtのモデルで、これはnormal4.txtと保存します。

JAGSコードができたので、実際にベイズ推定してみます。

RからJAGSを使うには、rjags、R2jags、runjagsがあります。

今回はR2jagsを使います。

まずR2jagsを呼びだします。

次に、以下のコードを書いて、JAGSでベイズ推定します。

もっと細かくいろいろ設定できるのですが、はじめて触るときはあんまり細かいところまで気にしないようにしたほうがいいように感じました。

1.normal1.txtの結果

平均は10.441、標準偏差は4.629でおおよそ最初の設定通りの結果になりました(当然といえば当然ですけど、何かうれしいです笑)。

Rhatも1.05以下で収束を示しています。

図は、ggmcmcで確認します。

NewImage

NewImage

NewImage

NewImage
他のモデルの結果はざーっと示しておきます。
当然のことながらほとんど同じ結果になっています。
2. normal2.txtの結果

3. normal3.txtの結果
半コーシー分布の書き方は、あまり自信がなかったのですが、結果を見る限りはさしあたりOKっぽいです。

4. normal4.txtの結果
これもsigmaを半コーシー分布にした(つもりの)モデルの結果です。

-研究

Copyright© 京極真の研究室 , 2018 All Rights Reserved.