R:円周率を力尽くで求めてみる（モンテカルロ法）～各務ゼミ生に送る～

おはこんばんにちわ。チャパティです。

この前りーやに「R何やってるかわからん！！」って言われたので、この記事を書き始めました。

基本的には各務ゼミ生が読むことを前提としていますが、神戸大学の経営学部で経営統計を受けてた人が理解できるくらいで書いていこうと思います。

Rとは

統計をするときに使うプログラミングです。

ビッグデータの解析に欠かせないRとは

もし、詳しく知りたい人がいたら、Twitter(@npokanzi)とかLINEで聞いてくれたらと思います。

それでは、本題に入りましょう。

円周率の計算の前に

普通に力ずくで求めるときは、外接正多角形と内接正多角形の辺の長さから求めるのが王道なんですが、今回はRを使うので、モンテカルロ法で攻めていきます。

モンテカルロ法
偶然現象の経過をシミュレーションする場合に、乱数を用いて数値計算を行い、問題の近似解を得る方法。コンピューターの発達によって広い分野で利用されている。
三省堂 大辞林より

要はn個の乱数を発生させて、試行回数をデータより莫大に増やすことによって、本来の答えに近い答えを求める方法ってことですね。

まだイメージできないかもなので、さいころの例で例えると

理論的にサイコロの出る目が等確率であることを証明するには、少なくとも以下の証明が必要でしょう。
・サイコロの各面に一切の歪みがないこと
・サイコロの重心が出目の彫に応じた位置になっていること
・サイコロを投げる場が水平で、かつ歪みがないこと
などなど、物理的な検証も必要なので面倒ですね。

しかし、実際に数万回サイコロを振ってみて、各目が1/6に近い確率で出ることが判明すれば、サイコロの出目が等確率であることの立派な証明となります。
https://yolo-kiyoshi.com/2017/10/29/post-391/

て、感じです。（めんどくさいからって、端折るな、、まぁゆるして）

ということでモンテカルロ法を使って円周率を求めて行きます。

円周率の計算

円周率の計算は以下の流れになります。

・半径ｒの円を考えると円の面積はπr²

・この円と接する正方形の面積は４r²

・ｒ×ｒの正方形の中に任意の乱数をｎ個生成したとき、その乱数が扇形の内側にある個数をｍとする

・この時
ｎ：ｍ＝４r²：πr²
→π＝４ｍ/ｎ

これをプログラムに書くと

set.seed(422)

iIter<- 10000000
z    <- 0
w    <-NULL
pdf( file = "円周率.pdf", height = 5, width = 7, family = "Times")

for( i in 1:iIter)
{
 x<-runif(1)
 y<-runif(1)
 if(x^2+y^2<1.0)
 {
   z<-z+1
 }
 if ( i %% 100 == 0)
 {
   v <- 4.0 * z / i
   w <- rbind( w, v )
   print( cbind(i,v))
 }
}

par( mar = c( 2, 2, 1.5, 0.5))
plot(w,type = "l",xlab = "",ylab = "", bty = "n")
##dev.copy2eps( file = "円周率.eps", height = 5, width = 7, family = "Times")
print(4.0*z/iIter)
dev.off()

これが10,000,000個乱数を生成したバージョンですね。

結果がこちら

i=100 v=3.2

i=1,000 v=3.108

i=10,000 v=3.1508

i=100,000 v=3.15024

…

i=10,000,000 v=3.142383

読み方（i→試行回数　v＝結果）

これをグラフにするとこうなります。

大体3.14で収束することがわかりますね！

まとめ

Rでモンテカルロ法を用いることで試行回数をめちゃめちゃ増やして、力業で近似値を求めることが出来ました。

ただ、10,000,000回やってもv=3.142383とまだまだ誤差があるので、効率のいい方法とは言えないかなぁと

他にもいろんな円周率の求め方があるので、

気になる方はこちらへ！

もう円周率で悩まない！πの求め方10選 - プロクラシスト

各務ゼミ生へ

今回の円周率を求めるコードを一行ずつ何してるかを書いたものが下の記事になります。