最近大学の授業で登場した、モンテカルロ・シミュレーションの当たり外れ法を用いて円周率を近似してみようと思います。Goで実装していますが、選んだ理由は最近ハマってるからってだけです。

##モンテカルロ法とは？
まず、モンテカルロ法ってなんだって話ですが、問題を確率を利用したシミュレーションで解く方法です。
今回は、正方形にランダムに点を落とし、内接する円の内側に落ちた点を数えることで円周率を近似する当たり外れ法を実装します。

###円周率を近似する際の考え方
簡単に言うと、半径1の円の面積は 1 × 1 × π = π になるので、円の面積がわかれば円周率もわかるということです。

そこで、図のように半径1の円を四分割したものを一辺の長さ1の正方形で囲い、その中のランダムな位置に N 個の点を落とします。落とした点のうち、円の内側に落ちた点の個数がそのまま1/4円の面積を近似しているため、 N 個のうちの円内の点の個数の割合を求め、 4 をかけて全円の面積にすることで、円周率を近似できます。

この図の例では、N = 10 とし、円内に7個の点が落ちているので、7 / 10 = 0.7 が1/4円の面積、 0.7 × 4 = 2.8 が全円の面積、すなわち円周率の近似となります。

お察しのとおり、この例ではNの値が小さすぎて皆さんご存知の π = 3.141592... とはかけ離れたものになってしまっているので、実際にプログラムを作成してシミュレーションする際にはもっと大きな値を指定します。

##Goによる実装
まず乱数についてですが、 math/rand パッケージの疑似乱数を使用します。rand パッケージの Intn 関数は、引数に指定した自然数 n に対し、 0 以上 n 未満の範囲の整数を rand.Seed 関数で指定されたシード値に基づいて一様に返します。

以下の関数ではシード値に 1 を指定し、 Intn 関数に 16777217 を渡すことで 0 以上 16777216 以下のいずれかの整数を取得します。取得した値を 16777216 で割ることで 0 以上 1 以下の小数を取得しています。
また、この 16777216 という数字ですが、除算において誤差の生じない 2 のべき乗の数となっています。コメントでご指摘いただきありがとうございました。

montecarlo.go

/*モンテカルロ法で円周率を求める関数*/
func CalcPi(n int) (int, float64) {
	var x, y, d float64
	count := 0											//円内に落ちた点を数える
	rand.Seed(1)										//シード値を設定

	for i := 0; i < n; i++ {
		x = float64(rand.Intn(16777217)) / 16777216		//x, y座標を乱数で決定
		y = float64(rand.Intn(16777217)) / 16777216
		d = x * x + y * y								//原点からの距離の2乗を計算

		if d <= 1.0 {count++}							//円内ならcount++
	}
	return count, 4.0 * float64(count) / float64(n)		//円周率を計算しreturn
}

二度乱数を取得、それぞれ x 座標、 y 座標とし、変数 d に原点(=円の中心)からの距離の2乗を計算して代入しています。 d の値が 1 以下ならば x, y ともに 1 以下、つまり、円内の点ということになるのでカウントを増やし、 N 回それを繰り返します。

最後に、円内に落ちた点の数を N で割った数に 4 をかけた値(=求めた円周率)と、円内に落ちた点の数を返しています。

###動作結果
実際にこの関数を使ってシミュレーションしてみました。

montecarlo.go

func main()  {
	n := 1000000
	start := time.Now()
	count, pi := CalcPi(n)
	end := time.Now()

	fmt.Println("\nPI = ", pi)
	fmt.Println(count, "points out of", n)
	fmt.Println("\nexecution time : ", end.Sub(start).Seconds(), "s")
}

実行結果は以下のようになりました。

N=1,000,000での実行結果.

PI =  3.141656
785414 points out of 1000000

execution time :  0.0560016 s

N = 1,000,000 で実行した結果、小数点以下 4 桁まで正しい値が出ましたね。
ちなみに N = 100,000,000 で実行すると、

N=100,000,000での実行結果.

PI =  3.14148656
78537164 points out of 100000000

execution time :  7.1651347 s

このようになり、大して変わりませんね。この方法の限界というものの存在がなんとなくわかります。また、実際にシードを変えて複数の乱数列でシミュレーションを行うと、意外とバラツキが大きいです。今回は試行回数を重ねて結果を集計などはしませんが、求めるものによっては信用に足る結果を出すために N の値を相当大きくしなければならないのだと思います。

しかし、小学生の頃から慣れ親しんでいる定数をこのような形で求められるのはなかなかおもしろいです。良かったら読んでいる方もプログラムを作成して試してみてはいかがでしょうか？

今回はシード値を 1 に固定して実行しましたが、毎回別の乱数列を使用したい場合は、 rand.Seed(time.Now().UnixNano()) のように現在時間をシード値に用いれば実行ごとに結果を変えることもできます。