# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add(["Distributions","Plots"])
import Random  # 乱数を扱うためのモジュールをインポート
using Distributions, Plots
ep = 0.001  # ε
p = 0.6
k = 10000
function f(n)
    Random.seed!(1)  # 乱数のシード値を固定
    Sn = rand(Binomial(n, p), k)  # B(n,p)に従う乱数をk個生成
    return count(abs.(Sn ./ n .- p) .<= ep) / k  # P(|Sn/n-p|<=ε)
end
m = 1000  # 凡例の横軸 = n / h
h = 1000  # 凡例の横軸最小単位
vf = zeros(m)
for i in 1:m
    vf[i] = f(i * h)
end
plot(1:m, vf, line=:line, xlab="n/1000", 
     ylab="P(|Sn/n-p|<=ε)", label="P(|Sn/n-p|<=ε)")
hline!([0.95], line=:dash, label="0.95")  # P=0.95を点線で追加
target_m = (1.96^2 * p * (1 - p)) / (h * ep^2)
# n = m * h の理論値を点線で追加
vline!([target_m], line=:dash, label="1.96^2*p*(1-p)/(10^3*ε^2)")