本文へスキップ
モンテカルロシミュレーション
2. 乱数はシュミレーションの基本−− 乱数の発生方法−−

 何の規則性もなしに、物事が起こることをコンピュータで実験するには、でたらめな数の列−−乱数−−が必要です。この章では、シミュレーションに必ず必要となる、乱数の発生方法と、その性質を理解します。

乱数とは
 乱数とは、なんらかの方法で、どの数も等しい確率で、でたらめになるように選んで、数字の列を作ったものです。例えば、さいころをふって出た目を記録すれば、それは1から6までの乱数になります。あるいは、1円のアルミの硬貨を 100枚用意しておき、その1枚、1枚に、00、01、02、・・・、99という番号を書いておき、それを十分かきまぜて、目をつぶって取り出せば、00から99までの乱数となります。


コンピュータによる乱数作り
 コンピュータで乱数を発生させるには、いくつかの方法がありますが、例えばある数字(これを“たね”と言います)をなんども掛けあわせると、0から9ま での数字がほぼ同数ずつ現われ、規則性も少ないので、乱数発生して使われます。このような方法で作られる乱数は、真の意味の乱数ではなく、疑似乱数と言い ます。
 コンピュータの疑似乱数を乱数として使用する場合は、次のような条件を満足しなければなりません。

@ 周期が長いこと
疑似乱数では、まったく同じ数列が繰り返し出現する場合があります。その周期が短いと乱数としての意味がなくなります。

A 統計的検定に耐えること
たとえ、周期が長くとも、
        0、 0、 0、 1、 1、 1、 2、 2、 2、 2
といった、ある種の法則に従って出現する数列は、乱数とは言えません。あとで述べる統計的検定という方法で、試験をします。

B 再現性があること
 シミュレーションに用いる場合、同一の乱数を使用して、異なったモデルでシミュレーションし、比較検討する場合があります。またプログラムのデバッグ時にも何度も同じ系列の乱数を必要とする場合があります。すなわち乱数に、再現性が必要されます。

C 乱数発生のスピードが早いこと
 シミュレーションでは、乱数を何度も繰り返し用いて実験を行い、多数回の実験結果から普遍性のある共通因子を求めます。このため乱数発生に要する計算スピードが早くなければなりません。

 次に述べるように,C言語には,rand()関数という乱数発生機構があり,これは,上記@〜Cの条件がよく満足されます。

rand()関数
 rand()関数は,Cのマニュアルによれば,標準ライブラリ関数のユーティリティ(stdlib.h)の中にあり,

int rand(void);


と定義されており,この関数が呼び出されるたびに,正の整数型の乱数が発生されます。整数型の範囲は,コンピュータあるいはコンパイラによって異なり,通常16ビットのパソコンでは最大32,767(すなわち215 - 1),32ビットのワークステーションでは最大2,147,483,647(すなわち232 - 1)となります。 rand()関数から出力される乱数は,0からこの整数の最大値の範囲です。また,システムによっては,stdlib.hの中で,最大値が

#define RANDMAX 32767


として,定義されており,RANDMAXで最大値を知ることができます。(Turbo Cでは誤っています。)
 本書においては,このrand()関数をすべてのプログラムで用いますので,詳しく調べてみることにしましょう。
 ためしに簡単なプログラム s0220.cをコンパイルして実行してみます。結果は,図2.2.1のとおりです。20個の0から整数型の最大値の範囲の乱数が出力されました。

 もう一度,このプログラムを実行してみても,出力は同じです。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;
}

図2.2.1 s0220.cの実行結果



0以上1未満の乱数発生−−RAND()マクロ−−
 シミュレーションで用いる乱数は,機種によらず,0以上で1未満であると便利です。このため,sO221.c のように,rand()関数の値を(最大値+1)で割って使うこととします。最大値は,16ビット系のパソコンの場合32767,32ビット系のワークステーションの場合2147483647です。s0221.cではこのことを,#defineで16ビット系の場合

#define RAND() ((double)rand()/(32767.0+1.0))

32ビット系の場合

#define RAND() ((double)rand()/(2147483647.0+1.0)

のように,マクロとして定義しています。(s0221.cでは,これらのうち使わない方を削除するか,コメントとして無効にしてからコンパイルして下さい。)このようにしておけば,シミュレーションのプログラムにおいて,プリプロッセッサが,すべてのRAND()という箇所で,この定義に従って置き換えてくれます。本書においては,大文字のRAND()が0以上1未満の乱数発生機構ということになります。なお,マクロ定義ではなく,0以上1未満の乱数発生を,関数として定義する方法もあります。

  s0221.c    0以上1未満の乱数
/*  s0221.c  0以上1未満の乱数  */
#include <stdio.h>                /* printf() */
#include <stdlib.h>               /* rand() */

/* 16ビット */
#define RAND()     ((double)rand()/(32767+1.0))
/* 32ビット */
/* #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;
}     

図2.2.2 s0221.cの実行結果


srand()関数
 次にrand()関数の“たね ”(seed)を設定するsrand()関数を使ってみます。srand()関数は,

void srand(unsigned seed);

と定義され,これにより,rand()関数をイニシャライズし,乱数列を変えることができます。プログラムs0222.cの行16のように,srand(2)を加えてみます。たしかに,先ほどとは異なった乱数列が出力されます。


  s0222.c    rand()を使ってみる
/*  s0222.c  0以上1未満の乱数(たねを設定) */
#include <stdio.h>                /* printf() */
#include <stdlib.h>               /* rand() */

/* 16ビット */
#define RAND()     ((double)rand()/(32767.0+1.0))
/* 32ビット */
/* #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;
}
     

図2.2.3  s0222.cの実行結果



実行するごとに異なった系列の乱数を発生させる
 rand()関数は,たねが同じであるかぎり,同一系列の乱数を発生します。これは,シミュレーションプログラムのデバッグ時には,実行するごとに同一の結果が得られ,大変都合がよろしい。しかし,場合によっては,実行するごとに別系列の乱数を発生させたい場合があります。これには,次のように srand()関数のたねをtime()関数で実行するごとに変えることにより実現できます。

srand((unsigned)time(NULL));

一様に分布しているか
 Cの乱数が,本当に一様分布しているかどうかを調べます。プログラム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 となるように10倍し,整数部分を取り出しています。行25では,それぞれのキザミに,いくつ落ちるかをf[i]としてカウントします。


理論分布との比較
 ところで,理論分布f0は10000 / 10 = 1000となるはずですが,プログラムs0230.cの結果が,果たして,一様乱数と言えるのか心配になります。一般に統計学では母集団をサンプリングして得られた分布と,既知の分布との間に差があるかどうかを判定することを,“適合度の検定”(以下では単に検定)と言います。例えば,全国民の年令別の構成がわかっている時,東京都民の年令構成が全国と差があるかどうかを,東京都民の数百人の年令構成を調査し,そのサンプリング・データから判断する場合などに,この検定が用いられます。
 さきほどの,一様乱数発生の例では,サンプリング・データはf(i),既知の分布はf0 = 1000に相当します。


検定
 検定は,次のように行ないます。
・・・・・ 式2.3.1

というものを考えてみます。これを通常 χ2 (カイ自乗と読む)と言います。この量は、サンプリング・データと、理論度数とのくい違いの程度を表わしています。もしサンプリング・データと理論度数のくい違いがまったくなければ、χ2 = 0 となり、くい違いが大きいほど χ2  は大きくなります。
 同じ母集団からのサンプリングのやり方は、何通りもあるので、 χ2 は”χ2分布”という分布に従うことが知られています。χ2分布は自由度と呼ばれるパラメータを持っていて、この例では9(キザミの数−1)です。χ2分布は図2.3.1のような形をしており、χ2分布表を表2.3.1に示します。この表は、  χα2  と  χα2  の右側の面積αの関係を表わしており、サンプリング・データから得られた  χ2  が χα2 より小さければ、”危険率αで理論度数とサンプリング・データが等しい”と言えることになります。


図2.3.1 χ2分布


頻度検定のプログラム
 プログラムs0230.cの 行32は xs として χ2 を計算しています。このプログラムの出力は図のように,χ2 = 9.17となります。一方、表2.3.1から危険率5%の χα2 は 16.919 ですから、xs はこの値より小さく、したがって5%の危険率で一様分布であると言ってよいことがわかります。すなわち、RAND()の頻度テストは合格したことになります。

  s0230.c   頻度検定プログラム
/*  s0230.c    頻度検定   (C) H.Ishikawa  1994 2018*/
#include <stdio.h>           /* printf() */
#include <stdlib.h>          /* rand() */

/* 16ビット */
#define RAND()     ((double)rand()/(32767.0+1.0))
/* 32ビット */ 
/* #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;
}

図2.3.2 s0230.cの実行結果


ポーカ検定とは
 頻度が一様だからと言って、ただちに一様乱数であるとは言えません。例えば、

    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, ・・・・

は頻度検定には合格しますが、繰り返しが規則的で乱数とは認められないでしょう。ポーカ検定は、無規則性の検定の手法で、整数の乱数の列を連続して5つず つまとめ、5桁の数字とみてグループを作り、トランプのポーカと同じように手を作って調べます。0から9の種類のカードを5枚くばったと考えると、次の” 場合”があります。

0  abcde:  5の数字が全部異なっている。
1  aabcd:  同じ数字2つからなるペアが一組ある。
2  aabbc:  ペアが二組ある。
3  aaabc:  同じ数字が3つある。
4  aaabb:  同じ数字が3つと、ペアが1つある。
5  aaaab:  同じ数字が4つある。
6  aaaaa:  全部の数字が同じである。 

 これらの起こる確率は理論的に求められ、その値を表2.4.1に示します。


ポーカ検定のプログラム
 この理論値と,実際にコンピュータで発生させた疑似乱数による数字のならびと比較してみましょう。プログラム s0240.cでは,行33から行35において,RAND()マクロにより5個の0から9までの整数のならびを作ります。それを n[0], n[1], n[2], n[3], n[4] とします。次に行37から行43において,図2.4.1のとおりの比較を行い,もし同じ数字であれば,countを1増加させていきます。全部同じ数字であればcount = 10となり,全部異なればcount = 0となります。上に述べた”場合”0〜6に対して,count表2.4.1の値をとりますから,行44から行50で,それを”場合”を表わす i に変換します。そして行51により,頻度f[i]をカウントし,行57において,χ2
・・・・・ 式2.4.1
によりもとめ、この χ2 を自由度6 のχ2 分布表を用いて検定します。


 このプログラムの実行結果を図2.4.2に示します。χ2の値は 2.602となり,CのRAND()関数は,ポーカテストにおいても合格します。


図2.4.1 ポーカーテストの方法


 s0240.c ポーカ検定
/* s0240.c   ポーカ検定    (C) H.Ishikawa  1994 2018*/
#include <stdio.h>              /* printf() */
#include <stdlib.h>             /* rand() */

/* 16ビット */
#define RAND()     ((double)rand()/(32767.0+1.0))
/* 32ビット */
/* #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];            /* 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;
}
  

図2.4.2 s0240.cの実行結果




 1. プログラム s0230.cで,たねの値を変えて実行し, χ2の値がどう変わるか,調べなさい。
  (解答省略)

 2. 1から6までの整数が等確率で出現するさいころのシミュレーションを行ないなさい。さいころをふるたびに,でたらめな数がでるようにしなさい。
  (解答) 
 a0220.c   2章2.の解答
/* a0220.c  2章2.の解答   
 *            サイコロ
 *            (C) H.Ishikawa  1994 2018
 */

#include <time.h>
#include <stdio.h>
#include <stdlib.h>

/* 16ビット */
#define RAND()     ((double)rand()/(32767.0+1.0))
/* 32ビット */
/* #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;
}


 3. 1桁の区間0から9までの乱数があったとき,同じ数が現われるまでの間隔を,ギャップといいます。その長さがrである確率は,

で表わされます(r = 0は,続いて同じ数字が現われた場合)。このことを用いて,Cの乱数のギャップ検定を行ないなさい。
  (解答)
 a0230.c   2章3.の解答
/* a0230.c  2章3.の解答 
 *           ギャップ検定  
 *        (C) H.Ishikawa  1994 2018
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

/* 16ビット */
#define RAND()     ((double)rand()/(32767.0+1.0))
/* 32ビット */
/* #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;
}
    


top