#_________________________________________
#Step 1: データを取得する
#-----------------------------------------
#作業ディレクトリを設定する
setwd("C:/Users/akita/OneDrive/Documents/52_10 Stats in Ling B 2021/StanPractice")
#ファイルを読み込む
grades = read.csv(file = "Grades.csv", header = T, stringsAsFactors = F, fileEncoding = "Shift-JIS")

#_________________________________________
#Step 2: Stanに読み込ませるデータを作る
#-----------------------------------------
data1 <- list(N = nrow(grades), y = grades[,"英語"])

#_________________________________________
#Step 3: Stanを実行する
#-----------------------------------------
library(rstan)
#-----------------------------------------
#MCMCを実施するときの設定
#-----------------------------------------
WARMUP = 500
ITER   = 1000
#-----------------------------------------
model1 <- stan_model(file = "Model01(intercept_prior).stan")
fit1   <- sampling(model1, data = data1, iter = ITER, warmup = WARMUP)

#_________________________________________
#Step 4: 結果を分析する
#-----------------------------------------
#事後分布への収斂がなされたかを調べ、
#事後分布に関する各種統計量を算出し解釈する
#-----------------------------------------
#事後分布からのサンプルを抽出
#-----------------------------------------
mcmc_sample = rstan::extract(fit1, permuted = FALSE)
#-----------------------------------------
#事後分布からのサンプルの確認
#-----------------------------------------
dim(mcmc_sample)
dimnames(mcmc_sample)
mcmc_sample[1,1,1]
mcmc_sample[1,"chain:1","mu"]
mcmc_sample[500,"chain:1","mu"]
#-----------------------------------------
#( 1 )チェーンの良し悪しの確認
#-----------------------------------------
library(bayesplot)
library(ggplot2)
#-----------------------------------------
#a. Traceplot
mcmc_trace(mcmc_sample, pars = c("mu", "sigma"))
#b. R-hat
Rhat(mcmc_sample[,,"mu"])
Rhat(mcmc_sample[,,"sigma"])
#c. Effective Sample Size
ess_bulk(mcmc_sample[,,"mu"])
#ess_bulk(mcmc_sample[,,"mu"])/ ((ITER-WARMUP)*4)
ess_tail(mcmc_sample[,,"mu"])
ess_bulk(mcmc_sample[,,"sigma"])
ess_tail(mcmc_sample[,,"sigma"])
#-----------------------------------------
#( 2 )事後分布の分析
#-----------------------------------------
mcmc_hist(mcmc_sample, pars = c("mu", "sigma"))
mcmc_dens(mcmc_sample, pars = c("mu", "sigma"))
fit1
summary(fit1)$summary


library(devEMF)
library(extraDistr)
emf(file = "Model01(intercept_prior)muPrior.emf", height = 3,width = 3)
curve(dnorm(x,50, 10), axes = F, xlab = "",ylab = "", xlim = c(70,85))
axis(1)
dev.off()
emf(file = "Model01(intercept_prior)sigmaPrior.emf", height = 3,width = 3)
curve(dhcauchy(x,10^-10, 10), axes = F, xlab = "",ylab = "", xlim = c(11,28))
axis(1, at = seq(12, 28 , by = 4))
dev.off()
emf(file = "Model01(intercept_prior)muPosterior.emf", height = 2,width = 2)
mcmc_dens(mcmc_sample, pars = c("mu"))
dev.off()
emf(file = "Model01(intercept_prior)sigmaPosterior.emf", height = 2,width = 2)
mcmc_dens(mcmc_sample, pars = c("sigma"))
dev.off()