Box-Muller transformation

예제 · 케이플러스 한성탁 · 29 Oct 2017

Box-Muller transformation으로 정규 난수 생성하기

박스-뮬러 변환은 가장 간단한 정규 난수 발생법 중 하나입니다. Box와 Muller에 의해 1958년 발표되었습니다.

$U_1, U_2 \sim \mathrm{U}(0,1)$ 일 때,

$X_1 = \sqrt{-2logU_1}cos(2\pi U_2) \sim \mathrm{N(0,1)}$

$X_2 = \sqrt{-2logU_1}sin(2\pi U_2) \sim \mathrm{N(0,1)}$

이 되어 서로 독립이고 같은 분포를 가지는(iid) 2개의 표준 정규 난수를 발생시키게 됩니다.

먼저 $X_1, X_2 \sim \mathrm{N(0,1)}$ 이라고 가정합니다. $X_1, X_2$ 를 좌표에 놓고 극좌표법으로 좌표를 구하면 $(X_1, X_2) = (R, \theta)$ 이고

$R^2 = X_1^2 + X_2^2$ $\theta = arctan \frac{X_2}{X_1}$

이 됩니다.

$R^2$을 보면 $X_1^2$과 $X_2^2$의 합으로 이루어져 있고, 표준 정규 난수의 제곱은 카이제곱 분포 $\chi^2(1)$를 따르고, 카이제곱 분포는 지수 분포로 표현할 수 있으며, 지수분포는 역변환법(Inverse transformation)에 의해 $\mathrm{U}(0,1)$의 꼴로 표현할 수 있습니다. 이를 식으로 표현하면 다음과 같습니다.

$R^2 = X_1^2 + X_2^2 \sim \chi^2(2) \sim \mathrm{Exp}\left(\frac{1}{2}\right) \sim -2log\mathrm{U}$

다른 방법으로 R의 누적분포함수를 직접 구하면 다음과 같습니다.

$\begin{align*} f \left( X_1, X_2 \right) & = \frac{1}{2\pi} \mathrm{exp} \left( \frac{-(X_1^2 + X_2^2)}{2} \right) \\ & = \frac{1}{\sqrt{2\pi}} \mathrm{exp} \left( \frac{-X_1^2}{2} \right)\frac{1}{\sqrt{2\pi}}\mathrm{exp} \left( \frac{-X_2^2}{2} \right) \\ & = f \left( X_1\right)f \left(X_2 \right) \end{align*}$ $\begin{align*} P \left(R \leq r_0 \right) & = \int_{x_1^2 + x_2^2 \leq r_0^2} \frac{1}{2\pi}\mathrm{exp} \left( \frac{-(x_1^2 + x_2^2)}{2} \right) dx_1dx_2 \\ & = \frac{1}{2\pi} \int_{0}^{2\pi} \int_{0}^{r_0} \mathrm{exp} \left( \frac{-r^2}{2} \right) rd\theta dr \\ & = \int_{0}^{r_0} \mathrm{exp} \left( \frac{-r^2}{2} \right) rdr \\ & = \int_{0}^{\frac{r_0^2}{2}} \mathrm{exp} \left( -z \right) dz \\ & = 1 - \mathrm{exp} \left( \frac{-r_0^2}{2} \right) \\ & = F\left( r_0^2 | \lambda = \frac{1}{2} \right) \end{align*}$

위의 식을 통해 $X_1, X_2$는 서로 독립이라는 것과 $P \left(R \leq r_0 \right)$ 은 지수분포의 누적분포함수라는 것을 확인 할 수 있습니다($r_0^2 >= 0$ 이므로). 이때 $ \lambda = \frac{1}{2}$ 입니다.

이제 위 식을 아래에서 위로 반대로 거슬러보면, $-2log \mathrm{U}$를 구하여서 이 값을 통해 두 표준 정규 난수의 제곱의 합을 알 수 있습니다. 그런데 $x_1, x_2$의 제곱의 합이 같은 경우는 무수히 많고(원에 대한 방정식 $x_1^2+x_2^2 = r^2$ 와 같습니다) 모두 정규분포를 따릅니다.

그러므로 이때 반지름이 정해진 경우 $\theta$는 $0 \leq \theta < 2\pi$ 범위에서 아무거나 공평하게 고르면 됩니다. 따라서 $ \theta = 2\pi U_2$, $ U_2 \sim \mathrm{U}(0,1) $ 로 정합니다.

따라서

$\begin{align*} X_1 = R \mathrm{cos}(\theta) = \sqrt{-2logU_1}\mathrm{cos}(2\pi U_2) \sim \mathrm{N(0,1)} \\ X_2 = R \mathrm{sin}(\theta) = \sqrt{-2logU_1}\mathrm{sin}(2\pi U_2) \sim \mathrm{N(0,1)} \end{align*}$

가 성립합니다.

아래는 박스-뮬러 변환법을 R 코드로 구현한 것입니다.

rnorm.boxmuller = function(n, mean=0, sd=1){
  U1 = runif(ceiling(n/2))
  U2 = runif(ceiling(n/2))
  Z1 = sqrt(-2*log(U1))*cos(2*pi*U2)
  Z2 = sqrt(-2*log(U1))*sin(2*pi*U2)
  c(Z1, Z2)[1:n] * sd + mean
}

rnorm.boxmuller(100)

꽤 간단하게 작성 가능합니다.