13. モンテカルロ法
13.1. Monte Carlo 法の基礎
多すぎて計算することができないアンサンブルから適切な方法によってサンプルを抽出し,抽 出されたサンプルの集合から対象のマクロな性質を明らかにしようとする手法をモンテカルロ法 という.
サンプルの抽出方法によって幾つかの方法に分類される.ランダムサンプリングは特に方針も なしにサンプリングを行うことをいう.マルコフ連鎖モンテカルロ(Markov Chain Monte Carlo)
は直前のサンプルを部分的に更新したものを新しいサンプルとする方法である.
13.2. 単純な例(ランダムサンプリング)
以下の式のような積分を計算することを考える.
d dr dx
x
0 4 1 0 1
0
1 2
第 2 項からわかるように,明らかに値はπ/4 である.従ってこの積分値を4倍するとπになる.
そこで,この積分値π/4が未知であると仮定して積分値を求めてみる.さらに,その値を 4 倍し た値を調べて,サンプリング数と値の精度について調べてみる.
この積分(の4倍)を計算するのに以下の方法を用いる.まず,和 Sn を導入する. 1 また は 0 であるような数を n 回足しあわせた結果である.初期的にはS0 0とする. Sn を求める ために,0 以上 1 未満の乱数 x と y を発生させる.y 1 x2 のとき,すなわち,x2 y2 1
の時は Sn Sn 1 とし,それ以外のときは Sn はそのままであるとする(すなわちゼロを加え るとする).更に 4Sn /n を求める.
13.3. プログラム例
13.4. 結果
3.07 3.08 3.09 3.1 3.11 3.12 3.13 3.14 3.15
0 20000 40000 60000 80000 100000
πの推定値
πの推定値 Sub monteCarlo()
'int_0^1 sqr(1-x^2)dx の積分値を求める.
Const max As Long = 100000 Const div As Long = 1000
Dim s As Double, xi As Double, eta As Double, cur As Integer
s = 0
ActiveCell.Value = "回数"
ActiveCell.Offset(0, 1) = "積分値"
ActiveCell.Offset(0, 2) = "πの推定値"
For i = 1 To max
If Rnd < Sqr(1 - Rnd ^ 2) Then s = s + 1#
End If
If i Mod div = 0 Then cur = i / div
ActiveCell.Offset(cur).Value = i
ActiveCell.Offset(cur, 1).Value = s / i ActiveCell.Offset(cur, 2).Value = 4 * s / i End If
Next i End Sub
13.5. 演習
A) プログラム例の計算結果を20個ほど並べて統計をとってみよう B) 別の積分を計算してみよう(下に行くほど工夫が必要)
(ア) 1 x xdx
0 2
(イ) 0sin xdx
(ウ) 0 2cos 2 xdx
(エ) 2ex dx
0
2
14. ランダムウォーク
14.1. ランダムウォーク
格子点上に各時間ステップに1格子間隔だけ動く粒子がある(これをランダムウォーカーと呼 ぶ).
すべての方向について,ランダムウォーカーがその方向へ行く確率が与えられている.n 時間 ステップ経過したときに,与えられた位置にランダムウォーカーいる確率はそれぞれいくらにな るか.
14.2. ランダムウォークの例
14.2.1. 1D(横軸は時間)
14.2.2. ランダムウォークの例(複数)
14.2.3. ランダムウォークの例(2D)
14.3. 1次元ランダムウォークの確率分布
水平線上に単位時間あたり1歩動くランダムウォーカーについて考える.右に行く確率を p
左に行く確率を (1 p) とすると, n 歩動いたあとで,スタート地点から右に r だけ離れた 場所にいる確率 p(n,r)は,
2 2
2
) 1 ( )
, (
r n r
n r
n n p p
r n
p C
1次元ランダムウォークの平均と分散を考える.一歩移動するときの平均(期待値)と分散は 以下のように求められる.
1 2 ) 1 ( ) 1 ( ) 1
(
p p p
m
2
2 2
2 2
) 1 2 ( 1
) 1 2 ( ) 1 ( ) 1 ( )
1 (
p
p p
p
n 歩移動するときの平均(期待値)と分散は以下のようになる.
) 1 2 ( ) , (
0
p n r n p r m
n
r
) 1 ( 4 ) , ( ) (
0
2
2 r m p n r np p
n
r
これは歩数が増えるたびに中心の位置が2p 1だけ移動し,バラツキの程度が4 p1 pだけ 大きくなることを示している.特に,
2
1
p のとき,中心位置は原点のままでバラツキの程度はn
になる.
14.4. 分布の例
14.4.1. p が一定で歩数が変わる p=0.5
右に行く確率と左に行く確率が同じなので y 軸について対称な分布になっている.
14.4.2. p が一定で歩数が変わる p=0.6
右に行く確率の方が高いので,歩数が大きくなると分布の中心が右にシフトしている.
分散も歩数が大きくなるほど大きくなっている.
14.4.3. 分布の例(同じ歩数で異なる p)
p を変えると分布がことなる.ピークの高さと分散の大きさも変化していることに注意.
14.5. プログラム例
14.6. 演習
A) 統計を取りなさい.
B) p=0.6 の場合のグラフを描きなさい.
C) p=0.6 の場合の統計を取りなさい.
Sub randomWalk1D() '1次元ランダムウォーク
Const max As Long = 10000 Const prob As Double = 0.5 Dim pos
ActiveCell.Value = "位置(p=" & prob & ")"
ActiveCell.Offset(1).Value = 0
For i = 2 To max + 1 If Rnd < prob Then pos = pos + 1 Else
pos = pos - 1 End If
ActiveCell.Offset(i).Value = pos Next i
End Sub