【C++11】<random>で学ぶ確率統計のお話

このネタは部誌に書こうと思ってたけど面倒になった部誌に収まるか分からないからこっちにあげます。平気そうなら部誌に…

C++11で<random>ライブラリが追加されました。このライブラリは乱数生成エンジンと乱数分布器によって構成されています。
・乱数生成エンジン … 一意にランダムな値を生成する
・乱数分布器 … 指定した分布に基づいて値を分布させる

乱数生成エンジンには以下の4種類が用意されています。
線形合同法
メルセンヌツイスター
・Lagged Fibonacci 法
・ハードウェアエントロピーソース
このうち、ハードウェアエントロピーソースは遅い代わりに予測不能な真の変数を生成します。線形合同法メルセンヌツイスター法などは、ハードウェアエントロピーソースよりは高速なもののある程度の規則性が存在し、初期シード(rand関数におけるsrand関数)の値によって値が変化します。
<cstdlib>(<stdlib.h>)のrand関数はランダム性が低く生成乱数の周期が短いなどの問題を抱えています*1が、これらの方法で比較的容易に回避することが可能です。

前述したようにハードウェアエントロピーソースは予測不可能ですが遅いです。よってこの乱数エンジンをメインに据えて使用することはあまり好ましくありません。よって、メルセンヌツイスター法などの初期シードの生成に1回使用するのが無難といえそうです。具体的には以下のようになります。

#include <iostream>
#include <random>
using namespace std;

int main()
{
	// random_deviceは予測不可能な乱数生成エンジン
	random_device Hes;
	// mt19937はメルセンヌツイスター法のdefineマクロ
	mt19937 Mt(Hes());
	
	for(int i=0; i < 100; i++){
		cout << Mt() << endl; // 出力
	}
	return 0;
}

上記コードのように

  1. random_deviceクラスで変数を確保
  2. メルセンヌツイスター法などの高速な乱数生成エンジンクラスで変数を確保し、初期値に1.で生成したクラスの乱数を使用
  3. 2.で生成したクラスの乱数を使用する

のが主な使用方法になると思います。勿論、速度を気にしないような状況であれば
random_deviceで生成した乱数をそのまま使ってしまってもいいと言えそうです。



乱数を分布させたいときに使用するのが乱数分布器です。
具体的には、1から3までのランダムな値が欲しいとか、1から100までの間の値で10が一番出るような値が欲しいといった場合に使用できます。
前者はともかく、後者の実装はなかなか大変だと思います。このような問題に簡単に対応できるようにできるのが乱数分布器というわけです。
では、実際の使用例を見てください。六面サイコロの目が出るように乱数を分布させた例です。

#include <iostream>
#include <random>
using namespace std;

int main()
{
	// random_deviceは予測不可能な乱数生成エンジン
	random_device Hes;
	// mt19937はメルセンヌツイスター法のdefineマクロ
	mt19937 Mt(Hes());
	// 整数分布なのでuniform_int_distributionを使用する
	// 実数全体の分布にしたいときはuniform_real_distribution<double>を使用する
	uniform_int_distribution<int> Dice(1, 6);

	for(int i=0; i < 100; i++){
		cout << Dice(Mt) << endl; // 出力
	}
	return 0;
}

見ていただければ大体はわかると思います。
ここで、uniform_int_distributionの初期値に渡しているのは一様分布の最小値と最大値(この場合は1と6)です。
そして、乱数を分布させる際に乱数生成エンジンの変数をカッコ無しで渡します。
他にも乱数分布器は多数存在します。その一例をこれから説明していきたいと思います。
詳細はこちら(Pseudo-random number generation - cppreference.com)を確認してください。



さて、ここからが本題です。randomライブラリで使用できる乱数生成エンジンと乱数分布器を使用して、確率の問題を考えていきましょう。

Q1. 六面サイコロを5回振った時、2回6の値が出る確率を求めよ。

二項分布の問題です。六面サイコロで6の値が出る確率は1/6です。よってp=1/6となります。
数式に当てはめることで確率を求めるのですが、今回は試行回数を重ねることでどの程度の確率になるかを算出します。
以下が算出を行うコードです。

#include <iostream>
#include <random>
#include <algorithm>
#include <vector>
using namespace std;

int main()
{
	unsigned int Trials = 10000; // 試行回数

	// random_deviceは予測不可能な乱数生成エンジン
	random_device Hes;
	// mt19937はメルセンヌツイスター法のdefineマクロ
	mt19937 Mt(Hes());
	// 試行回数5回, 成功確率1/6の二項分布
	binomial_distribution<int> Bin(5, (double)1/6);

	// vectorに入力
	vector<int> Vec;
	for(unsigned int i=0; i < Trials; i++){
		Vec.push_back(Bin(Mt));
	}

	// 要素数を出力
	cout << "試行回数: " << Trials << endl;
	for(int i = 0; i <= 5; i++){
		int Ct = count_if(Vec.begin(), Vec.end(), [i](int t){return i == t;});
		cout << i << ": " <<  Ct << " / " << Trials
		<< " = " << (double)Ct/Trials << endl;
	}
	return 0;
}

ちゃっかりラムダ式使ってごめんなさい
二項分布はbinomial_distributionを使用します。
初期化時の変数は試行回数と起こる確率です。使い方はサイコロの例とほぼ同じです。
ここで、初期化時の試行回数とコード内での試行回数とは別物であることに注意してください。
このコード内のTrials変数の値を書き換えることで試行回数を変えることができます。
試行回数10000回の時の結果は以下の通りになりました(実行するたびに結果は変わります)。

試行回数: 10000
0: 3902 / 10000 = 0.3902
1: 4112 / 10000 = 0.4112
2: 1614 / 10000 = 0.1614
3: 333 / 10000 = 0.0333
4: 38 / 10000 = 0.0038
5: 1 / 10000 = 0.0001

また、試行回数100000回の時の結果は以下の通りになりました(実行するたびに結果は変わります)。

試行回数: 100000
0: 40189 / 100000 = 0.40189
1: 40022 / 100000 = 0.40022
2: 16152 / 100000 = 0.16152
3: 3344 / 100000 = 0.03344
4: 283 / 100000 = 0.00283
5: 10 / 100000 = 0.0001

ここから、似たような値(0.161...)になっていることがわかります。
理論上はどのような値になるかを計算すると、p(2) = 0.1607...となります。
誤差は0.4%ほどですのでほぼ無いと言えるでしょう。
答えは16%ほどになります。


Q2. とある高専では毎年平均して40人が留年する。留年する確率はポアソン分布に従うとして
50人以上が留年する確率を求めよ。

ポアソン分布を使用するにはpoisson_distributionを使用します。引数には平均数を指定します。
他は基本的に同じです。結果の表示法を変えています。

#include <iostream>
#include <random>
#include <algorithm>
#include <vector>
using namespace std;

int main()
{
	unsigned int Trials = 10000; // 試行回数

	// random_deviceは予測不可能な乱数生成エンジン
	random_device Hes;
	// mt19937はメルセンヌツイスター法のdefineマクロ
	mt19937 Mt(Hes());
	// 平均数40のポアソン分布
	poisson_distribution<int> Pois(40);

	// vectorに入力
	vector<int> Vec;
	for(unsigned int i=0; i < Trials; i++){
		Vec.push_back(Pois(Mt));
	}

	// 要素数を出力
	cout << "試行回数: " << Trials << endl;
	int Ct = count_if(Vec.begin(), Vec.end(), [](int t){return t >= 50;});
	cout << "50以上: " << Ct << " / " << Trials << " = " << (double)Ct/Trials << endl;
	return 0;
}

実行結果はそれぞれ以下のようになります(実行する度に(ry

試行回数: 10000
50以上: 713 / 10000 = 0.0713
試行回数: 100000
50以上: 7027 / 100000 = 0.07027

この場合も似たような値(0.07...)となりました。理論値は0.0703...ですのでやはり殆ど違いはないと言えます。つまり答えは7%ほどになります。


Q3. 高専生200人が同じ数学の試験を受けたところ、平均点75点、標準偏差18点の正規分布に従った。
ある生徒が赤点である確率を求めなさい。赤点は60点未満とする。

正規分布はnormal_distributionを使用します。…あとは説明しなくても十分ですね?
コードは以下のようになります。

#include <iostream>
#include <random>
#include <algorithm>
#include <vector>
using namespace std;

int main()
{
	unsigned int Trials = 10000; // 試行回数

	// random_deviceは予測不可能な乱数生成エンジン
	random_device Hes;
	// mt19937はメルセンヌツイスター法のdefineマクロ
	mt19937 Mt(Hes());
	// 平均75, 標準偏差18の正規分布
	normal_distribution<double> Norm(75, 18);

	// vectorに入力
	vector<double> Vec;
	for(unsigned int i=0; i < Trials; i++){
		Vec.push_back(Norm(Mt));
	}

	// 要素数を出力
	cout << "試行回数: " << Trials << endl;
	int Ct = count_if(Vec.begin(), Vec.end(), [](int t){return t < 60;});
	cout << "赤点: " << Ct << " / " << Trials << " = " << (double)Ct/Trials << endl;
	return 0;
}

注意点として、normal_distribution<double>です。intにすると固まります。

実行結果はそれぞれ以下のようになります(実行(ry

試行回数: 10000
赤点: 2007 / 10000 = 0.2007
試行回数: 100000
赤点: 20183 / 100000 = 0.20183

それなりに近い値(0.20...)となりました。理論値は0.2023...ですのであまり違いはないと言えます。
答えは20%ほどです。意外と確率高い





// 雑記
ゲーセン行かないと暇すぎてしょうがない。
ゲーセンが生きる活力。アカンこれダメな人間だ

*1:このサイトに詳細が詳しく乗っています