[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 に入れています。


 

コメントをどうぞ

メールアドレスが公開されることはありません。