표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)

Bayesian과 Frequentist들의 가장 큰 차이점은 모수, parameter에 대한 관점의 차이이다.

 

Frequentist들은 모수 $\theta$를 Unknown constant로 가정하며,

Bayesian은 모수 $\theta$를 Random variable로 가정한다.

때문에, 대표적인 Frequentist의 방법인 MLE처럼 표본 $x$에 의해 값(ex, $\hat{\theta} = \cfrac{\sum_{i=1}^{N}x_i}{N}$)으로 계산되는 것과 달리, Bayesian에서 모수 $\theta$는 확률 분포를 갖게된다.

 

따라서, Bayesian 접근법에서는 모수를 추정하기 위해 모수의 사후 분포(posterior distribution)에서 표본을 생성하게 된다. 그러나 종종 모수의 사후 분포는 복잡한 형태를 갖게되며(잘 알려진 분포의 형태로 나타나지 않게 된다.), 이와 같은 경우 사후 분포에서 바로 표본을 샘플링 할 수 없게된다.

 

베이지안에서는 이러한 문제를 극복하기 위해, 여러가지 샘플링 알고리즘이 고안되었으며, 오늘은 그 중 첫번째로 기각 샘플링(rejection sampling)방법을 알아보겠다.

---

■ 기각 샘플링 - Rejection Sampling

기각 샘플링은 위에서 설명했던 것 처럼 표본을 생성하기 위해 사용되는 알고리즘이며, 'Acceptance-rejection algorithm'이라고 불리기도 한다.

 

기각 샘플링 방법은 우리가 알고싶은, 표본을 샘플링 하고싶은 분포를 $p(x)$ 타겟 분포(target density/distribution)가 있을 때,$p(x)$에서 직접적으로 표본 $x$를 샘플링 할 수 없는 경우,

표본을 쉽게 생성할 수 있는 $q(x)$라는 제안분포(proposal density/distribution)을 도입하여 일정한 조건을 통해 $\boldsymbol{q(x)}$에서 생성된 표본 $x_1,x_2,\cdots,x_N$ 을 $\boldsymbol{p(x)}$에서 생성된 표본으로 간주 할 수 있도록 하는 알고리즘이다.

 

일반적으로 $q(x)$는 표본 생성이 쉬운 Uniform distribution이나, Normal distribution같은 분포들이 많이 사용되며, 아래와 같은 조건이 만족하도록 상수 $M$이 설정된다.

 

$$p(x)\leq M * q(x)$$

 

예를 들어, 아래와 같은 복잡한 형태의 타겟 분포 $p(x)$가 있다고 하자.

 

$$p(x) = \begin{cases} 8x & 0\leq x < 0.25 \\ \cfrac{8}{3} - 8\cfrac{x}{3} & 0.25 \leq x \leq 1\\ 0 & \text{otherwise} \end{cases} $$

 

이 때, 쉬운 표본 생성을 위해 제안 분포 $q(x)$를 $Unif(0,1)$로 설정하고, $q(x)$이 $p(x)$를 커버하도록 $M=3$으로 설정할 수 있다.

orange line - q(x) , blue line - p(x)

그렇다면, 우리는 어떤 조건을 거쳐 $q(x)$에서 샘플링된 표본을 $p(x)$의 표본이라고 여길 수 있는 것일까?

기각 샘플링의 본격적인 의사결정과정을 살펴보자.

 

아래의 과정은 총 $T$개의 표본이 $p(x)$의 표본으로 accept될 때까지 반복 된다.

repeat until we get $x^{(1)},x^{(2)},\cdots,x^{(T)}$

 

1. $q(x)$에서 candidate $x^{*}$를 샘플링한다. : $x^{*}\sim q(x)$

   - $q(x)$에서 뽑힌 모든 $x$들이 모두 $p(x)$의 표본이 되는 것은 아니기 때문에, 여기서는 candidate이라는 표현을 쓴다.

 

2. candidate $x^{*}$에 대한 acceptance rate을 계산한다. : $\alpha = \cfrac{p(x^{*})}{Mq(x^{*})}$

   - 위에서 언급한 $p(x)\leq M * q(x)$ 이 조건에 의해 $0\leq\alpha\leq 1$을 만족하며, $\alpha$는 $x^{*}$가 $Mq(x^{*})$에 비해 $p(x^{*})$에 속할 정도,비율을 나타낸다.

 

만약 $q(x)=Unif(0,1)$에서 생성된 candidate $x^{*}$가 0.3이라고하자.

이 때, acceptance rate은 (하늘색 solid line/노란 dashed line)으로 정의되며, 두 선의 높이가 비슷할 수록 acceptance rate이 커진다.

 

3. Random variable $U$를 $Unif(0,1)$으로부터 뽑은 뒤, acceptance rate $\alpha$와 비교를 하게된다.

   - 만약 $U\leq \alpha$라면 candidate $x^{*}$를 $x^{(t)}$로 설정한다. 즉, $p(x)$의 표본으로 Accept하는 것이다.

   - 만약 $U>\alpha$라면, candidate $x^{*}$를 Reject하고 새로운 candidate $x^{*}$를 뽑아서, candidate이 accept될 때까지 1~3 단계를 반복한다.

 

기각 샘플링 알고리즘을 pseudo code로 간단히 다시 살펴보면 아래와 같다.

 

Algorithm1 : Rejection Sampling

Input : $p(x),q(x),M,T$   while $i\leq T$ :            repeat :               1. Draw candidate $x^{*}$ from $q(x)$               2. Compute acceptance rate $\alpha=\cfrac{p(x^{*})}{Mq(x^{*})}$ and draw $U$ from $Unif(0,1)$               If $U\leq \alpha$ : $x^{(i)}=x^{*}$ and break   end while   output : $\boldsymbol{x}=\{x^{(1)},\cdots,x^{(t)}\}$

 

위에서 살펴본 예로 기각 샘플링을 수행한 결과는 다음과 같다.

---

---

기각 샘플링 결과를 살펴보면, 더 많은 표본 수를 얻을 수록($T=10000) 타겟 분포와 더 유사한 근사 결과를 얻은 것을 확인 할 수 있다.

 

본문 내용 작성을 위해 사용된 python 코드는 아래 더보기 에서 확인 할 수 있다.

 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

####### Define density functions #######
def p(x) :
    if x>=0 and x<0.25 : 
        val = 8*x
    elif x>=0.25 and x<=1:
        val = 8/3 - 8*x/3
    else:
        val = 0
    return(val)
    
def q(x) :
    if x>=0 and x<=1 :
        val=1
    else : 
        val=0
    return(val)

x=np.arange(0,1,0.01)
y1=[];y2=[]
for i in range(0,len(x)):
    y1=np.append(y1,p(x[i]))
    y2=np.append(y2,q(x[i]))

####### Draw figure 1 #######
M=3

plt.plot(x,y1)
plt.plot(x,M*y2)
plt.ylim(0,3.5)

####### Draw figure 2 #######
M=3

plt.plot(x,y1)
plt.plot(x,M*y2)
plt.ylim(0,3.5)

# probability of p(x) for candidate
plt.vlines(0.3,0,p(0.3),color='skyblue')
plt.scatter(0.3,p(0.3),color='skyblue',s=50)

# probability of q(x) for candidate
# 살짝 빗겨나게 그림.
plt.vlines(0.31,0,M*q(0.30),color='y',linestyle='--')
plt.scatter(0.31,M*q(0.30),color='y',s=50)

####### Rejection sampling #######

T=10000
X=[]

while(len(X)<T):
    candi_x = np.random.uniform(size=1)
    alpha = p(candi_x)/(M*q(candi_x))
    
    U = np.random.uniform(size=1)
    if(U<=alpha):
        X=np.append(X,candi_x)
        
import seaborn as sns

plt.plot(x,y1,color='blue',linewidth=2)
plt.ylim(0,2.3)
sns.distplot(X,color='red')
plt.title('Sampling result : T=10000')