はじめに
E資格の次は統計検定の合格を！ということで、早速Pythonで統計について理解を深めていきたい。そこでまずは正規分布と戯れてみた。
正規分布とは
以下の確率密度関数で定義される分布である。
$$ f(x) = \frac{1}{\sqrt{2πσ^2}} \exp\left(-\frac{(x-u)^2}{2σ^2}\right) $$
確率密度関数をPythonで書いてみよう。
腐るほど、いろんな人がやったネタであるが、自分で書いてみよう。
確率変数Xと平均と分散を受け取って確率を返すというもの。
numpyなので複数のXに対して計算できちゃうぞ。
# 正規分布関数
def normal(x, u, v):
ret = 1 / np.sqrt(2 * np.pi * v) * np.exp(-(x-u)**2/(2*v))
return ret
#描画してみよう
先ほどの関数を使って平均0、分散1の標準正規分布を書いてみよう。
# matplotlibで描画
x = np.linspace(-5, 5,100)
y = normal(x, 0, 1)
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(1,1,1)
ax.plot(x, y)
するとこんな感じ。もう腐るほどみた図だと思う。
#モンテカルロ法により面積を求めてみよう
さて、確率密度関数というのは、-∞から∞まで積分すると1になるはずだ。
積分は大変なので、今回これをモンテカルロ法により求めてみよう。
モンテカルロ法とは、確率変数のサンプリングをコンピューター上で行うことで、数学的な問題を数値的に解く手法である。
今回は一様分布を用いる。
百聞は一見に如かずで、まずはXについて[-5, 5]、Yについて[0, 0.4]にの一様分布に従うX, Yを生成・プロットし、その上に先ほどの確率密度関数を重ねてみよう。
# -5～5までの間で一様にXを生成
randomX = -5 + np.random.rand(30000) * 10
# 0～0.4までの間で一様にYを生成
randomY = np.random.rand(30000) * 0.4
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1,1,1)
ax.scatter(randomX, randomY)
ax.plot(x, y, c="r")
そうするとこんなプロットが得られるはずだ
今、関数の下側の面積を求めたいわけだから、ランダムに領域内にプロットされた点の内、確率密度関数より下側にある点の割合を求め、領域全体面積(10 * 0.4)を掛ければよい。
※ -5～5の範囲外の確率はほぼ0近いので無視していることに注意してほしい。
コードは下の通りである。
upperCount = 0 #上にあるか
lowerCount = 0 #下にあるか
for i, (randomX_tmp, randomY_tmp) in enumerate(zip(randomX, randomY)):
normalY = normal(randomX_tmp, 0, 1)
if randomY_tmp < normalY:
lowerCount += 1
else:
upperCount += 1
print(lowerCount, upperCount, (lowerCount / len(randomX)) * 10 * 0.4)
計算すると1.008とほぼ1に近い値が得られた。
562 22439 1.0082666666666666
おわりに
次回からは様々な分布を取り扱っていきたい。
統計学の本は非常に抽象的に記載されているので、Pythonを使って具体的なイメージに落とし込んでいきたい。乞うご期待。