Python の Matplotlib で正規分布をσごとに塗りわける

正規分布関数

$$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) を使います。

In [1]:
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)
In [2]:
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()

In [3]:
x = np.linspace(-4, 4, 401)

# 横軸縦軸の表示範囲
plt.xlim(-4, 4)
plt.ylim(0, 0.41)

plt.plot(x, f(x), color='#0000ff');

σごとに塗りわける

In [4]:
# 塗りわける色
c = ['#0000ff','#1e90ff','#87ceeb','#add8e6']
In [5]:
# σごとに塗りわける
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]);

塗り分けた領域を白線で区切る

In [6]:
# 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]);

横軸の目盛ラベルを設定する

In [7]:
# 横軸の目盛の位置とラベルの設定
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]);

塗り分けた領域にラベルをつける

In [8]:
# 各領域のパーセントを表す文字列を配列として作成
# リスト sig の要素を全て 0 に
sig = [0]*4
for i in range(0, 4):
    sig[i] = '%5.2f' % (F(i,i+1)*100) + '%'
    print(sig[i])
34.13%
13.59%
 2.14%
 0.13%
In [9]:
# 領域にラベルをつける
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]);

枠と縦軸目盛を非表示に

In [10]:
# 枠と縦軸目盛を非表示に
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σ までを塗りつぶす

In [11]:
# 各領域のパーセントを表す文字列を配列として作成
# リスト sig の要素を全て 0 に
sig2 = [0]*4
for i in range(0, 4):
    sig2[i] = '%5.2f' % (F(-i-1,i+1)*100) + '%'
    print(sig2[i])
68.27%
95.45%
99.73%
99.99%
In [12]:
# 枠と縦軸目盛を非表示に
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]);

In [13]:
# 枠と縦軸目盛を非表示に
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]);

In [14]:
# 枠と縦軸目盛を非表示に
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]);

In [15]:
# 枠と縦軸目盛を非表示に
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]);