前回は実数の一様乱数を生成して、そのばらつき具合を統計的に検定しました。実数の一様乱数は整数の一様乱数から生成することができ、拙作のモジュール rand, mrand, lrand でもそのようにプログラムしています。今回は米国の NIST FIPS140-2 という規格の中にある乱数の統計検定 (以下 FIPS140-2 と記述) を参考に、整数の一様乱数をビット列 (bit sequence) として取り扱い、ビット単位で乱数のばらつき具合を検定してみましょう。
たとえば、ある範囲 (0 - 0xffffffff) で乱数が一様に分布しているならば、ビット 0 と 1 の個数はほぼ均等になるはずです。つまり、ビット単位で度数検定ができるわけです。この方法を FIPS140-2 では "The Monobit Test" と呼んでいます。このほかに、連続する複数のビットからひとつの数値を作り、それを用いて度数検定を行う方法があります。FIPS140-2 ではこれを "The Poker Test" と呼んでいます。
また、無規則性の検定として、同じビットが連続する長さ (連長) を計測する "The Runs Test" (NIST FIPS140-2) や、連続する n ビットの中で 1 が出現する回数を計測する「組み合わせテスト」などがあります。なお、組み合わせテストは FIPS140-2 には存在しません。ご注意ください。
今回は 4 つの方法について説明します。なお、今回の検定プログラムは FIPS140-2 の基本的な考え方を参考にしただけであり、規格を厳密にプログラムしたわけではありません。また、FIPS140-2 は 1 回のテストで 20000 ビットの乱数データを使用し、そのテストを複数回行うようです。これをそのまま Python でプログラムすると実行時間がかかるので、統計的な検定で必要と思われるデータ数しか入力していません。あしからずご了承くださいませ。
まず最初に、整数の乱数列をビット単位で取り出すクラスを作ります。
リスト : ビットストリーム
class BitStream:
def __init__(self, m):
self.prg = m.Random(1)
self.buff = self.prg.irand()
self.index = 31
def getbit(self):
if self.index < 0:
self.buff = self.prg.irand()
self.index = 31
bit = (self.buff >> self.index) ↦ 1
self.index -= 1
return bit
def getbits(self, n):
bits = 0
for _ in range(n):
bits <<= 1
bits |= self.getbit()
return bits
クラス名は BitStream としました。インスタンス変数 prg に乱数クラスのインスタンスを、buff には prg のメソッド irand() で生成した整数の一様乱数 (32 bit) をセットします。index は buff から取り出すビットの位置を管理します。
メソッド getbit() は buff から 1 ビット取り出します。もし、index が 0 よりも小さい場合、buff は空になったので新しい乱数を生成して buff にセットします。それから、index の位置にあるビットの値を bit にセットして、index の値を -1 します。最後に bit を返します。
メソッド getbits() はビットを n 個取り出して、数値に変換して返します。getbit() を呼び出してビットの取り出し、それを変数 bits にセットしていきます。最後に bits を返します。
等確率性の検定で一番簡単な方法は、先ほど説明した The Monobit Test です。FIPS140-2 では 20000 bit のデータ数に対してビット 1 の個数 x が次の条件を満たすことを求めています。
9725 < x < 10275
本ページでは FIPS140-2 の条件を使わないで、単純にカイ 2 乗検定を行うことにします。0 と 1 のビット数を計測してカイ 2 乗値を求めると、それは自由度 1 のカイ 2 乗分布に従うはずです。あとは危険率 5 % の棄却域に入っていない回数を求め、それが試行回数の 95 % 程度になればいいわけです。プログラムは次のようになります。
リスト : The Monobit Test
def test04(bs, n):
count = [0] * 2
for _ in range(n):
count[bs.getbit()] += 1
return get_chi2(count, n / 2.0)
# test04 の実行
def exec_test04():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test04(bs, 5000) < 3.841:
i += 1
print(m)
print(i)
関数 test04() の引数 bs は BitStream のインスタンス、n は検定で使用するビット数です。メソッド getbit() でビットストリームから 1 ビット取り出して、0 ならば count[0] の値を +1 し、1 ならば count[1] の値を +1 します。そして、関数 get_chi2() で count のカイ 2 乗値を求めて返します。
関数 exec_test04() は test04() を 10000 回呼び出して、自由度 1 のカイ 2 乗分布で危険率 5 % の棄却域 (3.841, ∞) に入らない回数を計測します。この値が 9500 に近い値であればいいわけです。入力するビット数は 5000 としました。ビット数が少ないと正確な検定はできません。ご注意くださいませ。
実行結果は次のようになりました。
表 : test04 の実行結果
(10000回の試行)
線形合同法 : 9733
M 系列乱数 : 9486
lagged fib : 9502
メルセンヌ : 9513
線形合同法は明らかに 9500 よりも多いですね。線形合同法の下位ビットはランダムではないので、危険率 5 % の棄却域に入る回数は少なくなるようです。つまり、乱数のばらつき具合が足りないということです。あとのアルゴリズムは 9500 に近い値となりました。
もちろん、test04 の返り値をカイ 2 乗検定することもできます。次のリストを見てください。
リスト : test04 のカイ 2 乗検定
# 自由度 1 のカイ 2 乗分布表
chi2_1 = [
2.706, 1.642, 1.074, 0.708, 0.455,
0.275, 0.148, 0.0642, 0.0158, 0.000,
]
def chi2_test04():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test04(bs, 5000)
for i in range(10):
if chi2 >= chi2_1[i]:
count[i] += 1
break
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
関数 chi2_test04() は、今までのカイ 2 乗検定とほとんど同じですが、乱数の生成に BitStream を使うことと、自由度 1 のカイ 2 乗分布表 chi2_1 との一致を検定することに注意してください。検定結果は次のようになりました。
表 : test04 の検定結果
(10000回の試行)
線形合同法 : 61.12 NG
M 系列乱数 : 8.56 OK
lagged fib : 7.65 OK
メルセンヌ : 12.23 OK
線形合同法は No Good で、あとのアルゴリズムは OK でした。ビット単位でみると、線形合同法の性能は良くないことがわかります。
次は連続する 4 ビットをひとつずつまとめて度数検定してみましょう。FIPS140-2 では、これを "The Poker Test" と呼んでいます。次の式を見てください
The Poker Test は 20000 bit / 4 bit = 5000 個の乱数データでテストを行い、上式の統計量 \(x\) を求めます。この値が条件 2.16 < x < 46.17 を満たしていないといけません。これをそのまま Python でプログラムすると実行時間がかかるので、ここでは単純なカイ 2 乗検定を行うことにします。
値は 0 から 15 までの 16 種類あり、ビット 0 と 1 が等確率で分布しているならば、0 から 15 までの出現頻度はほぼ等しくなるはずです。これをカイ 2 乗検定します。プログラムは次のようになります。
リスト : The Poker Test
def test05(bs, n):
count = [0] * 16
for _ in range(n):
count[bs.getbits(4)] += 1
return get_chi2(count, n / 16.0)
# test05 の実行
def exec_test05():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test05(bs, 16 * 100) < 24.996:
i += 1
print(m)
print(i)
関数 test05() の引数 n は 4 ビット単位で生成する乱数の個数です。0 - 15 の出現頻度は配列 count で計測します。メソッド getbits() でビットストリームから 4 ビット取り出して、それに対応する count の値を +1 します。そして、関数 get_chi2() で count のカイ 2 乗値を求めて返します。このとき、期待度数は n / 16 になります。
関数 exec_test05() は test05 を 10000 回呼び出して、自由度 15 のカイ 2 乗分布で危険率 5 % の棄却域 (24.996, ∞) に入らない回数を求めます。この値が 9500 に近い値であればいいわけです。生成する数値の個数は 1 区間に 100 個程度入るように 16 * 100 としました。
実行結果は次のようになりました。
表 : test05 の実行結果
(10000回の試行)
線形合同法 : 9946
M 系列乱数 : 9438
lagged fib : 9494
メルセンヌ : 9495
The Monobit Test で不合格だった線形合同法は、The Poker Test でも 9500 より多くなりました。他のアルゴリズムは 9500 に近い値となりました。
次は test05 の返り値をカイ 2 乗検定します。
リスト : test05 のカイ 2 乗検定
# 自由度 15 のカイ 2 乗分布表
chi2_15 = [
22.307, 19.311, 17.322, 15.733, 14.339,
13.030, 11.721, 10.307, 8.547, 0.000,
]
def chi2_test05():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test05(bs, 16 * 100)
for i in range(10):
if chi2 >= chi2_15[i]:
count[i] += 1
break
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
関数 chi2_test05() では、自由度 15 のカイ 2 乗分布表 chi2_15 を使います。あとは関数 chi2_test04() と同じです。
検定結果は次のようになりました。
表 : test05 の検定結果 線形合同法 : 863.81 NG M 系列乱数 : 6.95 OK lagged fib : 10.68 OK メルセンヌ : 9.31 OK
線形合同法は当然ですが No Good で、あとのアルゴリズムは OK でした。
次はビット列で無規則性 (独立性) の検定を行ってみましょう。ビット 0 と 1 が等確率で出現し、その出現する順序に規則性が無い場合、同じビットが n 個続く確率は \({\frac{1}{2}}^{n}\) になります。これを表にすると、次のようになります。
| 連長 | 確率 |
|---|---|
| 1 | 0.5 |
| 2 | 0.25 |
| 3 | 0.125 |
| 4 | 0.0625 |
| 5 | 0.03125 |
| 6+ | 0.03125 |
同じビットが続く長さ (連長) を n とすると、n が長くなるほど出現する確率は小さくなります。6+ は 6 以上の連長が現れる確率を表しています。この表から期待度数を求めることができるので、0 または 1 の連長の個数を実際に求めて度数検定すればいいわけです。これを「連の検定」といい、FIPS140-2 では "The Runs Test" と呼んでいます。FIPS140-2 の条件を表に示します。
| 連長 | 条件 |
|---|---|
| 1 | 2315～2685 |
| 2 | 1114～1386 |
| 3 | 527～723 |
| 4 | 240～384 |
| 5 | 103～209 |
| 6+ | 103～209 |
FIPS140-2 は 20000 bit の乱数データを使用するので、Python でプログラムすると実行時間がかかります。そこで、ここでも単純なカイ 2 乗検定を行うことにします。
プログラムは次のようになります。
リスト : The Runs Test
SIZE = 6
prob_table = [0.5, 0.25, 0.125, 0.0625, 0.03125, 0.03125]
# 連による検定
def test06(bs, n, bit):
count = [0] * SIZE
num = 0
while num < n:
# bit を探す
while bs.getbit() != bit: pass
# bit の長さをカウントする
c = 0
while bs.getbit() == bit: c += 1
if c >= SIZE: c = SIZE - 1
count[c] += 1
num += 1
# カイ 2 乗値を求める
chi2 = 0.0
for i in range(SIZE):
F = n * prob_table[i]
chi2 += (count[i] - F) ** 2 / F
return chi2
# test06 の実行
def exec_test06(bit):
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test06(bs, 3000, bit) < 11.070:
i += 1
print(m)
print(i)
関数 test06() の引数 bs は BitStream クラスのインスタンスで、n が計測する連長の総数、bit が計測するビットの種類 (0 or 1) です。ビットストリームから getbit() で 1 ビット取り出し、それが bit と同じ値であれば、その連長をカウントします。変数 c の値は "連長 - 1" になることに注意してください。c が SIZE 以上になった場合は SIZE - 1 にします。それから、count[c] の値を +1 しますこれで各連長の出現頻度を求めることができます。あとは、カイ 2 乗値を計算するだけです。
関数 chi2_test06() は test06() を 10000 回呼び出して、自由度 5 のカイ 2 乗分布の棄却域 (11.070, ∞) に入らない回数を求めます。最小の確率が 0.03125 なので、この区間のデータ数が 100 程度になるよう連長の総数を 3000 個に設定してテストすることにします。
実行結果は次のようになりました。
表 : test06 の実行結果
(10000回の試行)
| 1 | 0
-----------+------+------
線形合同法 | 9802 | 9828
M 系列乱数 | 9455 | 9474
lagged fib | 9506 | 9510
メルセンヌ | 9500 | 9473
線形合同法はビット 0, 1 ともに 9500 よりも多くなりました。これは The Monobit Test と同じです。他のアルゴリズムは、ビット 0, 1 ともに 9500 回に近い値になりました。
それでは test06() の返り値をカイ 2 乗検定してみましょう。次のリストを見てください。
リスト : test06 の検定
def chi2_test06(bit):
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test06(bs, 3000, bit)
for i in range(10):
if chi2 >= chi2_5[i]:
count[i] += 1
break
print(count, sum(count))
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
プログラムは chi2_test05() とほとんど同じですが、自由度 5 のカイ 2 乗分布 (chi2_5) と一致するか検定することに注意してください。検定結果は次のようになりました。
表 : test06 の検定結果
| 1 | 0
-----------+-----------+-----------
線形合同法 | 131.76 NG | 186.12 NG
M 系列乱数 | 6.89 OK | 16.90 OK
lagged fib | 6.95 OK | 7.61 OK
メルセンヌ | 7.70 OK | 6.89 OK
ビット 0, 1 ともに、線形合同法は当然ですが No Good で、あとのアルゴリズムは OK でした。
もうひとつ、簡単な方法を紹介しましょう。たとえば、連続するビットを 4 個取り出した場合、ビット 1 の個数は 0 から 4 まであります。ビット 0, 1 の出現確率が等しくて規則性が無いと仮定すると、ビット 1 の個数の確率は二項分布 \(B(4, 0.5)\) に従います。実際に値を求めると、次のようになります。
| 0 | 1 | 確率 |
|---|---|---|
| 4 | 0 | 0.0625 |
| 3 | 1 | 0.25 |
| 2 | 2 | 0.375 |
| 1 | 3 | 0.25 |
| 0 | 4 | 0.0625 |
あとは、0 から 4 までの出現頻度をカウントし、度数検定すればいいわけです。なお、ビット数を増やしてもかまいませんが、カイ 2 乗検定に必要なデータ数を確保するように注意してください。ビット数を増やすとそれだけ多くの乱数を生成する必要があるので、このページでは 4 ビットとしました。なお、このテストは参考文献『乱数の知識』に書かれている方法を参考にしたものであり、 FIPS140-2 の中に「組み合わせテスト」は存在しません。
プログラムは次のようになります。
リスト : 組み合わせテスト
COMB_SIZE = 5
comb_table = [0.0625, 0.25, 0.375, 0.25, 0.0625]
def test07(bs, n):
count = [0] * COMB_SIZE
for _ in range(n):
i = 0
for _ in range(COMB_SIZE - 1):
if bs.getbit() == 1:
i += 1
count[i] += 1
# カイ 2 乗値の計算
chi2 = 0.0
for i in range(COMB_SIZE):
F = n * comb_table[i]
chi2 += (count[i] - F) ** 2 / F
return chi2
# test07 の実行
def exec_test07():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test07(bs, 2000) < 9.488:
i += 1
print(m)
print(i)
関数 test07() の引数 bs は BitStream クラスのインスタンスで、n がテストする個数です。getbit() でビットストリームから 4 ビット取り出して、ビットが 1 の個数を変数 i でカウントします。そして、count[i] の値を +1 します。これで、各区間の観測度数を求めることができます。あとは count のカイ 2 乗値を求めて返すだけです。
関数 chi2_test07() は test07() を 10000 回呼び出して、自由度 4 のカイ 2 乗分布の棄却域 (9.488, ∞) に入らない回数を求めます。最低の確率が 0.0625 なので、この区間の値が 100 程度になるよう 2000 個でテストすることにします。
実行結果は次のようになりました。
表 : test07 の実行結果
(10000回の試行)
線形合同法 : 9810
M 系列乱数 : 9484
lagged fib : 9507
メルセンヌ : 9487
線形合同法は 9500 よりも多くなりました。これは The Monobit Test と同じです。他のアルゴリズムは 9500 に近い値になりました。
それでは test07() の返り値をカイ 2 乗検定してみましょう。次のリストを見てください。
リスト : test07 のカイ 2 乗検定
# 自由度 4 のカイ 2 乗分布表
chi2_4 = [
7.779, 5.989, 4.878, 4.045, 3.357,
2.753, 2.195, 1.649, 1.064, 0.000,
]
def chi2_test07():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test07(bs, 2000)
for i in range(10):
if chi2 >= chi2_4[i]:
count[i] += 1
break
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
プログラムは chi2_test06() とほとんど同じですが、自由度 4 のカイ 2 乗分布 (chi2_4) と一致するか検定することに注意してください。検定結果は次のようになりました。
表 : test07 の実行結果 線形合同法 : 205.68 NG M 系列乱数 : 4.91 OK lagged fib : 3.81 OK メルセンヌ : 8.32 OK
線形合同法は当然ですが No Good で、あとのアルゴリズムは OK でした。
2 回にわたって擬似乱数の簡単な検定方法を説明しましたが、すべての検定で良好な結果を残したのはメルセンヌ・ツイスターだけで、その次が M 系列乱数でした。簡単な検定しか行っていませんが、メルセンヌ・ツイスターが優れた擬似乱数生成アルゴリズムであることが確認できたと思います。擬似乱数の検定にはいろいろな方法があるので、興味のある方は調べてみてください。
上記の参考文献は『統計科学のための電子図書システム』で公開されています。優れた著作を公開されている脇本和昌様および関係各位に感謝いたします。
#
# randtest.py : 擬似乱数の検定
#
# Copyright (C) 2007-2022 Makoto Hiroi
#
import rand, mrand, lrand, mtrand
import math
# カイ 2 乗分布表 0.10 きざみ
chi2_1 = [
2.706, 1.642, 1.074, 0.708, 0.455,
0.275, 0.148, 0.0642, 0.0158, 0.000,
]
chi2_4 = [
7.779, 5.989, 4.878, 4.045, 3.357,
2.753, 2.195, 1.649, 1.064, 0.000,
]
chi2_5 = [
9.236, 7.289, 6.064, 5.132, 4.351,
3.656, 3.000, 2.343, 1.610, 0.000,
]
chi2_9 = [
14.684, 12.242, 10.656, 9.414, 8.343,
7.357, 6.393, 5.380, 4.168, 0.000,
]
chi2_15 = [
22.307, 19.311, 17.322, 15.733, 14.339,
13.030, 11.721, 10.307, 8.547, 0.000,
]
chi2_16 = [
23.542, 20.465, 18.418, 16.780, 15.339,
13.983, 12.624, 11.152, 9.312, 0.000,
]
chi2_19 = [
27.203, 23.900, 21.689, 19.910, 18.338,
16.850, 15.352, 13.716, 11.651, 0.000,
]
chi2_24 = [
33.196, 29.553, 27.096, 25.106, 23.337,
21.652, 19.943, 18.062, 15.659, 0.000,
]
#
# 標準正規分布
#
# 分割数
NORM_SIZE = 10
# 標準正規分布を 10 に分割
norm_table = [
1.29, 0.85, 0.53, 0.26, 0.00,
-0.26, -0.53, -0.85, -1.29, -1.0e301,
]
# 各区間の確率
norm_p = [
0.0985, 0.0992, 0.1004, 0.0993, 0.1026,
0.1026, 0.0993, 0.1004, 0.0992, 0.0985,
]
# カイ 2 乗値を求める
def get_chi2(count, F):
chi2 = 0.0
for x in count:
chi2 += (x - F) ** 2 / F
return chi2
#
# 度数検定
#
def test01(prg, n, k):
count = [0] * k
for _ in range(n):
count[int(prg.random() * k)] += 1
return get_chi2(count, float(n) / k)
# test01 の実行
def exec_test01():
for k, limit in [(10, 16.919), (20, 30.144), (25, 36.415)]:
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
n = k * 100 # 1 つの区間に 100 個
i = 0
for _ in range(10000):
if test01(prg, n, k) < limit:
i += 1
print(m)
print(k, i)
# test01 のカイ 2 乗検定
def chi2_test01(k, chi2_table):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
count = [0] * 10
for _ in range(2000):
chi2 = test01(prg, 100 * k, k)
for i in range(10):
if chi2 >= chi2_table[i]:
count[i] += 1
break
#
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
#
# 相関性のテスト
#
# おくれ k の系列相関係数
def test02(prg, n, k):
buff = [prg.random() for _ in range(k)]
sum = 0.0
for _ in range(n):
m = prg.random()
sum += buff[0] * m
del buff[0]
buff.append(m)
return math.sqrt(n) * (12.0 * sum / n - 3.0) / math.sqrt(13.0)
# test02 の実行
def exec_test02(k):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
i = 0
for _ in range(10000):
if -1.96 < test02(prg, 2500, k) < 1.96:
i += 1
print(m)
print(i)
# test02 の検定
def chi2_test02(k):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
count = [0] * 10
for _ in range(2000):
chi2 = test02(prg, 2500, k)
for i in range(NORM_SIZE):
if chi2 >= norm_table[i]:
count[i] += 1
break
chi2 = 0.0
for i in range(10):
F = 2000 * norm_p[i]
chi2 += (count[i] - F) ** 2 / F
# 自由度 9
print(m)
print(chi2, chi2 < 16.919)
#
# 分割表による検定
#
# 分割表のカイ 2 乗値を求める
def get_chi2_matrix(count):
total = 0.0
xn = len(count)
yn = len(count[0])
xcount = [0] * xn
ycount = [0] * yn
for x in range(xn):
for y in range(yn):
a = count[x][y]
total += a
xcount[x] += a
ycount[y] += a
# カイ 2 乗値を計算する
v = 0.0
for x in range(xn):
for y in range(yn):
F = xcount[x] * ycount[y] / total
a = count[x][y] - F
v += a * a / F
return v
# 分割表による無相関性の検定
def test03(prg, n, m, k):
count = [[0] * m for _ in range(m)]
buff = [prg.random() for _ in range(k)]
for _ in range(n):
a = prg.random()
x = int(buff[0] * m)
y = int(a * m)
count[x][y] += 1
del buff[0]
buff.append(a)
return get_chi2_matrix(count)
# test03 の実行
def exec_test03(k):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
i = 0
for _ in range(10000):
if test03(prg, 5 * 5 * 100, 5, k) < 26.296:
i += 1
print(m)
print(i)
# test03 の検定
def chi2_test03(k):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
count = [0] * 10
for _ in range(2000):
chi2 = test03(prg, 5 * 5 * 100, 5, k)
for i in range(10):
if chi2 >= chi2_16[i]:
count[i] += 1
break
chi2 = get_chi2(count, 200.0)
print(m)
print(count, sum(count))
print(chi2, chi2 < 16.919)
#
# Sum-Test
#
# 標準正規分布に従う乱数を生成
def nrand(prg, n):
v = 0.0
for _ in range(n):
v += prg.random()
v /= n
return (v - 0.5) * math.sqrt(12 * n)
# 標準正規分布に従っているか
def sum_test(n, k):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
i = 0
for _ in range(k):
if -1.96 < nrand(prg, n) < 1.96:
i += 1
print(m)
print(i)
# カイ 2 乗検定
def chi2_sum_test(n, k):
for m in [rand, mrand, lrand, mtrand]:
prg = m.Random(1)
count = [0] * NORM_SIZE
for _ in range(k):
x = nrand(prg, n)
for i in range(NORM_SIZE):
if x >= norm_table[i]:
count[i] += 1
break
#
v = 0.0
for i in range(NORM_SIZE):
F = k * norm_p[i]
v += (count[i] - F) ** 2 / F
# 自由度 9
print(m)
print(v, v < 16.919)
##### ビット単位での検定 #####
#
# ビット入力ストリーム
#
class BitStream:
def __init__(self, m):
self.prg = m.Random(1)
self.buff = self.prg.irand()
self.index = 31
def getbit(self):
if self.index < 0:
self.buff = self.prg.irand()
self.index = 31
bit = (self.buff >> self.index) & 1
self.index -= 1
return bit
def getbits(self, n):
bits = 0
for _ in range(n):
bits <<= 1
bits |= self.getbit()
return bits
#
# The Monobit Test
#
def test04(bs, n):
count = [0] * 2
for _ in range(n):
count[bs.getbit()] += 1
return get_chi2(count, n / 2.0)
# test04 の実行
def exec_test04():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test04(bs, 5000) < 3.841:
i += 1
print(m)
print(i)
# test04 の検定
def chi2_test04():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test04(bs, 5000)
for i in range(10):
if chi2 >= chi2_1[i]:
count[i] += 1
break
print(count, sum(count))
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
#
# The Poker Test
#
def test05(bs, n):
count = [0] * 16
for _ in range(n):
count[bs.getbits(4)] += 1
return get_chi2(count, n / 16.0)
# test05 の実行
def exec_test05():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test05(bs, 16 * 100) < 24.996:
i += 1
print(m)
print(i)
# test05 の検定
def chi2_test05():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test05(bs, 16 * 100)
for i in range(10):
if chi2 >= chi2_15[i]:
count[i] += 1
break
print(count, sum(count))
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
#
# The Runs Test
#
# 連の確率
SIZE = 6
prob_table = [0.5, 0.25, 0.125, 0.0625, 0.03125, 0.03125]
# 連による検定
def test06(bs, n, bit):
count = [0] * SIZE
num = 0
while num < n:
# bit を探す
while bs.getbit() != bit: pass
# bit の長さをカウントする
c = 0
while bs.getbit() == bit: c += 1
if c >= SIZE: c = SIZE - 1
count[c] += 1
num += 1
# カイ 2 乗値を求める
chi2 = 0.0
for i in range(SIZE):
F = n * prob_table[i]
chi2 += (count[i] - F) ** 2 / F
return chi2
# test06 の実行
def exec_test06(bit):
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test06(bs, 3000, bit) < 11.070:
i += 1
print(m)
print(i)
# test06 の検定
def chi2_test06(bit):
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test06(bs, 3000, bit)
for i in range(10):
if chi2 >= chi2_5[i]:
count[i] += 1
break
print(count, sum(count))
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)
#
# 組み合わせテスト
#
# 組み合わせの確率
COMB_SIZE = 5
comb_table = [0.0625, 0.25, 0.375, 0.25, 0.0625]
def test07(bs, n):
count = [0] * COMB_SIZE
for _ in range(n):
i = 0
for _ in range(COMB_SIZE - 1):
if bs.getbit() == 1:
i += 1
count[i] += 1
# カイ 2 乗値の計算
chi2 = 0.0
for i in range(COMB_SIZE):
F = n * comb_table[i]
chi2 += (count[i] - F) ** 2 / F
return chi2
# test07 の実行
def exec_test07():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
i = 0
for _ in range(10000):
if test07(bs, 2000) < 9.488:
i += 1
print(m)
print(i)
# test07 の検定
def chi2_test07():
for m in [rand, mrand, lrand, mtrand]:
bs = BitStream(m)
count = [0] * 10
for _ in range(2000):
chi2 = test07(bs, 2000)
for i in range(10):
if chi2 >= chi2_4[i]:
count[i] += 1
break
print(count, sum(count))
chi2 = get_chi2(count, 200.0)
print(m)
print(chi2, chi2 < 16.919)