スポンサーリンク

ヒストグラムと確率分布の同時プロット―フィッティング精度の検証用に

プログラミング

多数のデータ点が得られたとき

のようなフィッティングを行えば、それらのデータ点が従う分布を推定することができる。この場合、データのヒストグラムはその分布に類似した形状をとるはずであり、そのことを視覚的に確かめたい(または、その推定が正しいことを他者に主張したい)場合がある。

そのためには、データのヒストグラムとフィッティングされた確率分布(確率密度関数)を同時にプロットする必要があるが、それぞれ縦軸のスケールが異なるため、どちらかのスケールを変換しなければならない。ここではPythonによるプロットを例に、その手法について説明する(※プロット手法の概念を理解するためには、Pythonコードが読める必要はない)。

なお、この記事ではまず失敗例と成功例について示し、最後にその背景にある理論をサラッと説明する。

データの準備

ここでは、プロット例を作成するために、scikit-learnに標準で含まれている乳がんデータセット(breast_cancer)を用いる。その中の「良性」ラベルの「worst area」のみを抜き出し、ガウス分布へのフィッティングを行う。

import numpy as np
from sklearn import datasets
from scipy.stats import norm

data = datasets.load_breast_cancer() # 乳がんデータ読み込み
data_w = data["data"][:,np.where(data["feature_names"] == "worst area")[0][0]] # worst area のデータのみ抜き出す
data_wb = data_w[data["target"] == 1] # 良性のデータのみ用いる

mu, sigma = norm.fit(data_wb) # 正規分布にフィッティング
print("mu={}, sigma={}".format(mu, sigma)) # mu=558.8994397759103, sigma=163.3721295593988

その結果、\(\mu\fallingdotseq 558.9\), \(\sigma\fallingdotseq 163.4\) を得る。

データのプロット

ヒストグラムのスケールを変える

ヒストグラムの高さを調整し、それが確率密度関数と同じスケールになるようにする。

失敗例

ヒストグラムを確率密度関数に変換するため(?)、ヒストグラムの高さをすべてのデータ数で割る。

from matplotlib import pyplot as plt

# a: ヒストグラムの高さ, b: ビンの座標
a, b, _ = plt.hist(data_wb, bins=20)
plt.close() # a, b を取得して一度閉じる

n = len(data_wb) # データ点の数
b_width = b[1] - b[0] # ビンの間隔
a /= n # 高さをデータ数で割る
plt.bar(x=b[:-1], height=a, width=b_width, align="edge") # 再表示
plt.plot(b, norm.pdf(b, loc=mu, scale=sigma), color='#ff7f0e') # 確率分布をプロット
plt.title("hist -> pdf (Failed)")
plt.xlabel("worst area")
plt.ylabel("PDF?")
plt.savefig("hist_pdf_failed.png")
plt.show()

成功例

ヒストグラムを確率密度関数に変換するため、ヒストグラムの高さを、「すべてのデータ数 × ヒストグラムのビンの幅」で割る。

from matplotlib import pyplot as plt

# a: ヒストグラムの高さ, b: ビンの座標
a, b, _ = plt.hist(data_wb, bins=20)
plt.close() # a, b を取得して一度閉じる

n = len(data_wb) # データ点の数
b_width = b[1] - b[0] # ビンの間隔
a /= n * b_width # 高さを(データ数×ビンの幅)で割る
plt.bar(x=b[:-1], height=a, width=b_width, align="edge") # 再表示
plt.plot(b, norm.pdf(b, loc=mu, scale=sigma), color='#ff7f0e') # 確率分布をプロット
plt.title("hist -> pdf (Success)")
plt.xlabel("worst area")
plt.ylabel("PDF")
plt.savefig("hist_pdf_success1.png")
plt.show()

なお、この手法は plt.hist(density=True) とすることによっても適用できる。

from matplotlib import pyplot as plt

# ヒストグラムの下の面積を 1 にスケーリング
_, b, _ = plt.hist(data_wb, bins=20, density=True)
plt.plot(b, norm.pdf(b, loc=mu, scale=sigma))
plt.title("hist -> pdf (Success)")
plt.xlabel("worst area")
plt.ylabel("PDF")
plt.savefig("hist_pdf_success2.png")
plt.show()

確率分布のスケールを変える

確率分布に従って、全データ数と同じ数のサンプルを発生させたときの予測値をプロットする。このとき、プロットの x 座標にはヒストグラムの各ビンの中心を用いる。

失敗例

全データ数に、x 座標における確率密度関数の値をかける。

from matplotlib import pyplot as plt

_, b, _ = plt.hist(data_wb, bins=20)
b_mid = np.array([(b[i] + b[i+1]) / 2 for i in range(len(b)-1)]) # ビンの中心の座標を取得
n = len(data_wb) # データ点の数
# データ数に、確率密度関数の値をかける
y = norm.pdf(b_mid) * n
plt.plot(b_mid, y) # ビンの中心に、確率分布から予想される高さをプロット
plt.title("pdf -> hist (Failed)")
plt.xlabel("worst area")
plt.ylabel("number")
plt.savefig("pdf_hist_failed.png")
plt.show()

成功例

各ビンの幅における確率密度関数の積分を考える。これは、各ビンの両端における累積分布関数の差として求めることができる。

$$\int_{\alpha}^{\beta}p(x)dx=\int_{-\infty}^{\beta}p(x)dx-\int_{-\infty}^{\alpha}p(x)dx$$

この値を全データ数にかけ、ビンの中心にプロットする。

from matplotlib import pyplot as plt

_, b, _ = plt.hist(data_wb, bins=20)
b_mid = np.array([(b[i] + b[i+1]) / 2 for i in range(len(b)-1)]) # ビンの中心の座標を取得
cum = norm.cdf(b, loc=mu, scale=sigma) # 各ビンの端までの累積分布関数の値を計算
n = len(data_wb) # データ点の数
# データ数に、ビンの両端の累積分布関数の値の差をかける
y = [n * (cum[i+1] - cum[i]) for i in range(len(b_mid))]
plt.plot(b_mid, y) # ビンの中心に、確率分布から予想される高さをプロット
plt.title("pdf -> hist (Success)")
plt.xlabel("worst area")
plt.ylabel("number")
plt.savefig("pdf_hist_success.png")
plt.show()

解説

失敗例はともに、確率密度関数割合の関係を誤解したことが問題となっている。

確率密度関数 \(p(x=a)\) の値は、\(x=a\) となる確率を表したものではない。連続値をとる \(x\) がちょうど \(a\) となる確率は \(0\) である。

このあたりの議論は

を参照。

確率密度はあくまで「密度」であることに注意する。確率密度関数は、値の区間で積分をすることによって「確率変数がその区間に入る確率」を計算できる関数のことである。そのため、確率分布からヒストグラムを再現する際には、累積分布関数を考える必要がある。

また、ヒストグラムは「各ビンの範囲に入るデータ点の数を高さとした棒グラフ」のことである。したがって、これを確率「密度」関数に変換するためには、そのビンの幅で割る必要がある。

コメント