Box-Muller法による標準正規分布に従う乱数生成

概要
手法
1. サンプルの破棄を伴う場合
2. サンプルの破棄を伴わない場合
Pythonによる実装
1. サンプルの破棄を伴う場合
2. サンプルの破棄を伴わない場合
証明

概要

ボックス＝ミュラー法（Box-Muller's method）は、一様分布に従う確率変数から標準正規分布に従う確率変数を生成する手法である。この記事では、その手法と実装について紹介し、その後それらを導出する。

手法

サンプルの破棄を伴う場合

※こちらの手法を論ずる際には、確率変数に小文字を用いる。

標準一様分布から独立に $z_{1}$ と $z_{2}$ を生成し、$z_{1}\leftarrow 2z_{1}-1$, $z_{2}\leftarrow 2z_{2}-1$ の変換を行って、$(-1,1)$ の乱数列を得る。

そのうち、

$$z_{1}^{2}+z_{2}^{2}\leq 1\tag{1}$$

を満たさないものを破棄する。保持された $z_{1}$, $z_{2}$ に

$$y_{1}=z_{1}\left(\frac{-2\ln{r^{2}}}{r^{2}}\right)^{1/2}\tag{2}$$
$$y_{2}=z_{2}\left(\frac{-2\ln{r^{2}}}{r^{2}}\right)^{1/2}\tag{3}$$

ただし

$$r^{2}=z_{1}^{2}+z_{2}^{2}\tag{4}$$

の変換を行うと、$y_{1}$, $y_{2}$ はそれぞれ独立に平均 $0$, 分散 $1$ の標準正規分布に従う。

サンプルの破棄を伴わない場合

※こちらの手法を論ずる際には、確率変数に大文字を用いる。なお、小文字と大文字の間に対応関係はないため、注意されたい。

標準一様分布から独立に $X$ と $Y$ を生成し、

$$Z_{1}=\sqrt{-2\ln{X}}\cos{2\pi Y}\tag{5}$$
$$Z_{2}=\sqrt{-2\ln{X}}\sin{2\pi Y}\tag{6}$$

の変換を行うと、$Z_{1}$, $Z_{2}$ はそれぞれ独立に平均 $0$, 分散 $1$ の標準正規分布に従う。

サンプルの破棄を伴う場合、面積 $2\times 2=4$ の正方形と面積 $\pi\times 1\times 1=\pi$ の円の比率から、平均約25％のサンプルが脱落する。そのため、必要な数の乱数を生成するために、かなりの余分なサンプルを必要とする。その点、こちらの手法では、生成した $X$ と $Y$ の数だけ確実に乱数が得られる。

Pythonによる実装

サンプルの破棄を伴う場合

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

z1 = np.random.random(100000) * 2 - 1 # (-1, 1) の乱数生成
z2 = np.random.random(100000) * 2 - 1 # (-1, 1) の乱数生成

z1z2 = z1*z1 + z2*z2 # z1^2 + z2^2 を計算
z1 = z1[z1z2<=1] # z1^2 + z2^2 <= 1 を
z2 = z2[z1z2<=1] # みたすもののみ残す

rr = z1*z1 + z2*z2 # r^2 を計算
y1 = z1 * np.sqrt(-np.log(rr)*2/rr) # y1 に変換
y2 = z2 * np.sqrt(-np.log(rr)*2/rr) # y2 に変換

fig, ax = plt.subplots(1, 2, sharey="row") # ２分割キャンバスを作成
ax[0].hist(y1, bins=1000) # y1 のヒストグラム
ax[0].set_xlabel("y1")
ax[0].set_ylabel("number")
ax[1].hist(y2, bins=1000) # y2 のヒストグラム
ax[1].set_xlabel("y2")
fig.suptitle("Distribution of y1 & y2 (with rejection)")
plt.show()

サンプルの破棄を伴わない場合

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import norm
import math

X = np.random.random(100000) # (0, 1) の乱数生成
Y = np.random.random(100000) # (0, 1) の乱数生成

Z1 = np.sqrt(-np.log(X)*2)*np.cos(Y * 2 * math.pi) # Z1 に変換
Z2 = np.sqrt(-np.log(X)*2)*np.sin(Y * 2 * math.pi) # Z2 に変換

fig, ax = plt.subplots(1, 2, sharey="row") # ２分割キャンバスを作成
ax[0].hist(Z1, bins=1000) # Z1 のヒストグラム
ax[0].set_xlabel("Z1")
ax[0].set_ylabel("number")
ax[1].hist(Z2, bins=1000) # Z2 のヒストグラム
ax[1].set_xlabel("Z2")
fig.suptitle("Distribution of Z1 & Z2 (without rejection)")
plt.show()

証明

サンプルの破棄を伴う場合

式 $(1)$ の関係を満たすようにサンプリングされた $z_{1}$, $z_{2}$ は、中心 $0$ 、半径 $1$ の円の内部で一様分布する（$\because$ 式 $(1)$ を考慮しない状態で独立に得られた $z_{1}$, $z_{2}$ は四隅が $(1,1)$, $(1,-1)$, $(-1,-1)$, $(-1,1)$ で与えられる正方形の中で一様分布する。内部の円のみを切り出したとき、その中もまた一様分布）。

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

z1 = np.random.random(1000) * 2 - 1 # (-1, 1) の乱数生成
z2 = np.random.random(1000) * 2 - 1 # (-1, 1) の乱数生成

z1z2 = z1*z1 + z2*z2 # z1^2 + z2^2 を計算
z1 = z1[z1z2<=1] # z1^2 + z2^2 <= 1 を
z2 = z2[z1z2<=1] # みたすもののみ残す

fig = plt.figure() # 空のキャンバス作成
ax = fig.add_subplot(111) # 軸を１つ追加
ax.scatter(z1, z2, alpha=0.5) # z1, z2 の散布図
ax.set_aspect('equal', adjustable='box') # 縦横の軸を揃える
plt.title("Distribution of z1 & z2")
plt.xlabel("z1")
plt.ylabel("z2")
plt.show()

このとき、円の面積は $\pi\times 1\times 1=\pi$ なので、 $z_{1}$, $z_{2}$ の同時分布の確率密度関数は

$$p(z_{1},z_{2})=\frac{1}{\pi}\tag{7}$$

で与えられる。

$(z_{1}, z_{2})$ の同時分布の微小領域に入る値が、変換後は $(y_{1}, y_{2})$ の同時分布の微小領域に入るという対応関係を考えると、

$$p(y_{1},y_{2})=p(z_{1},z_{2})\left|\frac{\partial p(z_{1},z_{2})}{\partial p(y_{1},y_{2})}\right|\tag{8}$$

という関係が導かれる。この辺りの議論は

の式 $(1), (2)$ の導出に関する記述が参考になる。

上式のヤコビアンを計算するために、

$$z_{1}=r\cos{\theta}$$
$$z_{2}=r\sin{\theta}$$

という変換を行う。これは式 $(4)$ の関係をみたす。ここで

$$\left|\frac{\partial p(z_{1},z_{2})}{\partial p(r,\theta)}\right|=\cos{\theta}\cdot r\cos{\theta}-\sin{\theta}\cdot r(-\sin{\theta})$$
$$=r(\cos^2{\theta}+\sin^2{\theta})=r$$

である。また、変換式 $(2), (3)$ は

$$y_{1}=(-2\ln{r^{2}})^{1/2}\cos{\theta}\tag{2'}$$
$$y_{2}=(-2\ln{r^{2}})^{1/2}\sin{\theta}\tag{3'}$$

と書き換えられ、

$$\left|\frac{\partial p(y_{1},y_{2})}{\partial p(r,\theta)}\right|=-\frac{2}{r}(-2\ln{r^{2}})^{-1/2}\cos{\theta}\cdot(-2\ln{r^{2}})^{1/2}\cos{\theta}-(-2\ln{r^{2}})^{1/2}(-\sin{\theta})\cdot\left(-\frac{2}{r}(-2\ln{r^{2}})^{-1/2}\sin{\theta}\right)$$
$$=-\frac{2}{r}(\cos^2{\theta}+\sin^2{\theta})=-\frac{2}{r}$$

となる。したがって

$$\left|\frac{\partial p(z_{1},z_{2})}{\partial p(y_{1},y_{2})}\right|=\left|\frac{\partial p(z_{1},z_{2})}{\partial p(r,\theta)}\right|\left|\frac{\partial p(y_{1},y_{2})}{\partial p(r,\theta)}\right|^{-1}$$
$$=r\cdot\left(-\frac{r}{2}\right)=-\frac{r^{2}}{2}$$

が求められる。これと式 $(7)$ を、式 $(8)$ の両辺が非負となることに注意して代入すると

$$p(y_{1},y_{2})=\frac{1}{\pi}\frac{r^{2}}{2}$$
$$=\frac{1}{2\pi}r^{2}\tag{9}$$

となるが、ここで式 $(2'),(3')$ の辺々を２乗して足すと

$$y_{1}^{2}+y_{2}^{2}=(-2\ln{r^{2}})(\cos^{2}\theta+\sin^{2}\theta)$$
$$y_{1}^{2}+y_{2}^{2}=-2\ln{r^{2}}$$
$$r^{2}=\exp\left\{-\frac{1}{2}\left(y_{1}^{2}+y_{2}^{2}\right)\right\}$$
$$r^{2}=\exp\left(-\frac{1}{2}y_{1}^{2}\right)\exp\left(-\frac{1}{2}y_{2}^{2}\right)$$

より、式 $(9)$ は

$$p(y_{1},y_{2})=\frac{1}{2\pi}\exp\left(-\frac{1}{2}y_{1}^{2}\right)\exp\left(-\frac{1}{2}y_{2}^{2}\right)$$
$$=\left[\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}y_{1}^{2}\right)\right]\left[\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}y_{2}^{2}\right)\right]$$
$$=p(y_{1})p(y_{2})$$

より、確率変数 $y_{1}$, $y_{2}$ はそれぞれ独立に標準正規分布に従うことがわかる。

サンプルの破棄を伴わない場合

式 $(1)$ の関係をみたす $z_{1}, z_{2}$ について

$$X=z_{1}^{2}+z_{2}^{2}\tag{10}$$

と変数変換することを考える。ここからの議論は

の導出の流れに従う。

確率変数 $X$ の分布を $q(X)$ と置くと、これは $p(z_{1},z_{2})=\frac{1}{\pi}$ の同時分布を用いて

$$q(X)=\int_{0}^{1}\int_{0}^{\sqrt{1-z_{1}^{2}}}\delta\left(X-(z_{1}^{2}+z_{2}^{2})\right)p(z_{1},z_{2})dz_{1}dz_{2}$$
$$=\frac{1}{\pi}\int_{0}^{1}\int_{0}^{\sqrt{1-z_{1}^{2}}}\delta\left(X-(z_{1}^{2}+z_{2}^{2})\right)dz_{1}dz_{2}\tag{11}$$

と書ける。これを２次元単位球における積分（

を参照）に変換することを考えると、動径 $r$ には式 $(4)$ の関係をそのまま用いることができる。また、２次元単位球面の面素を $dS_{1,2}$ とおくと

$$dz_{1}dz_{2}=rdrdS_{1,2}$$

が成り立ち、式 $(11)$ の積分は

$$q(X)=\frac{1}{\pi}\int_{0}^{1}\int\delta\left(X-(z_{1}^{2}+z_{2}^{2})\right)rdrdS_{1,2}$$
$$=\frac{1}{\pi}\int_{0}^{1}\delta(X-r^{2})rdr\int dS_{1,2}\tag{11'}$$

と書き換えることができる。さらに、

$$\nu=r^{2}$$

の変数変換を行うと

$$d\nu=2rdr$$

より、これを代入して式 $(11')$ の積分を行うと

$$q(X)=\frac{1}{2\pi}\int_{0}^{1}\delta(X-\nu)d\nu\int dS_{1,2}$$
$$=\frac{1}{2\pi}\cdot 1\cdot\frac{(2\pi)^{2/2}}{\Gamma\left(\frac{2}{2}\right)}$$
$$=\frac{1}{2\pi}\frac{2\pi}{1}=1$$

より、確率変数 $X$ は標準一様分布に従う。

また、このとき式 $(2),(3)$ は

$$y_{1}=\frac{z_{1}}{r}(-2\ln{X})^{1/2}$$
$$y_{2}=\frac{z_{2}}{r}(-2\ln{X})^{1/2}$$

と書ける。ここで

$$\frac{z_{1}}{r}=\frac{z_{1}}{\sqrt{z_{1}^{2}+z_{2}^{2}}}=\cos{\theta}$$
$$\frac{z_{2}}{r}=\frac{z_{2}}{\sqrt{z_{1}^{2}+z_{2}^{2}}}=\sin{\theta}$$

とおけるが、このとき $\theta$ の値も区間 $(0,2\pi)$ で一様分布する。よって、標準一様分布に従う確率変数 $Y$ を用いて

$$\theta=2\pi Y$$

と書ける。以上より、$y_{1}\rightarrow Z_{1}$, $y_{2}\rightarrow Z_{2}$ と変数を書き換えると、「サンプルの破棄を伴う場合」の変換式 $(2), (3)$は

$$Z_{1}=\sqrt{-2\ln{X}}\cos{2\pi Y}\tag{5☆}$$
$$Z_{2}=\sqrt{-2\ln{X}}\sin{2\pi Y}\tag{6☆}$$

となる。

すなわち、２種類の確率変数 $z_{1}$, $z_{2}$ を生成して、それらを式 $(1)$ の条件で縛るのではなく、はじめからそれらの確率変数を式 $(5),(6)$ に代入すると、「サンプルの破棄を伴う場合」と同様に、標準正規分布に従う２種類の独立な確率変数が得られる。

参考：$X=z_{1}^{2}+z_{2}^{2}$ の分布

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

z1 = np.random.random(100000) * 2 - 1 # (-1, 1) の乱数生成
z2 = np.random.random(100000) * 2 - 1 # (-1, 1) の乱数生成

z1z2 = z1*z1 + z2*z2 # z1^2 + z2^2 を計算
z1 = z1[z1z2<=1] # z1^2 + z2^2 <= 1 を
z2 = z2[z1z2<=1] # みたすもののみ残す

X = z1*z1 + z2*z2 # X を計算

plt.hist(X, bins=100) # X のヒストグラム
plt.title("Distribution of X")
plt.xlabel("X")
plt.ylabel("number")
plt.show()