モンテカルロ法

TL;DR

「試行を十分な回数行えば、理論上の確率と同じものが求まるのではないか？」と予想し、実際に試してみると大体同じ数値が出てきます。

概要

表が$\frac{1}{2}$、裏が$\frac{1}{2}$の確率で出るコインがあるとします。このコインを1万回投げたとき、表が5000回、裏が5000回出るはずです。

逆に言えば、とあるコインを1万回投げて、表が5000回、裏が5000回出たとき、「表が$\frac{1}{2}$、裏が$\frac{1}{2}$の確率で出るコインなのではないか？」と推測できます。これがモンテカルロ法です。

例:円周率の計算

1辺の長さが$r$の正方形と、その内接円(半径$\frac{r}{2}$)を考えます。面積は、正方形が$r^2$、内接円が$\frac{r^2}{4}\pi$です。

ここで「正方形の面積に占める内接円の割合」を計算すると、
$$ \frac{\frac{r^2}{4}\pi}{r^2}=\frac{\pi}{4} $$

つまり、$r$に関係なく、この割合を4倍すれば円周率が求まるのです。

具体的な数字を出すために、次の方法を使います。

1. 正方形上のランダムな座標を100個指定する(この記事では「点を打つ」と表現)
2. 内接円の中に入っている点の割合を出す

「内接円の中に入っている点の割合」と「正方形の面積に占める内接円の割合」は同じはずなので、この割合を4倍すればいいのです。

計算手順

前述の通り、結果は$r$に依存しません。任意に決めていいのですが、この記事では計算を簡単にするため$r=2$とします。

具体的な手順は以下の通りです。

1. 0~2の乱数を2個発生させる
2. 片方を$x$座標、もう片方を$y$座標とみなす
3. $(x,y)$が内接円の中に入っているかどうか確認する
   • 円の中心は(1,1)、半径は1なので、$(x,y)$が内接円の中にある場合は$(x-1)^2+(y-1)^2 \leq 1$が成り立つ
4. これを100回繰り返し、内接円の中に入った割合を求めて4倍する

これを計算してみます。

表計算

ここでは Google Spreadsheet を使用しますが、Excel等でも大きくは変わりません。

A1セルに「**=rand()*2**」と入力することで02の乱数を発生させることが出来ます(randは01の乱数を発生させる関数で、これを2倍することで0~2の乱数にします)。これをA1:B100にコピーして100個の点を打ったとみなします。

C1セルに「**=(A1-1)^2+(B1-1)^2**」と入力して、計算式の左辺を算出します。これをC1:C100にオートフィルでコピーします。

D1セルに「**=if(C1<=1,1,0)**」と入力して、不等式が成立するなら1を、成立しない場合は0を出力させます。これをD1:D100にオートフィルでコピーします。

F1セルに「**=sum(D1:D100)/100**」と入力し、点が内接円の中に入っている割合を求めます。

F2セルに「**=F1*4**」と入力して4倍します。何回か更新しているうちに、およそ3.14だと知らなくても「円周率って大体これくらいなのでは？」と思えてくるはずです。

Python

プログラミング言語Pythonを用いてシミュレーションを行います。

      import random
a=0
for i in range(100):
  x=random.uniform(0,2) #0~2の乱数を生成してxに代入
  y=random.uniform(0,2) #0~2の乱数を生成してyに代入
  if((x-1)**2+(y-1)**2<=1):
    a+=1 #不等式が成立したらaを1増やす
a/100*4
    

R

プログラミング言語Rを用いてシミュレーションを行います。Rは統計処理に特化した言語です。

次のようなプログラムを実行することで、円周率を求めることが出来ます。

      a=0 #不等式が何回成立したかのカウンタ
for(i in 1:100){
  x=runif(1,min=0,max=2) #0~2の乱数を1個生成してxに代入
  y=runif(1,min=0,max=2) #0~2の乱数を1個生成してyに代入
  if((x-1)^2+(y-1)^2<=1){
    a=a+1 #不等式が成立したらaを1増やす
  }
}
a/100*4 #不等式が100回中何回成立したかの割合を4倍して円周率を出す
    

練習問題

各計算方法において、点を500個や1000個、それ以上に増やして下さい。

点が多ければ多いほど、正確な円周率が出る確率を高められます。