﻿from scipy.stats import chi2  # カイ二乗分布
import math  # 数学関係の機能をまとめたモジュールをインポート
import numpy as np  # モジュールNumpyをnpという名前でインポート
import matplotlib.pyplot as plt  # グラフ描画のためのモジュール
import pandas as pd  # pandasはデータを効率的に扱うためのモジュール

def hm(m, x):  # 関数hm(x)
    return (x/m)**(m/2) * math.exp((m - x)/2)

# 0<lamb<1に対してhm(x)=lambをみたすx(<m)を求める関数
def hm_inv_left(m, lamb):
    xmin = 0
    xmax = m
    xmid = (xmin + xmax) / 2
    for i in range(1, 50):  # 二分法
        if hm(m, xmid) < lamb:
            xmin = xmid
        else:
            xmax = xmid
        xmid = (xmin + xmax) / 2
    
    return xmid

# 0<lamb<1に対してhm(x)=lambをみたすx(>m)を求める関数
def 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
    xmid = (xmin + xmax) / 2
    for i in range(1, 50):  # 二分法
        if hm(m, xmid) < lamb:
            xmax = xmid
        else:
            xmin = xmid
        xmid = (xmin + xmax) / 2
    
    return xmid

#----------------------------------------------------------------
# αとmに対し，hm(xL)=hm(xR)，かつ
# ∫[0,xL] f^(m)_s0^2(x) dx + ∫[xR,∞] f^(m)_s0^2(x) dx = α
# をみたすxL = hmL_inv_lambda, xR = hmR_inv_lambdaを計算し
# 一覧表を出力する (m=1, 2, ..., m_max)
#----------------------------------------------------------------
alpha = 0.05  # 有意水準
m_max = 30  # m = 1, 2, ..., m_max.
m_vec = np.zeros(m_max)
xleft = np.zeros(m_max)
xright = np.zeros(m_max)
for m in range(1, m_max + 1):
    lambdamin = 0
    lambdamax = 1  
    for i in range(1, 50):  # 二分法
        lambdamid = (lambdamin + lambdamax) / 2
        xL = hm_inv_left(m, lambdamid)
        xR = hm_inv_right(m, lambdamid)
        IL = chi2.cdf(xL, df = m)  # IL = ∫[0,xL] f^(m)_s0^2(x) dx
        IR = chi2.sf(xR, df = m)  # IR = ∫[xR,∞] f^(m)_s0^2(x) dx     
        if IL + IR > alpha:
            lambdamax = lambdamid
        else:
            lambdamin = lambdamid
    
    m_vec[m-1] = m
    xleft[m-1] = xL
    xright[m-1] = xR

data = {'m':m_vec, 'xL':xleft, 'xR':xright}
df = pd.DataFrame(data)  # pandasの機能
print(df.to_string(index=False))
#------------------------------------------------------------------
# 参考：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 = np.zeros(N)
IR_lambda = np.zeros(N)
I_lambda = np.zeros(N)
for i in range(1, N+1):
    lamb = i/N
    xL = hm_inv_left(m, lamb)
    xR = hm_inv_right(m, lamb)
    I_L = chi2.cdf(xL, df = m)  # I_L = ∫[0,xL] f^(m)_s0^2(x) dx
    I_R = chi2.sf(xR, df = m)  # I_R = ∫[xR,∞] f^(m)_s0^2(x) dx
    IL_lambda[i-1] = I_L
    IR_lambda[i-1] = I_R
    I_lambda[i-1] = I_L + I_R

# 積分の値を順にプロットし，直線で結ぶ
x = np.linspace(0, 1, N)  # x=(0, 1/N, 2/N, ..., 1)
plt.title('m=3')
plt.xlabel('lambda')
plt.xlim(0,1)
plt.ylim(0,1)
plt.plot(x, I_lambda, label='I(lambda)', 
         linewidth = 1.8, color="blue")
plt.plot(x, IL_lambda, label='IL(lambda)', 
         linewidth = 1.8, ls=":", color="black")
plt.plot(x, IR_lambda, label='IR(lambda)', 
         linewidth = 1.8, ls="-.", color="black")
plt.legend()
plt.show()