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

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

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

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

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

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

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

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

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

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

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

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

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

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

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



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

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

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

充実した内容でとても勉強になります(感謝!)。
まず、サンプルサイズは100、平均は10、標準偏差は5というシミュレーションデータを生成します。
set.seed(123)x <- rnorm(100, 10, 5)
平均と標準偏差を確認すると、以下のようになります。
summary(x)
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## -1.546 7.531 10.310 10.450 13.460 20.940
sd(x)
## [1] 4.564079 
おおよそ平均は10、標準偏差は5のデータが作れました。

ヒストグラムは以下のような感じ。
hist(x)
NewImage
JAGSはStanと同様に、データをリストで渡す必要があるので、以下のコードを実行しておきます。
data <- list(x=x,N=100) 


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

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

さて今回のJAGSコードは以下になります。
model{
#prior
mu ~ dnorm(0,1.0e-4)
tau <- pow(sigma,-2)
sigma ~ dunif(0,1.0e+4)

#likelihood
for (i in 1:N) {
x[i] ~ dnorm(mu, tau)
}
} 
muは平均、tauは精度(分散の逆数)、sigmaは標準偏差で事前分布を指定しています。

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

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

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

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




以下のコードも同じく、お題に対応したものなので、結果を見比べるためにnormal2.txtと名づけて保存しておきましょう。
model {
#prior
mu ~ dnorm(0, 1.0e-4)
sigma ~ dunif(0, 1.0e+4)

#likelihood
for (i in 1:N) {
x[i] ~ dnorm(mu, 1/(sigma*sigma))
}
}
なお例では、sigmaの事前分布は一様分布ですが、半コーシー分布にするときはt分布(dt)を使うようです。
まずnormal1.txtのモデルの半コーシー分布版。これはnormal3.txtで保存します。
model{
#prior
mu ~ dnorm(0,1.0e-4)
tau <- pow(sigma,-2)
sigma ~ dt(0,1,1)T(0,)

#likelihood
for (i in 1:N) {
x[i] ~ dnorm(mu, tau)
}
} 
次にnormal2.txtのモデルで、これはnormal4.txtと保存します。
model {
#prior
mu ~ dnorm(0, 1.0e-4)
sigma ~ dt(0,1,1)T(0,)

#likelihood
for (i in 1:N) {
x[i] ~ dnorm(mu, 1/(sigma*sigma))
}
}
JAGSコードができたので、実際にベイズ推定してみます。

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

今回はR2jagsを使います。

まずR2jagsを呼びだします。
library(R2jags)
次に、以下のコードを書いて、JAGSでベイズ推定します。

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

1.normal1.txtの結果
fit1 <- jags(data=data, model.file="normal1.txt", parameters.to.save=c("mu", "sigma","tau"), n.chains=4, n.iter=2000, n.burnin=1000, DIC=TRUE)
print(fit1)
## Inference for Bugs model at "normal1.txt", fit using jags, 
## 4 chains, each with 2000 iterations (first 1000 discarded) 
## n.sims = 4000 iterations saved 
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat 
## mu 10.441 0.463 9.525 10.122 10.443 10.753 11.332 1.002 
## sigma 4.629 0.340 4.016 4.397 4.606 4.839 5.337 1.001 
## tau 0.047 0.007 0.035 0.043 0.047 0.052 0.062 1.001 
## deviance 588.500 2.072 586.471 587.041 587.854 589.310 593.946 1.002 
## n.eff 
## mu 1800 
## sigma 2900 
## tau 2900 
## deviance 4000 
## 
## For each parameter, n.eff is a crude measure of effective sample size, 
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## 
## DIC info (using the rule, pD = var(deviance)/2) 
## pD = 2.1 and DIC = 590.6 
## DIC is an estimate of expected predictive error (lower deviance is better).
平均は10.441、標準偏差は4.629でおおよそ最初の設定通りの結果になりました(当然といえば当然ですけど、何かうれしいです笑)。

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

図は、ggmcmcで確認します。
library(ggmcmc)
mcmc1 <- as.mcmc(fit1) 
ggs.result1 <- ggs(mcmc1)
ggs_histogram(ggs.result1)
NewImage
ggs_density(ggs.result1)
NewImage
ggs_traceplot(ggs.result1)
NewImage
ggs_running(ggs.result1)
NewImage
他のモデルの結果はざーっと示しておきます。
当然のことながらほとんど同じ結果になっています。
2. normal2.txtの結果
fit2 <- jags(data=data, model.file="normal2.txt", parameters.to.save=c("mu","sigma"), n.chains=4, n.iter=2000, n.burnin=1000, DIC=TRUE)
print(fit2)
## Inference for Bugs model at "normal2.txt", fit using jags, 
## 4 chains, each with 2000 iterations (first 1000 discarded) 
## n.sims = 4000 iterations saved ## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat 
## mu 10.461 0.464 9.512 10.154 10.465 10.785 11.334 1.001 
## sigma 4.631 0.328 4.034 4.405 4.614 4.845 5.319 1.002 
## deviance 588.446 1.980 586.482 587.021 587.872 589.271 593.722 1.003 
## n.eff 
## mu 4000 
## sigma 1600 
## deviance 1700 
## 
## For each parameter, n.eff is a crude measure of effective sample size, ## and Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## 
## DIC info (using the rule, pD = var(deviance)/2) 
## pD = 2.0 and DIC = 590.4 
## DIC is an estimate of expected predictive error (lower deviance is better).
3. normal3.txtの結果
半コーシー分布の書き方は、あまり自信がなかったのですが、結果を見る限りはさしあたりOKっぽいです。
fit3 <- jags(data=data, model.file="normal3.txt", parameters.to.save=c("mu", "sigma","tau"), n.chains=4, n.iter=2000, n.burnin=1000, DIC=TRUE)
print(fit3)
## Inference for Bugs model at "normal3.txt", fit using jags, 
## 4 chains, each with 2000 iterations (first 1000 discarded) 
## n.sims = 4000 iterations saved 
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat 
## mu 10.438 0.451 9.538 10.126 10.447 10.740 11.319 1.002 
## sigma 4.578 0.327 3.975 4.353 4.558 4.785 5.249 1.002 
## tau 0.048 0.007 0.036 0.044 0.048 0.053 0.063 1.002 
## deviance 588.392 1.986 586.478 586.986 587.807 589.142 593.516 1.002 
## n.eff 
## mu 1500 
## sigma 2200 
## tau 2200 
## deviance 2300 
## 
## For each parameter, n.eff is a crude measure of effective sample size, 
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## 
## DIC info (using the rule, pD = var(deviance)/2) 
## pD = 2.0 and DIC = 590.4 
## DIC is an estimate of expected predictive error (lower deviance is better).
4. normal4.txtの結果
これもsigmaを半コーシー分布にした(つもりの)モデルの結果です。
fit4 <- jags(data=data, model.file="normal4.txt", parameters.to.save=c("mu","sigma"), n.chains=4, n.iter=2000, n.burnin=1000, DIC=TRUE)
print(fit4)
## Inference for Bugs model at "normal4.txt", fit using jags, 
## 4 chains, each with 2000 iterations (first 1000 discarded) 
## n.sims = 4000 iterations saved 
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat 
## mu 10.459 0.463 9.526 10.157 10.459 10.773 11.354 1.002 
## sigma 4.581 0.321 4.003 4.355 4.566 4.795 5.244 1.003 
## deviance 588.409 1.949 586.474 587.022 587.834 589.166 593.729 1.001 
## n.eff 
## mu 2200 
## sigma 1100 
## deviance 4000 
## 
## For each parameter, n.eff is a crude measure of effective sample size, 
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## 
## DIC info (using the rule, pD = var(deviance)/2) 
## pD = 1.9 and DIC = 590.3 
## DIC is an estimate of expected predictive error (lower deviance is better).