「自然科学の統計学」第5章演習問題2の擬似標準正規乱数と標準正規分布の理論値の$\chi^2$適合度検定をPythonで実装してみた

import numpy as np
import pandas as pd
import itertools
from scipy.stats import norm

# カイ二乗値をLaTeX形式で表示するためにmatplotlibを使用
import matplotlib.pyplot as plt
#from matplotlib import rcParams
#rcParams['text.usetex'] = True
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
from scipy.stats import chi2

df = pd.read_csv('Goodness-of-fit_tests_for_pseudorandom_numbers.csv', index_col = 0)

# 生成した擬似乱数の数
N = 40000

# 標準正規分布を生成するための区間の配列を作成
base_arr = np.arange(-3,3.5, 0.5)

num_arr = [
    [-3],
    [ [base_arr[i], base_arr[i + 1]] for i  in range(len(base_arr) - 1)],
    [3]
]
num_arr =list( itertools.chain.from_iterable(num_arr))

# 確率を格納するための配列
prob_arr = []
for num_range in num_arr:
    if type(num_range) == list:
        # 区間の下限
        lower = num_range[0]
        #　区間の上限
        upper = num_range[1]
        # 区間の下限の上側確率
        lower_prob = norm.sf(x = lower, loc = 0, scale = 1.00)
        # 区間の上限の上側確率
        upper_prob = norm.sf(x = upper, loc = 0, scale = 1.00)
        # (区間の下限の上側確率 - 区間の上限の上側確率) = (区間の下限と上限の間の確率)
        prob = lower_prob - upper_prob
    else:
        if num_range < 0:
            # 下限値-3以下となる確率
            upper_prob = norm.sf(x = num_range, loc = 0, scale = 1.00)
            prob = 1 - upper_prob
        else:
            # 上限値3以上となる確率
            prob = norm.sf(x = num_range, loc = 0, scale = 1.00)
    prob_arr.append(prob)
    
df['確率'] = prob_arr
df['期待度数'] = df['確率'] * N
display(df)

def calc_chi_square(df):
    sum_chi = 0
    for index, row in df.iterrows():
        obj = df.at[index, '度数']
        exp = df.at[index, '期待度数']
        # カイ二乗値を計算
        sum_chi += (obj - exp) ** 2/exp
    return sum_chi

chi_squ = calc_chi_square(df)
txt = r"$\chi^2 = $" + str(chi_squ)

# matplotlibのLaTeXを利用して出力
fig, ax = plt.subplots(figsize = (0,0))
plt.text(0,0,txt, fontsize = 14)
ax = plt.gca()
ax.axis('off')
plt.show()

#  自由度 (度数の数 - 1)
v = len(df) - 1
print(f'自由度: {v}')
alpha = 0.05

chi2_pdf = chi2.ppf(1 - alpha , v)
print(f"自由度{v}の有意水準5%におけるカイ二乗値: {chi2_pdf}")

if chi_squ > chi2_pdf:
    print('擬似乱数は有意水準5%で適合しない')
else:
    print('擬似乱数は有意水準5%で適合する')

# 標準正規分布を生成するための区間の配列を作成
base_arr = np.arange(-3,3.5, 0.5)

num_arr = [
    [-3],
    [ [base_arr[i], base_arr[i + 1]] for i  in range(len(base_arr) - 1)],
    [3]
]
num_arr =list( itertools.chain.from_iterable(num_arr))

    if type(num_range) == list:
        # 区間の下限
        lower = num_range[0]
        #　区間の上限
        upper = num_range[1]
        # 区間の下限の上側確率
        lower_prob = norm.sf(x = lower, loc = 0, scale = 1.00)
        # 区間の上限の上側確率
        upper_prob = norm.sf(x = upper, loc = 0, scale = 1.00)
        # (区間の下限の上側確率 - 区間の上限の上側確率) = (区間の下限と上限の間の確率)
        prob = lower_prob - upper_prob

    else:
        if num_range < 0:
            # 下限値-3以下となる確率
            upper_prob = norm.sf(x = num_range, loc = 0, scale = 1.00)
            prob = 1 - upper_prob
        else:
            # 上限値3以上となる確率
            prob = norm.sf(x = num_range, loc = 0, scale = 1.00)

df['確率'] = prob_arr
df['期待度数'] = df['確率'] * N
display(df)

def calc_chi_square(df):
    sum_chi = 0
    for index, row in df.iterrows():
        obj = df.at[index, '度数']
        exp = df.at[index, '期待度数']
        # カイ二乗値を計算
        sum_chi += (obj - exp) ** 2/exp
    return sum_chi

chi_squ = calc_chi_square(df)
txt = r"$\chi^2 = $" + str(chi_squ)

# matplotlibのLaTeXを利用して出力
fig, ax = plt.subplots(figsize = (0,0))
plt.text(0,0,txt, fontsize = 14)
ax = plt.gca()
ax.axis('off')
plt.show()