2차원 Gaussian 분포 생성(Creating a 2d Gaussian Distribution)

Box-Muller Method는 균등분포(uniform distribution)를 따르는 두 개의 난수로부터 평균이 0이고 분산이 1인 표준정규분포(standard normal distribution) 난수를 생성하는 대표적인 방법이다. 이를 이용해서 2차원 가우시안 분포를 생성하는 알고리즘을 구현하자.

$$g(x, y; \sigma)= \frac{1}{2\pi \sigma^2}\exp\left( - \frac{x^2 +y^2 }{2\sigma^2 }\right)$$

 

먼저 \(\sigma=1\)인 표준 2차원 가우시안 분포를 생성한 뒤, 마지막에 \(\sigma\)를 곱하여 원하는 분산을 갖는 분포를 얻는다.

내장 함수 rand()를 이용하면 \([0,1]\) 구간에서 균등분포를 따르는 난수를 생성할 수 있다. \([0,1]\times[0,1]\) 영역에서 균등하게 분포하는 두 개의 독립적인 확률변수 \(x_1,x_2\)를 표준정규분포를 따르는 확률변수 \(y_1,y_2\)로 변환하려면 확률이 보존되어야 하므로

$$ dx_1 dx_2 = \frac{1}{\sqrt{2\pi}} e^{-y_1^2 /2 } dy_1 \frac{1}{\sqrt{2\pi}} e^{-y_2^2 /2 } dy_2 $$

을 만족시켜야 한다. 극좌표 $y_1 = \rho \cos \phi$, $y_2 = \rho\sin \phi$을 도입하면 

$$ dx_1 dx_2 = e^{-\rho^2/2} \rho d\rho \frac{d\phi}{2\pi}$$로 쓸 수 있으므로 미소 변환은

$$ dx_1 = e^{-\rho^2/2} \frac{d\rho^2}{2},\quad dx_2 = \frac{d\phi}{2\pi}$$ 로 선택할 수 있다. 따라서

$$ x_1 = \exp[-\rho^2/2] = \exp[- (y_1^2 + y_2^2 )/2],$$

$$ x_2 = \frac{1}{2\pi} \tan^{-1} (y_2/y_1).$$

그리고, 역변환은

$$ y_1 = \sqrt{-2 \ln (x_1) } \cos (2\pi x_2),\quad y_2 = \sqrt{-2 \ln(x_1) } \sin (2\pi x_2)$$

을 얻는데, 이를 Box-Muller 변환이라 한다.

삼각함수의 계산을 피하기 위해서 반지름이 1인 원 내부에서 정의되는 두 개의 균일한 랜덤 변수 $v_1, v_2$를 도입하면, 다음 변환에 의해서 반지름 1인 원 영역이 변의 길이가 1인 정사각형 영역으로 보내진다.

$$x_1 = R^2(=S) = v_1^2 + v_2^2, \quad 2\pi x_2 = \tan^{-1} (v_2/v_1).$$

따라서 반지름이 1인 원 내부의 랜덤 변수를 써서 2차원의 Gaussian 분포를 만들어주는 변환은

$$ y_1 = \sqrt{-2\ln S} \frac{v_1}{\sqrt{S}}, \quad y_2 = \sqrt{-2\ln S} \frac{v_2}{\sqrt{S}}$$

 

// 참고;  Numerical Receipe;
#define drand() ((double)rand() / RAND_MAX)   // [0,1]
// mean = 0; sigma = 1;
void normal_2D(double *x, double *y) {
    while (1) {
        double V1 = -1. + (2. * drand());
        double V2 = -1. + (2. * drand());
        double S = V1 * V1 + V2 * V2;
        if (S < 1.) {
            double rad = sqrt(-2.*log(S) / S);
            *x = V1 * rad;
            *y = V2 * rad;
            /* if (mean, sigma) were given,
            *x = meanx + sigma * (*x);
            *y = meany + sigma * (*y);
            */
            return;
        }
    }
};