帰ってきた「しっぽのさきっちょ」

モンテカルロ法による円周率の推定(その3 Gaussian) — プログラミング言語 Go

  1. モンテカルロ法による円周率の推定(その1)
  2. モンテカルロ法による円周率の推定(その2 CLI)
  3. モンテカルロ法による円周率の推定(その3 Gaussian) ← イマココ
  4. モンテカルロ法による円周率の推定(その4 PRNG)

推定結果の分布

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

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 で公開されている)

normal distribution from Wikimedia
normal distribution from Wikimedia

そこで $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 で描くことにした1

とりあえず結果だけ。

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

しまった。 ここまで 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 言語用に書き直した。

実際にはこんな感じ2

//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 言語に関するブックマーク集はこちら

参考図書

photo
プログラミング言語Go (ADDISON-WESLEY PROFESSIONAL COMPUTING SERIES)
Alan A.A. Donovan Brian W. Kernighan 柴田 芳樹
丸善出版 2016-06-20
評価

スターティングGo言語 (CodeZine BOOKS) Go言語によるWebアプリケーション開発 Kotlinスタートブック -新しいAndroidプログラミング Docker実戦活用ガイド グッド・マス ギークのための数・論理・計算機科学

著者のひとりは(あの「バイブル」とも呼ばれる)通称 “K&R” の K のほうである。

reviewed by Spiegel on 2016-07-13 (powered by G-Tools)

photo
数学ガール/乱択アルゴリズム
結城 浩
SBクリエイティブ株式会社 2014-02-14
評価

数学ガール/ゲーデルの不完全性定理 数学ガール/フェルマーの最終定理 数学ガール/ガロア理論 数学ガール 数学ガールの秘密ノート/式とグラフ

工学ガール,リサちゃん登場!

reviewed by Spiegel on 2015/04/19 (powered by G-Tools)

photo
数学ガールの秘密ノート/やさしい統計
結城 浩
SBクリエイティブ 2016-10-28
評価

数学ガールの秘密ノート/場合の数 数学ガールの秘密ノート/ベクトルの真実 はじめての深層学習(ディープラーニング)プログラミング ポートとソケットがわかればインターネットがわかる――TCP/IP・ネットワーク技術を学びたいあなたのために 数学ガールの秘密ノート/数列の広場 数学ガールの秘密ノート/微分を追いかけて 数学ガールの秘密ノート/整数で遊ぼう 数学ガールの秘密ノート/丸い三角関数 経済数学の直観的方法 確率・統計編 (ブルーバックス) 速習 Python 3 上: プログラミングの基礎編

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

reviewed by Spiegel on 2016-12-11 (powered by G-Tools)


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