「自然科学の統計学」第11章演習問題3-M系列による乱数発生プログラムの作成をPythonで実装してみた

import math
import numpy as np
import pandas as pd
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2

#  パラメータの設定
p = 127
q = 63

# 2進数の長さ = 最大数の定義
digits = 2 ** 3

# 生成する乱数のサンプルサイズ
N = 10000

# 一様分布との比較のためのビンのサイズ
bin_ = 0.1

# ランダムなバイナリa_nをa_pまで作成
ans = random.choices([0,1], k = p)

#  a_p以降はa_n = a_{n-q} + a_{n-p}で計算
for n in range(p, p * digits):
    an = ans[n - p] + ans[n - q]
    # 1+1=0で擬似的に排他的論理和を計算
    if an == 2:
        an = 0
    ans.append(an)

# 2進数の文字列同士から排他的論理和を計算する関数
def bitwise_addition_no_carry(bin1, bin2):
    # 長さを合わせるためにゼロ埋めする
    max_len = max(len(bin1), len(bin2))
    bin1 = bin1.zfill(max_len)
    bin2 = bin2.zfill(max_len)
    
    # 結果を保存するためのリスト
    result = []
    
    # 各ビットごとにXORを計算する
    for b1, b2 in zip(bin1, bin2):
        sum_bit = int(b1) ^ int(b2)  # XORを計算
        result.append(str(sum_bit))
    
    # 結果のリストを文字列に変換する
    return ''.join(result)
    
# 2進数を格納する配列
Xns_binary = []

# 10進数を格納する配列
Xns_decimal = []

# [0, 1)の一様乱数を格納する配列
pns = []

for n in range(p):
    # a_nからX_nを作成
    Xn = ans[n * digits: (n + 1) * digits]
    
    # 配列を文字列に変換
    Xn = ''.join(map(str, Xn))
    Xns_binary.append(Xn)
    
    # 2進数を10進数に変換
    Xn_decimal = int(Xn, 2)
    
    # 10進数を[0, 1)の一様乱数に変換
    pn = Xn_decimal / (2**8)
    pns.append(pn)
    Xns_decimal.append(Xn_decimal)
    
for n in range(p, N):
    # X_n = X_{n-p} + X_{n - q}(排他的論理和)
    Xn =bitwise_addition_no_carry(Xns_binary[n - p], Xns_binary[n - q])
    Xns_binary.append(Xn)
    
    # 2進数を10進数に計算
    Xn_decimal = int(Xn, 2)
    
    # 10進数を[0, 1)の一様乱数に変換
    pn = Xn_decimal / (2 ** digits - 1)
    pns.append(pn)
    Xns_decimal.append(Xn_decimal)
    
# 適合度検定のため、[0, 1)を等間隔でビニング
range_arr = np.arange(0, 1, bin_)
print(f'range_arr: {range_arr}')

df = pd.DataFrame(columns = ['p', '度数', '期待度数'])
df['p'] = range_arr
for i, f in enumerate(range_arr):
    # ビンの幅の中にある乱数の個数をカウント
    count = np.sum((f <= pns) & (pns < f + bin_))
    df.at[i, '度数'] = count
    # 期待度数 = 総数 / ビニング区間の数
    df.at[i, '期待度数'] = N / len(range_arr)
    df.at[i, 'p'] = f'{str(round(f,3))} ~ {str(round(f + bin_, 3))}'
    
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%で適合する')
print()
print()

import math
import numpy as np
import pandas as pd
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2

#  パラメータの設定
p = 127
q = 63

# 2進数の長さ = 最大数の定義
digits = 2 ** 3

# 生成する乱数のサンプルサイズ
N = 10000

# 一様分布との比較のためのビンのサイズ
bin_ = 0.1

# ランダムなバイナリa_nをa_pまで作成
ans = random.choices([0,1], k = p)

#  a_p以降はa_n = a_{n-q} + a_{n-p}で計算
for n in range(p, p * digits):
    an = ans[n - p] + ans[n - q]
    # 1+1=0で擬似的に排他的論理和を計算
    if an == 2:
        an = 0
    ans.append(an)

# 2進数を格納する配列
Xns_binary = []

# 10進数を格納する配列
Xns_decimal = []

# [0, 1)の一様乱数を格納する配列
pns = []

for n in range(p):
    # a_nからX_nを作成
    Xn = ans[n * digits: (n + 1) * digits]
    
    # 配列を文字列に変換
    Xn = ''.join(map(str, Xn))
    Xns_binary.append(Xn)
    
    # 2進数を10進数に変換
    Xn_decimal = int(Xn, 2)
    
    # 10進数を[0, 1)の一様乱数に変換
    pn = Xn_decimal / (2**8)
    pns.append(pn)
    Xns_decimal.append(Xn_decimal)

# 2進数の文字列同士から排他的論理和を計算する関数
def bitwise_addition_no_carry(bin1, bin2):
    # 長さを合わせるためにゼロ埋めする
    max_len = max(len(bin1), len(bin2))
    bin1 = bin1.zfill(max_len)
    bin2 = bin2.zfill(max_len)
    
    # 結果を保存するためのリスト
    result = []
    
    # 各ビットごとにXORを計算する
    for b1, b2 in zip(bin1, bin2):
        sum_bit = int(b1) ^ int(b2)  # XORを計算
        result.append(str(sum_bit))
    
    # 結果のリストを文字列に変換する
    return ''.join(result)

....

for n in range(p, N):
    # X_n = X_{n-p} + X_{n - q}(排他的論理和)
    Xn =bitwise_addition_no_carry(Xns_binary[n - p], Xns_binary[n - q])
    Xns_binary.append(Xn)
    
    # 2進数を10進数に計算
    Xn_decimal = int(Xn, 2)
    
    # 10進数を[0, 1)の一様乱数に変換
    pn = Xn_decimal / (2 ** digits - 1)
    pns.append(pn)
    Xns_decimal.append(Xn_decimal)

# 適合度検定のため、[0, 1)を等間隔でビニング
range_arr = np.arange(0, 1, bin_)
print(f'range_arr: {range_arr}')

df = pd.DataFrame(columns = ['p', '度数', '期待度数'])
df['p'] = range_arr
for i, f in enumerate(range_arr):
    # ビンの幅の中にある乱数の個数をカウント
    count = np.sum((f <= pns) & (pns < f + bin_))
    df.at[i, '度数'] = count
    # 期待度数 = 総数 / ビニング区間の数
    df.at[i, '期待度数'] = N / len(range_arr)
    df.at[i, 'p'] = f'{str(round(f,3))} ~ {str(round(f + bin_, 3))}'

# カイ二乗値を計算する関数
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%で適合する')
print()
print()