SymPy Plotting Backends で正規分布をσごとに塗りわける

正規分布関数

$$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}

誤差関数 erf(x) を使います。

In [1]:
from sympy import *
# 1文字変数の Symbol の定義が省略できる
from sympy.abc import *
# π
from sympy import pi
Pi = float(pi)

# SymPy Plotting Backends (SPB)
from spb import *

# NumPy
import numpy as np

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
In [2]:
def f(x):
    return exp(-x**2 / 2) / sqrt(2 * Pi)

def F(x1, x2):
    return 0.5*(erf(x2/sqrt(2))-erf(x1/sqrt(2)))

標準正規分布を plot()

In [3]:
plot(f(x), (x, -4, 4), grid = False, line_color='#0000ff',
     xlim = (-4, 4), ylim = (0, 0.41), size=(6.4, 2.4));

σごとに塗りわける

In [4]:
# 塗りわける色
c = ['#0000ff','#1e90ff','#87ceeb','#add8e6']
In [5]:
# 塗りわける x の範囲
x_range = []
for i in range(4):
    x_range.append(np.linspace(i, i+1, 101))
# 対応する y の値
y_range = []
for i in range(4):
    y_range.append(lambdify(x, f(x))(x_range[i]))

Fill = ([{'x': x_range[i], 'y1':y_range[i], 'color':c[i]} for i in range(4)] +
        [{'x':-x_range[i], 'y1':y_range[i], 'color':c[i]} for i in range(4)])
In [6]:
plot(f(x), (x, -4, 4), line_color = '#0000ff', grid = False, legend = False,
     xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
     fill=Fill);

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

In [7]:
p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         fill=Fill, show=False);

fig, ax = p.fig, p.ax

# vlines で白い縦線を引く
xsigma = np.linspace(-3, 3, 7)
fxsigma = lambdify(x, f(x))(xsigma)
ax.vlines(x=xsigma, ymin=0, ymax=fxsigma, color='white');

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

In [8]:
p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill=Fill, show=False);

fig, ax = p.fig, p.ax

xsigma = np.linspace(-3, 3, 7)
fxsigma = lambdify(x, f(x))(xsigma)
ax.vlines(x=xsigma, ymin=0, ymax=fxsigma, color='white');

# 横軸の目盛の位置とラベルの設定
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
            "μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"]);

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

In [9]:
# 各領域のパーセントを表す文字列を配列として作成
# リスト 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 [10]:
p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill=Fill, show=False);

fig, ax = p.fig, p.ax

xsigma = np.linspace(-3, 3, 7)
fxsigma = lambdify(x, f(x))(xsigma)
ax.vlines(x=xsigma, ymin=0, ymax=fxsigma, color='white');
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
            "μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"]);

# 領域にラベルをつける
ax.text(-0.5,0.20, sig[0], ha = 'center', color="white")
ax.text( 0.5,0.20, sig[0], ha = 'center', color="white")
ax.text(-1.5,0.04, sig[1], ha = 'center', color="white")
ax.text( 1.5,0.04, sig[1], ha = 'center', color="white")
ax.text(-2.5,0.05, sig[2], ha = 'center', color="black")
ax.text( 2.5,0.05, sig[2], ha = 'center', color="black")
ax.text(-3.5,0.02, sig[3], ha = 'center', color="black")
ax.text( 3.5,0.02, sig[3], ha = 'center', color="black");

枠と縦軸目盛を非表示に

In [11]:
p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill=Fill, show=False);

fig, ax = p.fig, p.ax

xsigma = np.linspace(-3, 3, 7)
fxsigma = lambdify(x, f(x))(xsigma)
ax.vlines(x=xsigma, ymin=0, ymax=fxsigma, color='white');
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
            "μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"]);

ax.text(-0.5,0.20, sig[0], ha = 'center', color="white")
ax.text( 0.5,0.20, sig[0], ha = 'center', color="white")
ax.text(-1.5,0.04, sig[1], ha = 'center', color="white")
ax.text( 1.5,0.04, sig[1], ha = 'center', color="white")
ax.text(-2.5,0.05, sig[2], ha = 'center', color="black")
ax.text( 2.5,0.05, sig[2], ha = 'center', color="black")
ax.text(-3.5,0.02, sig[3], ha = 'center', color="black")
ax.text( 3.5,0.02, sig[3], ha = 'center', color="black");

# 枠と縦軸目盛を非表示に
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False);

-nσ から +nσ までを塗りつぶす

In [12]:
# 各領域のパーセントを表す文字列を配列として作成
# リスト 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 [13]:
# -1σ から +1σ までを塗りつぶす
x0 = np.linspace(-1, 1, 101)
f0 = lambdify(x, f(x))(x0)

p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill={'x':x0, 'y1':f0, 'color':c[0]}, show=False);

fig, ax = p.fig, p.ax

ax.set_xticks([i for i in range(-4, 5)])
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["", "", "", "μ-σ", "μ",
            "", "", "", "μ+σ"]);

# 枠と縦軸目盛を非表示に
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False);

# 領域にラベルをつける
ax.text(0,0.15, sig2[0], ha = 'center', color="white");

In [14]:
# -2σ から +2σ までを塗りつぶす
x1 = np.linspace(-2, 2, 101)
f1 = lambdify(x, f(x))(x1)

p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill={'x':x1, 'y1':f1, 'color':c[1]}, show=False);

fig, ax = p.fig, p.ax

ax.set_xticks([i for i in range(-4, 5)])
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["", "", "μ-2σ", "μ-σ", "μ",
            "", "", "μ+2σ", "μ+σ"]);

# 枠と縦軸目盛を非表示に
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False);

# 領域にラベルをつける
ax.text(0,0.15, sig2[1], ha = 'center', color="white");

In [15]:
# -3σ から +3σ までを塗りつぶす
x2 = np.linspace(-3, 3, 101)
f2 = lambdify(x, f(x))(x2)

p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill={'x':x2, 'y1':f2, 'color':c[2]}, show=False);

fig, ax = p.fig, p.ax

ax.set_xticks([i for i in range(-4, 5)])
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["", "μ-3σ", "μ-2σ", "μ-σ", "μ",
            "", "μ+3σ", "μ+2σ", "μ+σ"]);

# 枠と縦軸目盛を非表示に
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False);

# 領域にラベルをつける
ax.text(0,0.15, sig2[2], ha = 'center', color="black");

In [16]:
# -4σ から +4σ までを塗りつぶす
x3 = np.linspace(-4, 4, 101)
f3 = lambdify(x, f(x))(x3)

p = plot(f(x), (x, -4, 4), line_color = '#0000ff', 
         grid = False, legend = False,
         xlim = (-4, 4), ylim = (0, 0.41), size = (6.4, 2.4), 
         xlabel = "", ylabel = "",
         fill={'x':x3, 'y1':f3, 'color':c[3]}, show=False);

fig, ax = p.fig, p.ax

ax.set_xticks([i for i in range(-4, 5)])
ax.set_xticks([-4, -3, -2, -1, 0, 
                4,  3,  2,  1])
ax.set_xticklabels(
  labels = ["μ-4σ", "μ-3σ", "μ-2σ", "μ-σ", "μ",
            "μ+4σ", "μ+3σ", "μ+2σ", "μ+σ"]);

# 枠と縦軸目盛を非表示に
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.axes.yaxis.set_visible(False);

# 領域にラベルをつける
ax.text(0,0.15, sig2[3], ha = 'center', color="black");