R: Stanをためす [統計]
久保本 10章の階層ベイズモデル。RStanだとなぜかRがおちるので、コマンドラインでコンパイルして実行する。
Stanのインストールは、Stan: Command Quick Startを参照。
kubo10.stanとして以下の内容を保存。
// Kubo Book Chapter 10
data {
int<lower=0> N; // sample size
int<lower=0> Y[N]; // response variable
}
parameters {
real beta;
real r[N];
real<lower=0> sigma;
}
transformed parameters {
real q[N];
for (i in 1:N) {
q[i] <- inv_logit(beta + r[i]); // 生存確率
}
}
model {
for (i in 1:N) {
Y[i] ~ binomial(8, q[i]); // 二項分布
}
beta ~ normal(0, 100); // 無情報事前分布
r ~ normal(0, sigma); // 階層事前分布
sigma ~ uniform(0, 1.0e+4); // 無情報事前分布
}
kubo10.stanのあるディレクトリで以下を実行。
library(coda)
library(parallel)
program <- "kubo10"
# read data
d <- read.csv(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/hbm/data7a.csv"))
N <- nrow(d)
Y <- d$y
# dump data
data.file <- paste(program, "Rdump", sep = ".")
dump(c("N", "Y"), data.file)
# Stan
stan.home <- ("~/Documents/src/stan-src-1.0.0")
cur.dir <- getwd()
setwd(stan.home)
system(paste("make ", cur.dir, "/", program, sep = ""))
setwd(cur.dir)
run <- function(chain, seed) {
cmd <- paste("./", program,
" --data=", data.file,
" --iter=10100 --warmup=100 --thin=10",
" --seed=", seed,
" --chain_id=", chain,
" --samples=", program, ".chain", chain, ".csv", sep = "")
system(cmd)
s <- read.csv(paste(program, ".chain", chain, ".csv", sep = ""), comment.char = "#")
mcmc(s, start = 101, thin = 10)
}
seeds <- c(1234, 12345, 123456)
samples <- as.mcmc.list(mclapply(1:3, function(i) run(i, seeds[i]), mc.cores = 4))
plot(samples[, c("beta", "sigma")])
summary(samples[, c("beta", "sigma")])
結果
> plot(samples[, c("beta", "sigma")])
> summary(samples[, c("beta", "sigma")])
Iterations = 101:10091
Thinning interval = 10
Number of chains = 3
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
beta 0.05743 0.3324 0.006069 0.009180
sigma 3.02590 0.3752 0.006850 0.008213
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
beta -0.5862 -0.1627 0.06333 0.2773 0.7223
sigma 2.3971 2.7603 2.99708 3.2526 3.8486
こんにちは。
久保本 10章をStanで勉強している者です。
こちらのコードを参考にしながら、動かしていて1つ分からないことがあり、お尋ねします。
// Kubo Book Chapter 10のmodelブロック
①
model {
for (i in 1:N) {
Y[i] ~ binomial(8, q[i]); // 二項分布
}
beta ~ normal(0, 100); // 無情報事前分布
r ~ normal(0, sigma); // 階層事前分布
sigma ~ uniform(0, 1.0e+4); // 無情報事前分布
}
これは
②
model {
for (i in 1:N) {
Y[i] ~ binomial(8, q[i]); // 二項分布
r[i] ~ normal(0, sigma); // 階層事前分布
}
beta ~ normal(0, 100); // 無情報事前分布
sigma ~ uniform(0, 1.0e+4); // 無情報事前分布
}
または
③
model {
for (i in 1:N) {
Y[i] ~ binomial(8, q[i]); // 二項分布
}
for (i in 1:N) {
r[i] ~ normal(0, sigma); // 階層事前分布
}
beta ~ normal(0, 100); // 無情報事前分布
sigma ~ uniform(0, 1.0e+4); // 無情報事前分布
}
ではいけないでしょうか。
久保本のp.229のBUGSコードによるとr[i]となっている点が気になりました。久保本に忠実に再現すると③のような気がしますが、②でもよいよいな。ただし、①と②と③は、rとbetaの推定値(mean)が異なりました。
質問の要点は
・r~は、rとするべきか、r[i]とするべきか
・r[i]は、y~binomのfor文に入れてよいか(②)、別なfor文に収めた方がよいか
ということです。
要領を得ない質問で恐縮ですが、ヒントなどいただけると幸いです。
よろしくお願いいたします。
by TK (2017-10-17 16:01)
このコードですが、今ならこうします。
data {
int<lower=0> N; // sample size
int<lower=0> Y[N]; // response variable
}
parameters {
real beta;
vector[N] r;
real<lower=0> sigma;
}
transformed parameters {
vector[N] logit_q = beta + r; // 生存確率のロジット
}
model {
Y ~ binomial_logit(8, logit_q); // 二項分布
beta ~ normal(0, 100); // 無情報事前分布
r ~ normal(0, sigma); // 階層事前分布
sigma ~ uniform(0, 1.0e+4); // 無情報事前分布
}
Stanのサンプリングはforループではなくベクトルにした方が高効率ですので、できるだけベクトルにすることが推奨されています。詳しくはStanのマニュアルあるいは松浦健太郎『StanとRでベイズ統計モデリング』をご覧ください。
また、質問とは直接関係ありませんが、パラメーターを逆ロジット変換するのではなく、bernoulli_logit分布であてはめた方が計算が安定するとのことです。
あと、事後平均の微妙な違いは乱数のせいですので気にしなくて大丈夫です。
by hiroki (2017-10-18 05:48)
おはようございます。
ご回答ありがとうございます。
forループではなく、ベクトル"vector[N] logit_q = beta + r;"ですね! あと、bernoulli_logitについても教えていただいてありがとうございます。inv_logitの方を使っていたので試してみます。事後平均の微妙な違いについても分かりました。
突然の質問に丁寧にお答えいただいて感謝です。
ありがとうございました!
by TK (2017-10-18 09:22)