﻿rm(list = ls())  # ワークスペースをクリア
ep <- 0.001  # ε
p <- 0.6 ; k <- 10000
f <- function(n) {
  set.seed(1)  # 乱数のシード値を固定
  # B(n,p)に従う乱数をk個生成し，Snに代入
  Sn <- rbinom(k, size = n, prob = p)
  return(sum(abs((Sn / n) - p) <= ep) / k)  # P(|Sn/n-p|<=ε)
}
m <- 1000  # 凡例の横軸 = n / h
h <- 1000  # 凡例の横軸最小単位
vf <- numeric(m)
for(i in 1:m) { 
  vf[i] <- f(i * h) 
}
plot(vf, xlim = c(0, m), ylim = c(0, 1), type = "l", 
     xlab = "n / 1000", ylab = "P(|Sn/n-p|<=ep)")
g <- function(m) { 0.95 }  # 定数関数
vg <- Vectorize(g)  # スカラー関数をベクトル化
curve(vg, lty = 2, add = TRUE)  # P=0.95を点線で追加
target_m <- (1.96)^2 * p * (1 - p) / (h * ep^2)
# n = m * h の理論値を点線で追加
lines(x = c(target_m, target_m), y = c(0, 1), type = "l", lty = 4)  
legend("bottomright", legend = c("P(|Sn/n-p|<=ep)", "0.95", 
       paste0("1.96^2*p*(1-p)/(10^3 ep^2)=", target_m)), 
       lty = c(1, 2, 4))  # 凡例の追加