またまた正規分布に従う乱数を作る方法です。前までの方法は三角関数を使ったり、中心極限定理を利用してたりしたので、C言語とか、Python、PHPみたいに、言語に関係なく使える方法になっていました(僕はC++が一番得意なのでC++で作っていましたが)。
ちなみに、前までの方法は下の記事のような感じです。気になったらご自由にどうぞ~。
すっごい大雑把な説明:正規分布に従う乱数の作成方法
方法その1(ボックス・ミュラー法):C++で正規分布に従う乱数を作る ~ソースコード付き~
方法その2(中心極限定理):正規分布に従う乱数を作る方法(ソースコードあり) ~中心極限定理を利用して~
ただ、C++を使う場合(より正確にはC++11以降のC++を使う場合です。まぁ、今はほとんどがC++11以降のコンパイラでしょうから問題ないかと思いますが)、色んな種類の乱数を作るためのライブラリ(random)が用意されていて、そいつを使った方が手っ取り早いんですな。
ってことで、今回はC++を使う場合に限りますが、正規分布に従う乱数はこんな方法でも作れるよーっていう紹介になります。
randomをインクルード!(ソースコードの節です)
まぁ、こちらのリファレンスを見てもらったらそれでいいような気もしますが(元も子もない笑)、説明していきます。
ソースコードは下のような感じです。説明はソースコードの下の方に書いています。
#include <iostream>
#include <fstream>
#include <random>
using namespace std;
const unsigned int RAND_NUM = 5000; // 作る疑似乱数の個数
int main()
{
random_device seed_gen; // 環境依存ではあるが、予測不能な乱数を作れる(らしい)。
mt19937 engine(seed_gen()); // メルセンヌ・ツイスター法による擬似乱数生成器初期化
ofstream fout("PSEUDO_RANDOM_NUMBER.csv"); // 出力ファイルを開く
if(!fout){ // ファイルをストリームにリンクできない場合、偽(=0)が返される。
cout << "出力ファイルを開けません\n";
return 1;
}
// 平均1.0、標準偏差0.5の正規分布で分布させる
normal_distribution<> dist(1.0, 0.5);
fout << "rand1" << "," << "rand2" << endl;
for (unsigned int i = 0; i < RAND_NUM; i++) {
// 予測不能な乱数を生成
unsigned int rand1 = seed_gen();
// 正規分布に従う乱数を生成
double rand2 = dist(engine);
fout << rand1 << "," << rand2 << endl;
}
return 0;
}
今回のプログラムでやっていることは前回までの乱数記事と同じで、5000個の乱数を発生させてファイルに出力しています。このプログラムを実行すると、実行ファイルと同じディレクトリに、PSEUDO_RANDOM_NUMBER.csvという名前のファイルが作られているかと思います(同じ名前のファイルが既にあった場合は上書きで更新されているかと思います)。
このライブラリを使うときに必要な手順は3つあって、
1.初期化
2.乱数の従う分布を決める
3.実際に乱数を取得する
という3ステップになります。詳しくは次の節で説明していきます。
ライブラリの使い方
というわけで、上の節でも軽く触れましたが、乱数の取得までは3つのステップに分解できるので、そのステップごとに説明していきます。
※いろいろ書いていますが、要するにこの節では、
random_device seed_gen;
mt19937 engine(seed_gen());
normal_distribution<> dist(1.0, 0.5);
double rand = dist(engine);
で正規分布に従う乱数を作れますよーってことが書いてあります。
1.初期化
main関数内の
random_device seed_gen; // 環境依存ではあるが、予測不能な乱数を作れる(らしい)。
mt19937 engine(seed_gen()); // メルセンヌ・ツイスター法による擬似乱数生成器初期化
という部分が初期化に当たります。
ここでやっていることはrandom_device型のseed_genというオブジェクトの作成と、mt19937型のengineというオブジェクトの作成の2つです。
random_device型というのは、予測不能な乱数生成器の型のことです。なので、
random_device seed_gen;
は、seed_genという名前の乱数生成器を作っているという意味になります。そして、乱数が欲しいときにはseed_gen()とすれば環境依存ではありますが、予測不能な乱数が得られます。なので、
mt19937 engine(seed_gen());
部分では、mt19937型のengineという名前のオブジェクトをseed_gen()で得られた乱数で初期化しているという意味になります。で、mt19937型というのは、メルセンヌ・ツイスターによる乱数を発生させる乱数生成器の型になります。※メルセンヌ・ツイスター法は乱数生成法の一種です。
つまり、まとめると、engineという名前のメルセンヌ・ツイスター法による乱数生成器を予測不能な乱数で初期化しているということになります。
ちなみに、もしも、engineの引数を定数にしていたら、毎回同じ数列が返されることになってしまいます。確認してみたい方は、上のプログラムで「mt19937 engine(seed_gen());」を「mt19937 engine(10);」みたいに変更して実行すると、何回実行しても、csvファイルの2列目はまったく変化しないはずです。
何回実行しても変わらないというのは乱数としては問題なので、それを防ぐためにseed_gen()で得られる再現性のない乱数を使ってるってわけですな。
ただし、「再現性のない乱数」とはいっても、このseed_gen()で得られる乱数には注意が必要で、環境依存なので、環境によって予測不能だったりそうでなかったりします。ちなみに、GCCの場合は9.2以降はまったく予測のできない乱数(=再現性のない乱数)が得られるんですが、それ以前の場合は再現性のある乱数(つまり、乱数としては欠陥品)になってしまいますのでご注意を。Visual Studioだと予測できない乱数になるっぽいです。
2.乱数の従う分布を決める
ステップ1までで乱数生成器を作ることはできた(engine)ものの、このままでは一体どんな確率分布に従うのかが定義できていません。なので、このステップで一様分布なのか、二項分布なのか、正規分布なのかというように、engineの作り出す乱数が従う確率分布を定義してやります。それをやってるのが、下に書いたソースコードになります。
normal_distribution<> dist1(1.0, 0.5);
今回は正規分布を作りたいのでこのようにしましたが、必要な分布に従って変える必要があります。これで、あとは乱数生成器を引数として渡してやれば正規分布(今回は平均値が1.0, 標準偏が0.5の正規分布)に従う乱数が得られます。
3.乱数を取得する
ステップ2でも軽く触れた通り、乱数生成器を渡してやれば乱数が得られます。具体的には次のような感じになります。
// 正規分布に従う乱数を生成
double rand2 = dist(engine);
rand2には正規分布に従う乱数が入っています。
続く
調子に乗って書いてたら、予想以上に分量が多くなってしまったので、プログラムの実行結果がどうなったかは次の記事で書きます。