 |
|
２． 乱数はシュミレーションの基本－－ 乱数の発生方法－－
|
何の規則性もなしに、物事が起こることをコンピュータで実験するには、でたらめな数の列－－乱数－－が必要です。この章では、シミュレーションに必ず必要となる、乱数の発生方法と、その性質を理解します。
乱数とは
乱数とは、なんらかの方法で、どの数も等しい確率で、でたらめになるように選んで、数字の列を作ったものです。例えば、さいころをふって出た目を記録すれば、それは１から６までの乱数になります。あるいは、１円のアルミの硬貨を
１００枚用意しておき、その１枚、１枚に、００、０１、０２、・・・、９９という番号を書いておき、それを十分かきまぜて、目をつぶって取り出せば、００から９９までの乱数となります。
コンピュータによる乱数作り
コンピュータで乱数を発生させるには、いくつかの方法がありますが、例えばある数字（これを“たね”と言います）をなんども掛けあわせると、０から９ま
での数字がほぼ同数ずつ現われ、規則性も少ないので、乱数発生して使われます。このような方法で作られる乱数は、真の意味の乱数ではなく、疑似乱数と言い
ます。
コンピュータの疑似乱数を乱数として使用する場合は、次のような条件を満足しなければなりません。
① 周期が長いこと
疑似乱数では、まったく同じ数列が繰り返し出現する場合があります。その周期が短いと乱数としての意味がなくなります。
② 統計的検定に耐えること
たとえ、周期が長くとも、
０、 ０、 ０、 １、 １、 １、 ２、 ２、 ２、 ２
といった、ある種の法則に従って出現する数列は、乱数とは言えません。あとで述べる統計的検定という方法で、試験をします。
③ 再現性があること
シミュレーションに用いる場合、同一の乱数を使用して、異なったモデルでシミュレーションし、比較検討する場合があります。またプログラムのデバッグ時にも何度も同じ系列の乱数を必要とする場合があります。すなわち乱数に、再現性が必要されます。
④ 乱数発生のスピードが早いこと
シミュレーションでは、乱数を何度も繰り返し用いて実験を行い、多数回の実験結果から普遍性のある共通因子を求めます。このため乱数発生に要する計算スピードが早くなければなりません。
次に述べるように，Ｃ言語には，rand()関数という乱数発生機構があり，これは，上記①～④の条件がよく満足されます。
rand()関数
rand()関数は，Ｃのマニュアルによれば，標準ライブラリ関数のユーティリティ（stdlib.h）の中にあり，
int rand(void);
と定義されており，この関数が呼び出されるたびに，正の整数型の乱数が発生されます。整数型の範囲は，コンピュータあるいはコンパイラによって異なり，通常１６ビットのパソコンでは最大32,767（すなわち215 - 1），３２ビットのワークステーションでは最大2,147,483,647（すなわち232 - 1）となります。 rand()関数から出力される乱数は，０からこの整数の最大値の範囲です。また，システムによっては，stdlib.hの中で，最大値が
#define RANDMAX 32767
として，定義されており，RANDMAXで最大値を知ることができます。（Ｔｕｒｂｏ Ｃでは誤っています。）
本書においては，このrand()関数をすべてのプログラムで用いますので，詳しく調べてみることにしましょう。
ためしに簡単なプログラム s0220.cをコンパイルして実行してみます。結果は，図２．２．１のとおりです。２０個の０から整数型の最大値の範囲の乱数が出力されました。
もう一度，このプログラムを実行してみても，出力は同じです。Ｃのrand()関数で出力される乱数は，再現性ががあります。
| s0220.c rand()を使ってみる |
/* s0220.c rand()関数 */
#include <stdio.h> /* printf() */
#include <stdlib.h> /* rand() */
#define NUMBER 20
int main(void)
{
int j,r;
for (j = 1; j <= NUMBER; j++){
r = rand();
printf("%10d\n",r);
}
return 0;
}
|
| 図２．２．１ s0220.cの実行結果 |
 |
０以上１未満の乱数発生－－RAND()マクロ－－
シミュレーションで用いる乱数は，機種によらず，０以上で１未満であると便利です。このため，sO221.c のように，rand()関数の値を（最大値＋１）で割って使うこととします。最大値は，１６ビット系のパソコンの場合32767，３２ビット系のワークステーションの場合2147483647です。s0221.cではこのことを，#defineで１６ビット系の場合
#define RAND() ((double)rand()/(32767.0+1.0))
３２ビット系の場合
#define RAND() ((double)rand()/(2147483647.0+1.0)
のように，マクロとして定義しています。（s0221.cでは，これらのうち使わない方を削除するか，コメントとして無効にしてからコンパイルして下さい。）このようにしておけば，シミュレーションのプログラムにおいて，プリプロッセッサが，すべてのRAND()という箇所で，この定義に従って置き換えてくれます。本書においては，大文字のRAND()が０以上１未満の乱数発生機構ということになります。なお，マクロ定義ではなく，０以上１未満の乱数発生を，関数として定義する方法もあります。
| s0221.c ０以上１未満の乱数 |
/* s0221.c ０以上１未満の乱数 */
#include <stdio.h> /* printf() */
#include <stdlib.h> /* rand() */
/* １６ビット */
#define RAND() ((double)rand()/(32767+1.0))
/* ３２ビット */
/* #define RAND() ((double)rand()/(2147483647.0+1.0)) */
#define NUMBER 20
int main(void)
{
int j;
double random;
for (j = 1; j <= NUMBER; j++){
random = RAND();
printf(" %1.10f\n", random);
}
return 0;
}
|
| 図２．２．２ s0221.cの実行結果 |
|
srand()関数
次にrand()関数の“たね ”（seed）を設定するsrand()関数を使ってみます。srand()関数は，
void srand(unsigned seed);
と定義され，これにより，rand()関数をイニシャライズし，乱数列を変えることができます。プログラムs0222.cの行16のように，srand(2)を加えてみます。たしかに，先ほどとは異なった乱数列が出力されます。
| s0222.c rand()を使ってみる |
/* s0222.c ０以上１未満の乱数（たねを設定） */
#include <stdio.h> /* printf() */
#include <stdlib.h> /* rand() */
/* １６ビット */
#define RAND() ((double)rand()/(32767.0+1.0))
/* ３２ビット */
/* #define RAND() ((double)rand()/(2147483647.0+1.0)) */
#define NUMBER 20
int main(void)
{
int j;
double random;
srand(2); /*たねを設定*/
for (j = 1; j <= NUMBER; j ++){
random = RAND();
printf(" %1.10f\n", random);
}
return 0;
}
|
| 図２．２．３ s0222.cの実行結果 |
|
実行するごとに異なった系列の乱数を発生させる
rand()関数は，たねが同じであるかぎり，同一系列の乱数を発生します。これは，シミュレーションプログラムのデバッグ時には，実行するごとに同一の結果が得られ，大変都合がよろしい。しかし，場合によっては，実行するごとに別系列の乱数を発生させたい場合があります。これには，次のように
srand()関数のたねをtime()関数で実行するごとに変えることにより実現できます。
srand((unsigned)time(NULL));
一様に分布しているか
Ｃの乱数が，本当に一様分布しているかどうかを調べます。プログラムs0230.cのように，RAND()マクロで乱数を次々発生させ，0.1ごとのキザミ幅，0～0.1, 0.1～0.2, 0.2～0.3, ･･･にそれぞれ何個入るかをカウントしてみます｡s0230.cの行24で乱数を発生させ，それぞれのキザミに対応し，i= 0, 1, 2, 3, ･･･,9 となるように１０倍し，整数部分を取り出しています。行25では，それぞれのキザミに，いくつ落ちるかをf[i]としてカウントします。
理論分布との比較
ところで，理論分布f0は10000 / 10 = 1000となるはずですが，プログラムs0230.cの結果が，果たして，一様乱数と言えるのか心配になります。一般に統計学では母集団をサンプリングして得られた分布と，既知の分布との間に差があるかどうかを判定することを，“適合度の検定”（以下では単に検定）と言います。例えば，全国民の年令別の構成がわかっている時，東京都民の年令構成が全国と差があるかどうかを，東京都民の数百人の年令構成を調査し，そのサンプリング・データから判断する場合などに，この検定が用いられます。
さきほどの，一様乱数発生の例では，サンプリング・データはf(i)，既知の分布はf0 = 1000に相当します。
検定
検定は，次のように行ないます。
 |
・・・・・ |
式２．３．１ |
というものを考えてみます。これを通常 χ2 （カイ自乗と読む）と言います。この量は、サンプリング・データと、理論度数とのくい違いの程度を表わしています。もしサンプリング・データと理論度数のくい違いがまったくなければ、χ2 = 0 となり、くい違いが大きいほど χ2 は大きくなります。
同じ母集団からのサンプリングのやり方は、何通りもあるので、 χ2 は”χ2分布”という分布に従うことが知られています。χ2分布は自由度と呼ばれるパラメータを持っていて、この例では９（キザミの数－１）です。χ2分布は図２．３．１のような形をしており、χ2分布表を表２．３．１に示します。この表は、 χα2 と χα2 の右側の面積αの関係を表わしており、サンプリング・データから得られた χ2 が χα2 より小さければ、”危険率αで理論度数とサンプリング・データが等しい”と言えることになります。
| 図２．３．１ χ2分布 |
 |
頻度検定のプログラム
プログラムs0230.cの 行32は xs として χ2 を計算しています。このプログラムの出力は図のように，χ2 = 9.17となります。一方、表２．３．１から危険率５％の χα2 は 16.919 ですから、xs はこの値より小さく、したがって５％の危険率で一様分布であると言ってよいことがわかります。すなわち、RAND()の頻度テストは合格したことになります。
| s0230.c 頻度検定プログラム |
/* s0230.c 頻度検定 (C) H.Ishikawa 1994 2018*/
#include <stdio.h> /* printf() */
#include <stdlib.h> /* rand() */
/* １６ビット */
#define RAND() ((double)rand()/(32767.0+1.0))
/* ３２ビット */
/* #define RAND() ((double)rand()/(2147483647.0+1.0))*/
#define NUMBER 10000 /* 発生させる乱数の数 */
int main(void)
{
int j; /* forのカウンタ */
int i; /* 0.1ごとのキザミの番号 */
int f[10] = {0}; /* キザミiに落ちる乱数の数 */
double f0; /* 理論度数 */
double xs = 0.0; /* カイ自乗の値 */
srand(15);
f0 = NUMBER/10.0;
for (j = 1; j <= NUMBER; j ++){
i = (int)(RAND() * 10);
f[i] = f[i] + 1;
}
printf(" i f[i]\n");
/*頻度の検定*/
for (i = 0; i <= 9; i ++){
printf("%10d %10d\n", i, f[i]);
xs = xs + (f[i] - f0) * (f[i] - f0) / f0;
}
printf("xs = %f\n", xs);
return 0;
}
|
| 図２．３．２ s0230.cの実行結果 |
|
ポーカ検定とは
頻度が一様だからと言って、ただちに一様乱数であるとは言えません。例えば、
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, ････
は頻度検定には合格しますが、繰り返しが規則的で乱数とは認められないでしょう。ポーカ検定は、無規則性の検定の手法で、整数の乱数の列を連続して５つず
つまとめ、５桁の数字とみてグループを作り、トランプのポーカと同じように手を作って調べます。０から９の種類のカードを５枚くばったと考えると、次の”
場合”があります。
０ abcde: 5の数字が全部異なっている。
１ aabcd: 同じ数字２つからなるペアが一組ある。
２ aabbc: ペアが二組ある。
３ aaabc: 同じ数字が３つある。
４ aaabb: 同じ数字が３つと、ペアが１つある。
５ aaaab: 同じ数字が４つある。
６ aaaaa: 全部の数字が同じである。
これらの起こる確率は理論的に求められ、その値を表２．４．１に示します。
ポーカ検定のプログラム
この理論値と，実際にコンピュータで発生させた疑似乱数による数字のならびと比較してみましょう。プログラム s0240.cでは，行33から行35において，RAND()マクロにより５個の０から９までの整数のならびを作ります。それを n[0], n[1], n[2], n[3], n[4] とします。次に行37から行43において，図２．４．１のとおりの比較を行い，もし同じ数字であれば，countを１増加させていきます。全部同じ数字であればcount = 10となり，全部異なればcount = 0となります。上に述べた”場合”０～６に対して，countは表２．４．１の値をとりますから，行44から行50で，それを”場合”を表わす i に変換します。そして行51により，頻度f[i]をカウントし，行57において，χ2を
 |
・・・・・ |
式２．４．１ |
によりもとめ、この χ2 を自由度６ のχ2 分布表を用いて検定します。
このプログラムの実行結果を図２．４．２に示します。χ2の値は 2.602となり，ＣのRAND()関数は，ポーカテストにおいても合格します。
| 図２．４．１ ポーカーテストの方法 |
 |
| s0240.c ポーカ検定 |
/* s0240.c ポーカ検定 (C) H.Ishikawa 1994 2018*/
#include <stdio.h> /* printf() */
#include <stdlib.h> /* rand() */
/* １６ビット */
#define RAND() ((double)rand()/(32767.0+1.0))
/* ３２ビット */
/* #define RAND() ((double)rand()/(2147483647.0+1.0)) */
#define NUMBER 10000 /* 発生させる乱数の数 */
int main(void)
{
int i; /* 表に示すi */
int count; /* 表に示すcount */
int j, j1, j2, k; /* forのカウンタ */
int f[7] = {0}; /* iごとの度数 */
int f0[7]; /* 理論度数 */
int n[5]; /* ５個の数字の並び */
double xs = 0.0; /* カイ自乗値 */
double diff; /* 理論値との差 */
srand(2);
f0[0] = NUMBER * 0.3024;
f0[1] = NUMBER * 0.5040;
f0[2] = NUMBER * 0.1080;
f0[3] = NUMBER * 0.0720;
f0[4] = NUMBER * 0.0090;
f0[5] = NUMBER * 0.0045;
f0[6] = NUMBER * 0.0001;
for (k = 1; k <= NUMBER; k ++){
for (j = 0; j < 5; j ++){
n[j] = (int)(RAND() * 10);
}
count = 0;
for (j1 = 0; j1 < 4; j1 ++){
for (j2 = j1 + 1; j2 < 5; j2 ++){
if (n[j1] == n[j2]) {
count = count + 1;
}
}
}
if (count == 10){
i = 6;
} else if (count == 6){
i = 5;
} else{
i = count;
}
f[i] = f[i] + 1;
}
printf(" i f[i] f0[i]\n");
for (i = 0; i <= 6; i++){
printf("%10d %10d %10d\n", i, f[i], f0[i]);
diff = f[i] - f0[i];
xs = xs + diff * diff / f0[i];
}
printf(" xs = %f\n", xs);
return 0;
}
|
| 図２．４．２ s0240.cの実行結果 |
|
１． プログラム s0230.cで，たねの値を変えて実行し， χ2の値がどう変わるか，調べなさい。
（解答省略）
２． １から６までの整数が等確率で出現するさいころのシミュレーションを行ないなさい。さいころをふるたびに，でたらめな数がでるようにしなさい。
（解答）
| a0220.c ２章２．の解答 |
/* a0220.c ２章２．の解答
* サイコロ
* (C) H.Ishikawa 1994 2018
*/
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
/* １６ビット */
#define RAND() ((double)rand()/(32767.0+1.0))
/* ３２ビット */
/* #define RAND() ((double)rand()/(2147483647.0+1.0)) */
#define END 20
int d;
int t;
int main(void)
{
srand((unsigned)time(NULL)); /* 乱数系列を変える */
for (t = 1; t <= END; t ++) {
d = (int)(RAND() * 6.0) + 1;
printf("%d\n",d);
}
returen 0;
}
|
３． １桁の区間０から９までの乱数があったとき，同じ数が現われるまでの間隔を，ギャップといいます。その長さがrである確率は，
で表わされます（r = 0は，続いて同じ数字が現われた場合）。このことを用いて，Ｃの乱数のギャップ検定を行ないなさい。
（解答）
| a0230.c ２章３．の解答 |
/* a0230.c ２章３．の解答
* ギャップ検定
* (C) H.Ishikawa 1994 2018
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* １６ビット */
#define RAND() ((double)rand()/(32767.0+1.0))
/* ３２ビット */
/* #define RAND() ((double)rand()/(2147483647.0+1.0)) */
#define NUMBER 10000
long r = 0; /* ギャップ */
long n; /* 乱数発生のfor */
long m; /* 同じ数が出た回数 */
long f[100] = {0}; /* ギャップの頻度 */
double f0; /* 理論値 */
double xs = 0.0; /* カイ自乗値 */
int main(void)
{
for (n = 1; n <= NUMBER; n ++) {
/* 0の出る間隔を調べる */
if ((long)(RAND() * 10.0) == 0) {
f[r] = f[r] + 1;
m = m +1;
r = 0;
} else {
r = r + 1;
}
}
for ( r = 0; r <= 19; r ++) {
f0 = pow((1.0 - 0.1),(double)r) * 0.1 * (NUMBER / 10);
printf("%6ld %6ld %6.4f\n", r, f[r], f0);
xs = xs + ((double)f[r] - f0) * ((double)f[r] - f0) / f0;
}
printf("xs = %f\n", xs);
return 0;
}
|
|
 |
