# パラメータ設定
set.seed(123) # 再現性のためのシード固定
N <- 10000    # 欲しいサンプル数

# 目的の関数 f(x) : Beta(2, 5) の密度関数
f <- function(x) dbeta(x, shape1 = 2, shape2 = 5)

# 提案分布 g(x) : Unif(0, 1) の密度関数（常に1）
g <- function(x) dunif(x, min = 0, max = 1)

# 定数 M (f(x) <= M * g(x) を満たす値)
M <- 2.5

# 結果を格納するベクトル
samples <- numeric(N)
count <- 0       # 受理された数
total_tries <- 0 # 試行総数

# 受理・棄却法のループ
while(count < N) {
  total_tries <- total_tries + 1
  
  # 1. 提案分布から候補 Y を生成
  Y <- runif(1, min = 0, max = 1)
  
  # 2. 一様乱数 U を生成
  U <- runif(1)
  
  # 3. 受理条件の判定
  if (U <= f(Y) / (M * g(Y))) {
    count <- count + 1
    samples[count] <- Y
  }
}

# 受理率の計算
accept_rate <- N / total_tries
cat("理論上の受理率 (1/M):", 1/M, "\n")
cat("実際のサンプリング受理率:", accept_rate, "\n")

# ヒストグラムの描画
hist(samples, breaks = 50, probability = TRUE, 
     main = "受理・棄却法によるBeta(2,5)のサンプリング", 
     xlab = "x", col = "lightblue", border = "white")

# 本物のBeta(2,5)の密度関数を重ねる
curve(dbeta(x, 2, 5), add = TRUE, col = "red", lwd = 2)