# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add(["DataFrames", "Distributions", "Plots"])
using DataFrames, Distributions, Plots
function hm(m, x)  # 関数hm(x)
    return (x/m)^(m/2) * exp((m - x)/2)
end
# 0<lamb<1に対してhm(x)=lambをみたすx(<m)を求める関数
function hm_inv_left(m, lamb)
    xmin = 0
    xmax = m
    xmid = (xmin + xmax) / 2
    for i in 1:50  # 二分法
        if hm(m, xmid) < lamb
            xmin = xmid
        else
            xmax = xmid
        end
        xmid = (xmin + xmax) / 2
    end
    return xmid
end
# 0<lamb<1に対してhm(x)=lambをみたすx(>m)を求める関数
function hm_inv_right(m, lamb)
    xmin = m
    xmax = 2 * m
    # hm(x)=lambをみたすxが含まれる区間[k*m,(k+1)*m]を探索
    while hm(m, xmax) > lamb
        xmin = xmax
        xmax = xmax + m
    end
    xmid = (xmin + xmax) / 2
    for i in 1:50  # 二分法
        if hm(m, xmid) < lamb
            xmax = xmid
        else
            xmin = xmid
        end
        xmid = (xmin + xmax) / 2
    end
    return xmid
end
#----------------------------------------------------------------
# αとmに対し，hm(xL)=hm(xR)，かつ
# ∫[0,xL] f^(m)_s0^2(x) dx + ∫[xR,∞] f^(m)_s0^2(x) dx = α
# をみたすxL, xRを計算し，一覧表を出力する (m=1, 2, ..., m_max)
#----------------------------------------------------------------
alpha = 0.05  # 有意水準
m_max = 30  # m = 1, 2,..., m_max
m_vec = zeros(m_max)
xleft = zeros(m_max)
xright = zeros(m_max)
for m in 1:m_max
    lambdamin = 0
    lambdamax = 1
    local xL, xR
    for i in 1:50  # 二分法
        lambdamid = (lambdamin + lambdamax) / 2
        xL = hm_inv_left(m, lambdamid)
        xR = hm_inv_right(m, lambdamid)
        IL = cdf(Chisq(m), xL)  # IL = ∫[0,xL] f^(m)_s0^2(x) dx
        IR = ccdf(Chisq(m), xR)  # IR = ∫[xR,∞] f^(m)_s0^2(x) dx
        if IL + IR > alpha
            lambdamax = lambdamid
        else
            lambdamin = lambdamid
        end
    end
    m_vec[m] = m
    xleft[m] = xL
    xright[m] = xR
end
df = DataFrame(:m => m_vec, :xL => xleft, :xR => xright)
println(df)
#------------------------------------------------------------------
# 参考：hm(xL)=hm(xR)=λをみたすxL<m<xRに対し
# I(λ) = ∫[0,xL] f^(m)_s0^2(x) dx + ∫[xR,∞] f^(m)_s0^2(x) dx
# のグラフを描画
#------------------------------------------------------------------
m = 3  # m = n - 1
N = 100  # [0,1]の分割数 (0<=λ<=1)
IL_lambda = zeros(N)
IR_lambda = zeros(N)
I_lambda = zeros(N)
for i in 1:N
    lamb = i / N
    xL = hm_inv_left(m, lamb)
    xR = hm_inv_right(m, lamb)
    IL_lambda[i] = cdf(Chisq(m), xL)  # I_L = ∫[0,xL] f^(m)_s0^2(x) dx
    IR_lambda[i] = ccdf(Chisq(m), xR)  # I_R = ∫[xR,∞] f^(m)_s0^2(x) dx
    I_lambda[i] = IL_lambda[i] + IR_lambda[i]
end
# 積分の値を順にプロットし，直線で結ぶ
x = collect(1/N:1/N:1)  # x=(1/N, 2/N, ..., 1)
plt = plot(x, I_lambda, label = "I(lambda)", linewidth = 1.8, color = :blue,
           title = "m=3", xlabel = "lambda", xlims = (0, 1), ylims = (0, 1))
plot!(plt, x, IL_lambda, label = "IL(lambda)", 
      linewidth = 1.8, linestyle = :dot, color = :black)
plot!(plt, x, IR_lambda, label = "IR(lambda)", 
      linewidth = 1.8, linestyle = :dashdot, color = :black)
display(plt)