﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
import matplotlib.pyplot as plt  # 図の出力に使用するモジュール
p = 1/6  # 1の目が出る確率
# (1) |Sn/n-1/6|<=1/60となる確率をモンテカルロ法で計算
def prb1(n, k):
    np.random.seed(1)  # 乱数のシード値を固定
    # 二項分布B(n,p)に従う乱数をk個生成
    Sn = np.random.binomial(n, p, k)
    prb = sum(abs(Sn/n - 1/6) <= 1/60) / k
    return prb

# モンテカルロ法で求めた確率を表示 (参考：理論値は0.9544)
print(prb1(2000, 10))
print(prb1(2000, 1000))
print(prb1(2000, 100000))

# (2) |Sn/n-1/6|<=0.03となる確率をモンテカルロ法で計算
def prb2(n, k):
    np.random.seed(1)  # 乱数のシード値を固定
    # 二項分布B(n,p)に従う乱数をk個生成
    Sn = np.random.binomial(n, p, k)  
    prb = sum(abs(Sn/n - 1/6) <= 0.03) / k
    return prb

vprb2 = np.zeros(1000)
for i in range(1,1001):
    vprb2[i-1] = prb2(i, 50000)

plt.xlabel('n')
plt.ylabel('P(|Sn/n-1/6|<=0.03)')
# prb2(n)のグラフを出力
plt.plot(vprb2, label = 'P(|Sn/n-1/6|<=0.03)')
# Prb=0.95のグラフを追加
plt.hlines(0.95, 0, 1000, label = '0.95', linestyle = 'dashed')
# n=593のグラフを追加
plt.vlines(593, 0, 1, label = 'n=593', linestyle = 'dashdot')  
plt.legend()  # グラフの凡例を追加
plt.show()