モンテカルロ法(乱択アルゴリズム)｜Excel VBAで学ぶ数学とアルゴリズム

今回は乱択アルゴリズム「モンテカルロ法」について解説していきます。
このモンテカルロ法、あえて処理の中で乱数を使ってランダム性を出すというアルゴリズムなのですが、ランダム性を出しているのにも関わらず最終的にはある1つの答えが浮かび上がってくるという、普通のアルゴリズムとは違った考え方を持つ面白いアルゴリズムです。

本ページではこの「モンテカルロ法」のアルゴリズムがどういった仕組みなのかを解説し、Excel VBAのコードとして実装していきます。

こういったアルゴリズムはPythonやC++などのがっつりしたプログラム言語では紹介されていることがあります。これは実際にその言語で使われるというのもそうですが、プログラムのコードとして見ることでアルゴリズムの中身をより深く理解することができるためです。

ここではそういった”がっつりとしたプログラミング言語”ではないVBAを使うことで、幅広い方に向けてモンテカルロ法のアルゴリズムの内容をご紹介していきます。

1. モンテカルロ法とは
2. モンテカルロ法で円周率を求める
- 2.1. サンプルコード
3. まとめ

モンテカルロ法とは

乱択アルゴリズム(Randomized Algorithm)とは

乱択アルゴリズムとはその名の通りでアルゴリズムの中に”乱数”、すなわちランダム性を取り入れたアルゴリズムです。「正しい答えを導き出すためのアルゴリズムなのにランダム性がある？」と疑問に思う方も多いと思いますが、場合によってはランダム性を入れたほうが処理時間の効率が良くなるケースもあるという考えが乱択アルゴリズムです。

$$A=Array(a,a,a,a,a,b,b,b,b,b)$$

たとえば上記の要素を持った配列$A$があったとします。
この配列は $a$ という要素を5つ、$b$ という要素を5つ持っていますがその順序は決まっていないものとします。このとき、ある1つの要素 $a$ を取り出すにはどのような処理を行うと効率が良いでしょうか？

最も単純な方法は配列の1番目の要素から順に確認していくという方法です。
しかしこの方法、すべての要素 $a$ が配列の後ろに固まっている場合を考えると処理時間は長くなってしまいます。これは逆から見た場合も、1つ飛ばしに見た場合にも同じことが言えます。

そもそものお話、配列内の要素の順番が決まっていない時点で、要素の組み合わせがどのような状態でも常に高速であると保証できるアルゴリズムは存在しません。

しかしこの確認する要素の順番をランダムにすることで、2つの要素 $a$ と $b$ の要素数が同じであれば、どのような並びの要素になっていたとしても確実に1/2の確率で要素 $a$ を取り出すことができます。これは5つの黒い球と5つの白い球の計10個の球が入った袋から黒い球を取り出す確率と同じであると考えるとイメージが付きやすいと思います。

このように「確率的に見たときに結果が良くなるものに対して処理の中にランダム性を持たせよう」というのが乱択アルゴリズムの考え方です。しかし、その確率において最悪の場合(たとえば1/2の確率なのに100回連続間違ってしまうケースなど)が起こってしまう可能性もあるという点だけは注意が必要です。
　　

モンテカルロ法(Monte Carlo method)とは

「モンテカルロ法」とは乱択アルゴリズムの一種で、数値計算やシミュレーションなどをあえて乱数(ランダム性)を使って行う手法のことです。この手法は主に確率を近似的に求める手法として使われることが多いです。$n$ 回シミュレーションを行って、ある事象が $m$ 回起これば、その事象の起こる確率は $m/n$ に近似されます。このとき試行回数が少なければ近似は荒く、試行回数が多ければより近似に近づいていきます。

具体例として「コインを投げたときに表が出る確率」を求めたいとします。
この例の場合、表と裏の2択のうちのいずれかが出るため、確率は1/2である「50%」と簡単に導き出すことができますが、ここではあえてモンテカルロ法を使ってこの確率の近似値を求めてみます。

この確率を調べるには実際にコインを投げて、下記のように計算すれば求めることができます。

コインの表が出る確率＝表が出た合計回数 ÷ コインを投げた総数

下図はVBAを使ってコイントスをシミュレーションしたもので、コインを投げる回数を「100回」「1000回」「10000回」と変化させていった時に、コインの表が出る確率を折れ線グラフ化したものです。横軸が「コインを投げた回数」、縦軸が「コインの表が出る確率」となっています。

このグラグからコインを投げる回数が増えるにしたがって、徐々にコインの表の出る確率が「0.5」、つまりは50%に収束していっていることがわかります。

このように入力される値が乱数であったとしても試行回数を増やしていけばある一定の数値に近似していきその確率を求めることができる、これがモンテカルロ法のアルゴリズムです。

ただ、注意しないといけない点はモンテカルロ法によって導き出された値が確実に正解とはいえない点です。たとえばこのコイントス、1000回投げて表が1回も出ないということもありえますが、こういったケースは極めて稀で基本的には正しい答えを導き出すことができるといわれています。ただ、アルゴリズムとして乱数を使っている以上”絶対はない”ということは頭に入れておきましょう。

上図のグラフ作成に使ったコイントスの確率を求めるVBAコードは下記の通りです。
Rnd関数を使っていことから、アルゴリズムとして乱数が導入されていることがわかります。

Const NUMBER_OF_TRIALS As Long = 100 'ｺｲﾝを投げる回数
Sub CoinTossProbability()

    Dim iRandNum As Long
    Dim iSum1 As Long

    With ThisWorkbook.ActiveSheet
        
        For i = 1 To NUMBER_OF_TRIALS
        
            'ﾗﾍﾞﾙ
            .Cells(1, i).Value = i
            
            '乱数 [0](裏) or [1](表)
            iRandNum = Int(2 * Rnd)
            
            '[1](表)が出た合計回数
            If iRandNum = 1 Then
                iSum1 = iSum1 + 1
            End If
            
            'これまでに[1](表)が出た確率
            .Cells(2, i).Value = iSum1 / i
        Next i
        
        '折れ線グラフ作成
        With ActiveSheet.Shapes.AddChart.Chart
            .ChartType = xlLine
            .SetSourceData Range(Cells(2, 1), Cells(2, NUMBER_OF_TRIALS))
            .Axes(xlValue, 1).MaximumScale = 1  'ｸﾞﾗﾌ最大値を[1]に設定
            .Axes(xlValue, 1).MajorUnit = 0.1   'ｸﾞﾗﾌ目盛りを[0.1]単位に設定
        End With
    
    End With
    
End Sub

Const NUMBER_OF_TRIALS As Long = 100 'ｺｲﾝを投げる回数

Sub CoinTossProbability()

Dim iRandNum As Long

Dim iSum1 As Long

With ThisWorkbook.ActiveSheet

For i = 1 To NUMBER_OF_TRIALS

'ﾗﾍﾞﾙ

.Cells(1, i).Value = i

'乱数 [0](裏) or [1](表)

iRandNum = Int(2 * Rnd)

'[1](表)が出た合計回数

If iRandNum = 1 Then

iSum1 = iSum1 + 1

End If

'これまでに[1](表)が出た確率

.Cells(2, i).Value = iSum1 / i

Next i

'折れ線グラフ作成

With ActiveSheet.Shapes.AddChart.Chart

.ChartType = xlLine

.SetSourceData Range(Cells(2, 1), Cells(2, NUMBER_OF_TRIALS))

.Axes(xlValue, 1).MaximumScale = 1 'ｸﾞﾗﾌ最大値を[1]に設定

.Axes(xlValue, 1).MajorUnit = 0.1 'ｸﾞﾗﾌ目盛りを[0.1]単位に設定

End With

End Sub

　
ちなみに名称の”モンテカルロ”についてですが、カジノで有名なモナコ公国の地区の名前が由来とされています。実はこのモンテカルロ法はギャンブラーの間では有名な手法で、かつてこの手法をつかったベット(賭け)方法でカジノを1つ破産させたというお話があるほど強力なアルゴリズムです。

ギャンブルで使われるモンテカルロ法はベットする額を数列と考え、勝ち負けの結果をもとにベットする額を調整していく手法です。かなりざっくりというのであれば「初めに$4賭けて負けても、次は$5賭けて勝てば収支は+$1でプラスだよね」というようなイメージです。
(※興味のある方は「ルーレットモンテカルロ法」などで検索してみて下さい)

モンテカルロ法で円周率を求める

直径に対する円周の長さの比率である円周率「$\pi$」は定数ですが、コイントスのときと同じように乱数を使ったモンテカルロ法でも求めることができます。$\pi$ を求める処理の流れは下記の通りです。

① 「1.0×1.0」の正方形内のランダムな位置に点を作成する　　
② 点が正方形内にピッタリ収まるような四分円の内側に作成されたかを判定する
　 (= 作成した点と原点の距離が1以下かを判定する)
　点が四分円の内側に作成された場合は点を配列 $p$ 内に入れる　
③ 上記の①と②を $n$ 回繰り返す
④ $p$ の要素数を4倍して生成した点の総数で割る($\pi=4p/n$)　　

Pi_30K.gif ©Nicoguaro (Licensed under CC BY 3.0)

この処理は上のgif画像のように「1×1」の正方形とそれピッタリ収まる四分円(1/4の円)を使います。
この「1×1」の正方形内にランダムで点を作成していき、点が四分円の内側にできる割合を求めることで円周率を求めることができます。そして作成する点の個数 $n$ が増えれば増えるほど $\pi$ の値が定数「3.141592…」の数値に収束していきます。

何故そうなるのか、まずは正方形と四分円の面積の割合を考えてみます。
面積の割合とは見方を変えると正方形内でランダムに点を作成した場合に四分円の内側に入る確率であると考えることもできます。
　
このとき「四分円の内側に入る確率」は「四分円の内側にはいった点の個数：$p$」を「全ての点の個数：$n$」で割った「$p/n$」で表すことができます。
そして、正方形と四分円の面積はそれぞれ下記のように表すことができます。

正方形の面積　：　$S_{1}={1}\times{1}=1$
四分円の面積　：　$S_{2}=({1}\times{1}\times{\pi})\div{4}=\pi/4$

つまり、上記の内容をまとめると下記のような関係で表すことができます。

面積と確率の関係：　$S_{1}:S_{2}=1:\pi/4=n：p$

この関係式を整えると最終的に「$\pi=4p/n$」と変形することができます。
よって、ランダム生成した点が四分円の中に入る確率を4倍したものが「$\pi$」であるといえます。

ここまで理解できればあとはコイントスの例と同じ要領です。
ランダム生成した点が四分円の中に入る確率は実際に点を作成して確認していきます。このとき、点は作れば作るほど”ある一定の数値”に収束していきます。その収束した数値に4をかければ「$\pi$」(の近似値)を求めることができるという訳です。

　
点が円の内側か外側かを判定する

四分円の内側と外側どちらにできた点なのかを判定するには作成した点と原点の距離が1以下かを調べることで判定することができます。このとき距離が1以下の場合は四分円の内側、1より大きい場合は外側となります。

「1×1」の正方形内に収まるような四分円とはつまり、単位円(半径1の円)の四分円です。
単位円は方程式で下図のように表すことができます。
これは直角三角形の辺の長さを求めるピタゴラスの定理からも導き出すことができます。

このことから点の座標と原点の距離である「r」を求めたときに、1以下であれば四分円の内側、1より大きければ四分円の外側にできた点であると判定することができます。

サンプルコード

円周率を求めるモンテカルロ法のアルゴリズムをVBAで書くと下記の通りです。
Excel VBAでの実装ということでShapeオブジェクトの「円」を点として扱います。
試行回数を増やせば増やすほど正しい円周率の値に近づいていきますが、その分処理の時間も増えていくので注意して下さい。(※上画像は試行回数10000回で処理時間は約10分近く)
また、視覚化するために点(円)を作成していますが、その処理を消すことで高速化が図れます。

Const NUMBER_OF_TRIALS As Long = 1000 '試行回数
Const MAGNIFICATION As Long = 100     '倍率(目安:試行回数の1/10) これを小さくするほど作成されるｴﾘｱが狭まっていく
Sub CoinTossProbability()

    Dim dX As Double
    Dim dY As Double
    Dim iSum As Long
    Dim oPt As Shape
    Dim dPI As Double
    Dim vR(0 To 2) As Long
    Dim vB(0 To 2) As Long
    
    '色設定
    vR(0) = 255: vR(1) = 0: vR(2) = 0
    vB(0) = 0: vB(1) = 0: vB(2) = 255

    With ThisWorkbook.ActiveSheet
        
        For i = 1 To NUMBER_OF_TRIALS
            DoEvents
            
            '乱数生成 [0～1]
            dX = Rnd
            dY = Rnd
            
        　　'<<円の内側と外側のどちらに含まれているかを判定>>
        
            '円の内側の場合
            If dX * dX + dY * dY <= 1 Then
            
                '円の内側に入った総数をｶｳﾝﾄ
                iSum = iSum + 1
                
                '青色の点を作成(点が大きい場合は高さと幅を調節 現状:5×5)
                Set oPt = CreateCircle(dX * MAGNIFICATION, dY * MAGNIFICATION, 5, 5, vB)
            
            '円の外側の場合
            Else

                '赤色の点を作成(点が大きい場合は高さと幅を調節 現状:5×5)
                Set oPt = CreateCircle(dX * MAGNIFICATION, dY * MAGNIFICATION, 5, 5, vR)
            End If
            
            
            '円周率の近似値
            dPI = 4 * iSum / i
        Next i

    End With
    
    
    MsgBox "π≒" & dPI
    
End Sub

'---------------------------------------------------------------------------------------------
'   円(ｼｪｲﾌﾟ)を作成する
'       dCoordX   ：X座標
'       dCoordY   ：Y座標
'       dWidth    ：幅
'       dHeight   ：高さ
'       vRGB()    ：色(R,G,B)
'       戻り値    ：作成した円
'---------------------------------------------------------------------------------------------
Function CreateCircle(dCoordX As Double, dCoordY As Double, _
                      dWidth As Double, dHeight As Double, vRGB() As Long) As Shape

    Dim shpCircle As Shape
    
    '円(ｼｪｲﾌﾟ)作成
    Set shpCircle = ThisWorkbook.ActiveSheet. _
                    Shapes.AddShape(msoShapeOval, dCoordX + (dWidth / 2), dCoordY + (dHeight / 2), dWidth, dHeight)

    '枠線非表示
    shpCircle.Line.Visible = msoFalse
    
    'ｶﾗｰ変更
    shpCircle.Fill.ForeColor.RGB = RGB(vRGB(0), vRGB(1), vRGB(2))
    
    Set CreateCircle = shpCircle

End Function