# 心理学データ解析応用/伴走サイトコード ----------------------------------------------------- # Programmed by kosugitti # Licence ; Creative Commons BY-SA license (CC BY-SA) version 4.0 ## Lesson 3. Seven scientist Rstan Version # 準備 ---------------------------------------------------------------------- rm(list = ls()) library(tidyverse) library(bayesplot) library(rstan) options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) ## MAP関数 map_estimation <- function(z) { density(z)$x[which.max(density(z)$y)] } ## MCMCサンプルをデータフレームにする関数 MCMCtoDF <- function(fit) { fit %>% rstan::extract() %>% as.data.frame() %>% tibble::as_tibble() %>% tibble::rowid_to_column("iter") %>% dplyr::select(-lp__) %>% tidyr::pivot_longer(-iter) -> MCMCsample return(MCMCsample) } ## MCMCデータフレームを要約する関数 MCMCsummary <- function(MCMCsample) { MCMCsample %>% dplyr::group_by(name) %>% dplyr::summarise( EAP = mean(value), MED = median(value), MAP = map_estimation(value), SD = sd(value), L95 = quantile(value, prob = 0.025), U95 = quantile(value, prob = 0.975), Rhat = posterior::rhat_basic(value) ) %>% mutate(across(where(is.numeric), ~ num(., digits = 3))) } # データ --------------------------------------------------------------------- x <- c(-27.020, 3.570, 8.191, 9.898, 9.603, 9.945, 10.056) # stanのコンパイル ## rstanの場合 library(rstan) options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) model <- rstan::stan_model("rstan/sevenScientist.stan") fit <- rstan::sampling(model, data = list(Y = x)) # トレースプロット --------------------------------------------------------------- traceplot(fit) # Stanのオプション指定 ------------------------------------------------------------ fit <- rstan::sampling(model, data = list(Y = x), chains = 4, iter = 5000, warmup = 1000 ) # 結果の出力 ------------------------------------------------------------------- fit %>% MCMCtoDF() %>% MCMCsummary()