[VBA] モンテカルロ法

 この記事では VBA で モンテカルロ法 (Monte Carlo method) を実装して、円周率の近似値や定積分を求める方法を解説します。

[VBA] モンテカルロ法による円周率の計算

 下の図のように半径 1 の四分円と、それに外接する正方形を考えます。

 モンテカルロ法による円周率の計算

 そして、この正方形の中に無作為に点を N 回発生させます。そのうちの何割かは四分円の中に入ることになります。その総数を Nc とします。Nc と N の比は、四分円の面積 Sc と正方形の面積 S の比にほぼ等しくなっているはずだと考えます。つまり

Nc/N = π/4

ですから、円周率は

π = 4(Nc/N)

と計算できるはずです。N を大きくするほど近似の精度は良くなります。とりあえず N を 100000 と設定しておきます。

'[VBA] モンテカルロ法による円周率の近似計算

Sub Monte_Carlo_Pi()

  Dim x As Double, y As Double
  Dim r As Double, p As Double
  Dim n As Long, nc As Long

  Randomize

  For n = 1 To 100000

    '無作為な点 (x, y) を1つ発生させます
    x = RND
    y = RND

    '原点からの距離の2乗を計算します
    r = x ^ 2 + y ^ 2

    '原点からの距離が 1 以内ならば
    '四分円に入ったカウント数に 1 を加えます
    If r <= 1 Then
      nc = nc + 1
    End If

  Next n

  '円周率の近似値を計算します
  p = 4 * nc / n

  Debug.Print p

End Sub

 Monte_Carlo_Pi() を実行するとイミディエイトウィンドウに

3.14356856431436

という円周率の近似値が表示されます。乱数を使っているので、実行する度に違った値が表示されます。

 Monte_Carlo_Pi() の構造について簡単に解説します。
 n は無作為な点の発生回数なので、For…Next ステートメントのカウンタ変数の上限値となっています。Monte_Carlo_Pi() では

For n = 1 To 100000

としてあるので、上限値を 100,000 としてありますが、より高い精度を求めるのであれば、この上限値を書き換えてください。x と y に 0 から 1 の乱数を入れていますが、これが発生した点の座標 (x, y) となっています。

r <= x ^ 2 + y ^ 2

という記述で (x, y) の原点からの距離の2乗を変数 r に入れます。

If r <= 1 Then

は r が 1 以内であること、つまり点が四分円の中にあるという条件で、それが満たされた場合に限って

nc = nc + 1

として四分円に入ったカウント数 nc に 1 を加えます。最後に nc と n の比率を使って

p = 4 * nc / n

として円周率の近似値を p に入れています。

[VBA] モンテカルロ法による定積分の計算

 関数 $f(x)=x^2$ の定積分
 
\[I=\int_{0}^{a}x^2dx\]
モンテカルロ法 で計算してみます。

 VBAモンテカルロ法による定積分

 図にある面積 $a\times a^2=a^3$ の長方形の中にランダムな点 $(x,\:y)$ を $N$ 個発生させて、そのうち黄色の網掛け部分、すなわち
 
\[y\leq f(x)\]
に入った回数が $n$ 回であったとすると、
 
\[I=\int_{0}^{a}x^2dx\simeq\frac{n}{N}a^3\]
というように近似計算できます。モンテカルロ法による定積分の計算プログラムの一例を載せておきます。

'[VBA] モンテカルロ法による定積分の計算

Sub MonteCarlo_Integral()

  Dim x As Double, y As Double
  Dim a As Double, itg As Double
  Dim n As Long, k As Long, ct As Long

  '積分の上限値
  a = 3

  Randomize

  '試行回数の設定
  n = 1000

  For k = 1 To n

    '無作為な点 (x, y) を1つ発生させます
    x = a * RND
    y = a ^ 2 * RND

    'y < f(x) ならば
    If y < x ^ 2 Then

      'カウント数に 1 を加えます
      ct = ct + 1

    End If

  Next k

  '積分近似値を計算します
  itg = ct * a ^ 3 / n

  Debug.Print itg

End Sub

 真値 (9) と比較しやすいように、上限を $a=3$ としてあります。試行回数は $N=10^3$ です。マクロを実行すると 9 の近傍値が現れるはずです。$N$ の値を大きくするほど精度は上がりますが、そのぶん処理に時間がかかります。$N=10^5$ でも 3 桁ほどの精度しかないので、定積分をわざわざこの手法で行うメリットはありませんが(シンプソンの公式など他にいくらでも有用な方法があります)、モンテカルロ法は色々な分野のシミュレーションに応用がきくので、そのアルゴリズムを知っておいて損はありません。

コメント

タイトルとURLをコピーしました