# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add(["Distributions","Plots"])
using Distributions, Plots
import Random  # 乱数を扱うためのモジュールをインポート
a = 0.5
# 整数 200, 400, 600, ..., 200000を格納した配列
n_values = collect(200:200:200000)
# 大数の強法則：m1
m1 = (exp(a) - exp(-a)) / (2 * a) 
function expXvars(n) 
    Random.seed!(1)  # 乱数のシード値を固定
    # U(-a,a)に従う乱数をn個生成し, Xに代入
    X = rand(Uniform(-a, a), n)
    return sum(exp.(X)) / n
end
v_expXvars = [expXvars(i) for i in n_values]  # スカラー関数をベクトル化
plt_m1 = plot(v_expXvars, xlab="n / 200", ylab="expXvars(n)", 
              label="expXvars", linestyle=:solid, linewidth=0.5)
plt_m1 = plot!(fill(m1, length(n_values)), label="m1", linestyle=:dash)
plt_m1 = plot!(legend=:topright)
# 大数の強法則：m2
m2 = a^2 / 3
function XSQvars(n)
    Random.seed!(1)  # 乱数のシード値を固定
    # U(-a,a)に従う乱数をn個生成し, Xに代入
    X = rand(Uniform(-a, a), n)
    return sum(X .^ 2) / n
end
v_XSQvars = [XSQvars(i) for i in n_values]  # スカラー関数をベクトル化  
plt_m2 = plot(v_XSQvars, xlab="n / 200", ylab="XSQvars(n)", 
              label="XSQvars", linestyle=:solid, linewidth=0.5)
plt_m2 = plot!(fill(m2, length(n_values)), label="m2", linestyle=:dash)
plt_m2 = plot!(legend=:topright)
# 大数の強法則：m3
m3 = 2 * a^3 / (3 * (exp(a) - exp(-a)))
function XSQvars_over_expXvars(n)
    Random.seed!(1)  # 乱数のシード値を固定
    # U(-a,a)に従う乱数をn個生成し, Xに代入
    x = rand(Uniform(-a, a), n)
    return sum(x .^ 2) / sum(exp.(x))
end
# スカラー関数をベクトル化
v_XSQvars_over_expXvars = [XSQvars_over_expXvars(i) for i in n_values]
plt_m3 = plot(v_XSQvars_over_expXvars, xlab="n / 200", 
              ylab="XSQvars_over_expXvars(n)", 
              label="XSQvars_over_expXvars", linestyle=:solid, linewidth=0.5)
plt_m3 = plot!(fill(m3, length(n_values)), label="m3", linestyle=:dash)
plt_m3 = plot!(legend=:topright)
# 図をまとめて表示
plts = plot(plt_m1, plt_m2, plt_m3, layout=(3, 1), size=(600, 600))