９．モンテカルロ法

1. ９．モンテカルロ法 元々は，決定論的な数学上の問題を乱数を用いて解くための方法。現在では，次のような問題に適用される。 ①確率的現象に対して乱数を用いてシミュレーション ②複雑すぎる問題の場合，想定できないパラメータを乱数と してみなしてシミュレーション。（ランダムウォーク）

2. ちょっと一息…モンテカルロって？ モナコ公国の観光・保養地の地名。カジノ（賭博場）や， 世界ラリー選手権が行われることで有名。 乱数を利用する方法が， カジノでの賭博に似ていることから， 乱数を使った方法はモンテカルロ法と呼ばれる。

3. ９．１　モンテカルロ法の計算方法（１）面積の計算９．１　モンテカルロ法の計算方法（１）面積の計算 ［方法］ ①領域を完全に含むことができる正方形を考え ②この範囲の乱数を発生させ， ③領域内に入った場合の回数をカウントして ④その割合を面積とする。 回が領域中にあった場合 回の試行中 正方形の１辺の長さを とすると 1 面積＝ 0 －1 1 －1

4. プログラム例 Sub ボタン1_Click() Randomize 159 T = 0 With Worksheets("Sheet1") For i = 1 To 10000 R1 = (Rnd - 0.5) * 2: R2 = (Rnd - 0.5) * 2 If (R1 * R1 + R2 * R2) < 1 Then T = T + 1 .Cells(i + 1, 1) = i - 1 .Cells(i + 1, 2) = T .Cells(i + 1, 3) = (T / i) * 4 ‘確認用 Next End With End Sub

5. 結果例 ４倍したものは，πに近づくはずである。

6. （２）乱数の発生方法①乗積合同法（multiplicative congruential method） を正の整数として，次の式により乱数を求める。 整数の乱数を求める方法であるが，これを実数化することで実数の乱数に直すことができる。 以下に示す方法は， に１から 30,000 までの任意の整数を与える。 ① （便宜上，VBの乱数を使って，これらを設定している） ②　乗積合同法により３系統の整数の乱数を求める。 ③ 法で除すことで 0 から 1 の実数に直し，和の小数部分をとる。

7. 乗積合同法を利用した実数の一様乱数（その１）乗積合同法を利用した実数の一様乱数（その１） Private Seed1 As Integer Private Seed2 As Integer Private Seed3 As Integer Private Function SetSeed() SetSeed = Rnd * 30000 + 1 ‘最初の３シード値はVBで設定 If SetSeed > 30000 Then SetSeed = 30000 End Function Public Sub initRandom() Randomize Seed1 = SetSeed(): Seed2 = SetSeed(): Seed3 = SetSeed() End Sub

8. 乗積合同法を利用した実数の一様乱数（その２）乗積合同法を利用した実数の一様乱数（その２） Public Function Random() As Double Dim R As Double ‘ 16ビット整数でも可能になるように記述 '以下２行 : Seed1 =171 *Seed1 mod 30269 と同じ Seed1 = (Seed1 Mod 177) * 171 - (Seed1 \ 177) * 2 If Seed1 < 0 Then Seed1 = Seed1 + 30269 '以下２行 : Seed2 =172 *Seed1 mod 30307 と同じ Seed2 = (Seed2 Mod 176) * 172 - (Seed2 \ 176) * 35 If Seed2 < 0 Then Seed2 = Seed2 + 30307 '以下２行 : Seed3 =170 *Seed3 mod 30323 と同じ Seed3 = (Seed3 Mod 178) * 170 - (Seed3 \ 178) * 63 If Seed3 < 0 Then Seed3 = Seed3 + 30323 '上記３系統の乱数を0～1の実数に変換 R = Seed1 / 30269# + Seed2 / 30307# + Seed3 / 30323# Random = R - Int(R) End Function

9. ②機械エプシロンを用いる方法（最小の（１＋ε）ー１＞の値）②機械エプシロンを用いる方法（最小の（１＋ε）ー１＞の値） 機械エプシロンを求め，これを最小単位として乱数を求める。 整数の乱数でも使用される引き算法（subtractive method）を使用している。 機械エプシロン：浮動小数点で最も小さな正の値。 ［ 定義］実数 がとりうる最小の正の値

10. 機械エプシロンを最小単位とする方法（その１）機械エプシロンを最小単位とする方法（その１） Private RandomIndex As Integer Private RandomTable(55) As Double Private Function MachineEP() As Double ‘機械エプシロンを求める Dim E As Double, ED As Double E = 1 Do ED = E: E = E / 2 Loop Until (E + 1) = 1 E = 1 + ED: MachineEP = E - 1 End Function

11. 機械エプシロンを最小単位とする方法（その２）機械エプシロンを最小単位とする方法（その２） Private Sub randomTableSet() ‘乱数表設定 Dim I As Integer, R As Double For I = 1 To 24 R = RandomTable(I) - RandomTable(I + 31) RandomTable(I) = R: If R < 0 Then RandomTable(I) = R + 1 Next For I = 25 To 55 R = RandomTable(I) - RandomTable(I - 24) RandomTable(I) = R: If R < 0 Then RandomTable(I) = R + 1 Next End Sub

12. 機械エプシロンを最小単位とする方法（その２）機械エプシロンを最小単位とする方法（その２） Public Sub initRandom2() ‘　乱数初期化 Dim I As Integer, J As Integer Dim R As Double Randomize: Seed = Rnd RandomTable(55) = Seed R = MachineEP() For I = 1 To 54 J = (21 * I) Mod 55 RandomTable(J) = R R = Seed - R: If R < 0 Then R = R + 1 Seed = RandomTable(J) Next randomTableSet‘以下３行はウォームアップ用 randomTableSet randomTableSet RandomIndex = 55 End Sub

13. 機械エプシロンを最小単位とする方法（その２）機械エプシロンを最小単位とする方法（その２） Public Function Random2() As Double RandomIndex = RandomIndex + 1 If RandomIndex > 55 Then randomTableSet RandomIndex = 1 End If Random2 = RandomTable(RandomIndex) End Function

14. （３）乱数の検定①頻度検定 ，分割された区間を とすると，区間 の期待値は， の全区間を分割し，乱数の総数を 実際の度数を として，次のような量を考える。 この値が大きければ，一様性が低いことになる。 検定（カイ二乗検定）という。

15. カイ二乗検定のプログラム（乱数によって赤字部分を変更する）カイ二乗検定のプログラム（乱数によって赤字部分を変更する） （この例は，VBに用意されたランダム関数を検定する方法である） Sub ボタン1_Click() Dim X(10) As Double Randomize M = 10: N = 100000 With Worksheets("Sheet1") For J = 1 To 100 For I = 0 To M - 1: X(I) = 0: Next For k = 1 To N ID = Int(Rnd * M): X(ID) = X(ID) + 1 Next E = N / M: Kai = 0: Kai = 0 For I = 0 To M - 1 EE = E - X(I): Kai = Kai + EE * EE / E Next .Cells(J + 1, 2) = Kai Next End With End Sub

16. 比較結果例 （ただし，実行によって大小関係は変わるので，この例だけで判断しないこと）

17. ②連検定 たとえば，０．５以上の値を＋，０．５未満を－として，同種符号が続く回数をカウントする。 この連の長さCは，以下の確率に従うので，この分布に近いかどうかを検定する。 で発生させた擬似乱数に対して増減の連続性等について調べる方法

18. 連検定のプログラム （乱数によって赤字部分を変更する） Sub Sheet4_ボタン1_Click() Dim X(100) As Double For I = 0 To 100 X(I) = 0 Next Randomize A = Rnd: NS = 10000 For I = 1 To NS N = 0 If A >= 0.5 Then Do While A >= 0.5: N = N + 1: A = Rnd: Loop Else Do While A < 0.5: N = N + 1: A = Rnd: Loop End If If N <= 100 Then X(N) = X(N) + 1 Next With Worksheets("Sheet4") For I = 1 To 100 .Cells(I + 1, 1) = I .Cells(I + 1, 2) = X(I) / NS Next End With End Sub

19. 検定結果例 理想とする分布の赤線に非常に接近しているので 連検定の結果は良好であるといえる。

20. （４）分布に従う乱数の発生①直接法の考え方（４）分布に従う乱数の発生①直接法の考え方 が確率密度 を持つ場合， 確率変数 は，区間［0, 1］で一様分布になるはずである。したがって，一様乱数 を区間［0, 1］で発生させ，上式を満たす を求めると，確率密度 を満たす乱数列になる。さらに，次の関係が成立する。 ［例］ として に一様乱数を用いると この分布を持つ乱数を得ることができる。

21. ②棄却法 を定義しておく。 をみたす定数 区間［0, 1］で一様乱数 区間［0, 1］で一様乱数 を発生させ， を採用し，それ以外の場合は， を満たせば， とし， 採用できるまで繰り返す方法。

22. ③合成法 のように，サンプリング可能な分布関数 における条件付密度関数 で表されている場合， を に対応する確率変数 に対して， を直接法，棄却法で求め，その にあうよう 乱数を発生させる。

23. （５）正規分布の乱数の発生方法（考え方）プログラムは，１０－２に示す（５）正規分布の乱数の発生方法（考え方）プログラムは，１０－２に示す 独立な乱数 が大きくなると，正規分布に近づく。 の和は， ①中心極限定理 とすると， 特に［0, 1] で一様分布する乱数を は，平均値０，分散１の正規分布に近づく。 程度でも，十分正規分布とみなす程度の結果が得られる。 ただし， すなわち，次のような単純な式で結果が得られる。

24. （考え方）プログラムは，１０－２に示す ②ボックス・マラー（Box-Muller）の方法 区間［0, 1] の一様乱数 を用いた の は平均０，分散１に従うことが知られているので，この方法 を用いる。