『東京大学の データサイエンティスト育成講座』という本の第2章の章末問題に、モンテカルロ法により円周率を求めてみようというものがありました。０から１までの範囲で一様な乱数を2個発生させ（x,y)、原点からの距離が1未満になる個数の割合から、円周率を求めるという原理です。4点(0,0),(1,0),(1,1),(0,1)がつくる正方形の面積は1で、4分の１の円の面積（半径＝１）は、円周率ｘ1×1/4 なので、上記の割合が面積の割合になることから、円周率が求まります。

numpy配列の要素のうち、ある条件を満たす要素のみを抽出する方法がわからなくてあれこれ調べてみたら、配列aからの要素の抽出は、

a[条件]

という書式でいけることがわかりました。たとえば、乱数を発生させて0.5未満のみ選びたければ、下のようになります。

x = np.random.uniform(0.0, 1.0, 10)
print(x)
x1=x[x<0.5]
print(x1)

実行結果の出力は、

[0.3476324  0.30311298 0.61869065 0.29137365 0.82916676 0.16149072
 0.02440774 0.36505252 0.2109926  0.72545332]
[0.3476324  0.30311298 0.29137365 0.16149072 0.02440774 0.36505252
 0.2109926 ]

という具合です。(x,y)の原点からの距離は、np.hypot(x,y)で求められます(npはnumpy）。math.hypot(x,y)を使った場合には、エラーになりました。

TypeError: only size-1 arrays can be converted to Python scalars

np.hypot(x,y)だとこのエラーは回避できました。（ｘ、ｙ）の点の描画は、点の形状が「小さな円」なら、’o’、「点」は、’.’（ピリオド）、「ピクセル」は’,’（カンマ））で指定して描けます。これらの知識を使って『東京大学の データサイエンティスト育成講座』第2章の章末問題やってみると、

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

n = 10000
x = np.random.uniform(0.0, 1.0, n)
y = np.random.uniform(0.0, 1.0, n)
d = np.hypot(x,y)
plt.figure(figsize = (5,5))

plt.plot(x[d<1], y[d<1],’,’, color=’red’)
plt.plot(x[d>=1],y[d>=1],’,’,color=’blue’)
print(‘円周率:’,4*sum(d<1)/n)

この方法で計算した円周率は、n = 10000 でも、例えば円周率: 3.1384などとなり、あまり正確ではありませんでした。もっと数を増やさないといけないようです。n = 10000000 にしたら、計算に15秒くらいかかりましたが、円周率: 3.1414896　と正確さが増しました。