プログラミング学習：モンテカルロ法で円周率を求める（GCI chapter2 python、pandas）

東大　GCIデータサイエンティスト育成講座演習コンテンツ Chapter 2の課題です。

乱数を発生させる方法を使って、円周率を求めるプログラムを作成する。

• 一様分布に従う乱数を2組発生させて、それぞれ10,000個の一様乱数を作る。
• x−y軸を使った中心(0,0)、半径1の円と、長さ1の正方形を考える。ここで先ほどのxとyの組み合わせの乱数10,000個のうち、円の内部に入る点は何組あるか。
• 半径1の1/4の円の面積と長さ1の長方形の面積の比は、π/4:1となる。これを利用して、円周率を求めよ。

モンテカルロ法については下記ページが参考になりました。

　モンテカルロ法と円周率の近似計算

1x1の正方形にランダムに点をN個打ち、そのうち半径1の円（の4分の1）内に入った点の数をXとすると、4X/Nがπの近似値になる。図を見るとイメージがつきやすいかと思います。

これをpandasを使ってやってみます。

・モジュールのインポート（numpy,pandas)

import numpy as np
import numpy.random as random
import pandas as pd

・乱数の発生

x=np.random.uniform(0, 1, 10000)
y=np.random.uniform(0, 1, 10000)
print(x,y)

[0.15469463 0.35501088 0.31071841 ... 0.67647336 0.10477222 0.77432113] [0.56009742 0.42374295 0.5874093  ... 0.62031945 0.66022409 0.92282591]

np.random.uniform(0, 1, 10000)で、0以上1未満の一様乱数を10,000個発生させ、それぞれの配列をx、yに格納します。

・x、yの配列データフレームに変換

montecarlo = pd.DataFrame({'x':x,'y':y})
montecarlo.head()

x
y

0
0.790832
0.774845
1
0.568365
0.950229
2
0.835321
0.895938
3
0.967860
0.286567
4
0.624064
0.281333

x,yを列名とし、それぞれの配列をデータフレームに変換しました。

・√x²+y²を計算する。

montecarlo["hypot"]=np.hypot(montecarlo["x"],montecarlo["y"])
montecarlo.head()

x
y
hypot

0
0.790832
0.774845
1.107159
1
0.568365
0.950229
1.107238
2
0.835321
0.895938
1.224935
3
0.967860
0.286567
1.009393
4
0.624064
0.281333
0.684547

行ごとにnp.hypotで√x²+y²を計算し、新しい列 hypotに値を格納します。

・半径1の円内に点が入るかどうかの判定。

montecarlo_small=montecarlo[montecarlo["hypot"]<=1]
montecarlo_large=montecarlo[montecarlo["hypot"]>1]
montecarlo_small.head()

x
y
hypot

4
0.624064
0.281333
0.684547
5
0.411707
0.040501
0.413694
6
0.278791
0.077373
0.289328
7
0.415050
0.306739
0.516096
9
0.799246
0.122930
0.808644

hypotの値が1以下なら円内、1より大きければ円外とし、それぞれmontecarlo_small、montecarlo_largeのdataframeに格納しました。

・グラフ描画

large=montecarlo_large.plot.scatter(x="x",y="y",figsize=(8,8),color="red")
montecarlo_small.plot.scatter(x="x",y="y",ax=large,color="blue")

円外の点を赤、円内に入る点を青でそれぞれ散布図で描画します。ポイントはfigsizeで図の大きさを正方形にすること。

1番目の図を変数largeに格納し、2番目の図でax=largeと指定することで1番目の図に上書きして2番目の図が描画されます。

・出力


円内の点が青、円外の点が赤で描かれました。所々空白がありますが、ちゃんと円が浮き上がって見えます。

・近似円周率の計算

(len(montecarlo_small)*4)/10000

3.1724

半径1の内に入っている点の数を数え、4倍して全体の点の数(10000)で割ると円周率の近似値となります。lenでmontecarlo_smallの行数を数えて、これを円内に入っている点の数としました。

円周率近似値の結果は、3.1724。もっと点の数を多くすると精度は上がりそうです。

グラフ描画などもっとスマートな方法がありそうな気がしますが、ひとまずこれで。

参考：
東京大学グローバル消費インテリジェンス寄付講座
http://gci.t.u-tokyo.ac.jp/
東京大学松尾研　GCIデータサイエンティスト育成講演習コンテンツ 公開ページ
http://weblab.t.u-tokyo.ac.jp/gci_contents/
プログラミング学習：東京大学のデータサイエンティスト育成講座（GCI）の演習コンテンツ
https://eneprog.blogspot.com/2018/04/python10gci-chapter1.html

モンテカルロ法と円周率の近似計算
https://mathtrain.jp/montecarlo

プログラミング学習:10までの素数を表示させる。(python GCI chapter1）
https://eneprog.blogspot.com/2018/04/python10gci-chapter1.html
プログラミング学習:numpyをつかって1から50までの自然数の和を計算する。(python GCI chapter2）
https://eneprog.blogspot.com/2018/04/numpy150python-gci-chapter2.html