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

Python
スポンサーリンク

ボックス=ミュラー法(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") # 2分割キャンバスを作成
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") # 2分割キャンバスを作成
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) # 軸を1つ追加
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変数関数の不定積分における置換積分からはじめて、最後は多変数関数の定積分における置換積分まで、概念図を交えながら、その公式と原理を解説する。

の式 \((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')\) の辺々を2乗して足すと

$$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}$$

と変数変換することを考える。

ここからの議論は

1次元正規変数の平方和がしたがう分布【カイ二乗分布】の導出
独立の正規分布にしたがう確率変数を複数用意し、それらを2乗して足し合わせて新たな確率変数を作ったとき、その変数はカイ二乗分布にしたがう。この記事では、確率変数の変換にともなう確率密度関数の変換公式を出発点にして、平方和とカイ二乗分布の関係を導出する。

の導出の流れにしたがう。

確率変数 \(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}$$

と書ける。

これを2次元単位球における積分(

標準正規分布の規格化条件から、M次元単位球の表面積を求める
正規分布$$\mathcal{N}(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\{-\frac{1}{2\sigma^2}(x-\mu)^2\}$$は様々な特徴を持つ重要な分布であるが、こ...

を参照)に変換することを考えると、動径 \(r\) には式 \((4)\) の関係をそのまま用いることができる。

また、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\) は標準一様分布にしたがう。

また、このとき式 \((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☆}$$

となる。

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

参考:\(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()

コメント