Algorithms with Python / 擬似乱数の検定

はじめに

今回は擬似乱数の統計的な検定について取り上げます。擬似乱数のアルゴリズムは、拙作のページ「乱数 [1] [2]」で説明しました。取り上げたアルゴリズムは線形合同法、M 系列乱数、lagged fibonacci 法の 3 つです。このときは簡単なテストしか行っていなかったので、今回は統計的な検定を試してみましょう。

●擬似乱数の性質

乱数には重要な性質が 2 つあります。ひとつは「等確率性 (等出現性)」です。たとえば、サイコロを考えてみましょう。サイコロの目はイカサマをしないかぎり、出る確率が 1/6 になります。したがって、サイコロを振る回数 (試行回数 n) を増やしていくと、どの目も出現する回数は理論値 (n/6) に限りなく近づいていきます。このような確率的な性質を「大数の法則」といいます。そして、ある範囲の整数値が等確率で現れる乱数を「整数の一様乱数」といいます。また、実数の一様乱数は、整数の一様乱数から生成することができます。

もうひとつの重要な性質が「無規則性 (独立性)」です。たとえば、サイコロには規則性がないので、目の出る順番はまったくでたらめになります。いま 2 が出たから次は 1 が出るとか、3, 4 と続いたから次は 5 が出るなどのような規則性はありません。ようするに、以前に出現した数によって次に出現する数が決まるようなことはないのです。これが無規則性で、出現した各数には相関性がないことから、無相関性とか独立性と呼ばれることもあります。

ただし、コンピュータで生成する擬似乱数の場合、無規則性を完全に実現することはできない、ということに注意してください。一般的な擬似乱数のアルゴリズムは、コンピュータのメモリに現在の状態 (数値) を記憶しておきます。乱数を生成するとき、メモリの内容をある手順 (A) で書き換えて新しい状態へ移行します。そして、その状態をある手順 (B) で数値に変換し、それを乱数として出力します。これを繰り返すことで乱数列を生成しています。

ここで、手順 (A) や (B) をどんなに複雑なものにしたとしても、あらかじめ「決められた手順 (アルゴリズム)」を実行していることにかわりはありません。したがって、擬似乱数の規則性をゼロにすることは原理的に不可能というわけです。また、コンピュータにあるメモリは有限なので、擬似乱数は周期性を持つこともわかるでしょう。

これに対し、何らかの物理的現象 (たとえば電子回路の熱雑音など) を利用して乱数を生成する方法があります。このような方法で生成された乱数を「物理乱数」と呼ぶ場合があります。物理乱数は規則性がまったくないので、まさに「真の乱数」といえるでしょう。そのかわり、同じ乱数列を再現するには、発生した乱数を記憶しておかない限り不可能です。つまり、物理乱数には「再現性」がまったくないのです。

擬似乱数は真の乱数ではありませんが、物理乱数にはない「再現性」があります。擬似乱数は、種 (シード) に同じ値を設定することで同じ乱数列を簡単に発生させることができます。シミュレーションやモンテカルロ法などの実験で乱数を用いる場合、この機能はとても役に立ちます。シードに別な値を設定すれば、異なる乱数列を簡単に発生させることができます。シードに同じ値を与えれば、追試も簡単に行うことができます。

このように、擬似乱数には再現性があるのでとても便利なのですが、真の乱数ではありません。このため、乱数を用いるアプリケーションを作成する場合、その乱数がアプリケーションに必要な性能を満たしているか確かめておく必要があります。乱数の性質は理論的に解析できるものがあります。一般に、理論的な解析は乱数列の一周期全体にわたって行われる場合がほとんどで、乱数列の一部分を利用する場合は、実際に乱数のばらつき具合を統計的に調べておいたほうがよいでしょう。

ところで、擬似乱数に求められる性能の基準はアプリケーションによって異なります。たとえば、簡単なゲームであれば線形合同法で十分ですが、複雑なシミュレーションになると線形合同法では力不足になる場合が多くなるでしょう。また、大量の乱数を消費する大規模なシミュレーションでは、擬似乱数の性能として長い周期性と高速性も求められています。このような場合、現時点ではメルセンヌ・ツイスターが最有力候補になると思います。

このほかに、擬似乱数は暗号でも用いられています。暗号で要求される擬似乱数の性質は等確率性と無規則性だけではなく、次の数が予測できないこと (予測不能性) が求められます。真の乱数は無規則であるため予測不能ですが、今まで説明した擬似乱数アルゴリズム、線形合同法、M 系列乱数、lagged fibonacci 法はすべて予測することが可能です。これはメルセンヌ・ツイスターも同じです。

これらのアルゴリズムをそのまま暗号で用いることはできません。Mersenne Twister Home Page - よくある質問によると、『Secure Hashing Algorithm (数ワードを圧縮して１ワードを生成する、非可逆的なアルゴリズム) と組み合わせて使う必要があります。』とのことです。本ページでは、暗号で使用する擬似乱数のことは考えないで、等確率性と無規則性を統計的に検定する簡単な方法を説明します。

●等確率性の検定

等確率性 (一様性) を検定する簡単な方法は、乱数の出現頻度をカイ 2 乗検定することです。これを「度数検定」といいます。基本的な考え方は簡単です。区間 [0, 1.0) の実数を生成する一様乱数を検定する場合、区間 [0, 1.0) を均等に k 分割します。乱数を n 個生成した場合、分割した各区間には n / k 個の乱数が入るはずです。各区間の観測度数を \(x_{i}\) 期待度数を \(F (= n / k)\) とすると、次に示す統計量は自由度 k - 1 のカイ 2 乗分布に従います。

\( \chi^{2} = \displaystyle \sum_{i=1}^k \dfrac{(x_{i} - F)^{2}}{F} \)

度数検定はプログラムも簡単です。次のリストを見てください。

リスト : 度数検定

# カイ 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() の引数 prg が擬似乱数クラスのインスタンス、n が生成する乱数の個数、k が分割数です。各区間の度数は count に格納します。あとは prg.random() で実数の一様乱数を生成し、各区間ごとに出現した乱数をカウントしてカイ 2 乗値を計算するだけです。

test01() は exec_test01() から呼び出します。検定を 10000 回行って、カイ 2 乗値が危険率 5 % の棄却域に入らない回数を求めます。擬似乱数が一様分布しているならば、その回数は 9500 回前後になるはずです。そうでなければ、その擬似乱数の等確率性に問題があるということです。

変数 k が分割数で、変数 limit が自由度 k - 1 のカイ 2 乗分布で危険率 5 % の値です。(limit, ∞) が棄却域になります。rand （プログラムリスト１) は線形合同法、mrand (プログラムリスト２) は M 系列乱数、lrand (プログラムリスト３) は lagged fibonacci 法のモジュールで、拙作のページ「乱数 [1] [2]」で作成したものと同じです。また、mtrand は Python の random モジュールのラッパークラスで、インターフェースを他のモジュールに合わせるために作成しました。次のリストを見てください。

リスト : random モジュール用のクラス

import random

class Random:
    def __init__(self, seed):
        random.seed(seed)

    def irand(self):
        return random.randint(0, 0xffffffff)

    def random(self):
        return random.random()

random モジュールのアルゴリズムにはメルセンヌ・ツイスターが使われています。生成する乱数の個数 n は、ひとつの区間に 100 個程度入るように k * 100 としました。

それでは実行結果を示します。

    表 : test01 の実行結果
         (10000回の試行)

           |  10  |  20  |  25  
-----------+------+------+------
線形合同法 | 9528 | 9478 | 9509
M 系列乱数 | 9515 | 9473 | 9493
lagged fib | 9489 | 9512 | 9499
メルセンヌ | 9504 | 9510 | 9489

どのアルゴリズムでも 9500 に近い値になりました。度数検定の結果は良好といえるでしょう。ところで、9500 に近い値になるはずといっても、どの程度までなら許されるのか基準が明確ではありませんね。そこで、関数 test01() の返り値がカイ 2 乗分布に従っているか、これをカイ 2 乗検定することにしましょう。次のリストを見てください。

リスト : test01 のカイ 2 乗検定

# カイ 2 乗分布表 (0.10 きざみ)
chi2_9 = [
  14.684, 12.242, 10.656, 9.414, 8.343,
   7.357,  6.393,  5.380, 4.168, 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,
]

# 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)

カイ 2 乗分布の区間 [0, ∞) を 10 分割します。このとき、各区間が等確率になるように分割します。その境界値を配列 chi2_9, chi2_19, chi2_24 に格納します。たとえば、chi2_9 において、[14.684, ∞) の確率が 0.10 で、[12.242, 14.684) の確率が 0.10 になります。

関数 chi2_test01() の引数 k は分割数で、chi2_table にはカイ 2 乗分布表を渡します。各区間の出現頻度は配列 count に格納します。各区間には 200 個の値が入るように 2000 回の試行を行います。そして、count のカイ 2 乗値を求めると、それは自由度 9 のカイ 2 乗分布に従うはずです。あとは、それを危険率 5 % の棄却域 (16.919, ∞) で判定すればいいわけです。

実行結果は次のようになりました。

        表 : test01 の検定結果

           |   10    |   20    |   25
-----------+---------+---------+---------
線形合同法 | 7.73 OK | 3.04 OK | 7.84 OK
M 系列乱数 | 6.74 OK |14.78 OK | 3.91 OK
lagged fib | 5.88 OK | 9.69 OK | 9.67 OK
メルセンヌ | 8.43 OK | 7.53 OK |12.89 OK

どのアルゴリズムでもカイ 2 乗値は棄却域に入っていません。test01() の返り値はカイ 2 乗分布に従っていて、擬似乱数は実数の区間 [0, 1) において一様に分布しているとみなしてよさそうです。

●無規則性の検定

●系列相関係数による検定

\(\begin{eqnarray} r_{xy} &=& \dfrac{1}{n} \dfrac{\sum_{i=1}^n (x_{i} - x_{m})(x_{i+k} - x_{m})}{\frac{1}{12}} \\ &=& \dfrac{12}{n} \left(\displaystyle \sum_{i=1}^n \left(x_{i}x_{i+k} - \dfrac{x_{i}}{2} - \dfrac{x_{i+k}}{2} + \dfrac{1}{4}\right) \right) \\ &=& \dfrac{12}{n} \left(\displaystyle \sum_{i=1}^n x_{i}x_{i+k} - \dfrac{n}{4} \right) \\ &=& \dfrac{12}{n} \displaystyle \sum_{i=1}^n x_{i}x_{i+k} - 3 \end{eqnarray}\)

前回説明しましたが、無規則性の検定には相関係数を用いることができます。擬似乱数列を \(x_{1}, x_{2}, \ldots, x_{n}\) として、対のデータ \((x_{i}, x_{i+k}) \ (k \geq 1)\) の相関係数を求めて検定することを考えます。実数の区間 [0, 1.0) の一様乱数は、平均が 1/2 で分散が 1/12 であることを用いると、相関係数は次の式で表すことができます。

参考文献『乱数の知識』によると、\(r_{xy}\) を「おくれ k の系列相関係数」といい、次式に示す統計量 \(Z\) は n が十分に大きいと正規分布 \(N(0, 1)\) に従うことが証明されているそうです。

\( Z = \dfrac{r_{xy} \sqrt n}{\sqrt 13} \)

標準正規分布において危険率 5 % の棄却域は (-∞, -1.96), (1.96, ∞) です。帰無仮説として「\(x_{i}\) と \(x_{i+k}\) に相関性はない」をたて、\(Z\) の値が棄却域に入っていれば、帰無仮説を棄却して \(x_{i}\) と \(x_{i+k}\) には相関性がある、という結論になります。

プログラムは簡単です。次のリストを見てください。

リスト : おくれ 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() の引数 prg は乱数クラスのインスタンス、n が生成する乱数の個数、k が「おくれ」を表す数値です。buff には直前に生成した k 個の乱数を格納しておきます。生成した乱数を m とすると、buff[0] * m が \(x_{i} * x_{i+k}\) に相当します。そのあとで buff を更新します。buff[0] を削除して、m を buff の後ろに追加するだけです。最後に統計量 Z を計算して返します。

関数 exec_test02() は test02() を 10000 回呼び出して、危険率 5 % の棄却域 (-∞, -1.96), (1.96, ∞) に入らない回数を求めます。test02() を呼び出すとき、生成する乱数の個数は、次に説明する「分割表によるカイ 2 乗検定」と同じにするため 2500 としました。

それでは実行結果を示します。

    表 : test02 の実行結果
         (10000回の試行)

           | k=1  | k=2  | k=3
-----------+------+------+------
線形合同法 | 9484 | 9513 | 9466
M 系列乱数 | 9469 | 9497 | 9468
lagged fib | 9498 | 9513 | 9494
メルセンヌ | 9535 | 9525 | 9508

どのアルゴリズムでも 9500 に近い値になりました。おくれ k の系列相関係数による検定の結果は、k = 1, 2, 3 では良好といえるでしょう。

ところで、度数検定と同様に test02() の返り値が標準正規分布に従っているか、カイ 2 乗検定することができます。次のリストを見てください。

リスト : test02 のカイ 2 乗検定

# 分割数
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,
]

# 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)

標準正規分布を 10 区間に分割し、その境界値を配列 norm_table にセットします。配列 norm_p は各区間の確率を表します。たとえば、[1.29, ∞) の確率は 0.0985 で、[0.85, 1.29) の確率は 0.0992 になります。

あとの処理は関数 chi2_test01() とほとんど同じです。test02() を 2000 回呼び出して各区間の観測度数を求めて、そのカイ 2 乗値を計算します。期待度数 F は 2000 * norm_p[i] で計算することに注意してください。

それでは検定結果を示します。

        表 : test02 の検定結果

           |  k = 1   |  k = 2   |  k = 3
-----------+----------+----------+----------
線形合同法 |  6.06 OK |  5.13 OK | 10.64 OK
M 系列乱数 |  6.47 OK |  8.02 OK |  3.02 OK
lagged fib |  8.86 OK |  7.96 OK |  8.31 OK
メルセンヌ |  7.98 OK | 12.81 OK |  5.43 OK

どのアルゴリズムでもカイ 2 乗値は棄却域に入っていません。test02() の返り値はカイ 2 乗分布に従っていて、おくれ k = 1, 2, 3 において相関性は無いと判断してもよさそうです。

●分割表による検定

無規則性 (独立性) の検定は、分割表を使ったカイ 2 乗検定でも行うことができます。次のリストを見てください。

リスト : 分割表による独立性の検定

# 分割表のカイ 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)

関数 get_chi2_matrix() は分割表のカイ 2 乗値を求めます。これは拙作のページ「統計学の基礎知識 [2]」で作成したプログラムと同じです。test03() は 5 行 5 列の分割表を使って独立性を検定します。引数 k は test02() と同じく「おくれ」を表す数値です。配列 buff の使い方は test02() と同じで、変数 x が \(x_{i}\) を、変数 y が \(x_{i+k}\) を表します。

関数 exec_test03() は test03() を 10000 回呼び出して、自由度 4 * 4 = 16 危険率 5 % の棄却域 (26.296, ∞) に入らない回数を求めます。test03() で生成する乱数の個数は、ひとつの区間に 100 個の乱数が入るように 2500 としました。これは test02() と同じ個数です。

それでは実行結果を示します。

    表 : test03 の実行結果
         (10000回の試行)

           | k=1  | k=2  | k=3
-----------+------+------+------
線形合同法 | 9481 | 9460 | 9506
M 系列乱数 | 9534 | 9519 | 9488
lagged fib | 9497 | 9481 | 9500
メルセンヌ | 9530 | 9484 | 9471

どのアルゴリズムでも 9500 に近い値になりました。分割表による独立性の検定でも、その結果は良好といえるでしょう。次に、test03() の返り値をカイ 2 乗検定してみましょう。次のリストを見てください。

リスト : test03 をカイ 2 乗検定

# カイ 2 乗分布表
chi2_16 = [
  23.542, 20.465, 18.418, 16.780, 15.339,
  13.983, 12.624, 11.152,  9.312,  0.000,
]

# 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(chi2, chi2 < 16.919)

関数 chi2_test03() は自由度 16 のカイ 2 乗分布表 chi2_16 を使います。あとの処理は今までのカイ 2 乗検定のプログラムと同じです。検定結果は次のようになりました。

        表 : test03 の検定結果

           |  k = 1   |  k = 2   |  k = 3
-----------+----------+----------+----------
線形合同法 | 10.07 OK |  6.02 OK |  7.94 OK
M 系列乱数 | 13.38 OK | 10.83 OK |  6.16 OK
lagged fib |  6.71 OK | 13.18 OK |  7.63 OK
メルセンヌ |  7.32 OK | 13.05 OK |  8.05 OK

どのアルゴリズムでもカイ 2 乗値は棄却域に入っていません。test03() の返り値は自由度 16 のカイ 2 乗分布に従っていて、おくれ k = 1, 2, 3 において相関性に問題は無いと判断してもよさそうです。

●Sum Test

拙作のページ「統計学の基礎知識」で説明しましたが、n がある程度大きい場合、母集団の分布がどのようなものであっても、標本平均 m の分布は正規分布で近似することができます。これを「中心極限定理」といいます。

●中心極限定理

平均 \(\mu\) 分散 \(\sigma^{2}\) の母集団から大きさ n (n >= 30) の標本を無作為抽出するとき、母集団がどのような分布であったとしても、その標本平均 \(m\) の分布は正規分布 \(N\left(\mu, \frac{\sigma^{2}}{n}\right)\) に従う。

区間 [0, 1.0) の一様乱数の場合、平均が 1/2 で分散が 1/12 になるので、n 個の乱数の平均値は正規分布 \(N(1/2, 1/12n)\) に従います。したがって、n 個の乱数の平均値を \(x_{m}\) とすると、次に示す統計量 \(Z\) は標準正規分布 \(N(0, 1)\) に従います。

\( Z = \dfrac{x_{m} - 0.5}{\sqrt{\frac{1}{12n}}} = (x_{m} - 0.5) \sqrt {12n} \)

統計量 \(Z\) が標準正規分布に従っているかテストすることで、擬似乱数の検定を行うことができます。このテストを "Sum Test" といいます。ここで、乱数の個数 n を大きくすると、それだけ正規分布への近似度が高くなります。そして、試行回数を増やしていくと、「大数の法則」によりその分布は限りなく正規分布へ近づいていきます。

したがって、擬似乱数で Sum Test を行う場合、n をある程度大きくして、試行回数を多くするほど、統計量 \(Z\) は正規分布に一致するようになるわけです。そうでなければ、そのアルゴリズムには統計的な欠点がある、ということになります。

プログラムは簡単です。次のリストを見てください。

リスト : 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)

関数 nrand() は標準正規分布 N(0, 1) に従う乱数を生成します。乱数を n 個加算して平均値を求め、統計量 Z を計算して返します。関数 sum_test() は nrand() を k 回呼び出して、危険率 5 % の棄却域 (-∞, -1.96), (1.96, ∞) に入らない回数を求めます。この回数が k の 95 % に近い値であればいいわけです。

実行結果は次のようになりました。

       表 : sum_test の実行結果

 n = 12    | 10000 | 100000 | 1000000
-----------+-------+--------+---------
線形合同法 |  9479 |  94885 |  950714
M 系列乱数 |  9521 |  95145 |  951205
lagged fib |  9555 |  94993 |  951124
メルセンヌ |  9486 |  95061 |  950854

 n = 50    | 10000 | 100000 | 1000000
-----------+-------+--------+---------
線形合同法 |  9487 |  95076 |  950436
M 系列乱数 |  9509 |  95139 |  950350
lagged fib |  9529 |  94947 |  949907
メルセンヌ |  9492 |  95107 |  950468

 n = 100   | 10000 | 100000 | 1000000
-----------+-------+--------+---------
線形合同法 |  9459 |  94949 |  950377
M 系列乱数 |  9499 |  95102 |  949949
lagged fib |  9546 |  95057 |  949665
メルセンヌ |  9504 |  94991 |  950153

加算する乱数の個数は 12, 50, 100 とし、試行回数は 1 万、10 万、100 万回としました。n = 12 の場合、正規分布への近似度はあまりよくありません。それでも、どのアルゴリズムも 95 % に近い値となりました。

次は nrand の返り値をカイ 2 乗検定してみましょう。

リスト : Sum Test のカイ 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
        # カイ 2 乗値を求める
        v = 0.0
        for i in range(NORM_SIZE):
            F = k * norm_p[i]
            v += (count[i] - F) ** 2 / F
        print(m)
        print(v, v < 16.919)

プログラムは、おくれ k の系列相関係数のカイ 2 乗検定 (chi2_test02) とほぼ同じです。nrand() を k 回呼び出して、その返り値の分布が標準正規分布と一致するか、カイ 2 乗検定しています。

検定結果は次のようになりました。

       表 : Sum Test の検定結果

 n = 12    |  10000   |  100000  | 1000000
-----------+----------+----------+-----------
線形合同法 |  1.74 OK | 29.67 NG | 113.11 NG
M 系列乱数 | 11.99 OK | 26.55 NG |  84.27 NG
lagged fib |  7.01 OK | 14.09 OK |  82.43 NG
メルセンヌ |  4.97 OK | 14.59 OK | 109.70 NG

 n = 50    |  10000   |  100000  | 1000000
-----------+----------+----------+----------
線形合同法 |  5.12 OK |  6.00 OK |  7.50 OK
M 系列乱数 | 19.14 NG |  6.10 OK |  8.07 OK
lagged fib |  7.05 OK | 13.93 OK | 36.26 NG
メルセンヌ |  8.79 OK |  4.11 OK | 14.37 OK

 n = 100   |  10000   |  100000  | 1000000
-----------+----------+----------+----------
線形合同法 |  5.70 OK | 17.69 NG |  7.34 OK
M 系列乱数 | 10.69 OK |  3.03 OK |  7.01 OK
lagged fib |  3.74 OK | 12.14 OK | 19.32 NG
メルセンヌ |  9.01 OK |  5.96 OK | 13.99 OK

n = 12 の場合、正規分布への近似度が低いため、試行回数を増やしていくと、どのアルゴリズムでも結果は No Good になります。n = 50 の場合、1 万回のテストで M 系列乱数が、100 万回のテストで lagged fibonacci 法が No Good になりました。M 系列乱数の場合、100 万回のテストでは OK になるので、たまたま悪い乱数列が生成されたのかもしれません。興味のある方はシードを変えて試してみてください。

n = 100 にすると、10 万回のテストで線形合同法が No Good に、100 万回のテストで lagged fibonacci 法が No Good になりました。乱数の個数 n と試行回数の両方を増やしても lagged fibonacci 法は No Good になるので、拙作のページ「乱数 [2]」で作成した lagged fibonacci 法のプログラムには統計的な欠陥があると考えられます。

また、参考文献『特集１乱数をきわめろ』によると、lagged fibonacci 法は UNIX 系 OS のＣ言語ライブラリ関数 random() で使われているそうです。この random() にも統計的な欠陥が見つかっていて、Sum Test を行うと正規分布から少し外れるそうです。今回の Sum Test でも同様の結果が得られたので、lagged fibonacci 法のアルゴリズムに何か問題点があるのかもしれません。

次回は乱数をビット単位で検定する方法を説明します。

●参考文献, URL

脇本和昌, 『身近なデータによる統計解析入門』,森北出版株式会社, 1973 年 (絶版)
脇本和昌, 『乱数の知識』,森北出版株式会社, 1970 年 (絶版)
和田維作, 宮坂電人『特集１乱数をきわめろ』, C MAGAZINE 2004 年 6 月号, ソフトバンククリエイティブ

上記の参考文献 1, 2 は『統計科学のための電子図書システム』で公開されています。優れた著作を公開されている脇本和昌様および関係各位に感謝いたします。

●プログラムリスト１

#
# rand.py : 線形合同法
#
#           Copyright (C) 2007-2022 Makoto Hiroi
#
# 定数
RAND_MAX = 0xffffffff

class Random:
    def __init__(self, seed = 1):
        self.seed = seed

    def irand(self):
        self.seed = (69069 * self.seed + 1) & RAND_MAX
        return self.seed

    def random1(self):
        return (1.0 / (RAND_MAX + 1.0)) * (self.irand() + 0.5)

    def random(self):
        return (1.0 / (RAND_MAX + 1.0)) * self.irand()

if __name__ == '__main__':
    a = Random()
    for i in range(8):
        print(a.irand())
    print('-----')
    a = Random()
    for i in range(8):
        print(a.random1())
    print('-----')
    a = Random()
    for i in range(8):
        print(a.random())

●プログラムリスト２

#
# mrand.py : M 系列乱数
#
#            Copyright (C) 2007-2022 Makoto Hiroi
#
# 参考文献
# 奥村晴彦, 『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
# M 系列乱数 (389, 390 ページ)
# 定数
P = 521
Q = 32
BITS = 32
RAND_MAX = 0xffffffff

class Msequence:
    def __init__(self, seed):
        self.buff = seed
        self.p = 0
        self.q = P - Q

    def irand(self):
        self.buff[self.p] ^= self.buff[self.q]
        value = self.buff[self.p]
        self.p += 1
        if self.p >= P: self.p = 0
        self.q += 1
        if self.q >= P: self.q = 0
        return value

class Random(Msequence):
    def __init__(self, seed = 1):
        self.seed = seed
        self.buff = []
        self.p = 0
        self.q = P - Q
        # 初期化
        # ビット単位でセットしているので遅い
        # オリジナルのプログラムは高速
        i = 0
        seed_bits = []
        for _ in range(P):
            n = 0
            for _ in range(BITS):
                if i < P:
                    # 線形合同法の最上位ビットを使う
                    b = self._rand() >> (BITS - 1)
                    seed_bits.append(b)
                    n = (n << 1) | b
                    i += 1
                    if i == P: a = Msequence(seed_bits)
                else:
                    n = (n << 1) | a.irand()
            self.buff.append(n)
        # warm up
        for _ in range(P * 3): self.irand()

    # 実数の一様乱数
    def random(self):
        return (1.0 / (RAND_MAX + 1.0)) * self.irand()

    # 線形合同法
    def _rand(self):
        self.seed = (1566083941 * self.seed + 1) & RAND_MAX
        return self.seed

# テスト
if __name__ == '__main__':
    a = Random()
    for _ in range(15): print(a.irand())

●プログラムリスト３

#
# lrand.py : lagged fibonacci 法
#
#            Copyright (C) 2007-2022 Makoto Hiroi
#
# 定数
RAND_MAX = 0xffffffff
P = 63
Q = 31

class Random:
    def __init__(self, seed = 1):
        self._init_buff(seed)    # _init_buff のほうがよい
        self.p = 0
        self.q = P - Q
        # warm up
        for _ in range(P * 5): self.irand()

    # 初期化
    def _init_buff(self, seed):
        buff = [0] * P
        buff[0] = seed & RAND_MAX
        for i in range(1, P):
            buff[i] = 1812433253 * (buff[i-1] ^ (buff[i-1] >> 30)) + i
            buff[i] &= RAND_MAX
        self.buff = buff

    def _bad_init_buff(self, seed):
        buff = [0] * P
        buff[0] = seed & RAND_MAX
        for i in range(1, P):
            buff[i] = (69069 * buff[i-1] + 1) & RAND_MAX
        self.buff = buff

    # 整数の一様乱数
    def irand(self):
        x = (self.buff[self.p] + self.buff[self.q]) & RAND_MAX
        self.buff[self.p] = x
        self.p += 1
        if self.p >= P: self.p = 0
        self.q += 1
        if self.q >= P: self.q = 0
        return x

    # 実数の一様乱数
    def random(self):
        return (1.0 / (RAND_MAX + 1.0)) * self.irand()

# test
if __name__ == '__main__':
    a = Random()
    for _ in range(16): print(a.irand())

Algorithms with Python

番外編 : 擬似乱数の検定

はじめに

●擬似乱数の性質

●等確率性の検定

●無規則性の検定

●系列相関係数による検定

●分割表による検定

●Sum Test

●参考文献, URL

●プログラムリスト１

●プログラムリスト２

●プログラムリスト３