ヒストグラムに折れ線グラフを追加して、近似曲線を描画するコード

Question

np.random.seed(123)
x=np.random.randn(1000)

fig,ax=plt.subplots()

counts,edges,patches=ax.hist(x,bins=25)
#近似曲線に用いる点を求める（ヒストグラムのビンの中点）
x_fit=(edges[:-1]+edges[1:])/2
#近似曲線をプロット
y=1000*np.diff(edges)*np.exp(-x_fit**2/2)/np.sqrt(2*np.pi)
ax.plot(x_fit,y)
plt.show()

python資格試験取得のために勉強しているのですが、ヒストグラムに折れ線グラフを追加して、近似曲線を描画する上記のコードで、下から3番目のコードで質問があります。yの値を求めようとしていること、np.diff、np.exp、np.sqrt、np.pi各々の意味は理解しているのですが、全体として式が複雑で、式の意味が理解出来ません。どなたかご解説いただくことは可能でしょうか。

Answer 1

np.exp(-(x_fit ** 2) / 2) / np.sqrt(2 * np.pi) の部分ですが、これは標準正規分布の確率密度関数の値を計算しています。

コード全体を順番に説明していくと、

1. x = np.random.randn(1000) で標準正規分布に従う観測値を1000個取得する

numpy.random.randn — NumPy v1.15 Manual

1. counts, edges, patches = ax.hist(x, bins=25) でヒストグラムを計算する。

2. np.exp(-(x_fit ** 2) / 2) / np.sqrt(2 * np.pi) で各ビンの中点における標準正規分布の確率密度関数の値を計算する。

3. ヒストグラムは規格化 (Σ ビンiの幅 * ビンiの度数 = 1) されていないので、確率密度関数の値を 1000 * np.diff(edges) を乗算してヒストグラムのスケールに合わせる。

ヒストグラムの度数 (count) を規格化するには、

規格化された度数 = 度数 / (ビンの幅 * サンプル数)

で計算されるので、逆に確率密度関数の値を元のヒストグラムのスケールに合わせるには、

確率密度関数の値 * (ビンの幅 * サンプル数)

としてスケールすればよい。

---

コードにコメントを入れて整理しました。

このコードの趣旨は、ある分布に従う値をいくつかサンプルしてヒストグラムを作り、そこにその分布の確率密度関数を一緒に描画したら、ヒストグラムでその分布を近似できていることが確認できるという内容だと思います。

import numpy as np
from matplotlib import pyplot as plt

# 乱数のシード初期化
np.random.seed(123)

# 標準正規分布に従う観測値を1000個取得する。
x = np.random.randn(1000)

fig, ax = plt.subplots()

# ヒストグラムを作成する。
counts, edges, patches = ax.hist(x, bins=25)

# 各ビンの中点
x_fit = (edges[:-1] + edges[1:]) / 2


def pdf(x):
    """標準正規分布の確率密度関数
    """
    return np.exp(-(x ** 2) / 2) / np.sqrt(2 * np.pi)


# 確率密度関数の値をヒストグラムに合わせてスケールする。
y = 1000 * np.diff(edges) * pdf(x_fit)

ax.plot(x_fit, y)
plt.show()

Answer 2

y=1000*np.diff(edges)*np.exp(-x_fit**2/2)/np.sqrt(2*np.pi)

1000: 図の大きさに合わせている
np.diff(edges): ヒストグラムの箱の幅を表している

np.exp(-x_fit**2/2)/np.sqrt(2np.pi):正規分布の確率密度関数を表している
ここ複雑ですが1つずつ追ってくとよく見る正規分布の確率密度関数の式ですね
数式だと前後逆転して
1/sqrt(2π)*exp(-x^2/2)
な感じかと思いますが

c.f. Wikipedia:正規分布