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

setwd(path)   # bayesian_inference.stan のあるディレクトリに移動する．適当なpathを入力
writeLines(readLines("bayesian_inference.stan"))

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)