keikoyone77’s blog

ベイズ推論とMCMCのフリーソフト　の勉強中に行った処理等のメモ

環境　macOS Sierra v10.12.6

サポートページのリンク
[特集] ベイズ推論とMCMCのフリーソフト - 岩波データサイエンス

●Rのインストール

　参考ページ　Rjpwiki RjpWiki - RjpWiki

Rのインストールのリンクから R-3.4.4.pkg をダウンロードしてダブルクリック

Rの起動 ::アプリケーション・フォルダにあるRのアイコンをクリック

●JAGSのインストール

　- サポートページ[特集] ベイズ推論とMCMCのフリーソフト - 岩波データサイエンス

　のリンク https://sourceforge.net/projects/mcmc-jags/files/からJAGS-4.3.5.dmgを

　ダウンロード

　- JAGS-4.3.0.dmg をダブルクリックしてJAGS-4.3.0のフォルダを開く

　-Readme.rtfを開いてインストラクションに従いインストール

--- Rのコンソールを開いて >.Platform$pkgType としてバージョン確認

　--- JAGS-4.3.0.mpkgを右クリックして開く　指示に従ってインストール

- {rjags}のインストール :: Rを開いて>install.packages("rjags")

依存対象のcodaもインストールされる

● データのダウンロード

　サポートページにデータ等へのリンクが示されている。

　ただp20のデータが上記に見つからなかったので、データフレームを作成する。

d <- data.frame(pref=c("A","B","C","D","E","F","G","H","I","J","A","B","C","D","E","F","G","H","I","J"), N=c(55,53,55,53,58,55,58,59,56,56,51,49,53,52,55,53,53,57,51,50), mean.Y=c(151.36,151.56,152.22,153.09,153.22,153.31,152.98,153.27,152.67,155.37,157.27,156.83,157.08,156.00,157.24,157.22,157.81,158.95,156.82,161.71), sd.Y=c(2.94,3.07,3.20,2.65,3.07,3.10,2.49,3.08,2.82,3.10,2.98,3.14,3.21,2.64,3.03,3.13,2.45,3.06,2.92,3.21), Age=c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1), X=c(1,1,0,1,1,0,0,0,1,0,1,1,0,1,1,0,0,0,1,0))

作成したデータを保存

save(d, file='FILE_PATH/data.Rdata')

●　JAGS によるベイズ推定
　サポートページの /kubo/jags1 にあるコードをダウンロードして実行しようと思ったが、Macでそのままでは動かなかったので修正。
　https://github.com/iwanami-datascience/vol1/tree/master/kubo/jags1

　1. working directoryを設定

setwd("WORK_DIR_PATH")

2. model.bug.txt をダウンロードしてWORK_DIR_PATHに保存
3. runjags.Rの後半部分を以下のように修正してWORK_DIR_PATHに保存

library(rjags)

load("/Users/yonezawakeiko/work/book/MCMC/data.Rdata") 
N.beta <- 3
list.data <- list(
	mean.Y = d$mean.Y,
	X = d$X, Age = d$Age,
	N = nrow(d), N.beta = N.beta
)
list.inits <- list(beta = rep(0, N.beta), sd = 1)

n.burnin <- 1000
n.chain <- 3
n.thin <- 10
n.iter <- n.thin * 3000

model <- jags.model(
	file = file.model, data = list.data,
	inits = list.inits, n.chain = n.chain
)
update(model, n.burnin) # burn in
post.mcmc.list <- coda.samples(
	model = model,
	#variable.names = c("mu", names(list.inits)),
	variable.names = names(list.inits),
	n.iter = n.iter,
	thin = n.thin
)
print(summary(post.mcmc.list))

4. 保存したファイルを実行

source('runjags2.R')

5. 計算結果が出力される

Iterations = 2010:32000
Thinning interval = 10 
Number of chains = 3 
Sample size per chain = 3000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean     SD Naive SE Time-series SE
beta[1] 152.900 0.4339 0.004574       0.004751
beta[2]   5.658 0.7615 0.008027       0.008530
beta[3]  -1.730 0.8849 0.009327       0.009518
sd        1.358 0.2575 0.002714       0.002714

2. Quantiles for each variable:

            2.5%     25%     50%     75%     97.5%
beta[1] 152.0378 152.618 152.902 153.173 153.78505
beta[2]   4.1515   5.163   5.662   6.161   7.17457
beta[3]  -3.4783  -2.296  -1.734  -1.150   0.02212
sd        0.9605   1.175   1.324   1.500   1.95791