Acceptance-Rejection Sampling (허용-기각 샘플링)

Acceptance-Rejection Sampling

Inverse Transform Method와 마찬가지로 확률분포에서 샘플링을 하기 위한 기법 중 하나이며 제안 분포 (Proposal Distribution)를 이용하여 생성된 샘플들의 분포가 목표 분포 (Target Distribution)를 따르도록 수정하는 과정을 포함한 샘플링입니다.

목표 분포의 특징은 직접 샘플을 생성하는 게 어렵거나 불가능하다는 특징이 있습니다. Acceptance-Rejection Sampling을 적용하여 목표 분포에 대해서 샘플링을 하기 위해서는 우선 확률 밀도 함수 (PDF)를 알고 있어야 합니다.

 

Motivation

Inverse Transform Method와 같이 CDF가 역함수가 가능한 상황에서 적용할 수 있는 기법과 다르게 Acceptance Rejection Sampling은 다양한 상황에서 적용이 가능한 방법입니다 (General Purpose Method).

 

Proposal Distribution

Rejection Sampling에서 제안분포는 쉽게 샘플을 생성할 수 있도록 하는 임의로 설정된 분포입니다. 목표 분포와 형태가 유사할수록 효과적인 샘플링이 가능해집니다. 보통 계산이 간단한 특징을 지닌 Uniform Distribution을 많이 사용하기도 합니다.

우선 하나의 예시로 설명을 드려보겠습니다. 아래 사진과 같이 검은색의 목표 분포가 있다고 하였을 때 우리는 샘플링을 위해서 Uniform Distribution을 가진 초록색 분포로 접근해 보겠습니다.

식으로 표현을 하면 다음과 같습니다.

 

Proposal Distribution을 제대로 이용하기 위해서는 상수배 (M)을 취해주어서 y값이 모두 Target Distribution을 모두 포함하도록 만들어줘야 합니다.

상수배를 10으로 해줬으며 1/10 * 10 = 1이므로 모든 영역을 다 커버할 수 있었습니다.

 

Sampling Process

다음은 샘플링 과정입니다. 가장 먼저 해줘야 할 일은 Proposal Distribution에서 샘플 하나 (x축 값, -4 ~ 6)를 추출하는 것입니다 (이때 샘플링을 위해서 파이썬 함수 uniform(-4,6) 구현할 수 있습니다)

추출이 완료되었으면 목표 분포 T와 상수배를 취한 제안 분포 G의 Likelihood를 비교해 주기 위해서 목표분포와 제안 분포에 각각 샘플값이 대입된 값을 계산하고 목표분포값 나누기 제안 분포값을 합니다.

A값이 1에 가깝다는 것은 타깃 분포 값이 제안 분포에 상수배를 한 값과 가까운 값일 것이며, 멀다는 것은 제안 분포에서 나온 값이 타겟 분포보다 훨씬 작다는 의미를 가집니다. 

A는 Uniform Distribution(0,1)에서 나온 u값과 비교하여 만약 크다면 Accept을 수행하며

작다면 Reject을 수행합니다.

 

Example

import numpy as np
import matplotlib.pyplot as plt

def target_distribution(x):
    return 0.3 * np.exp(-x**2/2) + 0.7 * np.exp(-(x-3)**2/2)

def proposal_distribution(x):
    if  x >= -4 and x <= 6:
        return 1/10
    else: 
        return 0
    #return np.exp(-x**2/2)

M = 10 # upper bound of the target distribution divided by the proposal distribution
N = 600 # number of samples

samples = []
sam_prob = []
rejected_samples = []
x_reject = []

while len(samples) < N:
    x = np.random.uniform(-4, 6)
    acceptance_prob = target_distribution(x) / (proposal_distribution(x)* M)

    u = np.random.uniform(0,1)
    if u < acceptance_prob:
        samples.append(x)
        sam_prob.append(u)
    else:
        rejected_samples.append(x)
        x_reject.append(u)
print("Acceptance rate: ", len(samples) / (len(samples) + len(rejected_samples)))

#plt.hist(samples, bins=50, density=True, alpha=0.9, label="Samples")
plt.scatter(samples, sam_prob, label="Accepted")
plt.scatter(rejected_samples, x_reject,color ='r', label="Rejected")

x = np.linspace(-4, 6, 100)
plt.plot(x, target_distribution(x), 'k-', label="Target distribution")
plt.plot(x, [1/10 * M] * len(x), 'g', label="Proposal distribution", linewidth=3)
plt.vlines(-4, 0,1, color='green', linewidth=3)
plt.vlines(6, 0, 1, color='green',  linewidth=3)
plt.hlines(0, -5, -4, color='green',  linewidth=3)
plt.hlines(0, 6, 7, color='green',  linewidth=3)
plt.legend()
plt.xlim(-5,7)
plt.ylim(-0.05,1.4)
plt.savefig("AR.png",dpi = 300)
plt.show()

 

Size of M

상수배를 수행하기 위해서 곱해지는 M은 얼마가 적당할까요? 

M은 곧 알고리즘의 효율성과 깊은 관련이 있으며 M이 작으면 Acceptace rate가 낮게 나오며 알고리즘이 비 효과적으로 수행됩니다. 만약 M이 너무 크다면 제안 분포가 너무나도 커져서 Acceptace rate가 매우 작아지게 됩니다.

M을 알아내기 위해서 한 가지 방법은 목표 분포의 upper bound를 제안분포로 나누어서 최대 비율을 찾는 방법이 있습니다. 또는 M을 작은 수부터 증가시켜서 점점 커지게 하여서 acceptance rate가 10~50% (합리적인 구간) 정도로 안착되었을 때 M을 선정하는 방법입니다.

acceptance rate와 M 간의 관계는 반비례이기 때문에 알고리즘의 속도를 빠르게 만들기 위해서는 목표 함수 <= M * 제안 분포를 만족하면서 가능한 작은 M값을 선정해야 합니다.