# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add(["Distributions","Plots"])
import Random  # 乱数を扱うためのモジュールをインポート
using Distributions, Plots
p = 1/6  # 1の目が出る確率
# (1) |Sn/n-1/6|<=1/60となる確率をモンテカルロ法で計算
function prb1(n, k)
    Random.seed!(1)  # 乱数のシード値を固定
    # 二項分布B(n,p)に従う乱数をk個生成
    Sn = rand(Binomial(n, p), k)
    prb = sum(abs.(Sn ./ n .- 1/6) .<= 1/60) / k
    return prb
end
# モンテカルロ法で求めた確率を表示 (参考：理論値は0.9544)
println(prb1(2000, 10))
println(prb1(2000, 1000))
println(prb1(2000, 100000))
# (2) |Sn/n-1/6|<=0.03となる確率をモンテカルロ法で計算
function prb2(n, k)
    Random.seed!(1)  # 乱数のシード値を固定
    # 二項分布B(n,p)に従う乱数をk個生成
    Sn = rand(Binomial(n, p), k)
    prb = sum(abs.(Sn ./ n .- 1/6) .<= 0.03) / k
    return prb
end
vprb2 = zeros(1000)
for i in 1:1000
    vprb2[i] = prb2(i, 50000)
end
# prb2(n)のグラフを出力
plot(vprb2, label = "P(|Sn/n-1/6|<=0.03)", 
     xlab = "n", ylab = "P(|Sn/n-1/6|<=0.03)")
# Prb=0.95のグラフを追加
hline!([0.95], linestyle = :dash, label = "0.95")
# n=593のグラフを追加
vline!([593], linestyle = :dot, label = "n=593")