기각 샘플링을 이용한 난수생성 알고리즘

이번 포스팅에서는 기각 샘플링 (rejection sampling)이라는 방법을 사용하여 임의의 확률밀도함수를 따르는 난수 (random number)를 생성하는 방법에 대해서 살펴보겠습니다. 추가로 C/C++ 프로그램에서 기각 샘플링을 구현한 예시도 소개하겠습니다.

 

반응형

 

여기서 다루는 내용은 연속적인 값을 가질 수 있는 변수가 특정한 확률 분포를 따르는 경우에 대한 것인데요. 이 때 해당 변수의 값이 특정 범위 내에 있을 확률은 확률밀도함수의 정적분으로 주어지는 것이 특징입니다. 확률밀도함수의 구체적인 정의와 특징에 대해서는 다음 포스팅에 더 자세하게 소개되어 있습니다.

 

수학 상식 : 확률 분포와 확률밀도함수

 

기각 샘플링의 원리

컴퓨터를 이용한 시뮬레이션을 하다 보면 주어진 확률 분포를 따르는 난수를 발생시켜야 할 때가 있는데요. 문제는 확률밀도함수의 형태에 따라 이 작업의 난이도가 매우 높아질 수 있다는 것입니다. 기각 샘플링은 이러한 문제를 해결하기 위한 방법 중의 하나로서, 우리가 원하는 확률밀도함수에 따라 직접적으로 난수를 추출하는 대신 이를 두 단계로 나눠서 샘플링을 하는 것입니다.

 

이를 위해서는 먼저 본래의 확률밀도함수 외에 또 다른 함수를 하나 더 상정해야 합니다. 이 새로운 함수를 덮개함수 (envelope function)라고도 하는데, 다음 두 가지 조건을 만족해야 합니다.

 

• 변수가 가질 수 있는 모든 값에 대해서 덮개 밀도함수는 확률밀도함수보다 커야 합니다. 이말인즉슨 덮개함수를 전체 범위에 대해서 적분하면 1보다 크다는 얘기가 되죠. 가장 단순한 예시로서 확률밀도함수의 최대값으로 주어지는 상수함수를 생각해볼 수 있는데요. 이것 역시 한 방법이긴 합니다만, 기각 샘플링의 효율성은 떨어지는 편입니다. 덮개함수의 형태가 본래 확률밀도함수와 유사할수록 좋습니다.
• 변수의 최소값에서부터 임의의 값까지의 범위에 대해 덮개함수를 정적분할 수 있어야 하고, 이렇게 주어지는 누적 덮개함수의 역함수를 쉽게 계산할 수 있어야 합니다. 확률밀도함수로부터 누적분포함수를 구하거나 누적분포함수의 역함수를 구하는 것이 어렵기 때문에 기각 샘플링을 사용하는 것이거든요. 그렇기 때문에 덮개함수는 다루기 쉬운 것일수록 좋습니다.

 

기각 샘플링의 첫번째 단계는 덮개함수를 확률분포로 하는 난수를 먼저 추출하는 것입니다. 그리고 추출된 값에서 실제 확률밀도함수와 덮개함수의 크기를 비교한 다음, 샘플을 받아들일지 혹은 기각할지의 여부를 결정하는 것이 두번째 단계입니다.

 

 

기각 샘플링을 수치적으로 구현하기 위해서 필요한 것은 연속균등분포 (continuous uniform distribution)에 따라 난수를 발생시키는 함수입니다. 예컨대 0에서 1 사이의 값을 가지는 실수형 변수를 리턴하는 함수를 반복적으로 호출하면서 균일한 확률 분포를 가지는 난수를 발생시킬 수 있으면 편리합니다.

 

마지막으로 한 가지 짚고 넘어갈 점은, 샘플을 받아들일지의 여부를 결정할 때 연속균등분포를 따르는 난수를 한번 더 추출해야 한다는 것입니다. 프로그램에서는 해당 함수를 한번 더 호출해야 한다는 얘기가 되겠습니다.

 

예시 : 가우스 함수

예시로서 하나의 피크를 가지는 확률밀도함수를 상정해서, 기각 샘플링을 통해 난수를 추출해 봅시다. 수학적으로는 피크에서의 위치를 평균값으로 하고 유한한 크기의 표준편차를 가지는 가우스 함수로 이를 구현할 수 있습니다. 가우스 함수는 정규 분포라는 확률 분포를 나타내는 확률밀도함수이기도 한데요. 가우스 함수의 정의와 특징에 대한 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

수학 상식 : 정규 분포와 가우스 함수

 

가우스 함수의 경우 누적분포함수를 해석적으로 구하는것도 어려울뿐더러, 누적분포함수의 역함수를 구하는 것은 더 어렵습니다. 그렇기 때문에 기각 샘플링의 좋은 예시가 될 수 있고, 덮개함수로서 아크탄젠트 함수의 미분을 생각해볼 수 있죠. 다시 말해서 덮개함수를 적분한 누적 덮개함수는 아크탄젠트 함수가 되고, 그 역함수는 탄젠트 함수가 되므로 매우 편리합니다.

 

정규 분포를 따르는 변수의 표본을 추출하기 위한 프로그램을 C언어로 작성해볼 수 있습니다. 연속균등분포를 따르는 난수를 발생시키기 위해서 수학 함수 라이브러리의 rand 함수를 사용했습니다. 다만 rand 함수는 0에서 RAND_MAX 사이의 정수를 리턴하는 관계로, rand 함수의 리턴값을 RAND_MAX로 나눠주면 0에서 1까지의 값을 가지고 연속균등분포를 따르는 난수를 얻을 수 있을 것입니다.

 

C언어 소스 코드입니다.

 

test1_rng_gauss.c [다운로드]

 

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>

// total number of try
unsigned int n_try_;
// number of accepted samples
unsigned int n_sample_acc_;
// number of rejected samples
unsigned int n_sample_rej_;

// mean value of Gaussian
double mean_;
// standard deviation
double sigma_;

/* function pointer to
 * the uniform random number generator */
double (*ptr_rng_uniform_)();

double fac_norm_gauss_;
double fac_int_env_;

/* function to initialize
 * random numger generator
 * of normal distribution 
 * mean_in : mean value
 * sigma_in : standard deviation
 * ptr_in_rng_uniform :
 *   function pointer
 *   for a uniform random number generator */
void init_rand_gauss(double mean_in,
                     double sigma_in,
         double (*ptr_in_rng_uniform)()) {
    n_try_ = 0;
    n_sample_acc_ = 0;
    n_sample_rej_ = 0;

    mean_ = mean_in;
    sigma_ = sigma_in;

    ptr_rng_uniform_ = ptr_in_rng_uniform;

    fac_norm_gauss_ = 1. / sqrt(2. * M_PI);
    fac_int_env_ = sqrt(M_PI);
}

double get_rand_double() {
    return (double)rand() /
           ((double)RAND_MAX);
}

/* function to sample
 * random numbers according to
 * the Gaussian distribution function
 * fout_accepted :
 *   file pointer to write accepted samples
 * fout_rejected :
 *   file pointer to write rejected samples */
double get_rand_gauss(FILE *fout_accepted,
                      FILE *fout_rejected) {

    int accepted = 0;
    double x_sample = 0.;
    double y_sample = 0.;

    while (accepted == 0) {
        double p_env = fac_int_env_ *
                       (*ptr_rng_uniform_)();

        y_sample = sqrt(2.) *
            tan(p_env * sqrt(M_PI) - 0.5 * M_PI);

        x_sample =
            mean_ + sigma_ * y_sample;

        double p_accept =
            (*ptr_rng_uniform_)();

        double f_pdf = fac_norm_gauss_ *
            exp(-0.5 * fabs(y_sample * y_sample));
        double f_env = fac_norm_gauss_ /
            (1. + 0.5 * fabs(y_sample * y_sample));

        n_try_ += 1;
        if (f_pdf > p_accept * f_env) {
            accepted = 1;
            n_sample_acc_ += 1;
        } else {
            accepted = 0;
            n_sample_rej_ += 1;
        }

        if (accepted == 0 &&
            fout_rejected != NULL) {
            fprintf(fout_rejected,
                "  %d  %d  %e\n", n_try_,
                n_sample_rej_, x_sample);
        }

        if (accepted != 0 &&
            fout_accepted != NULL) {
            fprintf(fout_accepted,
                "  %d  %d  %e\n", n_try_,
                n_sample_acc_, x_sample);
        }

    }

    return x_sample;
}

int main(int argc, char *argv[]) {

    double (*ptr_rng_uni)();

    ptr_rng_uni = &get_rand_double;
    srand((unsigned int)time(NULL));

    init_rand_gauss(0., 1., ptr_rng_uni);

    FILE *fout_acc;
    FILE *fout_rej;

    fout_acc =
        fopen("list_sample_acc.dat", "w");
    fout_rej =
        fopen("list_sample_rej.dat", "w");

    while (n_try_ < 100000) {
        double x =
            get_rand_gauss(fout_acc, fout_rej);
    }

    fclose(fout_acc);
    fclose(fout_rej);

    return 0;
}

 

기각 샘플링을 시도한 횟수, 가우스 함수의 평균과 표준편차 등을 전역변수로 가지고 있습니다. 그리고 ptr_rng_uniform_이라는 함수 포인터가 있는데요. 이는 0에서 1 사이의 연속균등분포에 따라 난수를 발생시키는 함수의 포인터입니다.

 

초기화를 위한 init_rand_gauss 함수를 호출하여 가우스 함수의 평균과 표준편차를 설정하고, 연속균등분포에 따라 난수를 발생시키는 함수의 포인터를 지정해줍니다. 이번 포스팅에서는 C언어 수학 함수 라이브러리에 있는 함수를 사용합니다만, 외부 라이브러리에 있는 난수발생 알고리즘을 사용하는 것 역시 가능합니다.

 

초기화가 되면, get_rand_gauss 함수를 반복적으로 호출하여 가우스 함수를 확률밀도함수로 하는 변수의 표본을 추출할 수 있게 됩니다. 여기서 get_rand_gauss 함수는 두 개의 파일 포인터를 매개변수로 받고 있는데, 받아들인 표본과 기각된 표본들을 각각의 파일에 기록하기 위한 용도입니다. 이게 필요없는 경우라면 NULL 포인터를 매개변수로 주면 되겠습니다.

 

 

프로그램을 빌드하고 실행시켜보면 정규 분포를 따르는 변수의 표본들을 얻을 수 있습니다. 평균값에서 멀리 떨어질수록 덮개함수화 확률밀도함수의 차이가 크기 때문에, 많은 샘플이 기각되는 것을 알 수 있는데요. 그럼에도 불구하고 평균값 근처에서는 확률밀도함수와 덮개함수가 매우 가깝기 때문에, 거의 대부분의 샘플을 받아들이게 되고 기각 샘플링의 효율이 그렇게 나쁘지 않은 것입니다.

 

---

 

앞에서 함수 포인터라는 것을 언급했었는데, 이를 사용하면 함수 자체를 다른 함수의 매개변수로 넘겨주는 것이 가능합니다. 구체적인 방법과 예시들은 다음 포스팅에 더 자세하게 소개되어 있습니다.

 

C/C++ 에서 함수를 매개변수로 사용하기