Rejection Sampling

확률 분포로부터 샘플 추출

uniform distribution이나 정규분포의 샘플은 컴퓨터를 이용하면 코드 한 줄이면 추출할 수 있다.

MATLAB의 경우는 아래와 같은 코드를 이용할 수 있다.

uniform_sample = rand(1,1)

normal_sample = randn(1,1)

그런데, 임의의 수식을 가진 확률밀도함수로부터 랜덤한 샘플을 추출하려면 어떻게 해야할까?

Rejection Sampling

대략적인 수식을 알고 있는 어떤 확률밀도함수가 있다고 하자.

가령 아래와 같은 함수 를 생각해볼 수 있다.


식 (1)

이 함수 는 정확히 말하면 확률밀도함수라고는 할 수 없다. 왜냐하면 부터 까지 이 함수를 적분했을 때의 전체 면적이 1이 아니기 때문이다. (정확히는 이다.)

하지만, 이와 같이 정확한 확률밀도 함수를 알기 곤란할 때나 확률밀도 함수의 수식은 있지만 해당 함수로부터 sample을 추출하기 어려운 경우 sampling 방법이 필요할 수 있다.

우리가 샘플을 추출하고자 하는 이 유사 확률 분포를 ‘타겟 분포(target distribution)’라고 이름 붙이고, 로 쓰도록 하자.

이 타겟 분포를 plot 해보면 다음과 같은 형태를 가지고 있다는 것을 알 수 있다.


그림 1. 수식을 알고 있는 임의의 유사-확률밀도 함수

만약에 우리가 이런 형태의 유사 확률밀도함수로부터 샘플을 추출하려고 하면 어떻게 샘플을 추출할 수 있을까?

이럴 때 사용할 수 있는 것이 rejection sampling이라고 할 수 있다.

제안 분포(혹은 유사 분포)

rejection sampling의 첫 단계는 제안 분포(proposal distribution)를 설정하는 것이다.

앞으로 제안 분포를 로 쓰도록 하자.

제안 분포 는 우리가 쉽게 샘플을 추출할 수 있는 분포를 이용해서 생성할 수 있다.

가령, uniform distribution이 생긴게 가장 단순하기 때문에 uniform distribution을 이용해서 제안 분포를 만들 수 있을 것이다.

교과서에서는 타겟 분포와 가장 유사한 형태의 제안 분포를 설정하는 것이 좋다고 되어 있다.

(하지만, uniform distribution을 쓰는게 제일 속편하다.)

이렇게 제안 분포를 일단 결정했으면 타겟 분포와 제안 분포를 일단 그려보자.

이 예시에서는 제안 분포로 결정한 uniform distribution의 정의역의 범위를 아래와 같이 설정하였다.

그 이유는 범위에서 가능한 타겟 분포를 모두 포함할 수 있도록 하기 위해서이다.


식 (2)

그러면, 우리의 제안분포는 다음과 같다.


식 (3)


그림 2. 타겟분포와 제안분포

그런 다음 제안 분포에 상수배()를 취해주어 축에서도 타겟 분포를 모두 포함할 수 있도록 해주자.


그림 3. 타겟분포와 제안분포에 상수배를 취해주어 타겟 분포를 모두 포함할 수 있도록 한 것.

샘플링 과정

위의 과정에서처럼 제안분포와 적절한 상수를 정했다면 rejection sampling을 수행해보도록 하자.

가장 먼저 해주어야 할 일은 제안분포 에서 샘플 하나()를 추출하는 것이다.


그림 4. 제안분포(여기서는 uniform distribution)에서 샘플 하나를 추출한다.

그런 다음, 타겟 분포 와 상수배를 취한 제안 분포 의 likelihood를 비교해준다.

즉, 제안분포 로부터 추출한 에 대해 의 함수 값을 비교해주는 것이다.

이 때, 크기를 비교할 때에는 두 함수값을 나눠줘서 비교하는 방식으로 수행한다.

즉, 아래와 같은 방법으로 비교를 수행하도록 하자.


식 (4)

식 (4)를 잘 생각해보면 타겟 분포의 크기가 만큼 높으면 1에 가까운 값이 나올 것이고 그렇지 않으면 작은 값이 나올 것이다.

이제 식 (4)의 값을 아래와 같은 정의역으로 정의되는 uniform distribution의 샘플 값과 비교할 수 있도록 하자.

즉, 이것을 그림으로 설명하면 다음과 같다.


그림 5. 주어진 에 대해 likelihood ratio를 비교하여 샘플의 accept/reject 여부를 결정한다.

위의 그림 5에서와 같이 likelihood ratio를 uniform distribution의 출력값과 비교하게 되면 타겟 분포 의 높이가 만큼 높은 곳일 수록 accept될 확률이 높다는 것을 알 수 있다.

특히, 비교를 위해 얻은 uniform distribution의 출력값은 0에서 1사이의 값이 랜덤하게 출력된다는 점에 주목해서 생각해보자.

알고리즘을 간략하게 쓰면 다음과 같다.

Set Repeat until

  1. Sample and .
  2. If , then accept and increment the counter by 1. Otherwise, reject.

Rejection sampling 결과

Rejection sample의 결과를 일부 표현하면 다음과 같다.


그림 6. accept 된 샘플들과 reject 된 샘플들

그리고 최종적으로 얻어진 샘플들을 원래의 구하고자 했던 타겟 분포 와 함께 히스토그램으로 그리면 다음과 같다.


그림 7. 타겟 분포 와 accept된 샘플들을 histogram으로 그린 것.
(다른 그림들과 축의 높이가 다른 것에 주목하자.)

이 결과를 얻기 위한 MATLAB 코드는 아래와 같다.

clear; close all; clc;

rng(1)
n = 500;
target = @(x) 0.3*exp(-0.2 * x.^2) + 0.7 * exp(-0.2 * (x - 10).^2);
xx = linspace(-10,20, 1000);

%% Uniform distribution을 proposal distribution으로 두고 sampling
% 성능에 크게 차이는 없으나 reject되는 sample 수가 조금 차이가 있음.

pseudo_dist2 = @(x) (x>=-10 * x<20) / 30;

figure;
plot(xx, target(xx));
hold on;
plot(xx, pseudo_dist2(xx)*21); % 21은 30 * 0.7을 얘기함. 30은 출력값이 1인 x의 범주, 0.7은 target의 최고 높이.

x_q = (rand(1, n) - 0.5) * 30 + 5; % -10에서 20사이의 uniform distribution

crits = target(x_q) ./ (pseudo_dist2(x_q) * 21);
coins = rand(1, length(crits));

x_p_uniform = x_q(coins<crits);

figure; h = histogram(x_p_uniform,'BinWidth',0.5, 'Normalization','pdf');
hold on; plot(xx, target(xx)/max(target(xx))*max(h.Values))

참고 문헌

  • An introduction to MCMC for Maching Learning / C. Andrieu et al., Machine Learning, 50, 5-43, 2003