トップ
新規
一覧
単語検索
最終更新
ヘルプ
ログイン
アールメカブ
R_Bayes_jags
をテンプレートにして作成
開始行:
[[Rの備忘録]]
授業用．
#contents
* 豊田秀樹『マルコフ連鎖モンテカルロ法』朝倉書店 [#f02738...
応用編を，BRugs パッケージではなく，[[JAGS:http://www-fis...
p.134の分散分析.
課題は，大教室で受けた講義の期末試験を同じ教室で受ける場...
rjags パッケージの読み込み
library(rjags)
まずデータは次のように用意．これを &ref(jagsData.txt); と...
y Lect Test LT
22 0 0 1
15 0 0 1
20 0 0 1
17 0 0 1
16 0 0 1
5 1 0 0
以下略
このデータを次のように読み込む．
d <- read.table("jagsData.txt", header = T)
初期値 list の準備
inits <-list(mu=1, a=0, b=0, c=0, tau.e=1)
モデルファイル&ref(jagsModel.txt);
model
{
for(n in 1:20){
y[n] ~ dnorm(theta[n], tau.e)
theta[n] <- mu + a*Lect[n] + b*Test[n] + c*LT[n]
}
mu~dnorm(0.0,1.0E-4)
a~dnorm(0.0,1.0E-4)
b~dnorm(0.0,1.0E-4)
c~dnorm(0.0,1.0E-4)
tau.e~dgamma(1.0E-3,1.0E-3); sigma.e <- 1.0/sqrt(tau.e)
}
モデルの定義
m <- jags.model(
file = "jagsModel.txt",
data = d,
inits = list(inits, inits, inits),
nchain = 3
)
//## リストとして読み込んだら data = list.data,
burn-in を行う
update(m, 1000)
MCMC 計算で事後分布からサンプリング，その結果をうけとる
x <- coda.samples(
m,
c("mu", "a", "b", "c", "tau.e"),
thin = 100, n.iter = 20000
)
class(x)
[1] "mcmc.list"
結果を見てみる
summary(x)
Iterations = 1100:21000
Thinning interval = 100
Number of chains = 3
Sample size per chain = 200
1. Empirical mean and standard deviation for each variab...
plus standard error of the mean:
Mean SD Naive SE Time-series SE
a -0.9706 1.51672 0.061920 0.062622
b -1.0420 1.50942 0.061622 0.072638
c 13.0763 1.48156 0.060485 0.058367
mu 5.0077 1.52229 0.062147 0.062471
tau.e 0.1066 0.03818 0.001559 0.001535
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
a -3.79705 -2.03311 -0.9784 0.0934 1.8935
b -4.02259 -2.00015 -1.0789 -0.1027 2.0586
c 10.27393 12.07503 13.0831 14.0117 15.8806
mu 1.87623 3.99396 5.1173 5.9961 7.8435
tau.e 0.04677 0.07727 0.1008 0.1297 0.2024
図を描く
par(ask=TRUE) # パラメータが5個あるので一枚では収まらない
plot(x)
上の実行例ではモデルファイルにデータ長を 20 をハードコー...
list.data <- list(
y = d$y,
Lect = d$Lect,
Test = d$Test,
LT = d$LT,
N = nrow(d)
)
自明だが，モデルファイルは次のように修正
model
{
for(n in 1:N){
y[n] ~ dnorm(theta[n], tau.e)
theta[n] <- mu + a*Lect[n] + b*Test[n] + c*LT[n]
}
mu~dnorm(0.0,1.0E-4)
a~dnorm(0.0,1.0E-4)
b~dnorm(0.0,1.0E-4)
c~dnorm(0.0,1.0E-4)
tau.e~dgamma(1.0E-3,1.0E-3); sigma.e <- 1.0/sqrt(tau.e)
}
このリストをdataに指定して実行
m <- jags.model(
file = "jagsModel.txt",
data = list.data,
inits = list(inits, inits, inits),
nchain = 3
)
* データブロックとモデルブロック [#ze9cc2a7]
JAGSのモデルブロックで，例えばデータブロックで初期化した ...
X[i, j, k, 1:L] ~ dmulti( p[i, j, k, 1:L], n[i, j, k] )
n[i, j, k] <- sum(X[i, j, k, ] )
Unable to resolve parameter n[1,1,1]
(one of its ancestors may be undefined)
と怒られる．2行目の処理はデータブロックに移さないといけな...
[[JAGSマニュアル:http://www-fis.iarc.fr/~martyn/software/...
//-- JAGSでは値を固定する母数の初期値として NA を指定する...
* varは予約語 [#d5840e35]
モデルにうっかり使うとエラーになるので，var0 などとする．
syntax error, unexpected VAR, expecting NAME or FUNC or ...
以下にエラー jags.model...
* DIC [#t798c6d2]
こんな感じで算出すればよい
x <- dic.samples(
m,
n.iter = 10000,
type = "pD"
)
## x$deviance
## x$penalty
x1 <- as.mcmc(x)
summary(x1)
* WinBUGSの[[structure()関数について>R_winBUGS_structure]...
終了行:
[[Rの備忘録]]
授業用．
#contents
* 豊田秀樹『マルコフ連鎖モンテカルロ法』朝倉書店 [#f02738...
応用編を，BRugs パッケージではなく，[[JAGS:http://www-fis...
p.134の分散分析.
課題は，大教室で受けた講義の期末試験を同じ教室で受ける場...
rjags パッケージの読み込み
library(rjags)
まずデータは次のように用意．これを &ref(jagsData.txt); と...
y Lect Test LT
22 0 0 1
15 0 0 1
20 0 0 1
17 0 0 1
16 0 0 1
5 1 0 0
以下略
このデータを次のように読み込む．
d <- read.table("jagsData.txt", header = T)
初期値 list の準備
inits <-list(mu=1, a=0, b=0, c=0, tau.e=1)
モデルファイル&ref(jagsModel.txt);
model
{
for(n in 1:20){
y[n] ~ dnorm(theta[n], tau.e)
theta[n] <- mu + a*Lect[n] + b*Test[n] + c*LT[n]
}
mu~dnorm(0.0,1.0E-4)
a~dnorm(0.0,1.0E-4)
b~dnorm(0.0,1.0E-4)
c~dnorm(0.0,1.0E-4)
tau.e~dgamma(1.0E-3,1.0E-3); sigma.e <- 1.0/sqrt(tau.e)
}
モデルの定義
m <- jags.model(
file = "jagsModel.txt",
data = d,
inits = list(inits, inits, inits),
nchain = 3
)
//## リストとして読み込んだら data = list.data,
burn-in を行う
update(m, 1000)
MCMC 計算で事後分布からサンプリング，その結果をうけとる
x <- coda.samples(
m,
c("mu", "a", "b", "c", "tau.e"),
thin = 100, n.iter = 20000
)
class(x)
[1] "mcmc.list"
結果を見てみる
summary(x)
Iterations = 1100:21000
Thinning interval = 100
Number of chains = 3
Sample size per chain = 200
1. Empirical mean and standard deviation for each variab...
plus standard error of the mean:
Mean SD Naive SE Time-series SE
a -0.9706 1.51672 0.061920 0.062622
b -1.0420 1.50942 0.061622 0.072638
c 13.0763 1.48156 0.060485 0.058367
mu 5.0077 1.52229 0.062147 0.062471
tau.e 0.1066 0.03818 0.001559 0.001535
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
a -3.79705 -2.03311 -0.9784 0.0934 1.8935
b -4.02259 -2.00015 -1.0789 -0.1027 2.0586
c 10.27393 12.07503 13.0831 14.0117 15.8806
mu 1.87623 3.99396 5.1173 5.9961 7.8435
tau.e 0.04677 0.07727 0.1008 0.1297 0.2024
図を描く
par(ask=TRUE) # パラメータが5個あるので一枚では収まらない
plot(x)
上の実行例ではモデルファイルにデータ長を 20 をハードコー...
list.data <- list(
y = d$y,
Lect = d$Lect,
Test = d$Test,
LT = d$LT,
N = nrow(d)
)
自明だが，モデルファイルは次のように修正
model
{
for(n in 1:N){
y[n] ~ dnorm(theta[n], tau.e)
theta[n] <- mu + a*Lect[n] + b*Test[n] + c*LT[n]
}
mu~dnorm(0.0,1.0E-4)
a~dnorm(0.0,1.0E-4)
b~dnorm(0.0,1.0E-4)
c~dnorm(0.0,1.0E-4)
tau.e~dgamma(1.0E-3,1.0E-3); sigma.e <- 1.0/sqrt(tau.e)
}
このリストをdataに指定して実行
m <- jags.model(
file = "jagsModel.txt",
data = list.data,
inits = list(inits, inits, inits),
nchain = 3
)
* データブロックとモデルブロック [#ze9cc2a7]
JAGSのモデルブロックで，例えばデータブロックで初期化した ...
X[i, j, k, 1:L] ~ dmulti( p[i, j, k, 1:L], n[i, j, k] )
n[i, j, k] <- sum(X[i, j, k, ] )
Unable to resolve parameter n[1,1,1]
(one of its ancestors may be undefined)
と怒られる．2行目の処理はデータブロックに移さないといけな...
[[JAGSマニュアル:http://www-fis.iarc.fr/~martyn/software/...
//-- JAGSでは値を固定する母数の初期値として NA を指定する...
* varは予約語 [#d5840e35]
モデルにうっかり使うとエラーになるので，var0 などとする．
syntax error, unexpected VAR, expecting NAME or FUNC or ...
以下にエラー jags.model...
* DIC [#t798c6d2]
こんな感じで算出すればよい
x <- dic.samples(
m,
n.iter = 10000,
type = "pD"
)
## x$deviance
## x$penalty
x1 <- as.mcmc(x)
summary(x1)
* WinBUGSの[[structure()関数について>R_winBUGS_structure]...
ページ名:
