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

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

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

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

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

 

モンテカルロ法とは

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

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

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

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

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

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

つまり逆に考えれば、要素を探索する順番もランダムにすれば、場合によっては最速で要素 \(a\) を取得することができる、というのが乱択アルゴリズムの考え方です。要は「運が良ければ最速の処理にもなるし、運が悪ければ遅くもなる!」というようなギャンブル的処理方法であるともいえます。
  

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

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

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

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

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

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

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

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

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

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

 
ちなみに名称の”モンテカルロ”についてですが、カジノで有名なモナコ公国の地区の名前が由来とされています。実はこのモンテカルロ法はギャンブラーの間では有名な手法で、かつてこの手法をつかったベット(賭け)方法でカジノを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\)」(の近似値)を求めることができるという訳です。

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

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

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

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

 

サンプルコード

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

 

まとめ

今回は乱択アルゴリズムの「モンテカルロ法」の解説でした。

特徴① アルゴリズム内で乱数を発生させてランダム性を与える
特徴② 主にある事象に対しての確率を求めることに使われるアルゴリズム
特徴③ 試行回数を増やせば増やすほど求めたい確率に収束していく
特徴④ この手法で導かれた答えは100%正しいとは言えない

主に確率を求めるために使われるアルゴリズムですが、円周率を求めたときのように考え方を工夫することで確率以外のものも調べることが可能です。たとえば素数判定のアルゴリズム「ミラーラビン素数判定法」でも同じ考えが利用されています。

ちなみにモンテカルロ法は特徴④の通りで、場合によっては失敗する可能性のあるアルゴリズムです。対して、処理時間を増やしてでも必ず正解を導き出す「ラスベガス法」というアルゴリズムも存在するのでぜひそちらのアルゴリズムもチェックしてみて下さい。

 
 icon-book オススメ書籍

Excel, VBA, アルゴリズム, 数学

Posted by Lic