Python：正規分布に関する関する関数（scipy.stats）のまとめ

16055 ワード
Python stats scipy 正規分布 Python テキストリンク

はじめに

統計学に関するプログラミングを勉強中です。正規分布を扱う時に便利な関数をまとめました。

紹介する関数

  • scipy.stats
  • scipy.stats.norm.pdf
  • scipy.stats.norm.rvs
  • scipy.stats.norm.cdf
  • scipy.stats.norm.ppf

正規分布の確率密度関数

$\mu$を標本平均、$\sigma^2$を標本分散（$\sigma$は標本偏差）とすると、正規分布の確率密度関数は以下になります。

N(x | \mu, \sigma^{2} ) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}

この式をPythonで直接コーディングするとすると、

sample01
>> import numpy as np

>> x = 1 
>> mu = 2
>> sigma = 0.4

>> 1 / (np.sqrt(2 * np.pi * sigma ** 2)) * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))

0.04382075123392136

このように長くなってしまいます。

Pythonのscipy.stats.pdfを使うことで、以下のように1文で記載ができます。

sample02
>> from scipy import stats
>> stats.norm.pdf(loc = 2, scale = 0.4, x = 1)

0.04382075123392136

正規分布は英語で、Normal Distributionといいます。また、確率密度関数は英語で、Probability Density Functionといいますので、この頭文字をとってpdfとなります。

ここで引数として、 loc = 平均、scale =　標準偏差　となります。

この関数さえ覚えておけば、面倒な正規分布の確率密度関数を改めて作成する必要がなく便利です。

また、以下のサンプルのようにstats.normをインスタンス化して、再利用することができます。

sample03
>> from scipy import stats
>> norm_dist = stats.norm(loc = 2, scale = 0.4)
>> V_x = np.arange(1, 11, 1)
>> norm_dist.pdf(x = V_x).round(3)

array([0.   , 0.022, 0.228, 0.499, 0.228, 0.022, 0.   , 0.   , 0.   ,
       0.   ])

matplotlibパッケージを利用して、正規分布のグラフを作成します。以下はJupyter Notebookを使った場合の例です。

sample04
#パッケージのインポート
import numpy as np
from scipy import stats
from matplotlib import pyplot as plt
plt.style.use('seaborn')
%matplotlib inline 

#変数の設定
x = np.arange(0, 4.1, 0.1)
mu = 2
sigma = 0.4

#グラフの作画
plt.plot(x, stats.norm.pdf(loc=mu, scale=sigma, x = x))

以下のグラフが作成されます。x = 2.0を中心（平均値）として左右対称の山のグラフが表示されます。

正規分布に従った母集団からのサンプリング

母集団から、無作為にサンプル抽出を行う関数としては、numpy.random.choiceがあります。

sample05
import numpy as np

#仮の母集団を作成（標準分布に従う）
population = np.random.standard_normal(100000) 

#仮の母集団から100個のサンプルを無作為抽出
sample = np.random.choice(population, size=100, replace=False) 

(実行結果)
array([-0.354,  1.001, -0.715,  1.493,  1.096, -0.044,  0.085,   
...(省略)...
        1.597,  0.027,  0.655, -0.72 ])

ここで引数のreplaceは、サンプルで重複を許すかを指定します。Falseだと重複せずに無作為抽出し、Trueだと重複を許します。(replaceは「引いたクジを元に戻す」のキーワードです）

母集団が正規分布に従っている場合は、stats.norm.rvsによって簡潔にコーディングできます。

sample06
import numpy as np
from scipy import stats

#母集団から100個のサンプルを無作為抽出
sample = stats.norm.rvs(loc=2.0, scale=0.4, size=100) 

(実行結果)
array([2.043, 1.577, 1.648, 2.167, 2.565, 2.579, 1.71 , 0.766,  
...(省略)...
       1.864])

rvsは、random variablesの略になります。

累積分布関数

$X$を確率変数、$P$を確率測度とするとき、任意の実数$x$に対して、

F(x) = P(X\leq x)

を累積分布関数または分布関数といいます。

サイコロを例にとると、確率変数$X$が取り得る値は、

1, 2, 3, 4, 5, 6

となり、この時のすべてのサイコロの目$x$で、

F(x) = P(X\leq x) = \frac{x}{6}

となります。

正規分布のような確率速度が連続型の関数の場合は、積分を用いて以下のように書きます。

F(x) = \int^x_{-\infty} \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}dx

Pythonには積分計算するパッケージがありますが、これをコーディングするのは煩雑になります。

scipy.stats.norm.cdfがこの計算を行います。

sample07
from scipy import stats

#１以下の累積分布関数の値
stats.norm.cdf(loc=2.0, scale=0.4, x=1)

(実行結果)
0.006209665325776132

cdfはCumulative Distribution Functionの略です。なお、このようにデータがある値以下となる確率を下側確率とよびます。

またこの逆で、下側確率が$y\%$となる%x%を求めるには、stats.norm.ppf関数を使います。

ppfは、Percent Point Functionの略です。

sample08
from scipy import stats
y = 0.25
stats.norm.ppf(loc=2, scale=0.4, q = y)

(実行結果)
1.7302040999215673

ちなみに正規分布は平均を中心として左右対称ですので、y=50%とした場合、xは標本平均と同値になります。

sample09
from scipy import stats
y = 0.5
stats.norm.ppf(loc=2, scale=0.4, q = y)

(実行結果)
2.0

まとめ

以上、正規分布を扱う場合に便利なPythonの関数を紹介しました。