「gnuplot で正規分布をσごとに塗りわける」の 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])
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])
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");
