이번 포스팅에서는 베이지안에서 연구되었던 샘플링 방법 중 가장 이상적인(?) 샘플링 방법인 깁스 샘플링에 대해서 살펴보겠다.
- 이전 포스팅-
2022.03.14 - [Statistics/Bayesian With Python] - 표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)
표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)
Bayesian과 Frequentist들의 가장 큰 차이점은 모수, parameter에 대한 관점의 차이이다. Frequentist들은 모수 $\theta$를 Unknown constant로 가정하며, Bayesian은 모수 $\theta$를 Random variable로 가정한다..
harang3418.tistory.com
표본 샘플링 방법(2) - 메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm)
저번 rejection sampling에 이어 이번 포스팅에서도 타겟 분포로부터 표본을 생성해낼 수 있는 샘플링 알고리즘 중 하나인 M-H Algorithm을 살펴보겠다. 2022.03.14 - [Statistics/Bayesian With Python] - 표본..
harang3418.tistory.com
■ 깁스 샘플링(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$를 먼저 생성하였다.
parameter $\mu$와 $\sigma$의 사전 분포는 각각 아래와같이 설정하였다.
$$p(\mu) = N(0,2) ~~~~~~~~~ p(\sigma^2)=Inv - Gamma(1,1)$$
충분한 수렴을 위해 총 iteration을 $T=10,000$으로 설정하고 깁스 샘플링을 수행해보자.
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):]))
'Statistics > Bayesian With Python' 카테고리의 다른 글
표본 샘플링 방법(2) - 메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm) (0) | 2022.03.16 |
---|---|
표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling) (0) | 2022.03.14 |
다양한 사전 분포 _ Prior distribution (0) | 2021.01.29 |
베이지안 모델링의 기본 (0) | 2021.01.13 |
댓글