あてはまりのよい確率分布を探したい

はじめに

データを眺めていると、ある分布に対してそれが正規分布に従うのか、対数正規分布か、それともガンマ分布の方が近いのか？、というようにどの分布の当てはまりがよいかが気になることがあると思います。

これを確認する方法を探してみたところ、scipy.statsを使えばできそうだったのと、fitterというライブラリもあったので、それらを試してみた結果を記述します。

実験

scipyを使う

実装はnumpy - Fitting empirical distribution to theoretical ones with Scipy (Python)? - Stack Overflowを少しだけ修正したものです。入力に対してscipy.statsに登録されているすべての確率分布のパラメータを最尤推定した結果の平均二乗誤差を比較することで最もあてはまりのよい分布を求めます。

scipyには80以上の確率分布が存在しているのですが、全部回してしまうと100行程度のデータでも数分かかってしまいますし、知らない分布が選ばれたところで仕方がないので、ある程度選択しておくのが現実的かなと思います。今回は正規分布、一様分布、ガンマ分布、レイリー分布の4つの分布で比較を行います。

from sklearn.datasets import load_iris
import pandas as pd
import scipy.stats as st
import matplotlib.pyplot as plt

#サンプルデータの読み込み
iris = load_iris()
data=pd.DataFrame(iris.data, columns=iris.feature_names)

#使用する分布
#正規分布、一様分布、ガンマ分布、レイリー分布
distributions=[st.norm,st.uniform,st.gamma,st.rayleigh]

result=[]
for distribution in distributions:
    y, x = np.histogram(data.iloc[:,1], bins=20, density=True)
    #xの値はヒストグラムの左端の値なので中心点に修正
    x = (x + np.roll(x, -1))[:-1] / 2.0
    
    params = distribution.fit(data.iloc[:,1])

    # Separate parts of parameters
    arg = params[:-2]
    loc = params[-2]
    scale = params[-1]

    # Calculate fitted PDF and error with fit in distribution
    pdf = distribution.pdf(x, loc=loc, scale=scale, *arg)
    sse = np.sum(np.power(y - pdf, 2.0))
    
    result.append((pdf,sse))

    df=pd.DataFrame()
    df["name"]=[i.name for i in distributions]
    df["sumsquare_error"]=[i[1] for i in result]
    df

今回の分布の中では正規分布のあてはまりが最も良いようです。

for i in range(len(result)):
    pd.Series(result[i][0], x).plot(label=distributions[i].name)
    plt.hist(data.iloc[:,1],density=True,alpha=0.4,bins=20)
    plt.legend()
    plt.show()

正規分布、ガンマ分布は元の分布に対してあてはまりが比較的良いことが確認できます。

fitterを使う

上に書いたものをライブラリ化したような感じです。

from fitter import Fitter
f = Fitter(data.iloc[:,1],distributions=['gamma', 'rayleigh', 'uniform','norm'],bins=20)
f.fit()
f.summary()

中でやっていることは同じなので結果は先ほどと同じになります。

気をつけたいこととしては、あてはまりの良い分布はbinの数によって変化することかなと思います。 fitterのbinのデフォルト値は100なのですが、irisデータはサンプルサイズが150と小さくデフォルトのbinを用いると 結果がおかしくなるため、今回の実験ではbinを20としています。

まとめ

任意のデータに対し、当てはまりの良い確率分布を導出する手法の調査を行いました。 確率分布の仮定が必要な統計モデリングを行う際や、 分布からデータ生成の仕組みなどを推測するのに役立つかもしれません。