正規分布に従う乱数を一様分布の乱数から生成する方法を
一様分布、正規分布、指数分布、ポアソン分布の乱数を生成する方法 (Scala) - Qiita
で書いた。

この方法で生成する乱数が正規分布に従うことを計算して示す。

（数式をいっぱい書く投稿は約6年ぶりである。→ 不偏分散はなぜ n – 1 で割るのか？）

一様分布の乱数から正規分布に従う乱数を生成する方法

一様分布の乱数から正規分布に従う乱数を生成する方法は次の通り。

$a$, $b$ を $0$ から $1$ までの一様な乱数とすると、次の式で計算した $x$, $y$ が正規分布に従った2つの乱数になる。 $x$, $y$ には相関はなく独立である。

\begin{align}
x =& \sqrt{-2\log{a}} \sin{2\pi b} \\
y =& \sqrt{-2\log{a}} \cos{2\pi b} \\
\end{align}

この方法はBox–Muller法というらしい。

対数関数や三角関数を組み合わせたこの式が正規分布の式とはほど遠く、なぜこの式で正規分布に従うのかがまったく想像つかなかったし、2つがペアで生成されるのに、その2つは独立だというのもすぐには理解できない。それを理解したかったのがこの記事のきっかけである。

厳密性は欠くかもしれないが、個人的に十分に納得感が得られる程度にかみ砕いている。

正規分布の乱数となる数学的な説明

正規分布に限らず一般に確率密度関数の変数を変換する方法を示してから、最後に上記式をあてはめてみて、正規分布の確率密度関数と一致することを示す。

確率密度関数を変数変換

確率変数 $x$ に対して確率密度関数を $p_x$ と書くことにする。

$p_x$ に微小な幅 $\Delta x$ を乗じた値は確率変数が $x$ から $x+\Delta x$ までの間を取る確率を表す。

ところで $x$ は別の変数 $a$ と関数 $f$ を使って次のように書けるものとする。

\begin{align}
x =& f(a) \\
x + \Delta x =& f(a + \Delta a) \\
\end{align}

以下の議論をシンプルにするため関数 $f$ が単調増加または単調減少だと仮定する。

$a$ を確率変数とみなしたときの確率密度関数 $p_a$ は $p_x$ とは次のような関係になる。

\begin{equation}
p_x\Delta x = p_a\Delta a \tag{1}
\end{equation}

右辺は確率変数が $a$ から $a+\Delta a$ までの間を取る確率を表す。

これは $x$ が $x$ から $x+\Delta x$ までの間を取る確率に等しいはずである。これは左辺の表現である。

$\Delta x$ は $\Delta a$ と $f$ の微分を使って次のように書ける。

\begin{align}
\Delta x =& (x + \Delta x) - x \\
=& f(a + \Delta a) - f(a) \\
=& f' \Delta a \\
=& \frac{dx}{da} \Delta a \tag{2} \\
\end{align}

(1)を変形して(2)を代入すると、

\begin{align}
p_x =& \frac{p_a\Delta a}{\Delta x} \\
=& \frac{p_a\Delta a}{\frac{dx}{da} \Delta a} \\
=& \frac{p_a}{\frac{dx}{da}} \tag{3} \\
\end{align}

2変数の確率密度関数を変数変換

次にこれを2変数に拡張した場合を考察する。

確率変数 $x$, $y$ に対して確率密度関数を $p_{xy}$ と書くことにする。

$p_{xy}$ に微小な領域 $\Delta(xy)$ を乗じた値は確率変数 $x$, $y$ が特定の微小な領域に収まる確率を表す。

ところで $x$, $y$ は別の変数 $a$, $b$ と2変数関数 $f$, $g$ を使って次のように書けるものとする。

\begin{align}
x =& f(a, b) \\
y =& g(a, b) \\
\end{align}

$a$, $b$ を確率変数とみなしたときの確率密度関数 $p_{ab}$ は $p_{xy}$ とは次のような関係になる。$\Delta(xy)$ は関数 $f$, $g$ を使って $\Delta(ab)$ を変換した領域である。

\begin{equation}
p_{xy}\Delta(xy) = p_{ab}\Delta(ab) \tag{4}
\end{equation}

この式は1変数の(1)式を2変数に拡張したものである。

右辺は確率変数 $a$, $b$ が微小な領域 $\Delta(ab)$ に収まる確率を表す。左辺は $x$, $y$ が微小な領域 $\Delta(xy)$ に収まる確率を表す。2つは言い換えているだけであるので等しい。

ここで、 $\Delta(ab)$ は $a$ が $a$ から $a + \Delta a$ 、$b$ が $b$ から $b + \Delta b$ の長方形の領域を表すものとする。

$\Delta(ab)$ を囲む長方形は以下の4つの頂点からなる。

$(a, b)$
$(a, b + \Delta b)$
$(a + \Delta a, b)$
$(a + \Delta a, b + \Delta b)$

その面積は次のとおり。

\Delta(ab) = \Delta a \Delta b

このとき $\Delta(xy)$ は長方形になるとは限らず、一般には平行四辺形で近似できる。 $\Delta(xy)$ の4つの頂点は以下となる。

$\left( x, y \right)$
$\left( x + \frac{\partial x}{\partial a}\Delta a, y + \frac{\partial y}{\partial a}\Delta a \right) $
$\left( x + \frac{\partial x}{\partial b}\Delta b, y + \frac{\partial y}{\partial b}\Delta b \right) $
$\left( x + \frac{\partial x}{\partial a}\Delta a + \frac{\partial x}{\partial b}\Delta b, y + \frac{\partial y}{\partial a}\Delta a + \frac{\partial y}{\partial b}\Delta b \right) $

一般に、平行四辺形の面積は2つの辺のベクトルの外積の大きさに等しい。したがって、次の4つの頂点からなる平行四辺形

$(0, 0)$
$(x_1, y_1)$
$(x_2, y_2)$
$(x_1 + x_2, y_1 + y_2)$

の面積は $\left| \overrightarrow{(x_1, y_1)} \times \overrightarrow{(x_2, y_2)} \right|$ つまり $|x_1 y_2 - y_1 x_2|$ となる。

これを踏まえて、微小領域の面積 $\Delta(xy)$ は、次のように書ける。

\begin{align}
\Delta(xy) =& \left|{\frac{\partial x}{\partial a}\Delta a \frac{\partial y}{\partial b}\Delta b - \frac{\partial y}{\partial a}\Delta a \frac{\partial x}{\partial b}\Delta b } \right| \\
=& \left|{\frac{\partial x}{\partial a} \frac{\partial y}{\partial b} - \frac{\partial y}{\partial a} \frac{\partial x}{\partial b} } \right| \Delta a \Delta b \\
=& \left|{\frac{\partial x}{\partial a} \frac{\partial y}{\partial b} - \frac{\partial y}{\partial a} \frac{\partial x}{\partial b} } \right| \Delta(ab) \tag{5}  \\
\end{align}

この式は1変数の(2)式を2変数に拡張したものである。

ここまでは準備である。

本題

ここから正規分布についての考察をする。

確率変数 $x$, $y$ は別の確率変数 $a$, $b$ と2変数関数 $f$, $g$ で以下のように表せるとする。これは最初に示した乱数生成式である。

\begin{align}
x =& f(a, b) = \sqrt{-2\log{a}} \sin{2\pi b} \\
y =& g(a, b) = \sqrt{-2\log{a}} \cos{2\pi b} \\
\end{align}

$a$, $b$ が一様分布に従うとき、 $x$, $y$ が正規分布に従うことを示したい。そのためには確率密度関数 $p_{xy}$ を計算し、正規分布の確率密度関数と一致することを示すのを目標とする。

計算しやすくするため

z = -2\log{a}

とおくと、

\begin{align}
x =& \sqrt{z} \sin{2\pi b} \\
y =& \sqrt{z} \cos{2\pi b} \\
\end{align}

となる。すると、次の計算ができる。

\begin{align}
x^2+y^2 =& (\sqrt{z} \sin{2\pi b})^2 + (\sqrt{z} \cos{2\pi b})^2 \\
=& z\left((\sin{2\pi b})^2 + (\cos{2\pi b})^2\right) \\
=& z \\
=& -2\log{a} \\
-\frac{x^2+y^2}{2} =& \log{a} \\
a =& e^{-\frac{x^2+y^2}{2}} \tag{6} \\
\end{align}

この結果はあとで利用する。

次に、偏微分 $\frac{\partial x}{\partial a}$, $\frac{\partial y}{\partial a}$, $\frac{\partial x}{\partial b}$, $\frac{\partial y}{\partial b}$ を計算する。

$\frac{\partial x}{\partial a}$ は次のとおり。

\begin{align}
\frac{\partial x}{\partial a}  =& \frac{\partial}{\partial a}\sqrt{z} \sin{2\pi b} \\
=& (\sin{2\pi b}) \frac{\partial}{\partial a}\sqrt{z} \\
=& (\sin{2\pi b}) \frac{\partial z}{\partial a} \frac{\partial}{\partial z}\sqrt{z} \\
=& (\sin{2\pi b}) \frac{\partial}{\partial a}(-2\log{a}) \frac{\partial}{\partial z}\sqrt{z} \\
=& (\sin{2\pi b}) \frac{-2}{a} \frac{\partial}{\partial z}\sqrt{z} \\
=& (\sin{2\pi b}) \frac{-2}{a} \frac{1}{2\sqrt{z}} \\
=& -\frac{\sin{2\pi b}}{a\sqrt{z}} \tag{7} \\
\end{align}

$\frac{\partial y}{\partial a}$ は $\sin$ が $\cos$ に代わるのみである。

\begin{align}
\frac{\partial y}{\partial a} =& -\frac{\cos{2\pi b}}{a\sqrt{z}} \tag{8} \\
\end{align}

$\frac{\partial x}{\partial b}$ は次のとおり。

\begin{align}
\frac{\partial x}{\partial b} =& \frac{\partial}{\partial b}\sqrt{z} \sin{2\pi b} \\
=& \sqrt{z} \frac{\partial}{\partial b} \sin{2\pi b} \\
=& 2\pi\sqrt{z} \cos{2\pi b} \tag{9} \\
\end{align}

$\frac{\partial y}{\partial b}$ は同様にして次のとおり。

\begin{align}
\frac{\partial y}{\partial b} =& \frac{\partial}{\partial b}\sqrt{z} \cos{2\pi b} \\
=& \sqrt{z} \frac{\partial}{\partial b} \cos{2\pi b} \\
=& -2\pi\sqrt{z} \sin{2\pi b} \tag{10} \\
\end{align}

これらの偏微分を(5)に代入すると次のとおりになる。

\begin{align}
\frac{\Delta(xy)}{\Delta(ab)} =& \left|{\frac{\partial x}{\partial a} \frac{\partial y}{\partial b} - \frac{\partial y}{\partial a} \frac{\partial x}{\partial b} } \right| \\
=& \left| \left(-\frac{\sin{2\pi b}}{a\sqrt{z}}\right) \left(-2\pi\sqrt{z} \sin{2\pi b}\right) - \left(-\frac{\cos{2\pi b}}{a\sqrt{z}}\right) \left(2\pi\sqrt{z} \cos{2\pi b}\right) \right| \\
=& \left| \left(\frac{\sin{2\pi b}}{a\sqrt{z}}\right) \left(2\pi\sqrt{z} \sin{2\pi b}\right) + \left(\frac{\cos{2\pi b}}{a\sqrt{z}}\right) \left(2\pi\sqrt{z} \cos{2\pi b}\right) \right| \\
=& \left| \left(\frac{(\sin{2\pi b})^2}{a\sqrt{z}}\right) \left(2\pi\sqrt{z}\right) + \left(\frac{(\cos{2\pi b})^2}{a\sqrt{z}}\right) \left(2\pi\sqrt{z}\right) \right| \\
=& \left| \frac{2\pi}{a} \left((\sin{2\pi b})^2 + (\cos{2\pi b})^2\right) \right| \\
=& \left| \frac{2\pi}{a} \right| \tag{11} \\
\end{align}

この結果と(4)、(6)を利用すれば次のとおりとなる。

\begin{align}
p_{xy} =& p_{ab} \frac{\Delta(ab)}{\Delta(xy)} \\
=& p_{ab} \left| \frac{a}{2\pi} \right| \\
=& p_{ab} \frac{1}{2\pi} e^{-\frac{x^2+y^2}{2}} \\
\end{align}

$a$, $b$ が $0$ から $1$ までの一様分布に従うとすると、その範囲内において $p_{ab}$ は定数の $1$ になるので、 $p_{xy}$ は最終的に以下となる。

\begin{align}
p_{xy} =& \frac{1}{2\pi} e^{-\frac{x^2+y^2}{2}} \\
=& \frac{1}{2\pi} e^{-\frac{x^2}{2}} e^{-\frac{y^2}{2}} \\
=& \left( \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}} \right) \left( \frac{1}{\sqrt{2\pi}} e^{-\frac{y^2}{2}} \right) \\
\end{align}

この式は確率変数 $x$ の正規分布に従う確率密度関数と確率変数 $y$ の正規分布に従う確率密度関数の積と一致する。このことから $x$, $y$ は相関がなく独立で、それぞれが正規分布に従うことがわかる。

以上。

補足

$f$, $g$ は $b$ について単調増加でも単調減少でもないが、 $b$ を $0$ から $1$ ではなく4等分した区間、 $0$ から $\frac{1}{4}$ 、 $\frac{1}{4}$ から $\frac{1}{2}$ 、、、に限定すればそれぞれの区間内では単調になる。それぞれの区間に限定しても $x$, $y$ は同じ正規分布に従うので、単にそれを重ね合わせたものと考えればよい。