9．モンテカルロ法

９．モンテカルロ法
元々は，決定論的な数学上の問題を乱数を用いて解くため
の方法。現在では，次のような問題に適用される。
①確率的現象に対して乱数を用いてシミュレーション
②複雑すぎる問題の場合，想定できないパラメータを乱数と
してみなしてシミュレーション。（ランダムウォーク）
ちょっと一息…
モンテカルロって？
モナコ公国の観光・保養地の地名。カジノ（賭博場）や，
世界ラリー選手権が行われることで有名。
乱数を利用する方法が，
カジノでの賭博に似ていることから，
乱数を使った方法はモンテカルロ法と呼ばれる。
９．１モンテカルロ法の計算方法
（１）面積の計算
［方法］
①領域を完全に含むことができる正方形を考え
②この範囲の乱数を発生させ，
③領域内に入った場合の回数をカウントして
④その割合を面積とする。
M 回が領域中にあった場合
正方形の１辺の長さを L とすると
N
回の試行中
y
ML2
面積＝
N
－1
1
0
－1
x
1
プログラム例
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
４倍したものは，πに近づくはずである。
結果例
（２）乱数の発生方法
①乗積合同法（multiplicative congruential method）
a, b, L を正の整数として，次の式により乱数を求める。
X 0  a, X 1  bX 0 mod L, , X n  bX n1 mod L
整数の乱数を求める方法であるが，これを実数化することで実数の乱数に
直すことができる。
以下に示す方法は，
①
Seed1, Seed 2, Seed 3 に１から 30,000 までの任意の整数を与える。
（便宜上，VBの乱数を使って，これらを設定している）
② 乗積合同法により３系統の整数の乱数を求める。
③ 法で除すことで 0 から 1 の実数に直し，和の小数部分をとる。
乗積合同法を利用した実数の一様乱数（その１）
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
乗積合同法を利用した実数の一様乱数（その２）
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
②機械エプシロンを用いる方法
（最小の（１＋ε）ー１＞の値）
機械エプシロンを求め，これを最小単位として乱数を求める。
整数の乱数でも使用される引き算法（subtractive method）を使用している。
X n   X n55  X n24  mod L
機械エプシロン：浮動小数点で最も小さな正の値。
［定義］実数 (1   )  1 がとりうる最小の正の値
機械エプシロンを最小単位とする方法（その１）
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
機械エプシロンを最小単位とする方法（その２）
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
機械エプシロンを最小単位とする方法（その２）
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
機械エプシロンを最小単位とする方法（その２）
Public Function Random2() As Double
RandomIndex = RandomIndex + 1
If RandomIndex > 55 Then
randomTableSet
RandomIndex = 1
End If
Random2 = RandomTable(RandomIndex)
End Function
（３）乱数の検定
①頻度検定
[0, 1] の全区間を分割し，乱数の総数を N ，分割された区間を [ai , bi ]
とすると，区間 [ai , bi ] の期待値は，
Ei  N  bi  ai 
実際の度数を
Ni
として，次のような量を考える。
2
(
E

N
)
i
2   i
Ei
i 1
n
この値が大きければ，一様性が低いことになる。
 2 検定（カイ二乗検定）という。
カイ二乗検定のプログラム
（乱数によって赤字部分を変更する）
（この例は，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
比較結果例
（ただし，実行によって大小関係は変わるので，この例だけで判断しないこ
と）
②連検定
[0, 1] で発生させた擬似乱数に対して増減の連続性等について調べる方法
たとえば，０．５以上の値を＋，０．５未満を－として，同種符号が続く回数を
カウントする。
この連の長さCは，以下の確率に従うので，この分布に近いかどうかを検定
する。
1
PC    
2
C
連検定のプログラム
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
検定結果例
理想とする分布の赤線に非常に接近しているので
連検定の結果は良好であるといえる。
（４）分布に従う乱数の発生
①直接法の考え方
確率変数
pi
が確率密度 G (x ) を持つ場合，
qi   Gx dx
pi

qi
を区間［0, 1］で発生させ，上式を満たす pi を求めると，確率密度 G (x )
は，区間［0, 1］で一様分布になるはずである。したがって，一様乱数
を満たす乱数列になる。さらに，次の関係が成立する。
qi  G  pi  
［例］
Gx   2 x
pi  G 1 qi 
( x  0)
log y
y  x log 2  x 
log 2
として
この分布を持つ乱数を得ることができる。
y に一様乱数を用いると
②棄却法
区間［0, 1］で一様乱数
G ( x)  D をみたす定数 D を定義しておく。
区間［0, 1］で一様乱数
pi , qi
を発生させ，
x  qi , y  Dpi
とし，
y  Dpi  G(x)
を満たせば， pi を採用し，それ以外の場合は，
採用できるまで繰り返す方法。
③合成法
G ( x)   f y ( x)dh( y )
のように，サンプリング可能な分布関数
h( y ) における条件付密度関数
f y (x) で表されている場合， h( y ) に対応する確率変数 y
を直接法，棄却法で求め，その
乱数を発生させる。
を
y に対して， f y (x) にあうよう
（５）正規分布の乱数の発生方法
（考え方）プログラムは，１０－２に示す
①中心極限定理
独立な乱数
a1 , a2 ,, an の和は，n
が大きくなると，正規分布に近づく。
特に［0, 1] で一様分布する乱数を a1 , a2 , , an とすると，
 n
n  12
x    a j  
2 n
 j 1
は，平均値０，分散１の正規分布に近づく。
ただし，n  12 程度でも，十分正規分布とみなす程度の結果が得られる。
すなわち，次のような単純な式で結果が得られる。
 12
12  12 12
x    a j  
 aj  6
2  12 j 1
 j 1
（考え方）プログラムは，１０－２に示す
②ボックス・マラー（Box-Muller）の方法
区間［0, 1] の一様乱数 a1 , a2 を用いた
x1   2 log a1  cos( 2a2 )
x2   2 log a1  sin( 2a2 )
の
x1 , x2
は平均０，分散１に従うことが知られているので，この方法
を用いる。

Download Report