﻿rm(list = ls())  # ワークスペースをクリア
p <- 1/6  # 1の目が出る確率
# (1) |Sn/n-1/6|<=1/60となる確率をモンテカルロ法で計算
prb1 <- function(n, k) {
  set.seed(1)  # 乱数のシード値を固定
  Sn <- rbinom(k, n, p)  # B(n,p)に従う乱数をk個生成し，Snに代入
  prb <- sum(abs((Sn / n) - (1 / 6)) <= 1 / 60) / k  
  return(prb)
}
# モンテカルロ法で求めた確率を表示 (参考：理論値は0.9544)
prb1(2000, 10) ; prb1(2000, 1000) ; prb1(2000, 100000)
# (2) |Sn/n-1/6|<=0.03となる確率をモンテカルロ法で計算
prb2 <- function(n, k) {
  set.seed(1)  # 乱数のシード値を固定
  Sn <- rbinom(k, n, p)  # B(n,p)に従う乱数をk個生成し，Snに代入
  prb <- sum(abs((Sn / n) - (1 / 6)) <= 0.03) / k
  return(prb)
}
vprb2 <- numeric(1000)
for(i in 1:1000) {
  vprb2[i] <- prb2(i, 50000)
}
# prb2(n)のグラフを出力
plot(vprb2, type = "l", xlab = "n", ylab = "P(|Sn/n-1/6|<=0.03)")  
# Prb=0.95のグラフを追加
lines(x = c(0,1000), y = c(0.95, 0.95), lty = 2)  
# n=593のグラフを追加
lines(x = c(593, 593), y = c(0,1), lty = 3)  
# グラフの凡例を追加
legend("bottomright", 
       legend = c("P(|Sn/n-1/6|<=0.03)", "0.95", "n=593"), lty = 1:3) 