一様乱数から正規乱数を生成する（ボックス=ミュラー法）

この記事では ボックス=ミュラー法 というものを紹介します。
Pythonで確認してから、数学的な証明を行います。

証明では前回の記事で紹介した変数変換を使用します。

ボックス=ミュラー法（Box-Muller法）

• ボックス=ミュラー法
  $[0,1]$上の一様分布上の独立な確率変数$U_1, U_2$があるとき、

  $$X = \sqrt{-2 \log{U_1}} \cos{2 \pi U_2} \\ Y = \sqrt{-2 \log{U_1}} \sin{2 \pi U_2}$$

  とすると$X,Y$は独立に標準正規分布${\it N}(0, 1)$に従う。

簡単に言うと 2つの独立な一様分布から2つの独立な標準正規分布を生成できる というものです。

これがあればもしPythonでnp.random.normal()が使えなくなってもnp.random.uniform()さえあれば正規分布に従う乱数を作れるということです。すごいですね！

Pythonで確認してみる

一様分布に上記の変換を施すと本当に正規分布になるのか確認してみましょう。

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

plt.style.use("seaborn")
font_cfg = {"fontsize": 14}


def box_muller(u_1, u_2):
    """ボックス=ミュラー法を適用する関数"""
    X = np.sqrt(-2 * np.log(u_1)) * np.cos(2 * np.pi * u_2)
    Y = np.sqrt(-2 * np.log(u_1)) * np.sin(2 * np.pi * u_2)  
    return X, Y


# 2つの独立な[0, 1]上の一様乱数を生成
U_1, U_2 = np.random.uniform(0, 1, (2, 100_000))

# ボックス=ミュラー法を適用
X, Y = box_muller(U_1, U_2)


# Plot
fig, axes = plt.subplots(2, 2, figsize=(9, 7), dpi=90)

# ヒストグラムの描画
for ax, d, s in zip(axes.ravel(), [U_1, U_2, X, Y], ["U_1", "U_2", "X", "Y"]):
    ax.hist(d, bins=100, density=True)
    ax.set_title(f"${s}$", **font_cfg)

# 一様分布と標準正規分布の確率密度の描画
arr = np.linspace(0, 1)
for ax in axes[0, :]:
    ax.plot(arr, stats.uniform.pdf(arr), c="r", label="$U(0, 1)$")
    ax.set_ylim(0, 1.3)
    ax.legend(**font_cfg)
    
arr = np.linspace(-4, 4)
for ax in axes[1, :]:
    ax.plot(arr, stats.norm.pdf(arr), c="r", label="$N(0, 1)$")
    ax.legend(**font_cfg)
    
plt.tight_layout()

$X,Y$は標準正規分布に従っていそうですね。
ここでは確認しませんが$X,Y$は独立です。（これは数式での証明で示します。）

これ、すごいですよね。対数やら三角関数を使うことで一様乱数が正規乱数に早変わりです。
初めて知ったときは驚きました。

Pythonで書くとこれだけの話なのですが、数学的に証明しようとすると少し大変です。
以降で実際に証明します。

ボックス=ミュラー法の証明

準備

証明に使う道具を簡単に紹介します。

• 2次元確率分布の変数変換
  前回記事で紹介した変数変換を使用します。
  　
  確率変数$\ (X, Y)$の同時確率（密度）関数を$\ f_{X,Y}(x, y)$ とし、
  $X = h_1(S, T), Y = h_2(S, T)$と表されるとき、
  $(S, T)$の同時確率（密度）関数$f_{S,T}(s, t)$は

  $$f_{S,T}(s,t) = f_{X,Y}(h_1(s, t), h_2(s, t)) |J(s,t)|$$

  となる、というものです。
  詳しくは前回の記事をご覧ください。
  
  

• 一様分布の確率密度関数
  $[a,b]$上の一様分布の確率密度関数は下記です。

  $$f_x(x) = \begin{cases} 1 / (b - a) & (a \leq x \leq b) \\ 0 & ({\it otherwise}) \end{cases}$$

  今回は$[0, 1]$上の一様分布を変数変換していきます。
  
  

• 標準正規分布の確率密度関数
  標準正規分布の確率密度関数は下記です。

  $$f_x(x) = \frac{1}{\sqrt{2 \pi}} e^{-x^2 / 2}$$

  変数変換の結果、この式を2つ掛けたものが得られれば証明となります。

---

計算

問題を確認します。

$[0,1]$上の一様分布からの独立な確率変数$U_1, U_2$があるとき、

$$X = \sqrt{-2 \log{U_1}} \cos{2 \pi U_2} \\ Y = \sqrt{-2 \log{U_1}} \sin{2 \pi U_2}$$

なる変数を施すと$X,Y$は独立に標準正規分布${\it N}(0, 1)$に従うことを示せ。

計算しやすくするために次のように2段階の変数変換に分けて考えます。

• ステップ1
  $U_1,U_2$ を

  $$r = \sqrt{-2 \log{U_1}} \\ \theta = 2 \pi U_2$$

  と変数変換する。

• ステップ2
  先に求めた$r,\theta$ を

  $$X = r \cos{\theta} \\ Y = r \sin{\theta}$$

  と変数変換する。

---

ステップ1

$U_1,U_2$ を

$$r = \sqrt{-2 \log{U_1}} \\ \theta = 2 \pi U_2$$

と変数変換する。

$U_1,U_2$ を $r,\theta$ で表すと、

$$U_1 = h_1(r) = e^{-r^2 /2} \\ U_2 = h_2(\theta) = \frac{1}{2 \pi} \theta$$

となります。よってヤコビアンは

$$\begin{aligned} J(r, \theta) &= det \left( \begin {array} {cc} (\partial/\partial r) e^{-r^2 / 2} & (\partial/\partial \theta) e^{-r^2 / 2} \\ (\partial/\partial r) \frac{1}{2 \pi} \theta & (\partial/\partial \theta) \frac{1}{2 \pi} \theta \end {array} \right) \\ &= det \left( \begin {array} {cc} -r e^{-r^2 / 2} & 0 \\ 0 & \frac{1}{2 \pi} \end {array} \right) \\ &= -r e^{r^2 / 2} \frac{1}{2 \pi} \end{aligned}$$

となります。これを変数変換の公式に入れて

$$\begin{aligned} f_{r,\theta}(r,\theta) &= f_{U_1,U_2}(h_1(r), h_2(\theta)) \ |J(r,\theta)| \\ &= f_{U_1}(h_1(r))f_{U_2}(h_2(\theta)) \ |J(r, \theta)| \\ &= f_{U_1}(e^{-r^2 /2}) f_{U_2}(\frac{1}{2 \pi} \theta) |-r e^{r^2 / 2} \frac{1}{2 \pi}| \\ &= \frac{1}{1-0} \frac{1}{1-0} r e^{r^2 / 2} \frac{1}{2 \pi} \\ &= r e^{r^2 / 2} \frac{1}{2 \pi} \\ \end{aligned}$$

となります。
（ややこしいですがヤコビアンの絶対値が残っただけです。）


これは $f_{r, \theta}(r, \theta) = f_r(r) f_{\theta}(\theta)$ の形なので $r,\theta$は独立でそれぞれ

$$f_r(r) = r e^{r^2 / 2} \quad 0 \leq r \leq \infty \\ f_{\theta}(\theta) = 1 / 2\pi \quad 0 \leq \theta \leq 2\pi$$

に従います。（$\theta$ の分布は $[0, 2\pi]$上の一様分布）

info:
$0 \leq r \leq \infty$ は $f_{U_1}(e^{-r^2 /2}) \neq 0$ より、$0 \leq re^{-r^2 /2} \leq 1$ となるので、
これを変形することで得られます。
同様にして$0 \leq x \leq 2\pi$ も得られます。

---

ステップ2

ステップ1で得られた $r, \theta$ を使って次の変数変換を行います。

先に求めた$r,\theta$ を

$$X = r \cos{\theta} \\ Y = r \sin{\theta}$$

と変数変換する。

まず、上の式を$r, \theta$ を$X,Y$ で表します。
$r$は

$$x^2 + y^2 = r^2 (\cos^2\theta + \cos^2\theta) = r^2$$

より

$$r = h_1(x,y) = \sqrt{x^2 + y^2}$$

です。$\theta$は

$$\frac{y}{x} = \frac{r\sin\theta}{r\cos\theta} = \tan\theta$$

より

$$\theta = h_2(x,y) = \tan^{-1}\left(\frac{y}{x} \right)$$

となります。

次に $r, \theta$ を使ってヤコビアンを求めます。
$(\tan^{-1}x)' = 1 / (1 + x^2)$ に注意すると

$$J(x,y) = det \left( \begin{array} {cc} x/\sqrt{x^2 + y^2} & y/\sqrt{x^2 + y^2} \\ -y / (x^2 + y^2) & x / (x^2 + y^2) \\ \end{array} \right) = \frac{1}{\sqrt{x^2 + y^2}}$$

となります。

最後に変数変換の公式にこれらを代入すると

$$\begin{aligned} f_{X,Y}(x,y) &= f_{r,\theta}(h_1(x,y),h_2(x,y)) |J(x,y)| \\ &= f_{r}(h_1(x,y)) f_{\theta}(h_2(x,y)) |J(x,y)| \\ &= f_{r}(\sqrt{x^2 + y^2}) f_{\theta}(\tan^{-1}(y/x)) |\frac{1}{\sqrt{x^2 + y^2}}| \\ &= \sqrt{x^2+y^2} e^{-(x^2+y^2)/2} \frac{1}{2\pi} \frac{1}{\sqrt{x^2 + y^2}} \\ &= \frac{1}{2\pi} e^{-(x^2+y^2)/2} \\ &= \frac{1}{\sqrt{2\pi}} e^{-x^2/2} \times \frac{1}{\sqrt{2\pi}} e^{-y^2/2} \\ \end{aligned}$$

となり、2つの標準正規分布のかけ算が得られました。


以上で $X,Y$が独立に標準正規分布に従うことが確認でき、ボックス=ミュラー法を数学的に証明することができました！

---

まとめ

この記事ではボックス=ミュラー法をPythonで確認し、数学的な証明も行いました。

統計検定1級取得のため勉強をしている中で出てきて面白かったことと、同じような勉強をしている方の参考になればと思い記事を書きました。

ちなみに統計検定1級の参考書としてよく使用される「現代数理統計学の基礎」の4章 演習問題 問8 もこのボックス=ミュラー法の証明です。

---

おまけ

「ステップ1」の計算の結果、 $r,\theta$は独立でそれぞれ

$$f_r(r) = r e^{r^2 / 2} \quad 0 \leq r \leq \infty \\ f_{\theta}(\theta) = 1 / 2\pi \quad 0 \leq \theta \leq 2\pi$$

に従うことが確認できました。

前述のように、$\theta$ の分布は$[0, 2\pi]$上の一様分布です。


一方、$r$ の方はこのままでは一般的な確率分布ではなさそうですが、$V=r^2$ なる変換を施すと$V$ は指数分布 ${\it Ex}(1/2)$ に従います。

つまり $r$ は指数分布 ${\it Ex}(1/2)$ 平方根を取ったものに従い

$$r \sim \sqrt{{\it Ex}(1/2)}$$

となります。

なので、ボックス=ミュラー法は下記のような変換を行っていると解釈できます。

1. 2つの独立な $[0,1]$上の一様分布から、指数分布と $[0,2\pi]$ 上の一様分布を生成する。
2. 指数分布の平方根と $[0,2\pi]$ 上の一様分布から2つの独立な標準正規分布を生成する。

---

参考文献

現代数理統計学の基礎 (久保川達也)