16회 · 베이즈의 귀환 · 노트

4 MCMC — 사후분포에서 뽑기

역사 1953년 Los Alamos의 Nicholas Metropolis, Arianna와 Marshall Rosenbluth, Augusta와 Edward Teller는 MANIAC 컴퓨터로 단단한 공 224개로 된 액체의 상태방정식을 계산했다. $224\times2$차원의 적분을 7회의 단순 몬테카를로로는 할 수 없었다. 대부분의 배치에서 공이 겹치기 때문이다. 그들은 공 하나를 조금 옮겨 보고, 에너지가 낮아지면 받아들이고 높아지면 확률적으로만 받아들였다. 이렇게 만든 마르코프 연쇄는 볼츠만 분포를 정상분포로 가졌다. 물리학의 도구였던 이 방법은 Hastings(1970)가 일반화하고, Geman 형제(1984)가 영상 복원에 쓴 Gibbs 표집을 거쳐, Gelfand와 Smith(1990)의 논문으로 통계학에 들어왔다. 정규화 상수 $m(\x)$를 몰라도 사후분포에서 표본을 뽑을 수 있다는 것이 열쇠였다. 1990년대 BUGS 소프트웨어와 함께 베이즈 추론은 손 계산의 한계에서 풀려났다.

4.1 마르코프 연쇄와 정상분포

정의 4.1 (전이핵과 정상분포) 상태공간 위의 수열 $X_0,X_1,\dots$에서 $X_{t+1}$의 조건부분포가 $X_t$에만 달려 있고 그 밀도가 $k(x,y)$(전이핵)이면 마르코프 연쇄라 한다. 밀도 $\pi$가

$$\int\pi(x)k(x,y)\,dx=\pi(y)\qquad\text{모든 }y$$

를 만족하면 $\pi$를 정상분포라 한다. $X_t\sim\pi$이면 $X_{t+1}\sim\pi$다. 전이핵이 받아들이지 않는 확률만큼 제자리에 머무는 덩어리를 가질 수 있으며, 그때 적분은 그 덩어리를 포함한다.

정의 4.2 (상세 균형) $\pi(x)k(x,y)=\pi(y)k(y,x)$가 모든 $x,y$에서 성립하면 연쇄가 $\pi$에 대해 상세 균형(가역성)을 만족한다고 한다.

정리 4.3 (상세 균형이면 정상분포) 상세 균형을 만족하는 $\pi$는 정상분포다.

증명 양변을 $x$로 적분하면 $\int\pi(x)k(x,\allowbreak y)dx=\pi(y)\int k(y,\allowbreak x)dx=\pi(y)$.

상세 균형은 “$x$에서 $y$로 가는 확률의 흐름과 $y$에서 $x$로 오는 흐름이 같다”는 뜻이다. 정상성은 들어오는 흐름의 합과 나가는 흐름의 합만 같으면 되므로, 상세 균형은 충분조건일 뿐 필요조건은 아니다(연습 4.2). 그러나 국소적인 조건이라 확인하기 쉽고, 다음의 알고리즘은 그것을 만족하도록 설계된다.

4.2 Metropolis–Hastings와 Gibbs

알고리즘 4.4 (Metropolis–Hastings) 목표 밀도 $\pi$(상수배까지만 알면 된다)와 제안 밀도 $q(x,\cdot)$가 주어진다. 현재 상태 $x$에서

1. 1.

   $y\sim q(x,\cdot)$를 제안한다.

2. 2.

   확률 $\alpha(x,\allowbreak y)=\min\Bigl\{1,\allowbreak \ \dfrac{\pi(y)q(y,x)}{\pi(x)q(x,y)}\Bigr\}$로 $y$로 옮기고, 아니면 $x$에 머문다.

$q$가 대칭($q(x,y)=q(y,x)$, 예컨대 $y=x+\sigma\epsilon$, $\epsilon\sim\N(0,1)$)이면 $\alpha=\min\{1,\allowbreak \pi(y)/\pi(x)\}$이고 이것이 1953년의 Metropolis 알고리즘이다.

정리 4.5 (Metropolis–Hastings의 정상분포) 알고리즘 4.4의 전이핵은 $\pi$에 대해 상세 균형을 만족한다. 따라서 $\pi$는 정상분포다.

증명 전략 상세 균형을 직접 확인한다. $\pi(x)$를 받아들임 확률의 $\min$ 안으로 넣으면 $x$와 $y$에 대칭인 식이 된다.

1. $x\neq y$이면 $y$로 옮기려면 $y$를 제안하고 받아들여야 하므로 $k(x,y)=q(x,y)\alpha(x,y)$이다.알고리즘 4.4

2. 양수를 $\min$ 안으로 곱하면

   $$\pi(x)k(x,y)=\pi(x)q(x,y)\min\Bigl\{1,\frac{\pi(y)q(y,x)}{\pi(x)q(x,y)}\Bigr\}=\min\bigl\{\pi(x)q(x,y),\ \pi(y)q(y,x)\bigr\} .$$

   $c\min\{a,b\}=\min\{ca,cb\}$, $c>0$

3. 오른쪽은 $x$와 $y$를 바꿔도 같으므로 $\pi(x)k(x,y)=\pi(y)k(y,x)$이다.대칭

4. $x=y$의 제자리 덩어리는 양변에 같은 항으로 나타나므로 상세 균형을 어기지 않는다.정의 4.2

5. 따라서 $\pi$는 정상분포다.정리 4.3

조건 확인 2는 $\pi(x)q(x,y)>0$일 때의 계산이다. 연쇄가 $\pi(x)>0$인 점에서 출발하면 그 뒤로도 $\pi>0$인 점에만 머문다($\pi(y)=0$인 제안은 $\alpha=0$으로 거절된다). $\pi(x)q(x,y)=0$이고 $\pi(y)q(y,x)>0$인 쌍에서는 $\alpha(y,x)=0$이라 양변이 모두 0이다. $\pi$를 상수배까지만 알아도 $\alpha$는 비만 쓰므로 계산이 바뀌지 않는다.

$\alpha$에는 $\pi(y)/\pi(x)$만 들어 있다. 사후분포 $\pi(\theta\mid\x)=\pi(\theta)f(\x\mid\theta)/m(\x)$에서 $m(\x)$는 비에서 지워진다. 1절에서 Laplace를 막았던 적분이 필요 없다.

알고리즘 4.6 (Gibbs 표집) $\theta=(\theta_1,\allowbreak \dots,\allowbreak \theta_d)$이고 각 완전 조건부분포 $\pi(\theta_j\mid\theta_{-j})$에서 뽑을 수 있다고 하자. 한 바퀴에서 $j=1,\dots,d$ 차례로 $\theta_j\sim\pi(\theta_j\mid\theta_{-j})$로 바꾼다($\theta_{-j}$는 가장 최근 값).

정리 4.7 (Gibbs는 받아들임 확률 1의 MH) 좌표 $j$를 완전 조건부분포로 바꾸는 단계는 제안 $q(\theta,\allowbreak \theta^\prime)=\pi(\theta^\prime_j\mid\theta_{-j})$($\theta^\prime_{-j}=\theta_{-j}$)의 MH이고 받아들임 확률은 1이다. 따라서 각 단계와 그 합성(한 바퀴)은 $\pi$를 정상분포로 갖는다.

증명 전략 결합밀도를 $\pi(\theta)=\pi(\theta_{-j})\pi(\theta_j\mid\theta_{-j})$로 쪼개 MH 비에 넣으면 모든 인수가 약분된다. 합성은 “$\pi$를 보존하는 핵을 차례로 쓰면 $\pi$가 보존된다”로 처리한다. 단계별 증명은 연습 4.6(유도형)에서 직접 쓴다. 완성된 증명은 그 연습의 풀이에 있다.

한 바퀴의 합성은 일반적으로 가역이 아니다. 그래도 정상분포는 지켜진다. 켤레 구조가 있는 계층모형에서는 완전 조건부분포가 정리 1.2의 꼴이 되어 Gibbs 표집이 자연스럽다(6절).

4.3 수렴과 진단

증명 없이 쓰는 사실 (마르코프 연쇄의 에르고딕 정리) 연쇄가 $\pi$를 정상분포로 갖고, $\pi$-기약(양의 확률의 어느 집합에나 언젠가 들어간다)이고 비주기적이면, 거의 모든 출발점에서 $\frac1N\sum_{t=1}^Ng(X_t)\to E_\pi g$ (확률 1)이다. 추가 조건 아래 중심극한정리도 성립하며

$$\mathrm{Var}\Bigl(\frac1N\sum_{t=1}^Ng(X_t)\Bigr)\approx\frac{\sigma_g^2}{N}\Bigl(1+2\sum_{k\ge1}\rho_k\Bigr),$$

$\rho_k$는 정상 상태에서 $g(X_t)$의 시차 $k$ 자기상관이다. 7회 정리 6.1의 $\sigma^2/N$에서 $N$을 $N_{\mathrm{eff}}=N/(1+2\sum\rho_k)$로 바꾼 것이며, $N_{\mathrm{eff}}$를 유효 표본 크기(ESS)라 한다. 3회 예제 6.6의 AR(1)에서 본 $n_{\text{eff}}\approx n(1-\rho)/(1+\rho)$가 그 특수한 경우다.

명제 4.8 (상관된 정규에서의 Gibbs) $(X,Y)$가 평균 0, 분산 1, 상관 $\rho$의 이변량 정규이면 Gibbs 표집 $X_t\sim\N(\rho Y_{t-1},\allowbreak 1-\rho^2)$, $Y_t\sim\N(\rho X_t,1-\rho^2)$의 $X$열은 계수 $\rho^2$의 AR(1)

$$X_t=\rho^2X_{t-1}+\eta_t,\qquad\eta_t\iid\N(0,1-\rho^4)$$

이다. 따라서 $\rho_k=\rho^{2k}$이고 $N_{\mathrm{eff}}/N=(1-\rho^2)/(1+\rho^2)$이다.

증명 전략 두 조건부 정규 걸음을 잡음 항으로 쓰고 하나를 다른 하나에 대입한다. 그러면 $X$열만의 재귀식이 나오고, AR(1)의 자기상관을 더한다(3회 예제 6.6). 단계별 증명은 연습 4.7(유도형)에서 직접 쓴다. 완성된 증명은 그 연습의 풀이에 있다.

$\rho=0.5,\ 0.9,\ 0.99$이면 $N_{\mathrm{eff}}/N=0.60,\allowbreak \ 0.105,\allowbreak \ 0.010$이다. 좌표축 방향으로만 움직이는 Gibbs는 대각선으로 길쭉한 분포에서 계단을 잘게 밟는다(그림 4.1). 모수의 상관을 줄이는 재모수화(중심화)가 MCMC의 효율을 좌우하는 이유다.

연쇄가 정상분포에 도달했는지는 한 연쇄만 보아서는 알 수 없다. Gelman과 Rubin(1992)은 흩어진 출발점에서 여러 연쇄를 돌려 연쇄 사이의 분산과 연쇄 안의 분산을 비교했다. 14회 분산분석의 논리 그대로다.

정의 4.9 ($\hat R$ (잠재적 축소 인자)) 길이 $n$의 연쇄 $m$개에서 연쇄 $j$의 평균을 $\bar\theta_j$, 분산을 $s_j^2$이라 하고 $W=\frac1m\sum s_j^2$, $B=\frac{n}{m-1}\sum(\bar\theta_j-\bar\theta)^2$이라 하자.

$$\hat V=\frac{n-1}nW+\frac1nB,\qquad\hat R=\sqrt{\hat V/W}.$$

모든 연쇄가 같은 분포에서 뽑히고 있으면 $\hat V$와 $W$는 모두 목표분포의 분산을 추정하므로 $\hat R\approx1$이다. 출발점의 영향이 남아 연쇄들이 서로 다른 곳을 돌고 있으면 $B$가 커서 $\hat R>1$이다. 실무에서는 $\hat R<1.01$을 기준으로 쓴다. $\hat R\approx1$은 필요조건일 뿐이다. 모든 연쇄가 같은 봉우리에 갇혀 있으면 $\hat R$도 속는다(7절).

그림 4.1 왼쪽: $\rho=0.95$인 이변량 정규에서 Gibbs 표집의 처음 40걸음. 좌표축 방향으로만 움직인다. 오른쪽: $X$열의 자기상관(점)과 명제 4.8의 $\rho^{2k}$(선). $\rho$가 1에 가까울수록 연쇄가 천천히 잊는다.

연습문제

연습 4.1 (두 상태 연쇄 ★) 전이행렬 $\bigl(\begin{smallmatrix}1-a&a\\b&1-b\end{smallmatrix}\bigr)$($0<a,b<1$)의 정상분포를 구하고 상세 균형을 확인하라. $P^n$이 정상분포로 가는 속도는 $|1-a-b|^n$임을 보여라.

풀이

$\pi P=\pi$에서 $\pi_1a=\pi_2b$이므로 $\pi=\bigl(\frac b{a+b},\allowbreak \frac a{a+b}\bigr)$이고, 이 식이 바로 상세 균형이다. $P$의 고윳값은 1과 $1-a-b$이므로 $P^n=\mathbf 1\pi+(1-a-b)^nR$ (고정된 행렬 $R$)이다.

연습 4.2 (가역이 아닌 정상분포 ★) 상태 $\{1,2,3\}$에서 $1\to2\to3\to1$을 확률 $0.9$로, 반대 방향을 $0.1$로 움직이는 연쇄의 정상분포가 균등임을 보이고, 상세 균형을 만족하지 않음을 보여라.

풀이

각 열의 합이 $0.9+0.1=1$이므로 전이행렬이 이중확률행렬이고 균등분포가 정상분포다. 그러나 $\pi_1k(1,\allowbreak 2)=0.3\neq\pi_2k(2,\allowbreak 1)=0.0\overline3$이다. 확률이 한 방향으로 순환하므로 가역이 아니다.

연습 4.3 (독립 제안 ★★) $q(x,y)=g(y)$(현재 상태와 무관)이면 MH 받아들임 확률이 $\min\{1,w(y)/w(x)\}$, $w=\pi/g$임을 보여라. 7회의 중요도 표집의 가중치와 비교하라. $g$의 꼬리가 $\pi$보다 가벼우면 어떤 일이 생기는가?

풀이

MH 비는 $\frac{\pi(y)g(x)}{\pi(x)g(y)}=\frac{w(y)}{w(x)}$이다. $w$는 7회의 중요도 가중치다. 중요도 표집이 가중치로 평균을 고친다면, 독립 MH는 가중치가 큰 점에 더 오래 머무는 방식으로 고친다. $g$의 꼬리가 가벼우면 $w$가 꼬리에서 한없이 커지고, 연쇄가 그런 점에 들어가면 매우 오래 갇힌다. 정상분포는 맞지만 수렴이 기하급수적이지 않다.

연습 4.4 ($\hat R$의 하한 ★) 연쇄 평균이 모두 같으면 $\hat R=\sqrt{(n-1)/n}<1$임을 보여라. $\hat R$이 1보다 큰 것은 무엇 때문인가?

풀이

$B=0$이면 $\hat V=\frac{n-1}nW$이고 $\hat R=\sqrt{(n-1)/n}$이다. $\hat R>1$이 되려면 $B/n>W/n$, 곧 연쇄 평균들의 흩어짐 $\frac1{m-1}\sum(\bar\theta_j-\bar\theta)^2$이 연쇄 안의 분산을 $n$으로 나눈 값 $W/n$보다 커야 한다. 연쇄들이 독립 표본이라면 연쇄 평균의 분산은 $W/n$ 정도이므로, 그보다 큰 흩어짐은 연쇄들이 아직 다른 곳에 있다는 증거다.

연습 4.5 (Gibbs의 무작위 순서 ★★) 매 단계 좌표 $j$를 균등하게 무작위로 골라 완전 조건부분포로 바꾸는 “무작위 순서” Gibbs의 전이핵은 $\pi$에 대해 상세 균형을 만족함을 보여라. (체계적 순서의 한 바퀴는 일반적으로 가역이 아니다.)

풀이

좌표 $j$의 갱신 핵 $K_j$는 받아들임 확률 1의 MH이므로 정리 4.5로 $\pi(\theta)K_j(\theta,\allowbreak \theta^\prime)=\pi(\theta^\prime)K_j(\theta^\prime,\allowbreak \theta)$이다. 무작위 순서의 핵은 $K=\frac1d\sum_jK_j$이고, 상세 균형은 선형이므로 합에도 성립한다. 체계적 순서 $K_1K_2\cdots K_d$의 수반은 $K_d\cdots K_1$이라 일반적으로 자기 자신과 다르다.

연습 4.6 (Gibbs는 MH: 증명 ★ 유도형) 정리 4.7를 다음 순서로 증명하라.

1. (a)

   $\theta^\prime_{-j}=\theta_{-j}$이면 $\pi(\theta)=\pi(\theta_{-j})\pi(\theta_j\mid\theta_{-j})$, $\pi(\theta^\prime)=\pi(\theta_{-j})\pi(\theta^\prime_j\mid\theta_{-j})$임을 조건부밀도의 정의로 보여라.

2. (b)

   제안 $q(\theta,\allowbreak \theta^\prime)=\pi(\theta^\prime_j\mid\theta_{-j})$와 $q(\theta^\prime,\allowbreak \theta)=\pi(\theta_j\mid\theta_{-j})$를 MH 비에 넣어 비가 1임을 보여라.

3. (c)

   정리 4.5로 각 좌표 단계가 $\pi$를 보존함을 결론짓고, 핵 $K_1$, $K_2$가 $\pi$를 보존하면 합성 $K_1K_2$도 $\pi$를 보존함을 적분의 순서를 바꿔 보여라.

풀이

(a) 조건부밀도의 정의 $\pi(\theta_j\mid\theta_{-j})=\pi(\theta)/\pi(\theta_{-j})$에서 $\pi(\theta)=\pi(\theta_{-j})\pi(\theta_j\mid\theta_{-j})$이다. $\theta^\prime$은 $\theta_{-j}$를 공유하므로 같은 식에서 $\pi(\theta^\prime)=\pi(\theta_{-j})\pi(\theta^\prime_j\mid\theta_{-j})$이다. (b) 제안은 좌표 $j$만 바꾸므로 $q(\theta,\allowbreak \theta^\prime)=\pi(\theta^\prime_j\mid\theta_{-j})$, 거꾸로 $q(\theta^\prime,\allowbreak \theta)=\pi(\theta_j\mid\theta_{-j})$이다. MH 비는

$$\frac{\pi(\theta^\prime)q(\theta^\prime,\theta)}{\pi(\theta)q(\theta,\theta^\prime)}=\frac{\pi(\theta_{-j})\pi(\theta^\prime_j\mid\theta_{-j})\cdot\pi(\theta_j\mid\theta_{-j})}{\pi(\theta_{-j})\pi(\theta_j\mid\theta_{-j})\cdot\pi(\theta^\prime_j\mid\theta_{-j})}=1$$

이므로 받아들임 확률은 $\min\{1,1\}=1$이다. (c) 정리 4.5로 각 $K_j$는 $\pi$에 대해 상세 균형을 만족하고, 정리 4.3로 $\pi K_j=\pi$이다. $\pi K_1=\pi$, $\pi K_2=\pi$이면 $\int\!\!\int\pi(x)K_1(x,\allowbreak y)K_2(y,\allowbreak z)\,dy\,dx=\int\Bigl[\int\pi(x)K_1(x,\allowbreak y)dx\Bigr]K_2(y,\allowbreak z)\,dy=\int\pi(y)K_2(y,\allowbreak z)\,dy=\pi(z)$이다. 피적분함수가 음이 아니므로 적분 순서를 바꿀 수 있다(Tonelli). 따라서 한 바퀴 $K_1\cdots K_d$도 $\pi$를 보존한다.

연습 4.7 (상관된 정규의 Gibbs: 증명 ★★ 유도형) 명제 4.8를 다음 순서로 증명하라.

1. (a)

   5회의 조건부 정규에서 $Y_{t-1}=\rho X_{t-1}+\sqrt{1-\rho^2}\,\epsilon_{t-1}$, $X_t=\rho Y_{t-1}+\sqrt{1-\rho^2}\,\epsilon^\prime_t$로 쓸 수 있음을 설명하라. 잡음 $\epsilon_{t-1},\allowbreak \epsilon^\prime_t$는 서로, 그리고 과거와 독립인 $\N(0,1)$이다.

2. (b)

   대입해 $X_t=\rho^2X_{t-1}+\eta_t$를 얻고 $\Var\eta_t=1-\rho^4$을 보여라.

3. (c)

   AR(1)의 자기상관이 $\rho^{2k}$임을 쓰고 $1+2\sum_{k\ge1}\rho^{2k}=\frac{1+\rho^2}{1-\rho^2}$을 기하급수로 계산하라.

풀이

(a) $(X,Y)$가 평균 0, 분산 1, 상관 $\rho$의 이변량 정규이면 $Y\mid X=x\sim\N(\rho x,\allowbreak 1-\rho^2)$이고 $X\mid Y=y\sim\N(\rho y,\allowbreak 1-\rho^2)$이다(5회 정리 3.13). 따라서 $Y_{t-1}$은 $\rho X_{t-1}$에 분산 $1-\rho^2$의 독립 잡음을 더한 것이고, $X_t$도 같다. 매 단계 새 표준정규를 뽑으므로 잡음들은 서로, 그리고 과거와 독립이다. (b) 대입하면 $X_t=\rho\bigl(\rho X_{t-1}+\sqrt{1-\rho^2}\epsilon_{t-1}\bigr)+\sqrt{1-\rho^2}\epsilon^\prime_t=\rho^2X_{t-1}+\eta_t$, $\eta_t=\rho\sqrt{1-\rho^2}\epsilon_{t-1}+\sqrt{1-\rho^2}\epsilon^\prime_t$이다. 두 잡음이 독립이므로 $\Var\eta_t=\rho^2(1-\rho^2)+(1-\rho^2)=(1-\rho^2)(1+\rho^2)=1-\rho^4$이다. (c) 정상 상태에서 $\Cov(X_t,\allowbreak X_{t-k})=\rho^{2k}\Var X$이므로 시차 $k$의 자기상관은 $\rho^{2k}$이다. $1+2\sum_{k\ge1}\rho^{2k}=1+\frac{2\rho^2}{1-\rho^2}=\frac{1+\rho^2}{1-\rho^2}$이고, $N_{\mathrm{eff}}/N$은 그 역수 $\frac{1-\rho^2}{1+\rho^2}$이다.

이 절을 PDF로 보기 (17쪽)노트북 4절

이 회차에서 읽은 절 0개 · 회차 차례