モンテカルロ法で円周率を求める(Go)

A Tour of Goを履修しました。練習として、モンテカルロ法で円周率を求めるシミュレーションをGoで実装したので、それをご紹介します。

モンテカルロ法による円周率の導出
シミュレーション
おわりに
参考文献

モンテカルロ法による円周率の導出

モンテカルロ法というと難しい印象を受けるかもしれませんが、やりたいことは「乱数を使って点をばらまき、円の中にある点を数える」だけです。

いったん、円の面積を求める公式を振り返ります。

$\displaystyle S = \pi r^2$

$S$ が円の面積、 $r$ が円の半径、 $\pi$ が円周率です。小学校では「円の面積=半径×半径×3.14」と習いますね。¹

ここでは後述するシミュレーションを簡単にするため、円全体ではなく円の4分の一（四分円）を考えます。

上図において、 $S_{1}$ を正方形の面積、 $S_{2}$ を四分円の面積とすると以下の比が成り立ちます。

$\displaystyle S_{1}:S_{2} = r^2 : \frac{\pi r^2}{4}$

これを次のように書き下します。

$\displaystyle \begin{align} \frac{S_{2}}{S_{1}} &= \frac{\frac{\pi r^2}{4}}{r^2} \\ \frac{S_{2}}{S_{1}} &= \frac{\pi}{4} \\ \end{align}$

$\pi$ について求めると次式になります。

$\displaystyle \pi = 4 \frac{S_{2}}{S_{1}}$

ところで、面積は点群の点の数とみなせます。

つまり上式から、 $S_{1}$ と $S_{2}$ それぞれの点の数を求めるだけで、円周率 $\pi$ が計算できることがわかります。

シミュレーション

今回はGoでシミュレーションを実装しました。

開発環境

シミュレーションプログラムを実行する開発環境は次の通りです。Goの標準パッケージだけで動きます。

項目	内容
OS	Windows 10
Go	1.20.3

サンプルプログラム

サンプルプログラム(pi.go)は以下の通りです。自由にコピペしていただいて構いません。ぜひ実際に動かして見てください。

package main

import (
    "fmt"
    "math"
    "math/rand"
)

func main() {
    var seed int64 = 23  // シード数
    var n1 int = 1000000 // 全ての点の数
    var n2 int           // 四分円の内部にある点の数
    var pi float64       // 円周率←これを計算する

    // シードを設定
    rand.Seed(seed)

    for i := 0; i < n1; i++ {
        // xy座標系に[0.0, 1.0)の区間で、点をランダムに配置
        x, y := rand.Float64(), rand.Float64()

        if math.Sqrt(x*x+y*y) < 1.0 {
            //四分円の内部にある場合はインクリメント
            n2++
        }
    }

    // 点の数の比から円周率を計算
    pi = float64(n2) / float64(n1) * 4.0
    fmt.Printf("pi is %v\n", pi)
}

実行結果は以下の通りです。

$ go run pi.go
pi is 3.141752

円周率に大体近い値が計算できました。

解説

サンプルプログラムは、「一辺の長さが1の正方形に点を1000000コばらまき、四分円の内部にある点を数えて、最後にその点の数から円周率を計算する」というアルゴリズムです。非常に簡単で行数も少ないので、1行ずつ読み解いてみてください。

なお、乱数のシードは設定しなくても計算できます。しかし、seedを設定しておくと実行するたびに値が変わってしまうことを防げます。全体の点数n1と合わせて、いろいろな値で実行してみてください。

おわりに

有名なシミュレーションである「モンテカルロ法による円周率の計算」を、これまでやったことがなかったので、この機会にできて満足です。

開発環境や言語仕様含めて、Goはサクッと使えるので結構気に入っています。メインはPythonですが、これからはGoも使っていきたいと思います。

参考になれば幸いです(^^)

参考文献

もしかして最近の小学校では、円周率は3.14でなく3かも？↩

ari23の研究ノート

メーカ勤務エンジニアの技術ブログです