正規分布に従う乱数を作る方法(ソースコードあり) ～中心極限定理を利用して～

前に、ボックス・ミュラー法を利用して正規分布に従う乱数を作るって記事を書きましたが、今回は中心極限定理を利用して作ってみます。

正規分布に従う乱数を作る方法として、このブログでは2通りの方法を紹介していて、その内の一つを実装してみたって記事が前のボックス・ミュラー記事になります。で、今回はもう一つの方法、中心極限定理を利用して乱数を作ろうって記事になります。

中心極限定理ってのは、同じ確率分布に従って、かつ独立な確率変数がn個あった場合（\( X_1, X_2, …, X_n \)みたいな感じで）、その平均値は正規分布に従いますよーっていう定理になります。そのときに従う正規分布っていうのは、平均値が真の平均値\( \mu \)で、分散は確率変数が従う確率分布の分散を確率変数の個数で割った値\( \sigma ^ 2 / n \)になるような正規分布になります。

この中心極限定理から考えてやると、同じ一様分布に従う独立な確率変数を複数個用意して、平均値を取ると、その平均値は正規分布に従うはずです。

そこで、正規分布に従う乱数が必要になった場合、乱数を複数個用意しておいて、その平均値を取れば欲しい乱数が得られるということになるはずです。

以上が、元々の考え方になります。というわけで、ソースコードを見ていきましょう。

実装
ソースコードの軽い解説を
注意点とか

実装

こちらがソースコードになります。5つの乱数の平均値を取って、その平均値（正規分布に従う乱数）は5000個発生させています。

#include <iostream>
#include <fstream>
#include <cstdlib>

using namespace std;

#define RAND_NUM 5000       // 作る疑似乱数の個数
#define ORIGINL_RAND_NUM 5  // 元となる乱数の個数(今回は、5つの乱数の平均を取る)

int main()
{
  int i, j;
  double average_rand = 0.0;

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

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

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

  for(i = 0; i < RAND_NUM; i++)
  {
    average_rand = 0.0;

    // 乱数を5回作ってそれぞれを順次足す（つまり、5つの乱数の合計値を計算する）
    for(j = 0; j < ORIGINL_RAND_NUM; j++)
    {
      average_rand += rand() % 10 + 1;
    }

    // ファイルに乱数の平均値を出力
    fout << average_rand / ORIGINL_RAND_NUM;
    fout << endl;
  }

  return 0;
}

このソースを実行してやるとcsvファイルが出力されるかと思います。そのcsvファイルの一番左の列に、正規分布に従う疑似乱数が入力されているはずです。ちなみに、僕の環境ではこんなcsvファイルが出力されたので、元データを知りたい方はそちらをご覧くださいませ～。

それで、下の画像がその実行結果をヒストグラムに直したものになります。グラフの見方は、ボックス・ミュラー記事のときと同じなので、そちらも読まれた方は、そのときと同じように見ていただければと思います。

各棒の上にある数字は横軸は、左から順に、1以上2未満、2以上3未満、…（以下同様）という風になっています。縦軸は出現回数を表しています。なので、各棒の内、一番左のものは1以上2未満の乱数が5個発生したということを表しています。

グラフを見た感じでは確かに正規分布に従ってるような感じがありますねー。

ソースコードの軽い解説を

まぁ、解説する必要もないかと思いますが一応。余裕で理解できたって方は次の節へどうぞ。

最初の方にある次の2文で、乱数の個数を決めています。

#define RAND_NUM 5000       // 作る疑似乱数の個数
#define ORIGINL_RAND_NUM 5  // 元となる乱数の個数(今回は、5つの乱数の平均を取る)

RAND_NUMが正規分布に従う乱数の個数を、ORIGINAL_RAND_NUMが独立な確率変数の個数を表すようにプログラムを組んでいます。なので、このソースでは、5つの確率変数の平均値を5000個作るということになります。ただ、実際にはconst int RAND_NUM = 5000みたいな感じにした方が良いんでしょうが。

で、乱数の平均値を計算したら出力しています。以下の部分です。

    // 乱数を5回作ってそれぞれを順次足す（つまり、5つの乱数の合計値を計算する）
    for(j = 0; j < ORIGINL_RAND_NUM; j++)
    {
      average_rand += rand() % 10 + 1;
    }

    // ファイルに乱数の平均値を出力
    fout << average_rand / ORIGINL_RAND_NUM;
    fout << endl;

average_rand += rand() % 10 + 1;の部分で、1から10までの一様乱数を取得して（0から9までの一様分布の乱数に1を足すので、1から10までの一様乱数になります）average_randに加算しています。