プログラマーは今こそアルゴリズムを書くべき!!~モンテカルロ法でπを計算してみよう(コードあり)~
常日頃思っていることとしてweb系だろうと業務系だろうとプログラマーたるものアルゴリズムを書く必要はあると思うんですね。
理由はいっぱいあるけど、数学の素養を持つということはプログラマーとしてかなり重要なことです。
というわけで今回はアルゴリズムについて。
特にモンテカルロ法について説明しようと思います。
このような流れで説明します。
アルゴリズムって?
アルゴリズムを自分の言葉で簡潔に説明すると、
「計算の手順で解があるもの」
です。
単純に「計算の手順」ならば『手続き』ですね。
他にも「計算の結果が有限時間内で出せる手続き」などいろんな言い方があると思うけど。
今回はアルゴリズムの一例としてモンテカルロ法を説明します。
モンテカルロ法
モンテカルロっていうギャンブルが盛んな街があるんですね。それから取った手法で、
ランダムな物を利用していろんな計算をしてみようという感じのもの。
定義をごちゃごちゃいってもあれなので具体的を説明します。
ここで質問ですが、みなさんはπ=3.14を出すためにどのような計算をしますか?
正多角形から計算する方法や三角関数をいじって出す方法があると思うのですが、モンテカルロ法を使うと簡単に出せます。早速説明していきます。
モンテカルロ法でπを計算する
まず次の図を考えます。
これは正方形と1/4円です。1/4円の中心は正方形の左下の角と一致します。
ここで定義として
- 縦、横ともに 1
- 1/4円の半径は 1
(単位は意味を持たないのでここでは気にしません)
とすると
- 正方形の面積は 1
- 1/4円の面積は π / 4
となります。
面積比は
1 : π / 4
ということになります。
次にこの正方形からはみ出さないようにランダムに点を打っていきます。
本当にランダムに。すると次の図のようになります(なんだか微妙にグロ…)。
すると1/4円の内部にある点と正方形の内部だけど、1/4円の外部にある点にわかれます(色分けしてるのでわかると思いますが)。
ここで点の数を数えて、
- 内部にある点の個数を InP (INにあるPoint)とする
- 外部にある点個数を OutP (OutにあるPoint)とする
で、ランダムに点がうたれているということは自然と比率が面積と一致するはず。
なので次のような式が得られます。
1 : π / 4 = (OutP + InP) : InP
これを整理すると次のようになります。
π = InP / ((OutP + InP) / 4)
というわけでπイコールの式になります。つまり、未知数である InP と OutP さえわかればπがわかるということ。
実装(ruby)
というわけで以上のアルゴリズムを実装すると次のようになります。
# 1 def square(n) return n*n end # 2 def monte(n) count = 0 for i in 1..n do if Math::sqrt(square(rand) + square(rand)) < 1 # 3 count += 1 end end return count end def pi(n) return monte(n).to_f/(n/4) # 4 end puts pi(1000) # 5
上から順位解説していきます。
#1. これはnの2乗を計算する関数です。あとで使います。
#2. これは内部の点の個数を返す関数です。つまり上記の式の InPを返します。
#3. ここが一番のポイントですが、n 回ループを回し、
Math::sqrt(square(rand) + square(rand))
が1より小さければカウントを大きければカウントしないという処理をしています。
sqrtの内部にrandが2つありますが、
- 最初のrandは x
- 次の randは y
と考えるとわかりやすいと思います。1 / 4の円の半径は1ということから原点と点の距離が1未満ならば内部の点に、1以上ならば外部の点になります。
またこの時
n = InP + OutP
となります。
#4. InP / ((InP + OutP) / 4)を返す。
つまりπを返す関数です。
monte(n).to_f
となっているのは結果を実数型で返すためです。
#5. ここで n を設定します。
以上がコードの簡単な説明でした。
いよいよ結果発表です。
結果
次のような結果になりました (n = 1000)
> 3.188
πは3.14なので少し違いますね。
次に n = 2000としましょう。すると、
> 3.152
少し近づきました
- n = 5000
> 3.1288
少し離れた笑
- n = 10000
> 3.1352
また近づいた!
- n = 50000
> 3.14032
今までで最高。ただし時間が…(0.5秒くらい遅い)
。。。。。
というわけで n の数をどんどん大きくするとより精度が増します。自明すぎるので特に解説はしません。
そんな感じで試しに上記のコードを写経して(コピペじゃなく)実際に試してみると理解がより深まるかと思います。
モンテカルロ法はアルゴリズムの超入門的な感じなのでファーストステップには是非!!
また、
「より綺麗なコーディングの仕方あるぜ!」
とか
「こんなコーディングも考えたのだが」
とかの意見があれば嬉しいです。
それではまたいつかアルゴリズムについて書いていこうと思います。
ではでは
