Rejection Sampling

절대온도 $T$일 때 질량 $m$인 이상기체의 속력은 Maxwell-Boltzmann 분포를 따른다.

$$P(v)dv = 4\pi \left(\frac{m}{2\pi kT}\right)^{3/2}v^2 e^{-mv^2/2kT}dv$$

속력 $v$는 0부터 $\infty$까지의 값을 갖는다. 그럼 이 속력분포를 가지는 $N$개의 속력 샘플을 만들고 싶은 경우 어떻게 할까? 우선 무한대 범위까지 모두 처리할 수 없으므로 적당한 속력 범위 내의 샘플만 만들자. 계산하기 쉽도록 하기 위해서 $x^2 = mv^2 /kT$인 차원없는 변수 $x$를 사용한다. 확률이 보존되어야 하므로 $P(v)dv= f(x)dx$에서 $x$에 대한 확률밀도함수는

$$ f(x) = \sqrt{\frac{2}{\pi}} x^2 e^{-x^2/2}$$

로 주어진다. 가능한 샘플 $x$는 다음의 단계로 얻어진다.

1. 주어진 $x$ 범위 $[x_\text{min}, x_\text{max}]$ 사이에서 랜덤하게 한 값 $x$를 취하고,

2. $[0,1]$ 사이에서 랜덤한 값 $p$을 취한 후 

3. $ p> f(x)$이면 $x$는 버리고 다시 과정 1로 돌아간다. 

 

다음 예는 $x$를 $[0,4]$사이에서 500 구간으로 나누어서 만든 50000개의 샘플의 분포와 Maxwell-Boltzmann 분포를 비교한 것이다.

그림의 붉은점은 accepted된 값이고 푸른점은 rejected된 값이다.

    const CSize sz = CSize(500, 500);
    raster.SetDimensions(sz, 24);
    srand(unsigned(time(0)));
    int hist[500] = {0};
    int required = 50000;
    const double xmin  = 0, xmax = 4;
    while (required) {
        double x = double(rand()) / RAND_MAX; //[0,1]
        double v = (xmax - xmin) * x + xmin;  //[xmin,xmax]
        double p = double(rand()) / RAND_MAX; //[0,1]
        double maxwell = MaxwellDist(v);
        CPoint q = CPoint(int(x * (sz.cx - 1)), int(p * (sz.cy - 1)));
        if (p <= maxwell) {
            raster.SetPixel(q.x, sz.cy - 1 - q.y, 0xFF0000);
            hist[q.x]++;
            required--;
        }
        else
            raster.SetPixel(q.x, sz.cy - 1 - q.y, 0xFF);	
    }
    int s = 0;
    for (int i = 0; i < sz.cx; i++) 
        s += hist[i];
    FILE *fp = fopen("maxwell.txt", "w");
    for (int i = 0; i < sz.cx; i++)
        fprintf(fp, "%f, %f\n", (xmax - xmin) * double(i) /(sz.cx - 1) + xmin, \
                                double(hist[i]) / s);
    fclose(fp);