# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add("StatsBase")
import Random  # 乱数を扱うためのモジュールをインポート
using StatsBase
# 表が出る確率がpのコインを繰り返し投げて
# 初めて表が出るまでに裏が出る回数を記録する関数
function Num_Tail(p)
    num_T = 0  # 初めて表が出るまでに裏が出る回数
    while true
        coin_HT = sample(["表", "裏"], ProbabilityWeights([p, 1-p]))
        if coin_HT == "表"  # 表が出たらwhile文を終了
            break
        end
        num_T += 1  # 裏が出たら加算 
    end
    return num_T  # 「初めて表が出るまでに裏が出る回数」を返す
end
function min_X1X2(n, p)
    X1 = zeros(n)  # 1回目の表が出るまでに裏が出る回数
    X2 = zeros(n)  # 表が1回出た後に2回目の表が出るまでに裏が出る回数
    Random.seed!(1)  # 乱数のシード値を固定
    for i in 1:n
        X1[i] = Num_Tail(p)
        X2[i] = Num_Tail(p)
    end
    # (min(X1[1],X2[1]),min(X1[2],X2[2]),...,min(X1[n],X2[n]))
    return min.(X1, X2)
end
n = 200000  # 試行回数
p0 = 1/6  # 帰無仮説のもとでコインの表が出る確率
p1 = 5/6  # 対立仮説のもとでコインの表が出る確率
v_min_X1X2_p0 = min_X1X2(n, p0)  # 帰無仮説のとき
v_min_X1X2_p1 = min_X1X2(n, p1)  # 対立仮説のとき
println("α(モンテカルロ法)：", mean(v_min_X1X2_p0 .<= 1))
println("α(参考：理論値)：", 1-(1-p0)*(1-p0)*(1-p0)*(1-p0))
println("β(p1) (モンテカルロ法)：", mean(v_min_X1X2_p1 .>= 2))
println("β(p1) (参考：理論値)：", (1-p1)*(1-p1)*(1-p1)*(1-p1))