正規分布

一番簡単な乱数は、ある範囲の数がまんべんなく出る一様乱数(uniform random numbers)です。乱数のところで説明したように、次のようにすると $0 \leqq x < 1$ の範囲の一様乱数を生成できます:

import numpy as np

rng = np.random.default_rng(12345)  # 12345は適当なタネ
rng.random()    # デフォルトは1個
0.22733602246716966
rng.random(10)  # 10個
array([0.31675834, 0.79736546, 0.67625467, 0.39110955, 0.33281393,
       0.59830875, 0.18673419, 0.67275604, 0.94180287, 0.24824571])

もっとたくさんの数、例えば百万個でやってみましょう。画面に書き出すとたいへんなので、ヒストグラムを描いてみましょう。

N = 1000000  # 百万

plt.hist(rng.random(N), edgecolor="black", bins=20, density=True)
一様乱数

ヒストグラム(度数分布図)を描く plt.hist() を使いました。オプション density = True は縦軸の目盛りを個数でなく密度(density)にすることを意味します。密度とは、棒グラフの棒の面積の和が 1 になるように付けた目盛りです。下で説明する確率密度関数の縦軸と同じものになります。

このヒストグラムを関数で表したものを確率密度関数(probability density function、略して p.d.f.)または単に密度関数といいます。数学的にいえば、任意の $a < b$ について、$a \leqq x < b$ の範囲に乱数が入る確率が $\int_a^b f(x) dx$ になるように選ばれた関数 $f(x)$ が確率密度関数です。rng.random() の確率密度関数は

\[ f(x) = \begin{cases} 1 & (0 \leqq x < 1) \\ 0 & (\text{$x < 0$, $1 \leqq x$}) \end{cases} \]

と表すことができます。

この rng.random() を二つ加えたらどうなるでしょうか。2倍するのではなく、独立な乱数を二つ加えるのです。

plt.hist(rng.random(N) + rng.random(N),
         bins=20, edgecolor="black", density=True)
一様乱数2個の和

三個ならどうなるでしょうか。

plt.hist(rng.random(N) + rng.random(N) + rng.random(N),
         bins=20, edgecolor="black", density=True)
一様乱数3個の和

だんだん釣り鐘の形になってきました。

rng.random() は平均値が 0.5 ですので、0.5 を引いて範囲を $-0.5 \leqq x < 0.5$ にして考えましょう。この分散は

\[ \int_{-0.5}^{0.5} x^2 dx = \frac{1}{12} \]

ですので、12個加えるとちょうど分散が 1 になります。

次の図が、rng.random() - 0.5 を12個加えた分布のヒストグラムです。オレンジの曲線は後で説明します。

plt.hist(rng.random(N) + rng.random(N) + rng.random(N) +
         rng.random(N) + rng.random(N) + rng.random(N) +
         rng.random(N) + rng.random(N) + rng.random(N) +
         rng.random(N) + rng.random(N) + rng.random(N) - 6,
         bins=np.arange(-4, 4.2, 0.2), edgecolor="black", density=True)

x = np.arange(-4, 4.1, 0.1)
plt.plot(x, np.exp(-x**2/2) / np.sqrt(2*np.pi))
一様乱数12個の和

オレンジの曲線は

\[ f(x) = \frac{1}{\sqrt{2\pi}} \exp\biggl(-\frac{x^2}{2}\biggr) \]

を描いたものです。ここで $\exp(z) = e^z$ は指数関数です。$e = 2.718\ldots$ は自然対数の底(てい)です。$\int_{-\infty}^{\infty} \exp(-x^2\!/2) dx = \sqrt{2\pi}$ ですので、これは確率密度関数です。この分布を標準正規分布(standard normal distribution)といいます。

乱数を足し合わせる上の実験は、一様分布の乱数 rng.random() から出発しましたが、(分散が有限な)どんな分布から出発しても、まったく同じことが成り立ちます。

数学的に言えば、平均 $\mu$、分散 $\sigma^2$ の確率変数 $X$ がどんな分布であっても、そこから引き出した数 $X_1$, $X_2$, ..., $X_n$ の平均値

\[ \bar{X} = \dfrac{X_1 + X_2 + \cdots + X_n}{n} \]

の分布は、平均値が元と同じ $\mu$ で、分散が $\sigma^2/n$ になりますので、

\[ \dfrac{\bar{X} - \mu}{\sqrt{\sigma^2 \! / n}} \]

の分布は平均 0、分散 1 になります。ここまでは $n$ の値にかかわらず言えることです。ここで、$n$ が十分大きくなると、この分布は標準正規分布に近づきます。このことを中心極限定理といいます。

$X$ が標準正規分布ならば、$\sigma X + \mu$ は平均 $\mu$、分散 $\sigma^2$ の正規分布(normal distribution)になります。平均 $\mu$、分散 $\sigma^2$ の正規分布を $N(\mu, \sigma^2)$ と表します。標準正規分布は平均 0、分散 1 の正規分布ですので $N(0, 1)$ と表されます。

なお、正規分布を研究した数学者ガウス(Carl Friedrich Gauss、1777-1855)に敬意を表して、正規分布のことをガウス分布(Gaussian distribution)ともいいます。

変数がある範囲に入る確率を求めるには、密度関数を積分しなければなりません。そこで、あらかじめ密度関数 $f(x)$ を積分した

\[ F(q) = \int_{-\infty}^q f(x) dx \]

を求めておくと便利です。この $F(q)$ を累積分布関数(cumulative distribution function、略して CDF)あるいは単に分布関数と呼びます。$F(q)$ がわかっていれば、確率変数 $X$ が $a \leqq X < b$ の範囲に入る確率は $F(b) - F(a)$ で求められます。

分布関数 $p = F(q)$ の逆関数 $q = F^{-1}(p)$ を分位関数または分位点関数(quantile function, percent point function)と呼びます。

Python の SciPy ライブラリから from scipy import stats でインポートした stats には、標準正規分布 $N(0,1)$ の

があります。また、正規分布の乱数(正規乱数)は NumPy で rng.normal(平均, 標準偏差, 個数) で生成できます。

物理現象の測定では、測定誤差がほぼ正規分布することがよくあります。これに対して、他の多くの分野では、観測値そのものが正規分布することはまずありませんが、中心極限定理のおかげで、観測値の平均値はほぼ正規分布します。

正規分布 $N(\mu, \sigma^2)$ に従う変数 $X$ が平均±標準偏差の範囲内に入る確率、つまり $\mu - \sigma < X < \mu + \sigma$ になる確率は、標準正規分布 $N(0,1)$ に従う変数 $Z$ が $-1 < Z < 1$ になる確率に等しく、Python で次のようにして求められます。

stats.norm.cdf(1) - stats.norm.cdf(-1)
0.6826894921370859

正規分布が左右対称であることを使えば、次のようにするほうが簡単です。

1 - 2 * stats.norm.cdf(-1)
0.6826894921370859

つまり、ほぼ68%の確率で正規分布は $\mu \pm \sigma$ の範囲に入ります。

同様にして、$\mu \pm 2\sigma$、$\mu \pm 3\sigma$ に入る確率を求めてください。

逆に、$\mu \pm k\sigma$ の範囲に入る確率が与えられたときに $k$ の値を求めるには stats.norm.ppf() を使います。例えば 95% の場合は

stats.norm.ppf(0.025)
-1.9599639845400545

ですので、ほぼ $\mu \pm 1.96\sigma$ です。ざっくり $\mu \pm 2\sigma$ で 95% と覚えておくと便利です。

Scrapboxのほうに書いた正規分布の何がいいの?もご参照ください。