「gnuplot で正規分布をσごとに塗りわける」の Python 版。
以下を参考に,標準正規分布がもつ確率密度関数のグラフを描いてみる。
正規分布関数
$$N(\mu, \sigma^2) \equiv \frac{1}{\sqrt{2\pi} \sigma} \exp\left\{ – \frac{(x-\mu)^2}{2\sigma^2}\right\}$$
簡単のために,$\mu = 0, \sigma = 1$ とした標準正規分布 $f(x)$ を使うことにする。
$$f(x) \equiv N(0, 1) = \frac{1}{\sqrt{2\pi}} \exp\left( – \frac{x^2}{2}\right)$$
誤差関数
$$\mbox{erf}(x) \equiv \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2}\, dt $$
正規分布の積分を誤差関数で表す
\begin{eqnarray}
F(x_1, x_2) &\equiv& \int_{x_1}^{x_2} f(t)\, dt \\
&=& \frac{1}{\sqrt{2\pi}} \int_{x_1}^{x_2}\exp\left( – \frac{t^2}{2}\right)\, dt \\
&& \quad\left(\frac{t}{\sqrt{2}} \equiv x, \ dt = \sqrt{2} dx \right) \\
&=& \frac{1}{\sqrt{\pi}} \int_{\frac{x_1}{\sqrt{2}}}^{\frac{x_2}{\sqrt{2}}} e^{-x^2} \,dx \\
&=& \frac{1}{2} \left(\mbox{erf}\left(\frac{x_2}{\sqrt{2}}\right)-\mbox{erf}\left(\frac{x_1}{\sqrt{2}}\right) \right)
\end{eqnarray}
誤差関数 scipy.special.erf(x)
を使います。
import matplotlib.pyplot as plt
import numpy as np
from scipy import special
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
# グラフのサイズ
plt.rcParams["figure.figsize"] = (6.4, 2.4)
def f(x):
return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)
def F(x1, x2):
return 0.5*(special.erf(x2/np.sqrt(2))-special.erf(x1/np.sqrt(2)))
標準正規分布を plt.plot()
x = np.linspace(-4, 4, 401)
# 横軸縦軸の表示範囲
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.plot(x, f(x), color='#0000ff');
σごとに塗りわける
# 塗りわける色
c = ['#0000ff','#1e90ff','#87ceeb','#add8e6']
# σごとに塗りわける
for i in range(4):
plt.fill_between(x, f(x), where=(-i-1<=x)&(x<= -i), fc=c[i])
plt.fill_between(x, f(x), where=( i<=x)&(x<= i+1), fc=c[i])
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.plot(x, f(x), color=c[0]);
塗り分けた領域を白線で区切る
# vlines で白い縦線を引く
xsigma = np.linspace(-3, 3, 7)
plt.vlines(xsigma, 0, f(xsigma), color='white')
for i in range(4):
plt.fill_between(x, f(x), where=(-i-1<=x)&(x<= -i), fc=c[i])
plt.fill_between(x, f(x), where=( i<=x)&(x<= i+1), fc=c[i])
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.plot(x, f(x), color=c[0]);
横軸の目盛ラベルを設定する
# 横軸の目盛の位置とラベルの設定
plt.xticks(
ticks = [-4, -3, -2, -1, 0, 4, 3, 2, 1],
labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
"μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"])
xsigma = np.linspace(-3, 3, 7)
plt.vlines(xsigma, 0, f(xsigma), color='white')
for i in range(4):
plt.fill_between(x, f(x), where=(-i-1<=x)&(x<= -i), fc=c[i])
plt.fill_between(x, f(x), where=( i<=x)&(x<= i+1), fc=c[i])
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.plot(x, f(x), color=c[0]);
塗り分けた領域にラベルをつける
# 各領域のパーセントを表す文字列を配列として作成
# リスト sig の要素を全て 0 に
sig = [0]*4
for i in range(0, 4):
sig[i] = '%5.2f' % (F(i,i+1)*100) + '%'
print(sig[i])
# 領域にラベルをつける
plt.text(-0.5,0.20, sig[0], ha = 'center', color="white")
plt.text( 0.5,0.20, sig[0], ha = 'center', color="white")
plt.text(-1.5,0.04, sig[1], ha = 'center', color="white")
plt.text( 1.5,0.04, sig[1], ha = 'center', color="white")
plt.text(-2.5,0.05, sig[2], ha = 'center', color="black")
plt.text( 2.5,0.05, sig[2], ha = 'center', color="black")
plt.text(-3.5,0.02, sig[3], ha = 'center', color="black")
plt.text( 3.5,0.02, sig[3], ha = 'center', color="black")
plt.xticks(
ticks = [-4, -3, -2, -1, 0, 4, 3, 2, 1],
labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
"μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"])
xsigma = np.linspace(-3, 3, 7)
plt.vlines(xsigma, 0, f(xsigma), color='white')
for i in range(4):
plt.fill_between(x, f(x), where=(-i-1<=x)&(x<= -i), fc=c[i])
plt.fill_between(x, f(x), where=( i<=x)&(x<= i+1), fc=c[i])
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.plot(x, f(x), color=c[0]);
枠と縦軸目盛を非表示に
# 枠と縦軸目盛を非表示に
fig, ax = plt.subplots()
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.text(-0.5,0.20, sig[0], ha = 'center', color="white")
plt.text( 0.5,0.20, sig[0], ha = 'center', color="white")
plt.text(-1.5,0.04, sig[1], ha = 'center', color="white")
plt.text( 1.5,0.04, sig[1], ha = 'center', color="white")
plt.text(-2.5,0.05, sig[2], ha = 'center', color="black")
plt.text( 2.5,0.05, sig[2], ha = 'center', color="black")
plt.text(-3.5,0.02, sig[3], ha = 'center', color="black")
plt.text( 3.5,0.02, sig[3], ha = 'center', color="black")
plt.xticks(
ticks = [-4, -3, -2, -1, 0, 4, 3, 2, 1],
labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
"μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"])
xsigma = np.linspace(-3, 3, 7)
plt.vlines(xsigma, 0, f(xsigma), color='white')
for i in range(4):
plt.fill_between(x, f(x), where=(-i-1<=x)&(x<= -i), fc=c[i])
plt.fill_between(x, f(x), where=( i<=x)&(x<= i+1), fc=c[i])
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.plot(x, f(x), color=c[0]);
-nσ から +nσ までを塗りつぶす
# 各領域のパーセントを表す文字列を配列として作成
# リスト sig の要素を全て 0 に
sig2 = [0]*4
for i in range(0, 4):
sig2[i] = '%5.2f' % (F(-i-1,i+1)*100) + '%'
print(sig2[i])
# 枠と縦軸目盛を非表示に
fig, ax = plt.subplots()
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.xticks(
ticks = [-1, 0, 1],
labels = ["μ-σ", "μ",
"μ+σ"])
plt.plot(x, f(x), color='gray')
# 領域にラベルをつける
plt.text(0,0.15, sig2[0], ha = 'center', color="white")
# -1σ から +1σ までを塗りつぶす
x1 = np.linspace(-1, 1, 101)
plt.fill_between(x1, f(x1), fc=c[0])
plt.plot(x1, f(x1), color=c[0]);
# 枠と縦軸目盛を非表示に
fig, ax = plt.subplots()
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.xticks(
ticks = [-2, -1, 0, 2, 1],
labels = ["μ-2σ", "μ-σ", "μ",
"μ+2σ", "μ+σ"])
plt.plot(x, f(x), color='gray')
# 領域にラベルをつける
plt.text(0,0.15, sig2[1], ha = 'center', color="white")
# -2σ から +2σ までを塗りつぶす
x2 = np.linspace(-2, 2, 101)
plt.fill_between(x2, f(x2), fc=c[1])
plt.plot(x2, f(x2), color=c[0]);
# 枠と縦軸目盛を非表示に
fig, ax = plt.subplots()
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.xticks(
ticks=[-3, -2, -1, 0, 3, 2, 1],
labels=["μ-3σ", "μ-2σ", "μ-σ", "μ",
"μ+3σ", "μ+2σ", "μ+σ"])
plt.plot(x, f(x), color='gray')
# 領域にラベルをつける
plt.text(0,0.15, sig2[2], ha = 'center', color="black")
# -3σ から +3σ までを塗りつぶす
x3 = np.linspace(-3, 3, 101)
plt.fill_between(x3, f(x3), fc=c[2])
plt.plot(x3, f(x3), color=c[0]);
# 枠と縦軸目盛を非表示に
fig, ax = plt.subplots()
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.xlim(-4, 4)
plt.ylim(0, 0.41)
plt.xticks(
ticks = [-4, -3, -2, -1, 0, 4, 3, 2, 1],
labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
"μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"])
# 領域にラベルをつける
plt.text(0,0.15, sig2[3], ha = 'center', color="black")
# -4σ から +4σ までを塗りつぶす
plt.fill_between(x, f(x), fc=c[3])
plt.plot(x, f(x), color=c[0]);