mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

C++ - 正規乱数(ボックス=ミューラー法)!

[ プログラミング, 数学 ] [ C言語 ]

こんばんは。

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

今回は、正規乱数を発生させて実際に正規分布になっているかを検証してみました。

まず、「正規乱数」とは「正規分布」を持つ「乱数」のことです。

「正規分布」は以下のようなグラフになり、

SEIKI_BUNPU_1

以下のような式で定義されます。

SEIKI_BUNPU_2

この「正規乱数」の生成方法として、今回は「ボックス=ミューラー法(Box-Muller transform」を使用します。 「ボックス=ミューラー法」は、0 より大きく 1 以下、すなわち (0,1] の一様乱数を正規乱数(正規分布を持つ乱数)に変換する方法で、計算式は以下のようになります。 (0, 1] の2個1組の一様乱数で2個の正規乱数を生成できます。

BOX_MULLER_1

※正規分布・正規乱数についての詳細は、別途お調べください。

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

記録

0. 前提条件

  • Cygwin 1.7.15
  • g++ (GCC) 4.7.1

1. C++ ソース作成

今回作成した C++ ソースは以下の通りです。 C++ なのでオブジェクト指向な作りにしています。 また、元となる一様乱数の生成は関数を使用しています。 以下の例では、平均を 10, 標準偏差を 2.5 とし、0 から 20 までの整数乱数に変換して検証しています。

【 ファイル名: RndnumBoxMuller.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
118
119
120
/*********************************************
 * ボックス=ミューラー法法による正規乱数生成
 *********************************************/
#include <iostream>  // for cout
#include <math.h>    // for cos, sin
#include <cstdlib>   // for rand()
#include <stdio.h>   // for printf()

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    // 各種定数
    static const int    m  = 10;            // 平均
    static const double s  = 2.5;           // 標準偏差
    static const int    n  = 10000;         // 発生させる乱数の個数
    static const double pi = 3.1415926535;  // 円周率
    // 各種変数
    double r1, r2;       // 一様乱数
    double x, y;         // 正規乱数
    int    hist[m*2+1];  // 件数格納用配列
    double scale;        // ヒストグラム用スケール
    int    i, j;         // ループインデックス

    public:
        // コンストラクタ
        Calc();
        // 正規乱数生成
        void generateRndnum();
        // 整数乱数
        void rnd(double s, double m, double *x, double *y);
        // 結果表示
        void display();
};

/*
 * コンストラクタ
 */
Calc::Calc()
{
    // 件数格納用配列初期化
    for (i = 0; i <= m * 2; i++) {
        hist[i] = 0;
    }
    // ヒストグラム用スケール
    scale = n / 100.0;
}
/*
 * 正規乱数生成
 */
void Calc::generateRndnum()
{
    // n 回乱数生成処理を繰り返す
    for (i = 0; i < n; i++) {
        // 整数乱数を2個生成
        rnd(s, m, &x, &y);
        // 生成された2個の整数乱数を件数格納用配列に格納
        hist[(int)x]++;
        hist[(int)y]++;
    }
}

/*
 * 整数乱数
 */
void Calc::rnd(double s, double m, double *x, double *y)
{
    // 一様乱数
    // ( 32Bitマシンでの (0,1] の実数乱数 )
    r1 = rand() / 2147483647.1;
    r2 = rand() / 2147483647.1;

    // 正規乱数計算
    *x = s * sqrt(-2 * log(r1)) * cos(2 * pi * r2) + m;
    *y = s * sqrt(-2 * log(r1)) * sin(2 * pi * r2) + m;
}

/*
 * 結果表示
 */
void Calc::display()
{
    // 0 ~ m * 2 を表示
    for (i = 0; i <= m * 2; i++) {
        // 件数表示
        printf("%3d:%4d | ", i, hist[i]);

        // ヒストグラム表示
        for (j = 1; j <= hist[i] / scale; j++)
            printf("*");
        printf("\n");
    }
}
/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 正規乱数生成
        objCalc.generateRndnum();

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

    // 正常終了
    return 0;
}

今回は32ビットOSでの作業でしたので、(0,1] の一様乱数を生成する際に、rand() 関数の値を int の「最大値+0.1」で除しています。 OSが32ビット以外なら int の最大値はこれとは異なります。

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

以下のコマンドで C++ ソースをコンパイルする。

1
$ g++ RndnumBoxMuller.cpp -o RndnumBoxMuller

何も出力されなければ成功です。

3. 実行

実際に実行して検証してみる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
$ ./RndnumBoxMuller
  0:   4 |
  1:  12 |
  2:  43 |
  3: 107 | *
  4: 256 | **
  5: 626 | ******
  6:1207 | ************
  7:1939 | *******************
  8:2597 | *************************
  9:3040 | ******************************
 10:3153 | *******************************
 11:2744 | ***************************
 12:1964 | *******************
 13:1201 | ************
 14: 617 | ******
 15: 320 | ***
 16: 121 | *
 17:  33 |
 18:   9 |
 19:   3 |
 20:   1 |

4. 判定

出力されたヒストグラムを確認してみると、綺麗な山になっているので正規乱数が正規分布になっていると言えるでしょう。


乱数生成回数をもっと増やしたりしてみてもよいでしょう。

以上。

Comments