正規分布関数

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

Maxima には組み込み関数として，誤差関数 erf(x) があります。

? erf;

-- Function: erf (<z>)
     The Error Function erf(z) (A&S 7.1.1)
     See also flag 'erfflag'.

In [1]:

f(x):= 1/sqrt(2*%pi) * exp(-x**2/2);
F(x1, x2):= float(1/2 * (erf(x2/sqrt(2)) - erf(x1/sqrt(2))));

Out[1]:

\[\tag{${\it \%o}_{1}$}f\left(x\right):=\frac{1}{\sqrt{2\,\pi}}\,\exp \left(\frac{-x^2}{2}\right)\]

Out[1]:

\[\tag{${\it \%o}_{2}$}F\left(x_{1} , x_{2}\right):={\it float}\left(\frac{1}{2}\,\left(\mathrm{erf}\left(\frac{x_{2}}{\sqrt{2}}\right)-\mathrm{erf}\left(\frac{x_{1}}{\sqrt{2}}\right)\right)\right)\]

標準正規分布を draw2d()

In [2]:

draw2d(
  dimensions=[640,240],
  xrange = [-4, 4], yrange = [0, 0.4],
  
  line_width = 2, color = blue,
  explicit(f(x), x, -4, 4)
)$

σごとに塗りわける

In [3]:

/* 塗りわける色 */
c: ["#0000ff","#1e90ff","#87ceeb","#add8e6"]$

In [4]:

draw2d(
  dimensions=[640,240],
  xrange = [-4, 4], yrange = [0, 0.4],

  /* σごとに塗りわけ; 塗りつぶし開始 */
  filled_func = 0, 
  makelist([fill_color = c[i], 
            explicit(f(x), x, -i, -i+1), 
            explicit(f(x), x, i-1, i)], i, 1, 3),
  /* 塗りつぶし終了 */
  filled_func = false, 
  
  line_width = 2, color = blue,
  explicit(f(x), x, -4, 4)
)$

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

In [5]:

draw2d(
  dimensions=[640,240],
  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  makelist([fill_color=c[i], 
            explicit(f(x), x,-i, -i+1), 
            explicit(f(x), x, i-1,i)], i, 1, 3),
  filled_func = false, 
  
  /* 縦線を引く; 端点には点を描かない */
  points_joined = true, point_size = 0, 
  line_width = 2, color = white, 
  makelist(points([[-i, 0], [-i, f(-i)]]), i, -3, 3),

  color = blue,
  explicit(f(x), x, -4, 4)
)$

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

In [6]:

draw2d(
  dimensions=[640,240],
  font = "Times", font_size = 16,
  
  /* 横軸の目盛ラベルを設定 */
  xtics = {["μ-4σ",-4], ["μ-3σ",-3], ["μ-2σ",-2], ["μ-σ",-1], ["μ",0], 
           ["μ+4σ", 4], ["μ+3σ", 3], ["μ+2σ", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  makelist([fill_color=c[i], 
            explicit(f(x), x,-i, -i+1), 
            explicit(f(x), x, i-1,i)], i, 1, 3),
  filled_func = false, 
  
  points_joined = true, point_size = 0, 
  line_width = 2, color = white, 
  makelist(points([[-i, 0], [-i, f(-i)]]), i, -3, 3),
  
  color = blue,
  explicit(f(x), x, -4, 4)
)$

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

In [7]:

/* 各領域のパーセントを表す文字列を配列として作成 */
for i:1 thru 4 do(
  sig[i]: printf(false, "~5,2f％", F(i-1,i)*100), 
  print(sig[i])
)$

34.13％

13.59％

2.14％

0.13％

In [8]:

draw2d(
  dimensions=[640,240],
  font = "Times", font_size = 16,
  xtics = {["μ-4σ",-4], ["μ-3σ",-3], ["μ-2σ",-2], ["μ-σ",-1], ["μ",0], 
           ["μ+4σ", 4], ["μ+3σ", 3], ["μ+2σ", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  makelist([fill_color=c[i], 
            explicit(f(x), x,-i, -i+1), 
            explicit(f(x), x, i-1,i)], i, 1, 3),
  filled_func = false, 
  
  points_joined = true, point_size = 0, 
  line_width = 2, color = white, 
  makelist(points([[-i, 0], [-i, f(-i)]]), i, -3, 3),
  
  /* 領域にラベルをつける*/
  color = white,
  label([sig[1], -0.5,0.20]), label([sig[1], 0.5,0.20]),
  label([sig[2], -1.5,0.05]), label([sig[2], 1.5,0.05]),
  color = black,
  label([sig[3], -2.5,0.06]), label([sig[3], 2.5,0.06]),
  label([sig[4], -3.5,0.02]), label([sig[4], 3.5,0.02]),
  
  color = blue,
  explicit(f(x), x, -4, 4)
)$

枠と縦軸目盛を非表示に

In [9]:

draw2d(
  /* user_preamble で直接 gnuplot 的設定 */
  ytics = false, user_preamble = "unset border;set xtics out;",
  xaxis = true, xaxis_type=solid, xaxis_width = 2, 
  
  dimensions=[640,240],
  font = "Times", font_size = 16,
  xtics = {["μ-4σ",-4], ["μ-3σ",-3], ["μ-2σ",-2], ["μ-σ",-1], ["μ",0], 
           ["μ+4σ", 4], ["μ+3σ", 3], ["μ+2σ", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  makelist([fill_color=c[i], 
            explicit(f(x), x,-i, -i+1), 
            explicit(f(x), x, i-1,i)], i, 1, 3),
  filled_func = false, 
  
  points_joined = true, point_size = 0, 
  line_width = 2, color = white, 
  makelist(points([[-i, 0], [-i, f(-i)]]), i, -3, 3),
  
  color = white,
  label([sig[1], -0.5,0.20]), label([sig[1], 0.5,0.20]),
  label([sig[2], -1.5,0.05]), label([sig[2], 1.5,0.05]),
  color = black,
  label([sig[3], -2.5,0.06]), label([sig[3], 2.5,0.06]),
  label([sig[4], -3.5,0.02]), label([sig[4], 3.5,0.02]),
  
  color = blue,
  explicit(f(x), x, -4, 4)
)$

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

In [10]:

/* 各領域のパーセントを表す文字列を配列として作成 */
for i:1 thru 4 do(
  sig2[i]: printf(false, "~5,2f％", F(-i,i)*100), 
  print(sig2[i])
)$

68.27％

95.45％

99.73％

99.99％

In [11]:

draw2d(
  dimensions=[640,240],
  ytics = false, 
  user_preamble = "unset border;set xtics out;",
  xaxis = true, xaxis_type=solid,
  
  font = "Times", font_size = 16,
  xtics = {["",-4], ["",-3], ["",-2], ["μ-σ",-1], ["μ",0], 
           ["", 4], ["", 3], ["", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  fill_color = "#0000ff", 
  explicit(f(x), x, -1, 1), 
  filled_func = false, 
    
  color = white,
  label([sig2[1], 0, 0.15]), 
  
  color = blue, line_width = 2, 
  explicit(f(x), x, -1, 1), 
  color = gray, 
  explicit(f(x), x, -4, -1), explicit(f(x), x, 1, 4)
)$

In [12]:

draw2d(
  dimensions=[640,240],
  ytics = false, 
  user_preamble = "unset border;set xtics out;",
  xaxis = true, xaxis_type=solid,
  font = "Times", font_size = 16,
  xtics = {["",-4], ["",-3], ["μ-2σ",-2], ["μ-σ",-1], ["μ",0], 
           ["", 4], ["", 3], ["μ+2σ", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  fill_color = "#1e90ff", 
  explicit(f(x), x, -2, 2), 
  filled_func = false, 
    
  color = white,
  label([sig2[2], 0,0.15]), 
  
  color = blue, line_width = 2, 
  explicit(f(x), x, -2, 2), 
  color = gray, 
  explicit(f(x), x, -4, -2), explicit(f(x), x, 2, 4)
)$

In [13]:

draw2d(
  dimensions=[640,240],
  ytics = false, 
  user_preamble = "unset border;set xtics out;",
  xaxis = true, xaxis_type=solid,
  font = "Times", font_size = 16,
  xtics = {["",-4], ["μ-3σ",-3], ["μ-2σ",-2], ["μ-σ",-1], ["μ",0], 
           ["", 4], ["μ+3σ", 3], ["μ+2σ", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  fill_color = "#87ceeb", 
  explicit(f(x), x, -3, 3), 
  filled_func = false, 
    
  color = black,
  label([sig2[3], 0,0.15]), 
  
  color = blue, line_width = 2, 
  explicit(f(x), x, -3, 3), 
  color = gray, 
  explicit(f(x), x, -4, -3), explicit(f(x), x, 3, 4)
)$

In [14]:

draw2d(
  dimensions=[640,240],
  ytics = false, 
  user_preamble = "unset border;set xtics out;",
  xaxis = true, xaxis_type=solid,
  font = "Times", font_size = 16,
  xtics = {["μ-4σ",-4], ["μ-3σ",-3], ["μ-2σ",-2], ["μ-σ",-1], ["μ",0], 
           ["μ+4σ", 4], ["μ+3σ", 3], ["μ+2σ", 2], ["μ+σ", 1]},

  xrange = [-4, 4], yrange = [0, 0.4],
  
  filled_func = 0, 
  fill_color = "#add8e6", 
  explicit(f(x), x, -4, 4), 
  filled_func = false, 
    
  color = black,
  label([sig2[4], 0,0.15]), 
  
  color = blue, line_width = 2, 
  explicit(f(x), x, -4, 4)
)$

Maxima で正規分布をσごとに塗りわける