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

推定結果の分布

さて，前回書いたコードを利用して，いよいよ円周率の推定結果を評価してみる。

spiegel-im-spiegel/pi: Estimate of Pi with Monte Carlo method.

CLI はこんな感じになった。

$ go run main.go estmt --help
Estimate of Pi with Monte Carlo method.

Usage:
  pi estmt [flags]

Flags:
  -e, --ecount int   Count of estimate (default 100)
  -p, --pcount int   Count of points (default 10000)

Global Flags:
      --config string   config file (default is $HOME/.pi.yaml)

まずは円周率の推定処理を10,000回繰り返してみる。また推定処理のためのランダムな点の数 $n$ を1,000個，10,000個，100,000個と変えて実行してみようか。

$ go run main.go estmt -e 10000 -p 1000 > estmt1k.dat
minimum value: 2.94800
maximum value: 3.31600
average value: 3.14199
standard deviation: 0.05104 (69.9%)

$ go run main.go estmt -e 10000 -p 10000 > estmt10k.dat
minimum value: 3.07240
maximum value: 3.20360
average value: 3.14178
standard deviation: 0.01654 (68.0%)

$ go run main.go estmt -e 10000 -p 100000 > estmt100k.dat
minimum value: 3.12360
maximum value: 3.16184
average value: 3.14163
standard deviation: 0.00518 (68.3%)

最後のはさすがに時間がかかるので，お茶でも飲みながら優雅に待ちましょう（笑）

標準エラー出力に最小値，最大値，平均値（$E$），標準偏差（$\sigma$）を出力してみた。標準偏差の後ろの括弧は $\left[ E-\sigma, E+\sigma \right]$ の範囲にある推定値の割合を示したものだ。

円周率の推定処理の試行回数が十分大きいなら推定値の分布は正規分布（またはガウス分布）に近似できる筈である。（以下の図は Wikimedia Commons のものを拝借した。 CC-BY-2.5 Generic で公開されている）

そこで $n=100,000$ のときの推定結果についてヒストグラムを描いてみる。幸いなことに gnuplot では簡単にヒストグラムを作図できる。こんな感じ（階級幅を0.001としている）。

gnuplot> filter(x,y)=int(x/y)*y
gnuplot> plot "estmt100k.dat" u (filter($1,0.001)):(1) smooth frequency with boxes lc rgb "black"

gnuplot の出力結果はこんな感じ。

んー。まぁ正規分布っぽい？

もうひとつ，正規確率の分布を調べてみよう。これも gnuplot で描こうと思ったけど，少し面倒そうなので，ズルして以下を参考に LibreOffice の Calc で描くことにした¹。

正規確率プロット（Normal Q-Q Plot）の作成 with Excel

とりあえず結果だけ。

プロットが直線状に並んでいれば正規分布であると言える。図から見る限り，概ね正規分布になっているようである。

しまった。ここまで Go 言語が全然出てこなかった。まぁ，いいや。多分あと1回続きます。

おまけ：誤差評価

モンテカルロ法を使ってどの程度の精度で円周率が求まるかの考察については以下が参考になる。

モンテカルロ法の誤差を考える

これも横着して結果だけを拝借すると， $n=100,000$ で推定を行った場合の値の分布は，中央値を $\pi$，99%信頼区間を $\frac{4.230}{\sqrt{100,000}} = 0.013$ として， $\left[ \pi - 0.013, \pi + 0.013 \right]$ の範囲になるようだ。

おまけの追記：正規確率の分布図について

正規確率の分布図（Q-Q プロット）を描くのに毎回 Excel や Calc を使うのもどうかという気がしたので，こちらのプログラム側であらかじめ計算して，結果のプロット・データを gnuplot に食わせるよう考えてみる。

まず qq サブコマンドを追加し，この qq サブコマンド時にデータファイルを読み込んで Q-Q プロットの計算を行うように CLI を変更する。（ついでに他のオプションも整理した）

$ go run main.go qq --help
make Q-Q plot data.

Usage:
  pi qq [data file] [flags]

Global Flags:
  -p, --pcount int   Count of points (default 10000)

実際の処理部分はこんな感じ。

//Execute output Q-Q plot data.
func Execute(cxt *Context) error {
    scanner := bufio.NewScanner(cxt.ui.Reader())
    pis := make([]float64, 0)
    for scanner.Scan() {
        pi, err := strconv.ParseFloat(scanner.Text(), 64)
        if err != nil {
            return errors.Wrap(err, "invalid data")
        }
        pis = append(pis, pi)
    }
    ecf := float64(len(pis))

    sort.Float64s(pis)
    for i, pi := range pis {
        rank := (float64(i+1) - 0.5) / ecf
        cxt.ui.Outputln(fmt.Sprintf("%v\t%v", qnorm(rank), pi))
    }

    return nil
}

qnorm() 関数は標準正規分布に対する累積分布関数の逆関数の値を返すのだが（Excel の NORM.S.INV 関数相当），Go 言語で書かれた適当なパッケージが見つからなかったので（もしあれば誰か教えて）以下のページのコードを Go 言語用に書き直した。

RangeVoting.org - Java normal distribution calculator II

実際にはこんな感じ²。

//qnorm function refers to http://rangevoting.org/Qnorm.html
// This function is licensed under GNU GPL version 3 or later.
func qnorm(p float64) float64 {
    const (
        split = 0.42
        a0    = 2.50662823884
        a1    = -18.61500062529
        a2    = 41.39119773534
        a3    = -25.44106049637
        b1    = -8.47351093090
        b2    = 23.08336743743
        b3    = -21.06224101826
        b4    = 3.13082909833
        c0    = -2.78718931138
        c1    = -2.29796479134
        c2    = 4.85014127135
        c3    = 2.32121276858
        d1    = 3.54388924762
        d2    = 1.63706781897
    )

    q := p - 0.5
    ppnd := float64(0)
    if math.Abs(q) <= split {
        r := q * q
        ppnd = q * (((a3*r+a2)*r+a1)*r + a0) / ((((b4*r+b3)*r+b2)*r+b1)*r + 1)
    } else {
        r := p
        if q > 0 {
            r = 1 - p
        }
        if r > 0 {
            r = math.Sqrt(-math.Log(r))
            ppnd = (((c3*r+c2)*r+c1)*r + c0) / ((d2*r+d1)*r + 1)
            if q < 0 {
                ppnd = -ppnd
            }
        }
    }
    return ppnd
}

では早速動かしてみよう。

$ go run main.go qq estmt100k.dat > qq100k.dat

生成した qq100k.dat ファイルを gnuplot に食わせる。こんな感じでいいだろう。

gnuplot> unset key
gnuplot> set style line 1 pointsize 0.1 pointtype 7 linecolor rgb "black"
gnuplot> plot "qq100k.dat" linestyle 1

結果はこんな感じ。

ついでにこの分布図にフィットする直線 $y=ax+b$ の $a, b$ 値を調べてみる。

gnuplot> f(x)=a*x+b
gnuplot> fit f(x) "qq100k.dat" via a,b
iter      chisq       delta/lim  lambda   a             b            
   0 5.5761271099e+04  0.00e+00 1.00e+00  1.000000e+00  1.000000e+00
   1 6.0252530865e-04 -9.25e+12 1.00e-01  5.277253e-03  3.141418e+00
   2 4.5071534394e-05 -1.24e+06 1.00e-02  5.177775e-03  3.141632e+00
   3 4.5071534393e-05 -1.23e-06 1.00e-03  5.177775e-03  3.141632e+00
iter      chisq       delta/lim  lambda   a             b            

After 3 iterations the fit converged.
final sum of squares of residuals : 4.50715e-005
rel. change during last iteration : -1.23364e-011

degrees of freedom    (FIT_NDF)                        : 9998
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 6.71421e-005
variance of residuals (reduced chisquare) = WSSR/ndf   : 4.50806e-009

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 0.00517777       +/- 6.715e-007   (0.01297%)
b               = 3.14163          +/- 6.714e-007   (2.137e-005%)

correlation matrix of the fit parameters:
                a      b      
a               1.000 
b               0.000  1.000

ここで $a$ が標準偏差， $b$ が平均値にマッチしている点に注目。分布図と上の直線を重ねあわせるとこうなる。

んー。こんなもんかな。

そうそう。 qq サブコマンドは，フィルタとしても機能するので

$ go run main.go estmt -e 100 -p 10000 | go run main.go qq > qq.dat
minimum value: 3.09640
maximum value: 3.18600
average value: 3.14158
standard deviation: 0.01654 (68.0%)

といった感じにパイプでつなぐこともできる。

ブックマーク

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

参考図書

プログラミング言語Go (ADDISON-WESLEY PROFESSIONAL COMPUTING SERIES): Alan A.A. Donovan (著), Brian W. Kernighan (著), 柴田芳樹 (翻訳); 丸善出版 2016-06-20; 単行本（ソフトカバー）; 4621300253 (ASIN), 9784621300251 (EAN), 4621300253 (ISBN), 9784621300251 (ISBN); 評価

著者のひとりは（あの「バイブル」とも呼ばれる）通称 “K&R” の K のほうである。この本は Go 言語の教科書と言ってもいいだろう。

reviewed by Spiegel on 2016-07-13 (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)

数学ガールの秘密ノート／やさしい統計: 結城浩 (著); SBクリエイティブ 2016-10-28 (Release 2016-11-10); Kindle版; B01MSJMKMW (ASIN); 評価

統計の本当に基礎の部分から。学業成績でよく聞く「偏差値」とは何を表していて何を意味しているのか。なんてなあたりから。

reviewed by Spiegel on 2016-12-11 (powered by PA-APIv5)

Calc でも Excel の関数がそのまま使えるようだ。助かる。 ↩︎
余談だが， Go 言語では untyped な定数を設定できる。型が評価されるのは，処理の中でその定数が使われた時点となる。数値の精度も使用時点で評価されるため，定義では大きい桁数の値を設定しても問題ない。（参考： Go の定数の話 - Qiita） ↩︎