ボックス＝ミュラー法（Box-Muller's method）は、一様分布にしたがう確率変数から、標準正規分布にしたがう確率変数を生成する手法です。
コンピュータで標準で生成できる乱数から、正規分布にしたがう乱数（正規乱数）を作るために使用します。
この記事では、正規乱数の生成する方法として、サンプルの破棄をともなう手法とボックス＝ミュラー法を紹介し、その後、それらの実装と導出を行います。
手法
正規乱数を生成するための手法を２種類紹介します。
それぞれを区別するため、確率変数の表現を以下のように使い分けます。
- サンプルの破棄をともなう手法：確率変数を小文字で書く
- ボックス＝ミュラー法（サンプルの破棄をともなわない手法）：確率変数を大文字で書く
それぞれの大文字と小文字の間に対応関係はないので、注意してください。
サンプルの破棄をともなう場合
標準一様分布（ \((0, 1)\) の値をとる一様分布）から独立に \(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}\) に
$$r^{2}=z_{1}^{2}+z_{2}^{2}\tag{2}$$
のもとで
$$y_{1}=z_{1}\left(\frac{-2\ln{r^{2}}}{r^{2}}\right)^{1/2}\tag{3}$$
$$y_{2}=z_{2}\left(\frac{-2\ln{r^{2}}}{r^{2}}\right)^{1/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\) の標準正規分布にしたがいます。
手法の比較
サンプルの破棄をともなう場合、範囲 \((-1, 1)\) をとる変数を２種類用意します。
これを２次元にプロットすると、各点は面積 \(2\times 2=4\) の正方形の中に均等に分布します。
しかし式 \((1)\) の条件より、面積 \(\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\) の円の内部で均等に分布します。
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})\) の同時分布の微小領域に入る値が、
式 \((3), (4)\) による変換後は \((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}$$
という関係が導かれます。
\(\left|\frac{\partial p(z_{1},z_{2})}{\partial p(y_{1},y_{2})}\right|\) はヤコビアンです。
このあたりの議論は
の式 \((1), (2)\) の導出を参考にしてください。
ヤコビアンを計算するために、式 \((2)\) の関係をみたすよう
$$z_{1}=r\cos{\theta}$$
$$z_{2}=r\sin{\theta}$$
という変換を行うと
$$\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$$
となります。
また、変換式 \((3), (4)\) は
$$y_{1}=(-2\ln{r^{2}})^{1/2}\cos{\theta}\tag{3'}$$
$$y_{2}=(-2\ln{r^{2}})^{1/2}\sin{\theta}\tag{4'}$$
と書き換えられるため、
$$\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}$$
となります。
ここで式 \((3'),(4')\) の辺々を２乗して足すと
$$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}^{2}+z_{2}^{2}\leq 1\tag{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\) には式 \((2)\) の関係
$$r^{2}=z_{1}^{2}+z_{2}^{2}\tag{2}$$
をそのまま用いることができます。
２次元単位球面の面素を \(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\) は標準一様分布にしたがいます。
また、このとき式 \((3),(4)\) は
$$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}$$
と変数を書き換えると、「サンプルの破棄をともなう場合」の変換式 \((3), (4)\) は
$$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()
Comments