bayesian_bootstrap_experiments.r

#ライブラリ
library(boot)
library(ggplot2)
library(tidyverse)
library(bayesboot)
library(simpleboot)

#関数
#正解率計算
calc_acc <- function(x, true){
  nrow(x[x$low <= true & x$high >= true, ]) / nrow(x)
}

#ガンマ分布でテスト
gamma_dist_bootstrap_test <- function(n_sample, p_shape, p_scale, n_loop){
  #結果格納
  bayes_boot <- 0
  normal_boot <- 0
  
  for(i in 1:n_loop){
    x <- rgamma(n_sample, shape = p_shape, scale = p_scale)
    
    nb <- one.boot(x, mean, 10000)
    bb <- bayesboot(x, mean, R = 10000, R2 = 1000)
    
    if(i == 1){
      tmp <- quantile(bb$V1, probs = c(0.025, 0.975))
      bayes_boot <- data.frame(low = tmp[1], high = tmp[2])
      
      tmp2 <- quantile(nb$t, probs = c(0.025, 0.975))
      normal_boot <- data.frame(low = tmp2[1], high = tmp2[2])
      
    }else{
      tmp <- quantile(bb$V1, probs = c(0.025, 0.975))
      
      bayes_boot <- rbind(bayes_boot, data.frame(low = tmp[1], high = tmp[2]))
      
      tmp2 <- quantile(nb$t, probs = c(0.025, 0.975))
      normal_boot <- rbind(normal_boot, data.frame(low = tmp2[1], high = tmp2[2]))
    }
  }
  
  bayes_boot$type <- "bayes"
  normal_boot$type <- "normal"
  
  return(rbind(bayes_boot, normal_boot))
  
}

#結果
x1 <- gamma_dist_bootstrap_test(5, 2, 1, 100)
x2 <- gamma_dist_bootstrap_test(10, 2, 1, 100)
x3 <- gamma_dist_bootstrap_test(20, 2, 1, 100)
x4 <- gamma_dist_bootstrap_test(50, 2, 1, 100)
x5 <- gamma_dist_bootstrap_test(100, 2, 1, 100)
x6 <- gamma_dist_bootstrap_test(1000, 2, 1, 100)



#精度
acc <- rbind(calc_acc(x1 %>% filter(type == "bayes"), 2), 
      calc_acc(x1 %>% filter(type == "normal"), 2),
      calc_acc(x2 %>% filter(type == "bayes"), 2),
      calc_acc(x2 %>% filter(type == "normal"), 2),
      calc_acc(x3 %>% filter(type == "bayes"), 2),
      calc_acc(x3 %>% filter(type == "normal"), 2),
      calc_acc(x4 %>% filter(type == "bayes"), 2),
      calc_acc(x4 %>% filter(type == "normal"), 2),
      calc_acc(x5 %>% filter(type == "bayes"), 2),
      calc_acc(x5 %>% filter(type == "normal"), 2),
      calc_acc(x6 %>% filter(type == "bayes"), 2),
      calc_acc(x6 %>% filter(type == "normal"), 2)
      )

#可視化
df_viz <- data.frame(acc = acc, 
      n = factor(c(5, 5, 10, 10, 20, 20, 50, 50, 100, 100, 1000, 1000)),
      type = rep(c("bayes", "normal"), 6))

#折れ線グラフ
g <- ggplot(df_viz, aes(x = n, y = acc, fill = type))
g <- g + geom_bar(stat = "identity", position = "dodge")
g