平均値と期待値
C++/Linuxシミュレーションプログラム
離散値
をとる確率変数
の期待値は
で与えられる。
サイコロ投げの場合、目の数の期待値は
となる。いま、N回サイコロを投げたときiの目が
回出たとすると平均値Mは
で与えられる。Nを大きくしていくと
となるので、Mは
に近づく。すなわち、
である。
式（１）を確認するシミュレーションプログラムを作成した。このプログラムを実行すると図１のようにサイコロ投げの回数（図１ではiで表されている）が大きくなるにつれて平均値meanは期待値3.5に近づくことが確認できる。ここで、PiはFi/N、すなわち目の数iの出た相対頻度である。
図１
プログラムのソースコードはリスト１に示すものである。プログラムのソースコードファイルは圧縮ファイルSimuExpect.zipとして用意している。ファイル名SimuExpect.zipをクリックしてダウンロードできる。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト１ MeanAndExpectation.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <iostream>
#include <iomanip>
#include "myrn.h" // 乱数生成用クラス型 WHrn のヘッダーファイル
using namespace std;
using namespace myrn;
//
// 0 < v < 1 の範囲の値を 0 ~ 5 の範囲の整数値に変換
//
int Cat( double v ){
return int(v * 6.0);
}
//
// Fiに格納された頻度とその比率の表示
// 平均値は、0 ~ 5 を 1 ~ 6 の値にシフトさせ、サイコロの目の数とみなして計算
//
void disp_data( int * Fi, int i ){
cout << endl << "i = " << i << endl;
cout << left << setw(5) << "Fi =";
for (int k = 0; k < 6; k++)
cout << right << setw(10) << Fi[k];
cout << endl << left << setw(5) << "Pi =";
for (int k = 0; k < 6; k++)
cout << right << setw(10) << fixed << setprecision(5) << (Fi[k] / double(i));
cout << endl;
double sum = 0.0;
for (int k = 0; k < 6; k++)
sum += (k+1.0)*Fi[k];
cout << "mean = " << fixed << setprecision(5) << (sum / i) << endl;
}
int main()
{
WHrn rn;
int Fi[] = {0, 0, 0, 0, 0, 0};
double Pi[6];
const int MaxN = 1000000;
for (int i = 1; i <= MaxN; i++){
Fi[Cat(rn.uni())]++; // 0 ~ 5 の整数をランダムに生成し、対応するFiの要素の値を １ 増やす
//
// 途中経過の表示
switch(i){
case 10: disp_data( Fi, i ); break;
case 20: disp_data( Fi, i ); break;
case 100: disp_data( Fi, i ); break;
case 1000: disp_data( Fi, i ); break;
case 10000: disp_data( Fi, i ); break;
default: break;
}
}
disp_data( Fi, MaxN );
return 0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
乱数生成のクラス型WHrnはリスト２のヘッダーファイルmyrn.hおよびリスト３のソースコードファイルmyrn.cppで与えられている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト２ myrn.h ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#ifndef MYRN_CK
#define MYRN_CK
namespace myrn {
// 32ビット符号無整数
#define myuint unsigned int
//
// Wichmann/Hill generator...周期 = 約1.0e12
//
class WHrn {
public:
myuint x;
myuint y;
myuint z;
double u;
WHrn(): x(11111), y(22222), z(33333) {}
virtual ~WHrn(){}
WHrn(myuint ix, myuint iy, myuint iz):
x(ix), y(iy), z(iz) {}
double uni(); // 0 < uni() < 1
};
}
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト３ myrn.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include "myrn.h"
#include <cmath>
using namespace std;
namespace myrn {
// 0 < uni() < 1
double WHrn::uni(){
do {
x = (171 * x) % 30269;
y = (172 * y) % 30307;
z = (170 * z) % 30323;
u = (x/30269.0) + (y/30307.0) + (z/30323.0);
u -= trunc(u);
} while (u <= 0.0);
return u;
}
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
リスト１のプログラムを実行すると図１の結果が得られるが、出力の上の方が画面からはみ出してしまって表示されていないので、表示出力をコピーしたものをリスト４に示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト４ 実行結果の出力 ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
okamoto@okamoto-XPS-8300:~/WorkSpace/MeanAndExpection/SimuExpect$ g++ -o prg MeanAndExpectation.cpp myrn.cpp
okamoto@okamoto-XPS-8300:~/WorkSpace/MeanAndExpection/SimuExpect$ ./prg
i = 10
Fi = 1 2 4 1 1 1
Pi = 0.10000 0.20000 0.40000 0.10000 0.10000 0.10000
mean = 3.20000
i = 20
Fi = 1 5 7 3 2 2
Pi = 0.05000 0.25000 0.35000 0.15000 0.10000 0.10000
mean = 3.30000
i = 100
Fi = 17 17 22 14 13 17
Pi = 0.17000 0.17000 0.22000 0.14000 0.13000 0.17000
mean = 3.40000
i = 1000
Fi = 171 165 195 159 147 163
Pi = 0.17100 0.16500 0.19500 0.15900 0.14700 0.16300
mean = 3.43500
i = 10000
Fi = 1666 1609 1763 1582 1649 1731
Pi = 0.16660 0.16090 0.17630 0.15820 0.16490 0.17310
mean = 3.51320
i = 1000000
Fi = 166408 166013 167012 167030 166812 166725
Pi = 0.16641 0.16601 0.16701 0.16703 0.16681 0.16673
mean = 3.50200
okamoto@okamoto-XPS-8300:~/WorkSpace/MeanAndExpection/SimuExpect$
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムで用いられている乱数生成のアルゴリズムはWhichmann/Hill generatorと呼ばれているものである（岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房、第8.2節「乱数」参照）。この乱数の周期は約
であるが、この周期分の乱数を呼び出す時間を計るプログラムをリスト５に示す。プログラムのソースコードファイルは圧縮ファイルTimeWHRN.zipとしてまとめて用意している。ファイル名TimeWHRN.zipをクリックしてダウンロードすることができる。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト５ TimeWHRn.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <iostream>
#include <iomanip>
#include <ctime>
#include <cmath>
#include "myrn.h" // 乱数生成用クラス型 WHrn のヘッダーファイル
using namespace std;
using namespace myrn;
int main()
{
cout << "The program started." << endl;
WHrn rn;
time_t t0;
time(&t0);
for (double i = 0; i < 1.0e9; i += 1.0){
rn.uni();
}
time_t t1;
time(&t1);
cout << t1 - t0 << "sec elapsed." << endl;
cout << "It will take about " << fixed << setprecision(2)
<< ((t1 - t0) * pow(10.0, 12) / (1.0e9 * 60.0 * 60.0))
<< " hours to call uni() 10**12 times." << endl;
return 0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
リスト５のプログラムの実行例を図２に示す。
図２
この実行例で使用されたPCの場合、１８時間内のシミュレーションであれば乱数の周期の問題はないことがわかる。より長い周期の乱数については別のウェブサイトで説明した。