﻿# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add(["Distributions","Statistics","Plots"])
using Distributions, Statistics, Plots
import Random  # 乱数を扱うためのモジュールをインポート
# Ynの漸近分布である正規分布N(μ,σ^2)のパラメータμ
function parammu(p)
    return 100 * p
end
# Ynの漸近分布である正規分布N(μ,σ^2)のパラメータσ
function paramsigma(n, p)
        return 100 * sqrt(p * (1-p) / n)
end 
# 受験者数m，設問数n，各設問の得点率pを引数として
# Zn，Ynのm人分の結果を配列Z，Yに格納する
function genZY(m, n, p)
    Z = zeros(m)  # 配列Zを初期化
    Y = zeros(m)  # 配列Yを初期化
    mu_p = parammu(p)  # 100p
    sigma_n_p = paramsigma(n,p)  # 100 * sqrt(p(1-p)/n)
    n_times_p = n * p  # np
    sqrt_npq = sqrt(n * p * (1 - p))  # sqrt(np(1-p))  
    Random.seed!(1) # 乱数のシード値を固定
    for i in 1:m
        # Be(p)に従うn個の乱数X1, X2, ..., Xnを生成
        X = rand(Binomial(1, p), n)
        Z[i] = (sum(X) - n_times_p) / sqrt_npq
        Y[i] = mu_p + sigma_n_p * Z[i]
    end
    return (Z, Y)
end
# 受験者数m，設問数n，各設問の得点率pを指定
m = 800000
n = 100
p = 0.6
# Zn，Ynのm人分の結果を配列Z，Yに格納する
(Z, Y) = genZY(m, n, p)
x1 = collect(-4:0.01:4)
vpdf1 = pdf.(Normal(0, 1), x1)
# Znのヒストグラム
histZ = histogram(Z, bins = 40, legend = false, 
                  title = "Histogram of Zn (n=100)", 
                  xlabel = "Zn", normed = true)
# N(0,1)の密度関数を重ねて表示
histZ = plot!(x1, vpdf1, color =:red)
histZ = xlims!(-4, 4)
# Ynのヒストグラム
histY = histogram(Y, bins = 40, legend = false, 
                  title = "Histogram of Yn (n=100)", 
                  xlabel = "Yn", normed = true)
x2 = collect(40:0.01:80)
vpdf2(k) = pdf.(Normal(parammu(p), paramsigma(n, p)), x2)
# N(μ,σ^2)の密度関数を重ねて表示
histY = plot!(x2, vpdf2(n), color =:red)
histY = xlims!(40, 80)
# 図をまとめて表示
plts = plot(histZ, histY, layout=(1, 2), size=(1000, 500))