C++ - 一様乱数の一様性検定（カイ２乗検定）！
Updated:
少し前に、線形合同法を使用して一様乱数を生成する C++ によるアルゴリズムについて紹介しました。
今回は、それらの生成した一様乱数が本当に一様かどうかを「カイ２乗検定」で検証してみました。
「カイ２乗検定」とは、今回のケースに合わせて簡単に言うと、
この値を検証してみるということになります。 「カイ２乗検定」の詳しい事は、別途サイトや統計関係の書籍をお調べください。
以下、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; // カイ２乗検定初期値
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");
// カイ２乗検定
e = e + (double)(hist[i] - f) * (hist[i] - f) / f;
}
// カイ２乗検定値表示
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. 判定
実行した結果が一様であるかどうかですが、ヒストグラムではそんなに大きなバラツキは確認できません。
そして、カイ２乗統計量は 6.72 という値になっています。
今回は 1 から 10 までの 10 個の整数で検証しましたので、カイ２乗検定でいうところの自由度が 9 ということになります。
統計関係書物等で「カイ２乗分布表」を調べてみると 、自由度 9、危険率（αパーセント点） 0.01 の値は 21.660 となっています。
明らかに 6.72 < 21.660 を満たしていますから、発生した乱数は危険率 1% で一様に分布していると判定できます。
乱数生成回数をもっと増やしたり、乱数生成時の定数を変更してみたりすると、もっと一様になるのではないでしょうか？
以上。
Comments