#_________________________________________
#Step 1: データを取得する
#-----------------------------------------
#作業ディレクトリを設定する
setwd("皆さんのディレクトリに変更してください")
#ファイルを読み込む
grades = read.csv(file = "Grades.csv", header = T, stringsAsFactors = F,fileEncoding = "Shift-JIS")
#ファイルの中身を確認する
head(grades)

#_________________________________________
#Step 2: Stanに読み込ませるデータを作る
#-----------------------------------------
#スクリプトの上では、Step 3より前に来るが
#実際に自らスクリプトを書く時は、Stanの
#ファイルを作成してからの方がスムーズ。
#-----------------------------------------
data1 <- list(
  N = nrow(grades)
  , y = grades[,"英語"]
  , x = grades[,"現代文"]
)

#_________________________________________
#Step 3: Stanを実行する
#-----------------------------------------
library(rstan)
#-----------------------------------------
#MCMCを実施するときの設定
#-----------------------------------------
WARMUP = 500
ITER   = 1000
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
#-----------------------------------------
model1 <- stan_model(file = "Model02(SimpleRegression).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","beta0"]
mcmc_sample[500,"chain:1","beta1"]
#-----------------------------------------
#( 1 )チェーンの良し悪しの確認
#-----------------------------------------
library(bayesplot)
library(ggplot2)
#-----------------------------------------
#a. Traceplot
mcmc_trace(mcmc_sample, pars = c("beta0", "beta1"))
#b. R-hat
Rhat(mcmc_sample[,,"beta0"])
Rhat(mcmc_sample[,,"beta1"])
#c. Effective Sample Size
ess_bulk(mcmc_sample[,,"beta0"])
#ess_bulk(mcmc_sample[,,"beta[1]"])/ ((ITER-WARMUP)*4)
ess_tail(mcmc_sample[,,"beta0"])
ess_bulk(mcmc_sample[,,"beta1"])
ess_tail(mcmc_sample[,,"beta1"])
#-----------------------------------------
#( 2 )事後分布の分析
#-----------------------------------------
mcmc_hist(mcmc_sample, pars = c("beta0", "beta1", "sigma"))
mcmc_dens(mcmc_sample, pars = c("beta0", "beta1", "sigma"))
fit1
summary(fit1)$summary