C++で正規分布に従う乱数を作る～ソースコード付き～

以前に、正規分布に従う疑似乱数を作る方法をご紹介しましたが、そういやソースコードは公開してなかったなと思いまして、C++で正規分布を作るプログラムを組んでみました（必要な部分は後でおさらいとして触れます）。

以前の記事では、正規分布に従う疑似乱数（面倒になったら、単に乱数って書きます。一様乱数なのか正規分布に従うのかは文脈から上手いこと判断してください）の作り方として、「三角関数を利用して計算する方法（ボックス・ミュラー法）」と「中心極限定理を利用する方法」の2通りを紹介していたのですが、今回は前者（ボックス・ミュラー法）についてのみ書いていきます。後者（中心極限定理を利用する方法）についてはまた今度書きます。

なんせ、記事がやたらと長くなっちゃいそうだったもんで。

というわけで、この記事ではボックス・ミュラー法で”正規分布に従う”疑似乱数を作るプログラムを書いていきます。

2つの方法
実装（ボックス・ミュラー法）
一応、ソースコードの説明
実行結果をヒストグラムにしてみた

2つの方法

どうやったら”正規分布に従う”疑似乱数が作れるかという、前回のおさらいみたいな節です。今回はボックス・ミュラー法のプログラムしか書いていませんが、一応、もう一つの方法も作り方の要点にだけ触れておきます。

この記事で扱っているボックス・ミュラー法での乱数の作り方です。こちらは、要するに独立な乱数を２つ（U1とU2）取っておいて、それを次の計算式に代入すれば”正規分布に従う”疑似乱数が作れるよーって方法です。

$$
X_1 = \sqrt{ -2 log_e U_1 } cos(2 \pi U_2) \\
X_2 = \sqrt{ -2 log_e U_1 } sin(2 \pi U_2)
$$

この場合はX1かX2か、どちらか片方だけを使っておけば良いかと思います。

2つ目は、中心極限定理を利用する方法です。乱数を複数個取っておいて、その平均を取ってやれば乱数が得られるって方法になります。例えば、乱数を5個（R1、R2、R3、R4、R5）取るということにしていたら、乱数Rは次のように表せるということになります。

$$ R = \frac{ R1 + R2 + R3 + R4 + R5 } { 5 } $$

さて、おさらいが終わったところで実装していきましょー。

実装（ボックス・ミュラー法）

ソースコードはこんな感じになります。

#define _USE_MATH_DEFINES
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>

using namespace std;

#define RAND_NUM 5000   // 作る乱数の個数

int main()
{
  int i;
  double rand1, rand2;   // 元となる乱数
  double pseudo_rand1, pseudo_rand2;

  ofstream fout("PSEUDO_RANDOM_NUMBER(BOX_MULLER).csv");    // 出力ファイルを開く

  if(!fout){                // ファイルをストリームにリンクできない場合、偽(=0)が返される。
    cout << "出力ファイルを開けません\n";
    return 1;
  }

  srand(10);  // 乱数のシード値を決める

  // ラベルを入力しておく
  fout << "X1,X2" << endl;

  // 乱数をRAND_NUM個作る
  for(i = 0; i < RAND_NUM; i++)
  {
    // まずは、乱数を2つ用意しておく
    // ただし、0以上, 1以下の実数を取るようにするために、乱数の最大値RAND_MAXで割る。
    rand1 = (double)rand() / RAND_MAX;    // rand()の返り値はint型なので、double型にキャストする
    rand2 = (double)rand() / RAND_MAX;    // でないと、RAND_MAXで割った瞬間に0になる。

    // 次に、本命の乱数(正規分布に従う疑似乱数)を計算する。
    pseudo_rand1 = sqrt(-2 * log(rand1)) * cos(2 * M_PI * rand2);
    pseudo_rand2 = sqrt(-2 * log(rand1)) * sin(2 * M_PI * rand2);

    fout << pseudo_rand1 << "," << pseudo_rand2 << endl;
  }

  return 0;
}

こちらのソースをビルドして作った実行ファイルを実行していただければ実行ファイルと同じディレクトリにPSEUDO_RANDOM_NUMBER.csvというファイルが作られます。そちらをExcelなどで開いていただくと、1行目のA列にはX1と、B列にはX2と入力されていて、それぞれの下に大量に数字が記録されているはずです。

その数字がボックス・ミュラー法で作り出した乱数になります。

A列に記録された数字が下に再掲した数式で言うところのX1、つまり、cosを使って作った乱数で、B列の数字がsinを使って作られたX2となっています。

$$
X_1 = \sqrt{ -2 log_e U_1 } cos(2 \pi U_2) \\
X_2 = \sqrt{ -2 log_e U_1 } sin(2 \pi U_2)
$$

一応、ソースコードの説明

やってることとしてはざっくりと5つで、

出力専用のファイル(.csv)を開く
乱数を2つ取得する
取得した乱数を上の式を使って正規分布に従う乱数に変換する
変換後の乱数をcsvファイルに出力する
2から4を5000回繰り返す

っていう5つになります。ソースを見れば分かるよーって方は次の”実行結果をヒストグラムにしてみた”の節まで飛んでいただければと思います。

この記事のソースでは、乱数を取得する（2.）で使ってる関数がC言語の標準ライブラリとして用意されてるrand()関数なので、本当に乱数になってるのかどうかは微妙なところですが。

正規分布に従う乱数に変換するという操作(3.)は

    pseudo_rand1 = sqrt(-2 * log(rand1)) * cos(2 * M_PI * rand2);
    pseudo_rand2 = sqrt(-2 * log(rand1)) * sin(2 * M_PI * rand2);

という部分でやっております。このpseudo_rand1、pseudo_rand2というのが、それぞれ上の式で言うところのX1、X2に当たって、rand1がU1に、rand2がU2に対応しています。

つまり、欲しい乱数はpseudo_rand1とpseudo_rand2に保存されてるってことですな。

で、その次の

    fout << pseudo_rand1 << "," << pseudo_rand2 << endl;

という部分で保存されてる乱数をファイルに保存しています。今回はcsvファイルに変換した乱数を保存していくので、pseudo_rand1とpseudo_rand2の間に,（コンマ）を挟むことで、Excelで開いたときにそれぞれの数値が別のセルに記録されているようにします。

というわけで、実行するとどうなったのかということを次の節で。

実行結果をヒストグラムにしてみた

上のプログラムの実行結果をそのまま表形式で張り付けても良かったんですが、分かりにくいだろうってことで、ブログにはヒストグラムだけを載せときます。

表形式ではどうなったのか興味があるって方は、一応こちらに出力されたcsvファイルを上げとくんで、ご自由にダウンロードしてご覧ください。

そのファイルのA列とB列の数値をヒストグラムに直してやると下のようなグラフになります。どちらも横軸に乱数の範囲を、縦軸にその範囲の乱数が発生した回数を取っています。

そして、ラベルが見づらくて申し訳ないんですが、左から順に、-4.0以上-3.5未満、-3.5以上-3.0未満というように、0.5刻みの範囲を示すようにしています。各棒の上にあるのが具体的な発生数になっています。

上の方のグラフ（図1）がX1（ソースコード中のpseudo_rand1）、下の方のグラフ（図2）がX2（pseudo_rand2）をヒストグラム化したものになります。

図1．cosを使って作った正規分布に従う疑似乱数X1(pseudo_rand1)

図2．sinを使って作った正規分布に従う疑似乱数X2(pseudo_rand2)

なんとなく、正規分布に従ってそうな感じがありますねー。階級の幅（刻み幅のこと）を0.5じゃなくて、0.25とか0.1みたいにもっと小さい値にすれば、もっと正規分布っぽくなりそうな気もしますが、とりあえず0.5でやめときました。

まぁ、上でもチラッと触れた通り、元の乱数がC言語の標準ライブラリに装備されたrand()関数なんで、本当に乱数になってるのかどうかはかなり怪しいところではあります。それに、今回の記事では周期とか連続性とかの問題についてはまったく検討してないんで、”この記事だけから”正規分布に従う乱数が欲しけりゃボックス・ミュラーだ！みたいな結論を導くのは危なそうっすねー。

ただまぁ、少なくとも正規分布に従う乱数「っぽく見える」乱数が欲しいときには使えるだろうってことは言えるだろうと思います。