﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
import matplotlib.pyplot as plt  # 図の出力に使用するモジュール
ep = 0.001  # ε
p = 0.6
k = 10000
def f(n):
    np.random.seed(1)  # 乱数のシード値を固定
    Sn = np.random.binomial(n, p, k)  # B(n,p)に従う乱数をk個生成
    return sum(abs(Sn/n - p) <= ep) / k  # P(|Sn/n-p|<=ε)

m = 1000  # 凡例の横軸 = n / h
h = 1000  # 凡例の横軸最小単位
vf = np.zeros(m)
for i in range(1, m+1):
    vf[i-1] = f(i * h)

plt.xlabel('n/1000')
plt.ylabel('P(|Sn/n-p|<=ε)')
plt.plot(vf, label = 'P(|Sn/n-p|<=ε)')
# P=0.95を追加
plt.axhline(y = 0.95, label = '0.95', linestyle = 'dashed')
target_m = 1.96**2 * p * (1 - p) / (h * ep**2)
# n = m * h の理論値を点線で追加
plt.axvline(x = target_m, 
            label = '1.96^2*p*(1-p)/(10^3ε^2)=%d'% target_m, 
            linestyle = 'dashdot')
plt.legend()  # 凡例の追加
plt.show()