표본 샘플링 방법(3) - 깁스 샘플링(Gibbs Sampling)

이번 포스팅에서는 베이지안에서 연구되었던 샘플링 방법 중 가장 이상적인(?) 샘플링 방법인 깁스 샘플링에 대해서 살펴보겠다.

---

- 이전 포스팅-

2022.03.14 - [Statistics/Bayesian With Python] - 표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)

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

2022.03.16 - [Statistics/Bayesian With Python] - 표본 샘플링 방법(2) - 메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm)

표본 샘플링 방법(2) - 메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm)

---

■ 깁스 샘플링(Gibbs Sampling)이란

 

깁스 샘플링 방법은 메트로폴리스-헤이스팅스 알고리즘처럼 MCMC에 기반한 샘플링 방법이다.

누차 강조했던 대로 베이지안 접근법에서는 모수 $\theta$에 대한 추론을 하기위해 사후 분포 $p(\theta|x)$의 형태를 아는 것이 중요하다. 하지만 대부분의 경우, 사후 분포가 잘 알려진 형태로 나타나지 않아 제안 분포(proposal distribution)을 사용하여 사후 분포를 추정 하였다.

 

이번에 살펴볼 깁스 샘플링의 경우 샘플링을 위한 제안 분포(proposal density)가 필요하지 않다는 점에서 큰 차이가 존재한다.

깁스 샘플링은 모수에 대한 사후 조건부 분포가 잘 알려진 형태로 나타날 경우, 사후 조건부 분포에서의 샘플링을 통해 사후 분포를 추정할 수 있게한다.

---

예를 들어 살펴보자.

표본 $\boldsymbol{X}$가 $N(\mu, \sigma^2)$ 분포에서 생성된 n개의 샘플이라고 가정하자.

$\mu$와 $\sigma^2$ 모두 Unknown parameter이고, $\mu$의 사전 분포를 $N(0,\tau^2)$를 설정하고, $\sigma^2$의 사전 분포를 $Inv-Gamma(\alpha,\beta)$라고 설정하였을 때, 모수의 사후 분포 $p(\mu,\sigma^2|\boldsymbol{X})$는 잘 알려진 분포의 형태가 아니다.

$$p(\mu,\sigma^2|\boldsymbol{X}) \propto p(\boldsymbol{X}|\mu,\sigma)p(\mu)p(\sigma^2)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$$

$$~~~~~~~~~~~~~~~~~~~~~~~\propto exp\{-\cfrac{n}{2\sigma^2}\Sigma_{i=1}^n(X_{i}-\mu)^2-\cfrac{\mu^2}{2\tau^2}\}(\sigma^2)^{-(\cfrac{n}{2}+\alpha)-1}exp\{-\cfrac{\beta}{\sigma^2}\}$$

 

하지만, $\mu$와 $\sigma^2$ 각각의 조건부 분포를 살펴보면 어떨까?

 

먼저 $\mu$의 조건부 분포 형태를 살펴 보자 ($p(\mu,\sigma^2|\boldsymbol{X})$ 에서 $\mu$제외 모두 상수 취급, drop out)

 

$p(\mu|\sigma^2,\boldsymbol{X})  \propto \exp\{-\cfrac{n}{2\sigma^2}(\mu^2-2n\mu\bar{X})-\cfrac{\mu^2}{2\tau^2}\}$

                         $\propto \exp\{-\cfrac{1}{2}[(\cfrac{n}{\sigma^2}+\cfrac{1}{\tau^2})\mu^2-2\mu\cfrac{n}{\sigma^2}(n\bar{X})]\}$

                         $\propto \exp\{-\cfrac{1}{2}(\cfrac{n}{\sigma^2}+\cfrac{1}{\tau^2})(\mu-\cfrac{\cfrac{n^2\bar{X}}{\sigma^2}}{(\cfrac{n}{\sigma^2}+\cfrac{1}{\tau^2})})^2\}$

 $p(\mu|\sigma^2\boldsymbol{X})$의 분포를 proportional하게 유도하면 정규분포의 exponential form이 나타난다.

따라서, $\mu$의 조건부 분포는 평균이 $\cfrac{\cfrac{n^2\bar{X}}{\sigma^2}}{(\cfrac{n}{\sigma^2}+\cfrac{1}{\tau^2})}$이고 분산이 $\cfrac{1}{(\cfrac{n}{\sigma^2}+\cfrac{1}{\tau^2})}$인 Noraml 분포인 것을 알 수 있다.

 

그렇다면 $\sigma^2$의 조건부 분포는 어떨까?

$p(\sigma^2| \mu,\boldsymbol{X})  \propto (\sigma^2)^{-(\cfrac{1}{2}+\alpha)-1}\exp\{-\cfrac{1}{\sigma}[\cfrac{n}{2}\Sigma_{i=1}^n(X_i - \mu)^2+\beta]\}$

 

$p(\sigma^2| \mu,\boldsymbol{X})$의 유도된 형태를 보면 바로 $Inverse~-~Gamma~distribution$의 형태임을 알 수 있다.

$$p(\sigma^2| \mu,\boldsymbol{X}) = Inv~-~Gamma((\cfrac{n}{2}+\alpha),\cfrac{1}{2}\Sigma_{i=1}^n(X_i - \mu)^2+\beta)$$

$$(\therefore x \sim Inv~-~Gamma(\alpha,\beta),~~p(x) = \cfrac{\beta^\alpha}{\Gamma(\alpha)}x^{-\alpha-1}exp\left(-\cfrac{\beta}{x}\right))$$

 

이렇듯 관심있는 모수의 Joint posterior distribution은 샘플링이 불가능하지만, 각각의 모수의 조건부 분포는 잘 알려진 분포로 수렴할 경우 사용할 수 있는 샘플링 방법이 Gibbs Sampling이다.

---

M-H Algorithm과 Gibbs Sampling 모두 베이지안 접근법에서 자주 사용되는 Sampling Algorithm이지만,

Gibbs Sampling의 경우 샘플링 되는 모수의 후보(candidate)가 서로 독립이라는 장점이 있다.

 

M-H Algorithm의 경우 제안 분포가 $p(\theta^{(t)}|\theta^{(t-1)})$의 꼴로, 이전 시점(t-1)의 후보가 바로 다음 시점(t) 표본을 추출할 때 관혀하기 때문에 가까운 시점의 표본은 독립일 수가 없다.

따라서, M-H Algorithm에서 독립 표본을 얻고 싶은 경우엔 burn - in period 이후에 일정한 텀을 두고 추출된 표본만을 사후 표본으로 고려하기도 한다.

 

하지만 Gibbs Sampling의 경우, (t) 시점의 모수는 이전 시점(t-1) 표본에 의지하는 것이 아닌 각각의 조건부 분포에서 iid하게 추출되기 때문에 그 자체로도 독립이 만족하게 된다.

---

■ Gibbs Sampling Algorithm

Step 1. 초기 값 $\mu^{(0)},\sigma^{2(0)}$을 설정한다.

 

    For iteration $t=1,\cdots,T$ : 

          Step 2. $\mu^{(t)}$를 $p(\mu|\sigma^{2(t-1)},X)$으로부터 샘플링한다.

          Step 3. $\sigma^{2(t)}를 $p( \sigma^2 | \mu^{(t)},X)$으로부터 샘플링한다.

 

Step 4. burn-in period $t_{0}$ 이전에 얻어진 후보들은 drop-out하고, $t_{0}$이후에 얻어진 후보들만 표본으로 간주한다.

---

■ 예제

위에서 보였던 예제를 사용하여 깁스 샘플링이 실제로도 잘 작동하는지 확인해보자.

우리가 추정하고자하는 unknow parameter $\mu$와 $\sigma^2$를 각각 3과 2로 설정하여 표본 $X_i,i=1,\cdots,10000$를 먼저 생성하였다.

$X_i \sim N(3,2)$

parameter $\mu$와 $\sigma$의 사전 분포는 각각 아래와같이 설정하였다.

$$p(\mu) = N(0,2) ~~~~~~~~~ p(\sigma^2)=Inv - Gamma(1,1)$$

 

충분한 수렴을 위해 총 iteration을 $T=10,000$으로 설정하고 깁스 샘플링을 수행해보자.

(left) Trace plot for $\mu$   (right) Trace plot for $\sigma$

Trace plot과 깁스 샘플링을 통해 추정한 posterior estimator값을 통해 샘플링의 결과가 안정적으로 수렴하면서 실제 값도 아주 근사하게 추정한 것을 확인 할 수 있다.

 

예제 코드 더보기

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as stats
from scipy.stats import invgamma

n = 10000
T = 10000
burn = T/2

## Sample X
X = np.random.normal(3,np.sqrt(2),n)

plt.hist(X,density=True,label = 'Histogram of X')

mu_candi = np.zeros(shape=T)
sigma_candi= np.zeros(shape=T)

# Set initial value
mu_candi[0]=0
sigma_candi[0]=1

# Prior distribution of parameters
tau = 2
alpha=1
beta=1

# running Gibbs Sampling
for t in range(1,T):
    post_var = 1/(n/sigma_candi[t-1]+1/tau)
    post_mean = post_var*(n*np.mean(X)/sigma_candi[t-1])
    mu_candi[t] = np.random.normal(post_mean,np.sqrt(post_var),1)
    
    shape = n/2+alpha
    scale = 1/2*np.sum((X-mu_candi[t])**2)+beta
    sigma_candi[t] = 1/np.random.gamma(shape=shape,scale=1/scale)
    
plt.plot(mu_candi[:])
plt.plot(sigma_candi[:])

print(r"Posterior estimator of mean :",np.mean(mu_candi[int(burn):]))
print(r"Posterior estimator of variance :",np.mean(sigma_candi[int(burn):]))