黄昏より暗きもの、血の流れより赤きもの

読者です 読者をやめる 読者になる 読者になる

黄昏より暗きもの、血の流れより赤きもの

自分の好きな事を好きなように書いて行きます。

【Python(matplotlib)】モンテカルロ法を使い、円周率の近似値を数値計算してみた(1.9 追記あり)

始めに

今日は筆者の勉強の為にモンテカルロ法と言うものを扱う。まず今日の筆者の目標を以下にまとめた。

確率統計:平均、標準偏差、信頼区間等を使い、実験においてどれくらいの精度の情報が欲しいかの目処立てができるようにする
Python(matplotlib):基本的なグラフの描画を覚える

統計学を上手く使うと、平均や標準偏差などの情報から実験結果の数値がどれだけの範囲に収まるかという目処が立てやすくなる。簡単に言えば実験における判断材料が増え、よりていねいに実験を進めたり、あるいは実験の計画を立てる事ができるようになるのだ。

モンテカルロ法を用いて円周率の近似値を得る流れ

さて以下の手順より、モンテカルロ法を用いて円周率の近似値を得る事を考える。この手法では円の情報などを座標で処理する。そこで後々プログラムに落とし込みやすくするために、「モンテカルロ法と円周率の近似計算|高校数学の美しい物語」に書かれている内容を座標の形に書き直してみた。(数式アレルギーの方ごめんなさい)


f:id:program_study:20160106195413p:plain:w300

xy平面上の正方形の領域をS=\{(x,y) | 0 \leq x \leq 1かつ 0 \leq y \leq 1\}
四分円の領域をT=\{(x,y) | 0 \leq x \leq 1 かつ 0 \leq y \leq1 かつ x^2 + y^2 \leq 1\}とする。このとき以下の操作を行なう。
[STEP1]:四分円の的Tめかけて、ダーツをn回投げる。
[STEP2]:n回の試行のうち、四分円の的Tの領域(図の赤の領域)に当たった回数をカウントし、この値をxとする。
[STEP3]:四分円Tの領域に当たったかどうかは、k回目の試行においてダーツの当たった位置をP_k(x_k,y_k)とし、原点からの距離が1以下(D_k=\sqrt{{x_k}^2+{y_k}^2} \leq 1)で判定する。(上の図の赤い点)
[STEP4]:\frac{4x}{n}の値を円周率の近似値として採用する(図を見ても分かる通り、4分円なので\frac{x}{n}を4倍する)

ここで\frac{x}{n}とは、n回ダーツを投げ、的Tに入る事象がどれだけ観測されるかの確率である。後述する大数の法則よりnを限りなく大きくすればこの値が\frac{\pi}{4}に近づく性質を利用し、πの近似値を求めていく。

P_k(x_k,y_k)の値に、区間[0,1]における一様乱数[numpy.random.rand()]*1を用いる。これは理論上、0から1の乱数が同じ割合で登場するものである。

ダーツが四分円の的Tに当たる確率

ここでダーツの正方形が四分円の的Tに当たる確率について説明したい。面積とは、ある領域内を、限りなく細かい正方形で敷き詰めたものと解釈される。

領域Tの面積は、Tの部分に含まれる正方形の個数となり、(Tの面積) = 1 * 1 * π * 1/4 = π/4 となる。このエリアにダーツが当たれば良いので、ダーツが四分円の的Tに当たる確率p= \frac{\pi}{4}となる。

円周率の近似値を予測する

統計学を活用し、実験の計画立てを行なう

今度は統計学的な手法により、モンテカルロ法を行なうとどのような円周率の値となるのか?その実験結果の予想立てを行なう。

試行回数は多いに越す事はない(大数の法則)

まず実験に必要な背景として大数の法則というものを紹介する。大数の法則とは、

大数の法則

同じ分布における確率変数X_1,X_2,X_3,……,X_nにおいて、そのn個の確率変数の算術平均\overline{X} = \frac{X_1+X_2+…+X_n}{n}は、試行回数nを限りなく大きくすれば真の平均μに近づく

良く「実験結果は複数回計測しその平均値を取れ」とか「スロットを打てば打つ程、BIGを引く確率が元の値に近づく*2」と言われる。これは大数の法則より多く試行を繰り返せば繰り返す程、平均値が元の平均値に落ち着くと言う事が根拠にあたる。

95%の信頼区間の情報から、試行回数の目処を立てる

まずn回の試行のうち、ダーツが的Tに入る場合とダーツが的Tに入らない場合という成功と失敗からなる試行とみなす事ができる。このことから成功回数Xの分布は二項分布に従うとみなすことができ、公式より平均\overline{X}=np標準偏差 \sqrt{V(X)}=\sqrt{np(1-p)}と計算出来る。

さらに二項分布はnを限りなく大きくすると正規分布に近づくことが分かっている。この事を上手く利用し、円周率の近似値の95%が収まる区間(信頼区間)を考えると、

P(\overline{X} -1.96 * \sqrt{V(X)}  \leq X \leq \overline{X} + 1.96 * \sqrt{V(X)} ) = 0.95
P(\frac{4\overline{X} - 4 * 1.96 * \sqrt{V(X)}}{n} \leq \frac{4X}{n} \leq \frac{4\overline{X} + 4 * 1.96 * \sqrt{V(X)}}{n}) = 0.95

π=3.14と仮定して見積もると、分布の95%が収まる区間は[3.14-\frac{0.805}{\sqrt{n}}  ,3.14+\frac{0.805}{\sqrt{n}}]となる。ここでnに適当な値を入力し、Excelで計算した結果を以下に示す。

  • n=100のとき、円周率の値の95%が収まる区間は[3.05947,3.22052]
  • n=1000のとき、円周率の値の95%が収まる区間は[3.11453,3.16546]
  • n=10000のとき、円周率の値の95%が収まる区間は[3.13194,3.14805]

ここで、95%の確率で円周率の誤差を0.001以内に抑えるためには\frac{0.805}{\sqrt{n}} \leq 0.001n \geq 648,025(回) 以上計測しなければならない。しかしながら使用言語がPython(インタプリタ言語)の関係で、実行時間がかかるのが面倒である。その為試行回数をn=10000までとして実験を行なう。

実験1:n=10000までの円周率の近似値の変化に注目し、大数の法則が成り立つか?

実験1では今回の状況により、大数の法則がどれだけ成り立つかグラフの上でチェックする。まず、n=500回目までの円周率の近似値の変化について調べてみる。このソースコードにて調査をしたところ、以下のような結果が得られた。


f:id:program_study:20160106194211p:plain

100回前後までは数値の振動が激しく、最低でも500回はダーツを投げないと使い物にならない事がわかる。

今度はn=10000回までの円周率の近似値の変化を見てみたい。


f:id:program_study:20160106194234p:plain

さて、円周率の値の95%が収まる区間は[3.13194,3.14805]であった。ともすれば95%の確率でδ=3.14805 - 3.13194 = 0.0161の範囲で数値が振動する事が予想できる。

グラフを見てもn=8,000回以降で数値の振動が少なくなる事がわかり、誤差を抑えて円周率の近似値を得やすい事が分かる。又直感的に大数の法則が成り立っていることも、上の図より読み取れる。

実験2:n=10,000までの円周率の近似値の変化を調べる

実験2ではダーツを投げる試行回数をn=100,1000,10,000の場合の計測を3回行ない、それぞれの場合における円周率の近似値をチェックする。検証に使ったPCはmacbook air(Mac OS X 10.9.5,1.6Ghz Intel Core i5,2GB)、Pythonのバージョンは2.7。実験に用いるソースコードは以下となる。

モンテカルロ法ソースコード(Python2.7+numpy+matplotlib)

n=100の場合の実行結果

図(左から1,2,3回目)

f:id:program_study:20160106194255p:plainf:id:program_study:20160106194303p:plainf:id:program_study:20160106194311p:plain

回数 1回目 2回目 3回目 標本平均 標本標準偏差
x(赤い点の数)*3 73 80 78 77 3.605
円周率の近似値 2.92 3.20 3.12 3.08 0.144
実行時間 0.201 0.225 0.221 0.22 0.012

n=1,000の場合

図(左から1,2,3回目)

f:id:program_study:20160106194324p:plainf:id:program_study:20160106194330p:plainf:id:program_study:20160106194337p:plain

計測結果
回数 1回目 2回目 3回目 標本平均 標本標準偏差
x(赤い点の数) 791 774 800 788.3 13.20
円周率の近似値 3.164 3.096 3.200 3.153 0.052
実行時間 1.681 1.827 1.731 1.747 0.074

n=10,000の場合

図(左から1,2,3回目)

f:id:program_study:20160106194354p:plainf:id:program_study:20160106194402p:plainf:id:program_study:20160106194408p:plain

計測結果
回数 1回目 2回目 3回目 標本平均 標本標準偏差
x(赤い点の数) 7870 7838 7885 7864 24.00
円周率の近似値 3.148 3.135 3.154 3.145 0.096
実行時間 17.606 17.279 17.179 17.35 0.2233

最後に:実験2の講評

乱数の個数に応じ、信頼区間から円周率の近似値を割り出し、実際に計測してみた。すると信頼区間どおりに行かない事が多かった。信頼区間の計算の仕方、あるいは話の立て方が違っていたのか?この辺も反省するべき点だと考えている。

今日の記事を通し統計的な知識を引っ張って来て、実験の予測立てができたような気がする。今後はどういうタイミングで検定すべきか?と言う事に気を払っていきたい。以下、実験2を通して思った事を書いて行く。

状況に応じて試行回数nを決めよう

複数回モンテカルロ法を実行すると、試行回数nに依らず各回において結果にばらつきがある事がわかる。このような事から複数回実行する上での数値のばらつきを考慮した上で、モンテカルロ法を使う必要があると感じた。

プログラムの実行時間と数値の精度のどちらを優先するか?状況に応じて試行回数nを決定する必要がある事が分かった。

試行回数を10倍に上げたときの実行時間は?

今度はプログラムの実行時間について講評したい。まずプログラムのループ部分に注目すると、n回ループしているので計算量のオーダーはO(n)である事が分かる。次に試行回数をk倍すると、ループの回数もk倍となる。

最後にループ回数をnとしたときのプログラムの実行時間をT(n)とおき、T(kn) = kT(n) (比例)と実行時間がk倍となる事を予想する。

そこで実際にn=1,000とn=10,000のときの標本平均値を比較してみる。試行回数が10倍となったとき、17.35÷1.747 = 9.33(倍)となり、実際に約10倍実行時間がかかってしまった。

Pythonで大きいnを扱うのは不向き?

統計的な理由により、精度を上げて円周率の近似値を取る為にn=1,000,000回程度試行しなくてはならない。しかし現実にはn=10,000を処理するのに17.35[秒]かかる。このためn=1,000,000を処理するのに約1,735[秒]と約30分かかると見積もることができる。

流石にアニメ1本見ながら実験結果待ちをするわけにはいくまい。計算量を減らすアルゴリズムが他にあるのか?あるいはC言語のようなコンパイラ型の言語を使って実装するべきなのか?この辺も今後の課題である。

お世話になったWebサイト

*1:Numpyによる乱数生成まとめ:Qiita(yubais氏)より

*2:この界隈で言う「確率の収束」と言うものに当たる

*3:プログラムで計測するのを忘れた為、\pi=\frac{4x}{n}の式から個数を再計算した。