C++ - 一様乱数の一様性検定(カイ2乗検定)!

Updated:


少し前に、線形合同法を使用して一様乱数を生成する C++ によるアルゴリズムについて紹介しました。

今回は、それらの生成した一様乱数が本当に一様かどうかを「カイ2乗検定」で検証してみました。

「カイ2乗検定」とは、今回のケースに合わせて簡単に言うと、

CHI2_KENTEI

この値を検証してみるということになります。 「カイ2乗検定」の詳しい事は、別途サイトや統計関係の書籍をお調べください。

以下、C++ によるサンプルソースです。

0. 前提条件

  • Cygwin 1.7.15
  • g++ (GCC) 4.7.1

1. C++ ソース作成

今回作成した C++ ソースは以下の通りです。
(C++ なのでオブジェクト指向な作りにしている)

File: chi_2_rndnum.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
/*********************************************
 * 線形合同法による一様乱数の一様性検証
 *********************************************/
#include <iostream>  // for cout
#include <math.h>    // for pow
#include <stdio.h>   // for printf

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    // 宣言
    static const int a = 1103515245;  // 乗数
    static const int c = 12345;       // 加数
    static const unsigned int m = pow(2, 31);  // 法(2の31乗)
    static const int n = 1000;        // 発生させる乱数の個数
    int r = 12345;                    // 乱数の種の初期値
    static const int m_max = 10;      // 整数乱数の範囲
    double f = n / m_max;             // 期待値
    double s = 40.0 / f;              // ヒストグラム用スケール
    int rank;                         // 整数乱数
    int hist[m_max+1];                // 件数格納用配列
    double e = 0.0;                   // カイ2乗検定初期値
    int i, j;                         // ループインデックス

    public:
        // コンストラクタ
        Calc();
        // 一様乱数生成
        void generateRndnum();
        // 1 ~ 10 の整数乱数
        int rnd();
        // 結果表示
        void display();
};

/*
 * コンストラクタ
 */
Calc::Calc()
{
    // 件数格納用配列初期化
    for (i = 1; i <= m_max; i++) {
        hist[i] = 0;
    }
}
/*
 * 一様乱数生成
 */
void Calc::generateRndnum()
{
    for (i = 0; i < n; i++) {
        rank = rnd();  // 1 ~ 10 の整数乱数
        hist[rank]++;  // 1 ~ 10 別にカウント
    }
}

/*
 * 1 ~ 10 の整数乱数
 */
int Calc::rnd()
{
    // 0 ~ 2の31乗 未満の整数乱数
    r = (a * r + c) % m;

    // 0 ~ 1 未満の実数乱数に 10 を乗じて 1 を加えることで
    // 1 ~ 10 の整数乱数にする
    return m_max * (r / (m - 0.9)) + 1;
}

/*
 * 結果表示
 */
void Calc::display()
{
    for (i = 1; i <= m_max; i++) {
        // 件数表示
        printf("%3d:%3d ", i, hist[i]);

        // ヒストグラム表示
        for (j = 0; j < hist[i] * s; j++)
            printf("*");
        printf("\n");

        // カイ2乗検定
        e = e + (double)(hist[i] - f) * (hist[i] - f) / f;
    }
    // カイ2乗検定値表示
    printf("χ2 = %f\n", e);
}
/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 一様乱数生成
        objCalc.generateRndnum();

        // 結果表示
        objCalc.display();
    }
    catch (...) {
        cout << "例外発生!" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

2. C++ ソースコンパイル

$ g++ chi_2_rndnum.cpp -o chi_2_rndnum -std=c++11

(”-std=c++11” は警告を出力しない為のおまじない)
何も出力されなければ成功です。

3. 実行

$ ./chi_2_rndnum
  1:105 *******************************************
  2: 89 ************************************
  3:109 ********************************************
  4: 99 ****************************************
  5:100 *****************************************
  6: 97 ***************************************
  7:103 ******************************************
  8:116 ***********************************************
  9: 89 ************************************
 10: 93 **************************************
χ2 = 6.720000

4. 判定

実行した結果が一様であるかどうかですが、ヒストグラムではそんなに大きなバラツキは確認できません。
そして、カイ2乗統計量は 6.72 という値になっています。

今回は 1 から 10 までの 10 個の整数で検証しましたので、カイ2乗検定でいうところの自由度が 9 ということになります。
統計関係書物等で「カイ2乗分布表」を調べてみると 、自由度 9、危険率(αパーセント点) 0.01 の値は 21.660 となっています。
明らかに 6.72 < 21.660 を満たしていますから、発生した乱数は危険率 1% で一様に分布していると判定できます。


乱数生成回数をもっと増やしたり、乱数生成時の定数を変更してみたりすると、もっと一様になるのではないでしょうか?

以上。





 

Sponsored Link

 

Comments