情報基礎A 「Cプログラミング」（ステップ７・乱数と中心極限定理）

このページでは、確率や統計の分野で極めて重要な中心極限定理について考える。（このページはまだ書きかけ）

1. 擬似乱数

コンピュータを使って、サイコロ投げなどのランダムな事象を模倣する際に、乱数を自動発生させる機能があると便利だ。 さもないと、「人間がサイコロを振って、出た目を、パソコンに入力する」といった作業を強いられることになる。 そのため、標準的なCのライブラリ関数には、rand( )という関数が用意されている。

使い方は、まず宣言部に

#include <stdlib.h>

と書いてから、引数無しで rand() と呼ぶと、0からRAND_MAXまでの間の整数が 等しい確率で、ランダムに返される。つまり、呼び出す度に、値が変わる。 ここで、RAND_MAXという記号には、通常、int型整数の最大値(2147483647)が設定されている。 （つまり、プログラムの中では 2147483647 と書いても、RAND_MAXと書いても同じこと）。

その他の関数

コンピュータが生成する乱数は、ある種の計算の結果出力された擬似的な乱数であって、本当の意味の乱数ではない。 そして、上で紹介した関数rand()は、乱数としての性質があまり優れていないと言われている。 標準的なUnix環境では、異なる乱数生成アルゴリズムに基づく random() 関数や lrand48() 関数が 組み込まれているので、より本格的なシミュレーションをする際には、チェックしておくと良い。 さらに、（標準では組み込まれていないが）優れた擬似乱数生成アルゴリズムがネット上で公開されている（例えば Mersenne Twister, ISAAC 等)。

擬似乱数を使って、「サイコロ」を10万回振り、出た目（0〜5）の数と割合を出力するプログラム例を以下に示す：

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

main()
{
  int dice[6] = {0,0,0,0,0,0} ;
  int i,cnt ;

  for (cnt=0; cnt<100000; cnt=cnt+1) {
    i = rand()%6 ;
    dice[i] = dice[i] + 1 ;
  }

  for (i=0; i<6; i=i+1) {
    printf("%d  %d  %f\n",i,dice[i],(float)dice[i]/cnt) ;
  }
}
実数値の乱数

物理シミュレーション等では、実数の乱数が必要になる場合が多い。上記のrand()関数の値をその最大値で割って

(float)は、その右の値を float と見なして計算せよ、という指示である。

x = (float)rand()/RAND_MAX ;

のように記述すれば、0から1の間に一様に分布する実数の乱数（一様乱数）を得ることができる。

icon-pc　練習：モンテカルロ法

乱数を用いた算法は、一般にモンテカルロ法と呼ばれている。

以下の考え方に基づいて、シミュレーションによって円周率の値を推定するプログラムを作成してみなさい。

1: カウンター c をリセットする
2: 4までをn回反復する
3: 　　0から1の間で一様に分布する二つの乱数 x, y を生成する
4: 　　もし x*x + y*y が1より小さければ（(x,y)が半径1の円の内部ならば）カウンター c をひとつ増やす
5: nが十分大きければ、この時点で、c/nは半径1の円の面積の1/4 (すなわちπ/4)になっているはずである。
icon-hint　ヒント

サンプル数を$n$とすると、一般に、統計的な揺らぎは$1/\sqrt{n}$程度で減少するので、この種の計算では、誤差を1/10にするには、100倍の反復が必要となる。


2. 乱数の和

$X$を、サイコロの目など、ランダムに生じる数を表すとしよう。そして、同様の、しかし互いに独立な$n$回の試行の 1回目を$X_1$, 2回目を$X_n$, ... とし、それらの和を回数で割った量 $$ \bar{X}_n = \frac{X_1 + X_2 + \cdots + X_n}{n} $$ を考えると、大数の法則によって、$n \to \infty$で$\bar{X}_n$は$X$の平均値$\mu$に収束する。 1から6の目が同じ確率で出る（いかさまではない）サイコロの場合は、各目の出る確率は1/6なので、 $$ \bar{X}_n \to 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + \cdots + 6 \times \frac{1}{6} = 3.5 $$ となるわけだ。

では、次に、$n$で割らないで、単純に加えた数 $X_1 + X_2 + \cdots + X_n$ は$n$を大きくすると、どのように振る舞うであろうか。 サイコロの例を念頭に考えると、$n$を増やすとともに、その和もまた増大することは明らかであるが、それに加えて、総和した値の変動幅も大きくなるのではないだろうか？

そこで、実際に、コンピュータ上でサイコロを振って、この様子を確かめてみた。 下のグラフは、『サイコロを$n$回振って、その目の合計を調べる』作業を10万回繰り返し、目の合計値がどのように分布するかをプロットしたものである。 ここで、$n$は、3, 5, 10の場合について調べた。

c-7-array-statistics-clt-dice

すると、よく見慣れたベル型の分布が現れるのが分かる。 この様子を数学的に表現したのが中心極限定理である。 すなわち、 $X_i$の平均を$\mu$、分散を$\sigma^2$とすると、それらの$n$項分の和の分布は、平均$n \mu$, 分散$n \sigma^2$の正規分布に漸近する、 すなわち、$n \to \infty$ で $$ X_1 + X_2 + \cdots + X_n \sim N(n \mu, n \sigma^2) $$ となる。 上のグラフからも分かるように、元々の$X_i$の分布は（平均と分散が有限ならば）どんな形でもよく、しかも、$n$が5程度でも、非常に正規分布に近い分布形が得られる点に注目したい。 このことはまた、少数であっても、何らかのランダムな要因が加算的に作用すると、その結果として、正規分布が現れる可能性を示唆しており、正規分布が様々な現象に顔を出すひとつの理由と考えられる。

上記の性質を使うと、一様乱数$X_i$を$n$個加えた量を使って、 $$ Y = \frac{X_1 + X_2 + \cdots + X_n - n \mu}{\sqrt{n}\sigma} $$ を計算すれば($n$が10程度でも）正規分布$N(0,1)$に従う乱数を生成できることになる（ 計算速度を求めなければ、この方法は十分に実用的である。 ただし、一様乱数から正規分布する乱数を生成するアルゴリズムとしては Box-Mullerの方法が一般には良く使われている。）

icon-pc　練習：正規乱数の生成

上記のアルゴリズムに従い、1〜6の目を持つサイコロの値を加えることで、正規分布$N(0,1)$に従う乱数を発生させ、得られた頻度分布と正規分布 $$ f(y) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{y^2}{2}\right) $$ とを比較するプログラムを作成し、$n$の値を変えて、分布形の変化（両者の差異）を確認せよ。

icon-hint　ヒント

サイコロの場合、平均$\mu=3.5$, 分散 $\sigma^2 = \frac{(1-3.5)^2 + (2-3.5)^2 + \cdots + (6-3.5)^2}{6}=2.9166666\cdots$である。

$X_1 + X_2 + \cdots + X_n$の値が$x$であるような頻度が$p(x)$であったとしよう（ここで、$\sum_x p(x) = 1$とする）。 すると、変数変換 $y = (x - n \mu)/(\sqrt{n}\sigma)$ に対して、分布関数は $$ p(y)dy = p(y(x)) \left| \frac{dx}{dy} \right| dy $$ のように変換されるので、$f(y)$と比較すべきは $p(y(x)) \left| \frac{dx}{dy} \right|$ であることに注意 （ここで、$\frac{dx}{dy} = \sqrt{n} \sigma$）。