﻿rm(list = ls())  # ワークスペースをクリア
# 表が出る確率がpのコインを繰り返し投げて
# 初めて表が出るまでに裏が出る回数を記録する関数
Num_Tail <- function(p) {
  num_T <- 0  # 初めて表が出るまでに裏が出る回数
  while (TRUE) {
    coin_HT <- sample(c("表", "裏"), 1, replace=TRUE, prob=c(p, 1-p))
    if (coin_HT == "表") {  # 表が出たらwhile文を終了
      break
    }
    num_T <- num_T + 1  # 裏が出たら加算
  }
  return(num_T)  # 初めて表が出るまでに裏が出る回数を返す
}
min_X1X2 <- function(n, p) {
  X1 <- numeric(n)  # 1回目の表が出るまでに裏が出る回数
  X2 <- numeric(n)  # 表が1回出た後に2回目の表が出るまでに裏が出る回数
  set.seed(1)  # 乱数のシード値を固定  
  for (i in 1:n) {
    X1[i] <- Num_Tail(p)
    X2[i] <- Num_Tail(p)
  }
  # pmin(X1,X2) = (min(X1[1],X2[1]),min(X1[2],X2[2]),...,min(X1[n],X2[n]))
  return(pmin(X1, X2))  
}
n <- 200000  # 試行回数
p0 <- 1/6  # 帰無仮説のもとでコインの表が出る確率
p1 <- 5/6  # 対立仮説のもとでコインの表が出る確率
v_min_X1X2_p0 <- min_X1X2(n, p0)  # 帰無仮説のとき
v_min_X1X2_p1 <- min_X1X2(n, p1)  # 対立仮説のとき
sum(v_min_X1X2_p0 <= 1) / n  # α(モンテカルロ法)
1 - (1 - p0) * (1 - p0) * (1 - p0) * (1 - p0)  # α(参考：理論値)
sum(v_min_X1X2_p1 >= 2) / n  # β(p1) (モンテカルロ法)
(1 - p1) * (1 - p1) * (1 - p1) * (1 - p1)  # β(p1) (参考：理論値)