モンテカルロ法による円周率の推定（その1）

乱数（random number）についていい例題がないかなぁ，と考えて「円周率をモンテカルロ法で求めるやつやればいいぢゃん」と思いついた。ので早速試してみる。ちなみに「モンテカルロ法」というのは数値計算やシミュレーションに乱数を用いる手法をさす。

モンテカルロ法による円周率の推定

では乱数を使ってどうやって円周率を求めるのか。まずは以下のように原点を中心とした半径 $1$ の円を考える。ただしここでは第一象限のみを対象とする。

そして $0 \le y \le 1$ および $0 \le y \le 1$ の範囲でランダムに点をプロットしていく。（以下の図は Wikimedia Commons のものを拝借した。 CC-BY-3.0 で公開されている）

全ての点 $n$ が領域内に均等にプロットされていれば，円の内側に入る点の数 $m$ は，面積比から，以下の式のようになると期待できる。

\[ m \simeq \frac{1}{4}{\pi}n \]

この式を $\pi$ を求める形に変形すると

\[ \pi \simeq \frac{4m}{n} \]

となる。プロットした点が円の内側かどうかは原点からの距離で判定できる。すなわち

\[ \sqrt{x^2 + y^2} \le 1 \]

を満たしていればよい。

math/rand パッケージ

Go 言語にはコア・パッケージとして math/rand が用意されていて，このパッケージを使って擬似乱数を発生させることができる。今回は $0 \le r \le 1.0$ の範囲で乱数を発生させればいいのだが，生憎そのものズバリな関数が用意されていない。たとえば rand.Float64() が吐く値の範囲は $0 \le r \lt 1.0$ なのでそのままでは使えないのだ。

そこで，こんなコードを考えてみる。

package main

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

func main() {
    rand.Seed(time.Now().UnixNano())
    for i := 0; i < 10; i++ {
        fmt.Println(float64(rand.Int63n(10000001)) / float64(10000000))
    }
}

rand.Int63n(n) 関数は $0 \le i \lt n$ の範囲で整数を吐く。 $n=10,000,001$ なら $0 \le i \le 10,000,000$ の範囲になる。これを $10,000,000$ で割って $0 \le r \le 1.0$ の範囲の乱数を作るのである。

実際には2次元座標なので複素数（complex）表現にして

package main

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

func main() {
    rand.Seed(time.Now().UnixNano())
    for i := 0; i < 10; i++ {
        c := complex(float64(rand.Int63n(10000001))/float64(10000000), float64(rand.Int63n(10000001))/float64(10000000))
        fmt.Printf("%v (%v)\n", c, cmplx.Abs(c))
    }
}

とすればよい。ちなみに cmplx.Abs() 関数は複素数の絶対値を取るもので， $\sqrt{x^2 + y^2}$ と同じである。

では，以上を踏まえてランダムな点を生成する gencmplx パッケージを作ってみよう。 channel と goroutine を使ってこんな感じにするかな。

package gencmplx

import "math/rand"

//New returns generator of random complex number
func New(s rand.Source, count int64) <-chan complex128 {
    ch := make(chan complex128)
    r := rand.New(s)
    go func(r *rand.Rand, count int64) {
        for i := int64(0); i < count; i++ {
            ch <- complex(float64(r.Int63n(10000001))/float64(10000000), float64(r.Int63n(10000001))/float64(10000000))
        }
        close(ch)
    }(r, count)
    return ch
}

後々のことを考えて，乱数の rand.Source¹ と生成する点の個数は引数で指定するようにした。

いっぽう， gencmplx パッケージの呼び出し側はこんな感じになる。

package main

import (
    "fmt"
    "math/rand"
    "time"

    "github.com/spiegel-im-spiegel/pi/gencmplx"
)

func main() {
    c := gencmplx.New(rand.NewSource(time.Now().UnixNano()), int64(10000))
    for p := range c {
        fmt.Printf("%v\t%v\n", real(p), imag(p))
    }
}

ここでは少なめに1万個の点を生成している。 channel c からの値の取り出しは for-range 構文を使うと記述がシンプルになり c が close() するまでループしてくれる。

早速これを動かしてみよう。

$ go run main.go > plot-1.dat

これで1万個の点を plot-1.dat に格納したことになる。 plot-1.dat を gnuplot に食わせてみるとこんな感じ。

うーん。均等？なのかなぁ。まぁ，この辺の評価については「その3」以降で。

最後に，生成した点の集合から円周率を推定するところまでやってみよう。 main() 関数はこのように変える。

package main

import (
    "fmt"
    "math/cmplx"
    "math/rand"
    "time"

    "github.com/spiegel-im-spiegel/pi/gencmplx"
)

func main() {
    c := gencmplx.New(rand.NewSource(time.Now().UnixNano()), int64(100000))
    n := int64(0) // total
    m := int64(0) // plot in circle
    for p := range c {
        n++
        if cmplx.Abs(p) <= float64(1) {
            m++
        }
    }
    fmt.Printf("n = %v, m = %v, 4m/n = %v\n", n, m, float64(4*m)/float64(n))
}

点の数は10万個まで増やしている。実行結果は

$ go run main2.go
n = 100000, m = 78397, 4m/n = 3.13588

と，まぁそれっぽい値が出てきた。

今回はここまで。

ブックマーク

Go 言語に関するブックマーク集はこちら。

参考図書

プログラミング言語Go: アラン・ドノバン (著), ブライアン・カーニハン (著), 柴田芳樹 (著); 丸善出版 2016-06-20 (Release 2021-07-13); Kindle版; B099928SJD (ASIN); 評価

Kindle 版出た！一部内容が古びてしまったが，この本は Go 言語の教科書と言ってもいいだろう。感想はこちら。

reviewed by Spiegel on 2021-05-22 (powered by PA-APIv5)

数学ガール／乱択アルゴリズム: 結城浩 (著); SBクリエイティブ 2011-02-25 (Release 2014-03-12); Kindle版; B00I8AT1FO (ASIN); 評価

工学ガール，リサちゃん登場！

reviewed by Spiegel on 2015-04-19 (powered by PA-APIv5)

rand.Source は interface として定義されていて，これを満たす type であれば他の擬似乱数アルゴリズムも使えるようになっている。 ↩︎