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

いま一辺の長さが2の正方形に内接するように半径1の円があったとします。 この正方形の上に上から均等に砂粒を振りまきます。 ある程度たくさんまいたところで円の内側に 入った砂粒の数と振りまいた全部の砂粒の数を比べれば円の面積と正方形の 面積のだいたいの比がわかるという寸法です。 円の面積の公式S=pi*r^2はわかっているものとすれば円周率は求まります。

まず、本当は必要ないのですが、見てわかりやすくするために右の図を 用意しましょう。正方形に見えなかったり真円に見えないのはこの表示の 中でx軸とy軸の尺度の取り方が違っているためです。

> x <- seq(-1,1,by=0.001)

> plot(c(-1,1),c(-1,1),type="n")

> par(pch=".")

> abline(v=-1)

> abline(v=1)

> abline(1,0)

> abline(-1,0)

> points(x,sqrt(1-x^2),col="red")

> points(x,-sqrt(1-x^2),col="red")

> xr <- runif(100000)

> yr <- runif(100000)

> xr <- xr*2-1

> yr <- yr*2-1

> length((1:100000)[xr^2+yr^2<=1])

[1] 78524

> 78524*4

[1] 314096

> points(xr,yr,col="yellow")

それでは上から黄色の砂を振りまきましょう。「runif({n})」は [0,1]区間にn個の一様乱数を発生させます。これを2倍して1引くと [-1,1]区間の一様乱数になります。これでx座標、y座標をそれぞれ 決定します。あとは円内に入った砂粒がいくつか数えます。

5行目の説明をします。まず「length({vector})」は中のベクトルの 次元がいくつか(つまり要素数）を数えます。

次に「(1:100000)[xr^2+yr^2<=1]」において大括弧[　]は限定の意味で使います。 ここでは(xr,yr)の組の1〜100000の番号のうちxr^2+yr^2<=1を満たすものです。 ここで「<=」は以下という意味です。この限定の使い方は他にも例えば 「xr[5]」と書いてxrベクトルの5番目の要素と指定することができます。

そこで円内に入る砂粒の数がカウントできたので、円と正方形の面積の比を 求めるとここでは78524:100000となりました。正方形の面積は2*2=4なので 円の面積はpi*1*1=3.14096とわかり、円周率が求まります。






いま求めた方法は確率的なシミュレーションなので円周率といっても 円周率に近いある確率変数の実現値にすぎません。 このようなシミュレーションでこの確率変数がどんな感じに 分布するのか調べてみましょう。今度は実験を10000回繰り返します。 砂粒の数は多ければ多いほどいいし実験回数も多ければ多いほどいいのですが 計算時間が長くかかりすぎるのでこのくらいにします。 これでもPentium2で数十分はかかるでしょう。10000回の実験で 各回に求めた円周率をベクトルzに書き込みます。

> z <- 1:10000

これはただ単に次数10000のベクトルを用意しただけです。

>for (i in 1:10000){

+ xr <- runif(100000)

+ yr <- runif(100000)

+ xr <- xr*2-1

+ yr <- yr*2-1

+ z[i] <- length((1:100000)[xr^2+yr^2<=1])*4

+ }

実験を10000回繰り返しています。長い命令はこのように書きかけでEnterを押すと Rは+記号によりその続きを要求してきますので、複数行にわたって書くことができます。

> hist(z)

ここでヒストグラムを表示します。

> hist(z,nclass=20)

ヒストグラムの階級数を増やしてみました。

今回モンテカルロ法を用いたのは単にモンテカルロ法がどんな物かの 紹介をするためだけで、本来円周率を求めるためには適しません。 この方法は計算回数を多くしてもそれほど精度がよくなるわけではなく、 円周率の計算ならもっと少ない計算量でもっとよい精度を出すことができます。 それではモンテカルロ法の何がいいかというとそのような方法が確立していない 複雑で困難な問題にも容易に適用できてある種の見通しを得ることが できることです。

ホームヘ