﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
# 表が出る確率がpのコインを繰り返し投げて
# 初めて表が出るまでに裏が出る回数を記録する関数
def Num_Tail(p):
    num_T = 0  # 初めて表が出るまでに裏が出る回数
    while(True):
        coin_HT = np.random.choice(['表', '裏'], size = 1, p = [p,1-p])
        if coin_HT == '表':  # 表が出たらwhile文を終了
            break
        num_T += 1  # 裏が出たら加算   
    return num_T  # 初めて表が出るまでに裏が出る回数を返す

def min_X1X2(n, p):
    X1 = np.zeros(n)  # 1回目の表が出るまでに裏が出る回数
    X2 = np.zeros(n)  # 表が1回出た後に2回目の表が出るまでに裏が出る回数
    np.random.seed(1)  # 乱数のシード値を固定
    for i in range(1, n+1):
        X1[i-1] = Num_Tail(p)
        X2[i-1] = Num_Tail(p)
    
    # (min(X1[0],X2[0]),min(X1[1],X2[1]),...,min(X1[n-1],X2[n-1]))
    return np.minimum(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)  # 対立仮説のとき
print('α(モンテカルロ法)：', sum(v_min_X1X2_p0 <= 1) / n)
print('α(参考：理論値)：', 1-(1-p0)*(1-p0)*(1-p0)*(1-p0))
print('β(p1) (モンテカルロ法)：', sum(v_min_X1X2_p1 >= 2) / n)
print('β(p1) (参考：理論値)：', (1-p1)*(1-p1)*(1-p1)*(1-p1))