2017年12月6日 更新

# ベイズ統計学を用いた平均値の差の推測

4,818 view 0
setwd(path)   # bayesian_inference.stan のあるディレクトリに移動する．適当なpathを入力

library(rstan)
rstan_options(auto_write = TRUE)

x <- c(49, 66, 69, 55, 54, 72, 51, 76, 40, 62, 66, 51, 59, 68, 66, 57, 53, 66, 58, 57)
y <- c(41, 55, 21, 49, 53, 50, 52, 67, 54, 69, 57, 48, 31, 52, 56, 50, 46, 38, 62, 59)
N_x <- length(x)
N_y <- length(y)
data <-list(N_x = N_x, N_y = N_y, x = x, y = y)

model <- stan_model("bayesian_inference.stan")
fit <- sampling(model, iter = 90000, chains = 4, data = data)
print(fit)

diff <- rstan::extract(fit)$diff (p <- sum(ifelse(diff > 0, 1, 0)) / length(diff)) plot(density(diff), xlim=c(-1, 20), ylim=c(0, 0.15), xlab="Difference between mean values of class X and class Y.", ylab="Density", col="black", main="", lwd=2) legend("topleft", legend="X組とY組の平均値の差の分布", lty=1, col="black", cex=0.8) mu_x <- rstan::extract(fit)$mu_x
mu_y <- rstan::extract(fit)\$mu_y

plot(density(mu_x),
xlim=c(40, 70), ylim=c(0, 0.2),
xlab="Mean values in class X and class Y", ylab="Density",
col="red", main="", lwd=2)

par(new=T)

plot(density(mu_y),
xlim=c(40, 70), ylim=c(0, 0.2),
xlab="", ylab="",
col="blue", main="", lwd=2)

legend("topleft",
legend=c("X組の平均値の分布", "Y組の平均値の分布"),
lty=c(1,1),
col=c("red", "blue"),
cex=0.8)


bayesian_inference.R
16 件