疑似乱数

講義

はじめに

この講義のメインテーマは「確率微分方程式の数値計算」である． 確率微分方程式とは不確定な要素（「ノイズ」という）を含む 微分方程式である．この不確定な要素を導入するために， 不確定要素を含まない微分方程式の数値計算方法と異なる操作を しなければならない．この講義では，不確定要素を含むとなぜ 異なる数値計算法を用いなければならないのかを学習してゆく． 第一回は「擬似乱数」について考えてゆく．乱数とはランダムな数列， すなわちそれぞれが関係を持たない数の並びである．擬似乱数とは コンピュータが作る乱数である． 擬似乱数には様々な種類があるが，それを作る方法は一様乱数を 作ることに集約される．他の乱数は一様乱数を変数変換によって 求められる． 今回は擬似乱数を求めるアルゴリズム，擬似乱数ライブラリを使う方法， そして擬似乱数の応用例として数値積分について考えてゆく．

一様乱数

はじめに， $0\le \xi < 1$ の範囲で一様分布する乱数を作る方法を 考える．これを拡張して $\alpha \le \eta < \beta$ である一様乱数$\eta$を 作りたい場合には，上記の$\xi$を作ったあとで， 変数変換 $\eta = (\beta-\alpha) \xi +\alpha$ を行えばよい． 一様乱数は「線形合同法」と呼ばれる方法によって作られる．

$$x_{i+1}\equiv a x_{i} +b \quad ({\rm mod }m)$$ (1.1)

この方法によって，ランダムな数の並び $x_i$ が作られる．ただし， この $0 \le x_i < m$ である整数であることに注意する． この $x_i$ を $m$ で割ることによって，我々がほしい $0\le \xi < 1$ の範囲で一様分布する乱数 $\xi_i$ が作られる．

標準正規乱数

平均 0 ，分散 1 で正規分布する乱数を標準正規乱数という．ここで， 平均 $\mu$ 分散 $\sigma^2$ の正規分布の確率密度関数 $P(x;\mu, \sigma^2)$は，

$$P(x;\mu, \sigma^2)=\frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$ (1.2)

である． 標準正規乱数の生成には Box-Muller-Wiener アルゴリズムを用いる． これは，ふたつの$[0,1)$の一様乱数 $\xi_1, \xi_2$ を用いて， ふたつの標準正規乱数 $\eta_1, \eta_2$ を作るものである．ここで， $\xi_1, \xi_2$ と $\eta_1, \eta_2$ の関係は，

$$\eta_1= \sqrt{-2\log{\xi_1}} \sin (2\pi \xi_2)$$ (1.3)

$$\eta_2= \sqrt{-2\log{\xi_1}} \cos (2\pi \xi_2)$$ (1.4)

である． 平均 $\mu$ 分散 $\sigma^2$ で正規分布する乱数 $\eta$ を作りたい場合には， 標準正規乱数 $\xi$ を発生させたあとで， $\eta = \sigma \xi + \mu$ と 変数変換すればよい．

数値積分

擬似乱数の応用例として有名なのは数値積分である．例えば数値積分によって 単位円の面積を求めるためにはどうすればよいのかを考えてみよう． 下の絵は４分の１円を表している．ゼロから１までの一様乱数を24個 発生させて，２個を一組として座標をプロットしたのが，図中の○と×で ある．○は円の内部，×は円の外部にプロットされたことを示している． 全体の点の個数12個のうち，単位円の内部にあったものは９個．すなわち， この４分の１円の面積は 9/12=3/4 であり，これを４倍したものが 単位円の面積の近似値，すなわち 3 である．一般に，乱数によって 作られる点の数を２桁増やすと，近似値の精度は１桁増えると言われて いる．

上の過程をプログラムで表すにはどうしたらよいのかを考えてみよう． $i$ 番目に発生させたランダムな点を $(\xi_i, \eta_i)$ とする． 境界の方程式は $y=\sqrt{1-x^2}$ なので， $\sqrt{1-\xi_i^2}$ と $\eta_i$ の大きさを比較する．計算値の方が大きいときには， ○の個数をひとつ増やす．これを続けていき，最終的に○の個数と 全体の個数の比を取ればよい．

課題

課題1-1 乱数列の生成

線形合同法

$$x_{i+1}\equiv a x_{i} +b \quad ({\rm mod }m)$$ (1.5)

を用いて20個の「区間 [0,1]の一様乱数」を発生させるプログラムを作成しなさい． ただし $a=12869$， $b=6925$， $m=2^{14}$ として $x_0$ は自分で決定する ものとする．

課題1-2

drand48 ライブラリや rand ライブラリを用いて20個の「区間 [0,1]の一様乱数」を 発生させるプログラムを作成しなさい． ライブラリの使い方は ``man drand48'' コマンドなどを実行することによって 調べることができる．

課題1-3

drand48 ライブラリや rand ライブラリを用いて20個の標準正規乱数を発生させる プログラムを作成しなさい． $\xi_1, \xi_2$ と $\eta_1, \eta_2$ の関係は，

$$\eta_1= \sqrt{-2\log{\xi_1}} \sin (2\pi \xi_2) \eta_2= \sqrt{-2\log{\xi_1}} \cos (2\pi \xi_2)$$ (1.6)

である．

課題1-4

講義で示した手順で単位円の面積を求めるプログラムを作ろう． 発生させる点の数を多くするとπの値に近づいていくだろうか？

課題1-5

drand48 ライブラリや rand ライブラリを用いて

$$\int_0^1 x \exp (x-1) dx$$ (1.7)

の値を数値積分によって求めなさい．ちなみに非積分関数の不定積分の 結果は

$$(x-1) \exp (x-1)$$ (1.8)

であり，解析的な積分の結果は1/e=0.367879 となる．

プログラムの例

課題1-1

課題1-2

課題1-3

課題1-4

結果

課題1-1

 1 0.208142594
 2 0.846181989
 3 0.274430811
 4 0.857291102
 5 0.228407502
 6 0.619361520
 7 0.499786377
 8 0.780870438
 9 0.830739200
 10 0.552707076
 11 0.775926292
 12 0.208325699
 13 0.202527016
 14 0.583714843
 15 0.790331423
 16 0.577244699
 17 0.531404495
 18 0.649880946
 19 0.230482817
 20 0.325032055

課題1-2

 
 1 0.574890494 0.574890494
 2 0.310593992 0.310593992
 3 5.08008860e-02 5.08008860e-02
 4 0.370710492 0.370710492
 5 0.128306121 0.128306121
 6 0.294693470 0.294693470
 7 0.974656940 0.974656940
 8 5.69514558e-02 5.69514558e-02
 9 0.931047559 0.931047559
 10 0.934973478 0.934973478
 11 0.591664016 0.591664016
 12 0.996473730 0.996473730
 13 1.51172932e-02 1.51172932e-02
 14 0.983691454 0.983691454
 15 0.127344206 0.127344206
 16 0.772093415 0.772093415
 17 0.847713649 0.847713649
 18 0.998692393 0.998692393
 19 0.723982871 0.723982871
 20 0.216916546 0.216916546

課題1-3

 
 1 0.976870060
 2 -0.390993834
 3 1.77211404
 4 -1.67907572
 5 1.94711351
 6 -0.561624587
 7 7.93600827e-02
 8 0.212229669
 9 -0.150182709
 10 0.346893370
 11 -2.26973817e-02
 12 1.02426445
 13 -0.296179235
 14 2.88029671
 15 -2.01067281
 16 0.280923039
 17 -4.72256029e-03
 18 0.574806452
 19 0.786423862
 20 0.165869489

課題1-4