24. Markov chain Monte Carlo (MCMC) inference

24.1. Introduction

앞 장에서 배웠던 몬테 카를로 법은 고차원에 대해 잘 작동하지 않는다는 점이었다. 고차원 분포에서 샘플링하는 가장 널리 쓰이는 방법은 마르코프 연쇄 몬테 카를로(MCMC)이다. MCMC의 기본 아이디어는 목표하는 확률분포를 정지 분포로 갖는 마르코프 연쇄를 만든 뒤 상태 공간에서 무작위 걸음을 해서 각 상태 \mathbf{x}에서 소비하는 시간이 p^{\ast}(\mathbf{x})에 비례하도록 하는 것이다. 이 연쇄로부터 (상호 연관된) 표본들을 추출함으로써, p^{\ast}에 대한 몬테 카를로 적분을 수행할 수 있게 된다.

변분 추론과 비교하자면, 변분 추론의 이점들은 작은 문제에 대해 더 빠르다는 것, 결정론적이라는 것, 멈출 지점을 판단하기 쉽다는 것, 로그 가능도에 대한 하한을 제공한다는 것이 있다. 마르코프 연쇄 몬테 카를로 샘플링의 이점들은 구현하기 쉽다는 것, 더 넓은 범위의 모델에 적용가능하다는 것, 큰 모델이나 데이터셋에 대해 더 빠르다는 것 등이 있다.

24.2. Gibbs sampling

유명한 마르코프 연쇄 몬테 카를로법으로 깁스 샘플링(깁스 동역학, 열 목욕)이 있다. 이는 좌표 하강법의 마르코프 연쇄 몬테 카를로 버전이다.

24.2.1. Basic idea

기본 아이디어는 각 단계마다 하나의 변수를 다른 변수들은 분포에 존재한다는 조건하에 샘플링을 하는 것이다. 관측가능한 변수면 샘플링하지 않는다.

이 때 p(x_{i} | \mathbf{x}_{-i})은 변수 i에 대한 전체 조건부함수라 불린다. 일반적으로 x_{i}는 변수들 중 일부에만 의존할 수 있는데, p(\mathbf{x})를 그래프 모델로 나타낸다면 i의 마르코프 담요들을 살펴봄으로써 의존성을 추론할 수 있다. 즉, x_{i}를 알기 위해서는 그의 근방만 알면 충분하다. 이런 관점에서 깁스 샘플링은 분산 알고리즘이라 할 수 있다. 하지만 병렬 알고리즘은 아닌데, 샘플들이 순차적으로 추출되어야 하기 때문이다.

그러므로 마르코프 연쇄가 발화되기 전까지는, 즉 정지 분포에 도달하기 전까지는 초기 표본들을 버리는 것이 필요하다.

24.2.2. Example: Gibbs sampling for the Ising model

쌍별 마르코프 무작위 장/조건부 무작위 장은 다음의 형태를 갖는다.

p(x_{t} | \mathbf{x}_{-t}, \mathbf{\theta}) = \prod_{s \in \mathrm{nbr}(t)} \phi_{st}(x_{s}, x_{t})

간선 잠재 함수가 \phi(x_{s}, x_{t}) = e^{J x_{s} x_{t}}인 경우에는 전체 조건부함수를 구하면 p(x_{t} = +1 | \mathbf{x}_{-t}, \mathbf{\theta}) = \mathrm{sigm}(2 J \eta_{t})이 된다. (\eta_{t} = x_{t}(a_{t} - d_{t}), 즉 근방 노드 중 부호가 같은 노드의 수에서 부호가 다른 노드의 수를 뺀 것에 현재 부호를 곱한 것)

이는 국소적 증거 항 \phi_{t}와 결합 가능하다. 예를 들어, 가우시안 관측 모델로 \phi_{t} = \mathcal{N}(y_{t} | x_{t}, \sigma^{2})인 경우에는 전체 조건부함수는 p(x_{t} = +1 | \mathbf{x}_{-t}, \mathbf{y}_{t}, \mathbf{\theta}) = \mathrm{sigm}(2 J \eta_{t} - \log \frac{\phi_{t}(+1)}{\phi_{t}(-1)})이 된다. 즉, x_{t}가 각 상태에 진입할 확률은 그 근방과의 유사도와 데이터에 대한 유사도를 모두 고려하게 되는 것이다.

24.2.3. Example: Gibbs sampling for inferring the parameters of a GMM

가우시안 혼합 모델의 매개변수를 깁스 샘플링을 이용해 추론해 보자. 반켤레사전분포를 사용할 때 전체 결합분포는 다음과 같다.

p(\mathbf{x}, \mathbf{z}, \mathbf{\mu}, \mathbf{\Sigma}, \mathbf{\pi}) = p(\mathbf{x} | \mathbf{z}, \mathbf{\mu}, \mathbf{\Sigma}) p(\mathbf{z} | \mathbf{\pi}) p(\mathbf{\pi}) \prod_{k=1}^{K} p(\mathbf{\mu}_{k}) p(\mathbf{\Sigma}_{k}) = (\prod_{i=1}^{N} \prod_{k=1}^{K} (\pi_{k} \mathcal{N}(\mathbf{x}_{i} | \mathbf{\mu}_{k}, \mathbf{\Sigma}_{k}))^{\mathbf{1}_{z_{i} = k}}) \mathrm{Dir}(\mathbf{\pi} | \mathbf{\alpha}) \prod_{k=1}^{K} \mathcal{N}(\mathbf{\mu}_{k} | \mathbf{m}_{0}, \mathbf{V}_{0}) \mathrm{IW}(\mathbf{\Sigma}_{k} | \mathbf{S}_{0}, \nu_{0})

혼합 분포 각각에 대한 사전분포는 통일시킨다. 이 때 전체 조건부 분포는 다음과 같다.

p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\mu}, \mathbf{\Sigma}, \mathbf{\pi}) \propto \pi_{k} \mathcal{N}(\mathbf{x}_{i} | \mathbf{\mu}_{k}, \mathbf{\Sigma}_{k})

p(\mathbf{\pi} | \mathbf{z}) = \mathrm{Dir}(\{ \alpha_{k} + \sum_{i=1}^{N} \mathbf{1}_{z_{i} = k} \}_{k=1}^{K})

p(\mathbf{\mu}_{k} | \mathbf{\Sigma}_{k}, \mathbf{z}, \mathbf{x}) = \mathcal{N}(\mathbf{\mu}_{k} | \mathbf{m}_{k}, \mathbf{V}_{k})

\mathbf{V}_{k}^{-1} = \mathbf{V}_{0}^{-1} + N_{k} \mathbf{\Sigma}_{k}^{-1}

\mathbf{m}_{k} = \mathbf{V}_{k}(\mathbf{\Sigma}_{k} N_{k} \bar{\mathbf{m}}_{k} + \mathbf{V}_{0}^{-1} \mathbf{m}_{0})

N_{k} = \sum_{i=1}^{N} \mathbf{1}_{z_{i} = k}

\bar{\mathbf{x}}_{k} = \frac{\sum_{i=1}^{N} \mathbf{1}_{z_{i} = k} \mathbf{x}_{i}}{N_{k}}

p(\mathbf{\Sigma}_{k} | \mathbf{\mu}_{k}, \mathbf{z}, \mathbf{x}) = \mathrm{IW}(\mathbf{\Sigma}_{k} | \mathbf{S}_{k}, \nu_{k})

\mathbf{S}_{k} = \mathbf{S}_{0} + \sum_{i=1}^{N} \mathbf{1}_{z_{i} = k} (\mathbf{x} - \mathbf{\mu}_{k})(\mathbf{x} - \mathbf{\mu}_{k})^{T}

\nu_{k} = \nu_{0} + N_{k}

24.2.3.1. Label switching

혼합 모델에서의 깁스 샘플링은 근본적 약점이 있다. 모델의 매개변수 \mathbf{\theta}와 지시함수 \mathbf{z}를 특정할 수 없다는 것이다. 같은 가능도를 도출하는 서로 다른 여러 은닉 라벨의 순열이 존재하기 때문이다. 즉, 표본들의 몬테 카를로 평균을 취하는 것만으로는 사후평균을 계산할 수 없다. 그 대신 모든 최빈값에 대해 평균을 내서 \mathbb{E}[\mathbf{\mu}_{k} | \mathcal{D}]이 모든 k에 대해 같아지는 값을 찾아야 한다. 이를 라벨 변환 문제라 한다.

이 문제는 기대값 최대화나 변분 기대값 최대화 문제에서는 등장하지 않는다. 단일 최빈값에 분포를 가두기 때문이다. 하지만 최빈값이 여러 개인 분포에서는 항상 등장할 수 있는 문제이다. 1차원 문제에서는 매개변수에 대해 조건을 걸어서 특정을 할 수 있다. 하지만 가능도가 사전분포를 압도하게 될 수 있으므로 항상 활용가능한 방법은 아니다. 또한, 고차원으로 확장가능한 방법도 아니다. 다른 방법은 표본들을 후처리해 손실 함수를 최소화하는 순열을 찾는 것이다. 그러나 이 방법은 느리다.

아마도 최선의 방법은 유일하게 결정 가능하지 않은 해답을 갖는 질문을 아예 하지 않는 것이다. 예를 들면 데이터 i가 클러스터 k에 속하는지가 아니라 데이터 i와 데이터 j가 같은 클러스터에 속하는지를 알아보는 방법이 있다. 후자의 문제는 라벨링에 불변이고, 관측가능한 통계량에만 의존한다. 이는 혼합 모델에 대해 확장가능하다는 이점도 있다. 이 경우 은닉 클러스터는 잘 정의되지 않지만 분할은 잘 정의되기 때문이다.

24.2.4. Collapsed Gibbs sampling

일부 경우에 대해서는 미지 변수 중 일부를 적분해 빼내고 나머지만 샘플링하는 방법도 있다. 이를 붕괴 깁스 샘플링이라 한다. 이는 더 저차원에서 샘플링을 하므로 훨씬 더 효과적이다.

정확하게는, \mathbf{z}를 샘플링하고 \mathbf{\theta}를 적분해 빼낸다고 하자. \mathbf{\theta}는 마르코프 연쇄에 참가하지 않는다. 그러므로 조건부 독립인 표본 \mathbf{\theta}_{s} \sim p(\mathbf{\theta} | \mathbf{z}_{s}, \mathcal{D})들을 추출할 수 있다. 이는 단순 결합분포에서 추출한 샘플들에 비해 훨씬 더 낮은 분산을 갖는다. 이를 라오-블랙웰화라 하며, 다음 정리와 연관이 있다.

Theorem. (Rao-Blackwell) \mathbf{z}\mathbf{\theta}가 연관 확률변수고 f(\mathbf{z}, \mathbf{\theta})가 스칼라 함수라 하자. 그러면 \mathrm{var}_{\mathbf{z}, \mathbf{\theta}} [f(\mathbf{z}, \mathbf{\theta}] \geq \mathrm{var}_{\mathbf{z}} [\mathbb{E}_{\mathbf{\theta}} [f(\mathbf{z}, \mathbf{\theta}] | \mathbf{z}]]이 성립한다.

이 정리는 \mathbf{\theta}를 적분해 빼낸 추정으로부터의 분산이 항상 직접 몬테 카를로 추정을 한 분산 이하임을 보증한다.

24.2.4.1. Example: collapsed Gibbs for fitting a GMM

완전 켤레사전분포를 사용하는 가우시안 혼합 모델을 고려해 보자. 이 경우 \mathbf{\mu}_{k}, \mathbf{\Sigma}_{k}, \mathbf{\pi}를 적분해 빼낼 수 있고, \mathbf{z}만 샘플링하면 충분하다. 이 때 전체 조건부함수는 다음과 같아진다.

p(z_{i} = k | \mathbf{z}_{-i}, \mathbf{x}, \mathbf{\alpha}, \mathbf{\beta}) \propto p(z_{i} = k | \mathbf{z}_{-i}, \mathbf{\alpha}) p(\mathbf{x}_{i} | \mathbf{x}_{-i}, z_{i} = k, \mathbf{z}_{-i}, \mathbf{\beta}) (\mathbf{\beta} = (\mathbf{m}_{0}, \mathbf{V}_{0}, \mathbf{S}_{0}, \nu_{0})은 클래스 조건분포의 초매개변수)

첫 번째 항은 다음과 같다.

p(z_{i} = k | \mathbf{z}_{-i}, \mathbf{\alpha}) = \frac{N_{k, -i} + \frac{\alpha}{K}}{N + \alpha - 1} (N_{k, -i} = \sum_{n \neq i} \mathbb{1}_{z_{n} = k} = N_{k} - 1)

두 번째 항은 다음과 같다.

p(\mathbf{x}_{i} | \mathbf{x}_{-i}, \mathbf{z}_{-i}, z_{i} = k, \mathbf{\beta}) = p(\mathbf{x}_{i} | \mathcal{D}_{-i, k}) (\mathcal{D}_{-i, k} = \{\mathbf{x}_{j} : z_{j} = k, j \neq i \})

이 때 \mathbf{\theta}_{k}에 대한 켤레사전분포를 사용한다면 위의 항을 닫힌 형태로 계산할 수가 있다. 또한, 각 클러스터에 대한 충족 통계량을 구함으로써 예측가능도도 업데이트할 수 있다. 현재 클러스터 z_{i}에서 \mathbf{x}_{i}의 통계량을 제거한 뒤 각 클러스터의 사후예측분포에서 \mathbf{x}_{i}을 계산하면 된다. 그 방법으로 새 클러스터를 계산한 뒤에, 그 새 클러스터에 \mathbf{x}_{i}의 통계량을 추가하면 된다. 순차적으로 p(z_{i} | \mathbf{z}_{1 : i - 1}, \mathbf{x}_{1 : i})을 샘플링함으로써 샘플을 초기화할 수 있다. 가우시안 혼합 모델의 경우에는 단계마다 O(NKD) 시간이 걸린다.

붕괴 깁스 샘플링과 일반 깁스 샘플링을 비교하면 붕괴 샘플링이 일반적으로 더 잘 동작함을 알 수 있다. 하지만, 둘 다 국소적 해에 갇힐 수 있다.

24.2.5. Gibbs sampling for hierarchical GLMs

데이터가 상호 연관된 복수의 소스를 가지고 어떤 소스의 신뢰도가 다른 것들보다 높다면 모든 데이터를 동시에 모델링해서 통계적 세기를 빌려오는 것이 좋을 것이다. 이러한 문제를 푸는 자연스러운 방법은 계층적 베이지안 모델링, 다단계 모델링을 사용하는 것이다. 이에 깁스 샘플링을 사용해 보자.

y_{ij} = x_{ij}^{T} \mathbf{w}_{j} + \epsilon_{ij}

\mathbf{w}_{j} \sim \mathcal{N}(\mathbf{\mu}_{w}, \mathbf{\Sigma}_{w})

\mathbf{\mu}_{w} \sim \mathcal{N}(\mathbf{\mu}_{0}, \mathbf{V}_{0})

\mathbf{\Sigma}_{w} \sim \mathrm{IW}(\eta_{0}, \mathbf{S}_{0}^{-1})

\sigma^{2} \sim \mathrm{IG}(\frac{\nu_{0}}{2}, \frac{\nu_{0}\sigma_{0}^{2}}{2})

의 형태를 가질 때, 깁스 샘플링의 전체 조건부 함수는 다음과 같다.

p(\mathbf{w}_{j} | \mathcal{D}_{j}, \mathbf{\theta}) = \mathcal{N}(\mathbf{w}_{j} | \mathbf{\mu}_{j}, \mathbf{\Sigma}_{j})

\mathbf{\Sigma}_{j}^{-1} = \mathbf{\Sigma}^{-1} + \frac{1}{\sigma^{2}} \mathbf{X}_{j}^{T} \mathbf{X}_{j}

\mathbf{\mu}_{j} = \mathbf{\Sigma}_{j}(\mathbf{\Sigma}^{-1} \mathbf{\mu} + \frac{1}{\sigma^{2}} \mathbf{X}_{j}^{T} \mathbf{y}_{j})

p(\mathbf{\mu}_{w} | \mathbf{W}_{1 : J}, \mathbf{\Sigma}_{w}) = \mathcal{N}(\mathbf{\mu} | \mathbf{\mu}_{N}, \mathbf{\Sigma}_{N})

\mathbf{\Sigma}_{N}^{-1} = \mathbf{V}_{0}^{-1} + J \mathbf{\Sigma}^{-1}

\mathbf{\mu}_{N} = \mathbf{\Sigma}_{N}(\mathbf{V}_{0}^{-1} \mathbf{\mu}_{0} + J \mathbf{\Sigma}^{-1} \bar{\mathbf{w}}) (\bar{\mathbf{w}} = \frac{1}{J} \sum_{j} \mathbf{w}_{j}

p(\mathbf{\Sigma}_{w} | \mathbf{\mu}_{w}, \mathbf{w}_{1 : J}) = \mathrm{IW}((\mathbf{S}_{0} + \mathbf{S}_{\mu})^{-1}, \nu_{0} + J)

\mathbf{S}_{\mu} = \sum_{j} (\mathbf{w}_{j} - \mathbf{\mu}_{w})(\mathbf{w}_{j} - \mathbf{\mu}_{w})^{T}

p(\sigma^{2} | \mathcal{D}, \mathbf{w}_{1 : J}) = \mathrm{IG}(\frac{\nu_{0} + N}{2}, \frac{\nu_{0} \sigma_{0}^{2} + \mathrm{SSR}(\mathbf{w}_{1:J})}{2})

\mathrm{SSR}(\mathbf{w}_{1 : J}) = \sum_{j=1}^{J} \sum_{i=1}^{N_{j}} (y_{ij} - \mathbf{w}_{j}^{T} \mathbf{x}_{ij})^{2}

\mathbb{E}[y_{j} | \mathbf{x}_{ij}] = \mathbf{x}_{ij}^{T} \hat{\mathbf{w}}_{j}

\hat{\mathbf{w}}_{j} = \mathbb{E}[\mathbf{w}_{j} | \mathcal{D}] \simeq \frac{1}{S} \sum_{s=1}^{S} \mathbf{w}_{j, s}

이 방법은 정규화를 잘 수행한다.

다층 선형 회귀 데모.

24.2.6. BUGS and JAGS

깁스 샘플링의 장점 중 하나는 거의 모든 모델에 적용할 수 있는 소프트웨어를 설계할 수 있다는 것이다. 모델 특정 (그래프 형태)와 전체 조건부함수로부터 샘플링하는 라이브러리만 있으면 된다. 이러한 패키지의 예제는 깁스 샘플링을 사용한 베이지안 업데이트 (BUGS), 또 다른 깁스 샘플러(JAGS) 등이 있다. 그러나 수작업으로 최적화한 코드보다 느리다는 단점이 있다.

24.2.7. The Imputation Posterior (IP) algorithm

전가 사후분포(IP) 알고리즘은 깁스 샘플링의 특수한 경우로, 변수를 두 클래스로 나눈다. 은닉 변수 \mathbf{z}와 매개변수 \mathbf{\theta}. 이는 기대값 최대화의 마르코프 연쇄 몬테 카를로 버전으로, E 단계가 I 단계로, M 단계가 P 단계로 대체된 것이다. 이를 일반화한 것이 데이터 증강이다.

24.2.8. Blocking Gibbs sampling

깁스 샘플링은 한 번에 하나의 변수만 업데이트하므로 (단일 사이트 업데이트) 느릴 수 있다. 변수들의 상호 연관성이 크다면, 현 상태로부터 이동하는 데 오랜 시간이 걸린다. 이동의 크기는 조건부 분포의 분산에 의존한다. x_{1} 방향의 분산이 l이고 지지집합의 크기가 L이면 독립적인 표본을 얻는 데 O(\frac{L}{l}^{2}) 단계가 걸린다.

일부 경우에 대해서는 변수의 그룹을 한 번에 업데이트할 수 있다. 이를 블록화 깁스 샘플링 또는 블록 깁스 샘플링이라 불린다.

쏠린 2D 가우시안에 대해 깁스 샘플링의 속도 예제.

24.3. Metropolis Hastings algorithm

깁스 샘플링은 간단하지만 적용가능한 모델이 제한되어 있으며 느리다는 단점이 있다. 그래서 이를 일반화한 메트로폴리스 헤이스팅스(MH) 알고리즘이 있다.

24.3.1. Basic idea

메트로폴리스 헤이스팅스의 기본 발상은 각 단계에서 제안 분포 q(\mathbf{x}^{\prime} | \mathbf{x})의 확률로 현 상태 \mathbf{x}에서 다음 상태 \mathbf{x}^{\prime}으로 이동하는 것이다. 제안 분포는 특정 조건만 만족한다면 어떻게 골라도 상관없다. 널리 사용되는 분포는 q(\mathbf{x}^{\prime} | \mathbf{x}) = \mathcal{N}(\mathbf{x}^{\prime} | \mathbf{x}, \mathbf{\Sigma})가 있다. 이를 무작위 걸음 메트로폴리스 알고리즘이라 한다. q(\mathbf{x}^{\prime} | \mathbf{x}) = q(\mathbf{x}^{\prime})으로 이전 상태에 의존하지 않는 제안 분포를 만든다면 이를 독립 추출기라 한다. 이는 중요성 추출과 비슷하다.

\mathbf{x}^{\prime}으로의 이동이 제안된 상태에서, 각 상태에서 보내는 시간의 비율이 p^{\ast}(\mathbf{x})에 비례하는지를 체크해서 이 제안을 채택할지 결정할 수 있다. 제안 분포가 대칭적(q(\mathbf{x}^{\prime} | \mathbf{x}) = q(\mathbf{x} | \mathbf{x}^{\prime}))이면, 채택 확률은 다음과 같다.

r = \min(1, \frac{p^{\ast}(\mathbf{x}^{\prime})}{p^{\ast}(\mathbf{x})})

이는 탐욕적 이동뿐만 아니라 하향 이동도 가능하게 한다. 제안 분포가 대칭적이지 않으면 채택 확률은 다음의 헤이스팅스 보정으로 주어진다.

r = \min(1, \alpha)

\alpha = \frac{p^{\ast}(\mathbf{x}^{\prime})q(\mathbf{x} | \mathbf{x}^{\prime})}{p^{\ast}(\mathbf{x}) q(\mathbf{x}^{\prime} | \mathbf{x})}

메트로폴리스 헤이스팅스 알고리즘을 유용하게 만드는 중요한 이유 중 하나는 \alpha를 계산할 때 그저 대상 분포만 알면 된다는 것이다. 즉 표준화 상수를 알지 못해도 p^{\ast}로부터 표본을 추출할 수 있다. \tilde{p}(\mathbf{x}) = Z p^{\ast}(\mathbf{x})를 점별로 계산하는 것으로 충분하다.

24.3.2. Gibbs sampling is a special case of MH

깁스 샘플링은 메트로폴리스 헤이스팅스의 특수 케이스이다. 정확히는, q(\mathbf{x}^{\prime} | \mathbf{x}) = p(x_{i}^{\prime} | \mathbf{x}_{-i}) \mathbf{1}_{\mathbf{x}_{-i}^{\prime} = \mathbf{x}_{-i}})의 제안 분포를 사용할 때의 메트로폴리스 헤이스팅스가 바로 깁스 샘플링이다. 이 때 채택 확률은 100%가 되지만 그것이 깁스 샘플링의 수렴성을 보장하는 것은 아니다. 한 번에 하나의 좌표만 업데이트하기 때문이다. 그래서 다른 제안 분포들이 존재한다.

24.3.3. Proposal distribution

대상 분포 p^{\ast}에 대해, 제안 분포 q가 허용 가능하거나 적절하려면 지지집합이 대상 분포보다 커야 한다. 물론, 제안 분포의 모양도 매우 중요하다. 제안 분포로 가우시안 혼합 모델을 사용하는 경우 분산이 너무 작다면 연쇄는 하나의 최빈값에 갇힌다. 분산이 너무 크다면 대부분의 이동이 기각되고 연쇄는 초기값에 매우 달라붙을 것이다.

깁스 샘플링의 큰 이점은 제안 분포를 따로 고를 필요가 없다는 것이다. 채택률 자체는 크게 중요하지는 않다. 최빈값을 평균으로 갖고 분산이 0인 분포도 채택률은 100%지만 쓸모가 없기 때문이다. 가우시안 커널의 분산을 증가시킴으로써 적절한 채택률을 고를 수가 있는데, 가우시안 대상 분포에 대해서는 25%-40%가 최적인 것으로 알려져 있다. 이렇게 제안 분포를 튜닝하기 위한 짧은 최초의 샘플링들을 파일럿 수행이라 한다.

가우시안 혼합 모델에 대한 마르코프 연쇄 몬테 카를로 vs 메트로폴리스 헤이스팅스.

24.3.3.1. Gaussian proposals

상태 공간이 연속이라면 가우시안 제안 분포의 공분산에는 헤시안 \mathbf{H}를 국소적 최빈값 \hat{\mathbf{w}}에서 구한 값을 쓴다. 이는 헤시안이 국소적 곡률과 각 차원의 길이 스케일을 모델링한다는 점을 이용해서 깁스 샘플링의 느리다는 단점을 보완한다.

두 가지 자연스러운 접근법이 있다. 독립적 제안 q(\mathbf{w}^{\prime} | \mathbf{w}) = \mathcal{N}(\mathbf{w}^{\prime} | \hat{\mathbf{w}}, \mathbf{H}^{-1})을 사용하는 것과 무작위 걸음 제안 q(\mathbf{w}^{\prime} | \mathbf{w}) = \mathcal{N}(\mathbf{w}^{\prime} | \mathbf{w}, s^{2} \mathbf{H}^{-1})을 사용하는 것. 이 때 s^{2} = \frac{2.38^{2}}{D}가 점근적으로 최적임이 알려져 있다. 이 때 채택률은 0.234가 된다.

1D 로지스틱 회귀에 대한 메트로폴리스 헤이스팅스.

24.3.3.2. Mixture proposals

여러 제안 분포를 혼합한 혼합 제안도 있다. q(\mathbf{x}^{\prime} | \mathbf{x}) = \sum_{k=1}^{K} w_{k} q_{k}(\mathbf{x}^{\prime} | \mathbf{x})의 형태를 갖는다. (w_{k} > 0, q_{k}는 올바른 제안 분포)

24.3.3.3. Data-driven MCMC

최적의 제안 분포는 단지 이전 은닉 상태뿐만 아니라 관측 가능한 데이터에도 의존하는 q(\mathbf{x}^{\prime} | \mathbf{x}, \mathcal{D})의 형태를 가진다. 이를 데이터 기반 마르코프 연쇄 몬테 카를로라 한다. 이러한 제안 분포를 만들기 위해서는 전방 모델로부터 (\mathbf{x}, \mathcal{D}) 쌍을 샘플링하고 구별적 분류기를 학습해 p(\mathbf{x} | f(\mathcal{D}))을 예측하는 것이다. (f(\mathcal{D})는 관측 가능한 데이터에서 추출한 특성)

보통 \mathbf{x}는 고차원이기 때문에 전체 상태 벡터 p(\mathbf{x} | f(\mathcal{D})) 를 예측하기는 어렵다. 그 대신에 각 상태에 대한 벡터 p(x_{k} | f_{k}(\mathcal{D}))을 학습한다. 이후 q(\mathbf{x}^{\prime} | \mathbf{x}, \mathcal{D}) = \pi_{0} q_{0}(\mathbf{x}^{\prime} | \mathbf{x}) + \sum_{k} \pi_{k} q_{k} (x_{k}^{\prime} | f_{k}(\mathcal{D}))의 제안을 이용한다. 여기서 q_{0}은 데이터에 독립적인 제안 분포이고 (ex. 무작위 걸음) q_{k}는 상태 공간의 k번째 부분을 업데이트한다. 복수 변수에 대한 결합 변화를 제안할 수도 있지만 이는 쉽지 않다.

이 과정은 생성 후 테스트의 형태를 가진다. 구별적 제안 분포 q(\mathbf{x}^{\prime} | \mathbf{x})는 새 가설들을 생성하고, 이는 사후분포 비율 \frac{p(\mathbf{x}^{\prime} | \mathcal{D})}{p(\mathbf{x} | \mathcal{D})}을 계산함으로써 테스트된다. 담금질 단계를 추가함으로써, 이 알고리즘은 사후 최빈값을 찾는 데 쓰일 수 있다. 이를 시뮬레이션화된 담금질이라 한다. 이는 제안 분포가 가역적일 필요가 없다는 장점이 있다.

24.3.4. Adaptive MCMC

알고리즘이 진행되는 도중 제안 분포의 매개변수를 변화시킬 수 있다. 이를 적응적 마르코프 연쇄 몬테 카를로라 한다. 이는 처음에는 넓은 공분산으로 시작한 뒤 최빈값을 찾고 나면 공분산을 줄이는 식의 알고리즘을 가능케 한다. 하지만, 마르코프 특성을 깨면 안 된다는 것에 주의할 필요가 있다. 즉, 제안 분포의 매개변수가 연쇄의 전체 기록에 의존하면 안 된다. 이를 보장하는 방법은 매개변수의 적응을 시간에 따라 흐릿하게 하는 것이다.

24.3.5. Initialization and mode hopping

마르코프 연쇄 몬테 카를로를 초기화할 때에는 0이 아닌 확률에서 초기화할 필요가 있다. 모델이 결정론적 제한 조건을 갖고 있다면 이는 어렵기 때문에, 초기화할 때는 대개 국소적 최빈값에서 초기값을 갖도록 하고는 한다. 상태 공간이 이산적이라면 최적화를 여러 번 돌려서 여러 최빈값에서 평균을 내는 것이 효율적이고, 상태 공간이 연속적이라면 최빈값의 확률밀도가 낮으므로 각 최빈값에서 국소적으로 탐색을 하는 것이 필요하다.

24.3.6. Why MH works

메트로폴리스 헤이스팅스 알고리즘은 왜 p^{\ast}로부터 표본을 추출하는가? 이 알고리즘은 다음의 전이 행렬에 대한 마르코프 연쇄를 정의한다.

p(\mathbf{x}^{\prime} | \mathbf{x}) = q(\mathbf{x}^{\prime} | \mathbf{x}) r(\mathbf{x}^{\prime} | \mathbf{x}) if \mathbf{x}^{\prime} \neq \mathbf{x}, 이외엔 q(\mathbf{x} | \mathbf{x}) + \sum_{\mathbf{x}^{\prime} \neq \mathbf{x}} q(\mathbf{x}^{\prime} | \mathbf{x}) (1 - r(\mathbf{x}^{\prime} | \mathbf{x}))

이를 분석해 보자. 먼저 연쇄가 p(\mathbf{x}^{\prime} | \mathbf{x}) p^{\ast}(\mathbf{x}) = p(\mathbf{x} | \mathbf{x}^{\prime}) p^{\ast} (\mathbf{x}^{\prime})을 만족할 경우 세부 균형을 만족한다고 한다. 이 경우 p^{\ast}는 그 자체의 저이 분포가 된다. 그리고 다음의 정리에 의해 메트로폴리스 헤이스팅스 알고리즘의 정당성이 보장된다.

Theorem. 메트로폴리스 헤이스팅스 알고리즘에 의해 정의된 전이 행렬이 에르고딕이고 기약이면 p^{\ast}는 그 자신의 유일한 제한 분포이다.

24.3.7. Reversible jump (trans-dimensional) MCMC

매개변수의 수가 서로 다른 모델들이 있다고 하다. 모델을 m, 미지수를 \mathbf{x}_{m} \in \mathcal{X}_{m}이라 하자. 서로 다른 차원을 갖는 공간에서의 샘플링을 차원 변화 마르코프 연쇄 몬테 카를로라 한다. 먼저 모델 지시자 m \{ 1, \cdots, M \}을 샘플링한다. 그 뒤 모든 매개변수를 \prod_{m=1}^{M} \mathcal{X}_{m}로부터 추출할 수도 있지만 이는 매우 비효율적이고, \mathcal{X} = \cup_{m=1}^{M} \{m\} \times \mathcal{X}_{m}로부터 샘플링하는 것이 효율적이다.

이 접근법의 난점은 서로 다른 차원에서 모델을 이동할 때 발생한다. 메트로폴리스 헤이스팅스 채택율을 계산할 때 서로 다른 공간에서의 확률분포를 비교하게 되므로 의미가 없어지기 때문이다. 그래서 저차원 공간을 추가 확률변수로 증강해 차원을 맞추는 가역 점프 마르코프 연쇄 몬테 카를로(RJMCMC) 등의 방법이 있다. 이는 실제로는 상당히 어렵지만 연속적 매개변수를 적분해 빼낼 수 있다면 이산적 상태 공간만 남기 때문에 훨씬 더 쉬워진다. 이를 붕괴 가역 점프 마르코프 연쇄 몬테 카를로라 한다. 이를 이용해 연관 벡터 기계나 보조 벡터 기계에 대한 베이지안적 대안을 제공하는 것도 가능하다.

24.4. Speed and accuracy of MCMC

마르코프 연쇄 몬테 카를로와 관련된 이론적, 실제적 이슈를 알아보자.

24.4.1. The burn-in phase

정지 분포에 도달하기 이전 추출된 표본들은 버려진다. 이를 발화 단계라 한다. 이 단계를 결정하는 것은 어려운데, 이것은 마르코프 연쇄 몬테 카를로법의 약점이기도 하다.

무작위 걸음의 예제.

24.4.2. Mixing rates of Markov chains

마르코프 연쇄가 수렴하는 데 걸리는 시간을 혼합 시간이라 부른다. 즉, 상태 x_{0}으로부터 혼합 시간은 모든 \epsilon > 0에 대해 \tau_{\epsilon}(x_{0}) = \min \{ t : \lVert \delta_{x_{0}}(x)T^{t} - p^{\ast} \rVert_{1} \leq \epsilon \}이 되는 최소 시간이다. (T는 전이 행렬) 이 때 혼합 시간은 \tau_{\epsilon} = \max_{x_{0}} \tau_{\epsilon}(x_{0})으로 정의된다. \gamma를 전이 행렬의 첫 번째 고유값과 두 번째 고유값의 차라고 할 때, \tau_{\epsilon} \leq O(\frac{1}{\gamma} \log \frac{n}{\epsilon})임이 알려져 있다.

다른 접근법은 상태 공간의 모양을 조사해 보는 것이다. 연쇄의 전도도\phi = \min_{S : 0 \leq p^{\ast}(S) \leq 0.5} \frac{\sum_{x \in S, x^{\prime} \in S^{c}} T(x \to x^{\prime})}{p^{\ast}(S)}로 정의할 때, 즉 해당 집합에서 다른 집합으로 전이할 확률의 최소치로 정의할 때 \tau_{\epsilon} \leq O(\frac{1}{\phi^{2}} \log \frac{n}{\epsilon})임이 알려져 있다. 즉, 전도도가 낮을 수록 혼합 시간은 길어진다. 최빈값이 잘 분리된 분포라면 전도도가 낮고 혼합 시간은 길어질 것이다.

24.4.3. Practical converge diagnostics

연쇄의 혼합 시간을 계산하는 것은 일반적으로 매우 어렵다. 전이 행렬의 계산이 어렵기 때문이다. 실제로는 수렴을 진단하는 휴리스틱이 많이 쓰이는데, 실제로 수렴성을 정확히 판단해 주지는 못한다. 간단한 방법 중 하나는 과분산된 시작점으로부터 여러 연쇄를 시작시킨 뒤 추적 플롯을 찍어봐서 같은 지점으로 수렴했는지 알아보는 것이다.

24.4.3.1. Estimated potential scale reduction (EPSR)

수렴의 정도를 정량적으로 알아볼 수 있는 방법은 다음과 같다. 기본 발상은 각각 연쇄에서의 분산과 연쇄들간의 분산을 비교하는 것이다. 발화 이후 D개의 변수를 가진 C개의 연쇄로부터 각각 S개의 샘플 x_{isc}를 추출하고 y_{sc}\mathbf{x}_{1:D, s, c}로부터 유도되는 관심 있는 스칼라값이라 하면 수열 평균과 전체 평균을 \bar{y}_{.c} = \frac{1}{S} \sum_{s=1}^{S} y_{sc}, \bar{y}_{..} = \frac{1}{C} \sum_{c=1}^{C} \bar{y}_{.c}로 정의하자. 수열간 분산과 수열 내 분산을 B = \frac{S}{C-1}\sum_{c=1}^{C} (\bar{y}_{.c} - \bar{y}_{..})^{2}, W= \frac{1}{C} \sum_{c=1}^{C} [\frac{1}{S-1} \sum_{s=1}^{S} (y_{sc} - \bar{y}_{.c})^{2}]로 정의하면, \hat{V} = \frac{S-1}{S}W + \frac{1}{S} B이라 할 때 W와 \hat{V}\mathrm{var}[y]를 각각 과소추정, 과대추정한다. 이로부터, 추정 잠재 크기 감소량(EPSR) \hat{R} = \sqrt{\frac{\hat{V}}{W}}로 정의하면 계속 샘플링을 했을 때 사후분산이 감소할지 아닐지를 판단할 수 있다. \hat{R} \simeq 1이라면 이 추정은 신뢰할만 하다고 할 수 있다.

24.4.4. Accuracy of MCMC

마르코프 연쇄 몬테 카를로로부터 추출되는 표본들은 자동으로 상호 연관되므로 독립적인 표본들이 가지는 정보량을 감소시킨다. 이를 정량화하면 다음과 같다. 어떤 함수 f에 대해 X \sim p()일 때 f(X)의 평균을 구하고 싶다고 하자. 평균의 참값을 f^{\ast} = \mathbb{E}[f(X)]이라고 하고 이에 대한 몬테 카를로 추정을 \bar{f} = \frac{1}{S} \sum_{s=1}^{S} f_{s}이라 하자. 이 추정의 분산에 대한 마르코프 연쇄 몬테 카를로 추정은 다음과 같다.

\mathrm{var}_{\mathrm{MCMC}}[\bar{f}] = \mathbb{E}[(\bar{f} - f^{\ast})^{2}] = \mathrm{var}_{\mathrm{MC}}(\bar{f}) + \frac{1}{S^{2}} \sum_{s \neq t} \mathbb{E}[(f_{s} - f^{\ast})(f_{t} - f^{\ast})]

이는 자동 상호 연관 함수(ACF) \rho_{t} = \frac{\frac{1}{S-t} \sum_{s=1}^{S-t} (f_{s} - \bar{f})(f_{s+t} - \bar{f})}{\frac{1}{S-1} \sum_{s=1}^{S} (f_{s} - \bar{f})^{2}}을 정의해 측정할 수 있다. 이는 n번째마다 샘플을 보존하는 솎아내기로 계산을 쉽게 할 수 있다. 유효 표본수(ESS) S_{\mathrm{eff}} = \frac{\mathrm{var}_{\mathrm{MC}}(f)}{\mathrm{var}_{\mathrm{MCMC}}(\bar{f})   }을 정의하면 표본의 집합의 정보량을 이로부터 추정할 수 있다.

24.4.5. How many chains?

그래서 몇 개의 연쇄를 수행해야 할까? 보통은 3개 정도의 연쇄를 각각 길이 100000 정도로 수행하고 표본의 첫 절반을 버린 뒤 표본을 취한다. 국소적 최빈값으로 초기값을 잡으면 발화를 기다리느라 표본의 첫 절반을 버리지 않아도 될 것이다.

24.5. Auxiliary variable MCMC

가끔 더미 보조 변수를 써서 원본 변수들 사이의 상호 연관을 줄여 샘플링의 효율성을 비약적으로 높일 수 있다. 원본 변수를 \mathbf{x}, 보조 변수를 \mathbf{z}라고 할 때 \sum_{\mathbf{z}} p(\mathbf{x}, \mathbf{z}) = p(\mathbf{x})여야 하고 p(\mathbf{x}, \mathbf{z})p(\mathbf{x})보다 더 샘플링하기 쉬울 수 있기 때문이다. 이후에는 그냥 \mathbf{z}를 버리면 된다.

24.5.1. Auxillary variable sampling for logistic regression

프로빗 회귀의 잠재 변수 표현은 다음과 같다.

z_{i} = \mathbf{w}^{T} \mathbf{x}_{i} + \epsilon_{i}

\epsilon_{i} \sim \mathcal{N}(0, 1)

p(y_{i} = 1) = \mathbf{1}_{z_{i} \geq 0}

이는 보조 변수 깁스 샘플링으로 치환하기 쉽다.

로지스틱 회귀에 대해서는 어떻게 할까? \epsilon_{i}로지스틱 분포 p_{\mathrm{logistic}}(\epsilon) = \frac{e^{-\epsilon}}{(1 + e^{-\epsilon})^{2}}를 따른다고 하면 p(y_{i} = 1 | \mathbf{x}_{i}, \mathbf{w}) = \int_{-\mathbf{w}^{T} \mathbf{x}_{i}}^{\infty} f(\epsilon) d \epsilon = \mathrm{sigm}(\mathbf{w}^{T} \mathbf{x}_{i})이 된다. 이 때 p(\mathbf{z} | \mathbf{w}, \mathcal{D}), p(\mathbf{w} | \mathbf{z}, \mathcal{D})로부터 샘플링을 함으로써 보조 변수 깁스 샘플링을 할 수 있다. p(\mathbf{w} | \mathbf{z}, \mathcal{D})로부터 직접 샘플링하는 것은 불가능하기 때문에, 한 가지 방법은 \lambda_{i} = (2 \phi_{i})^{2}, \phi_{i} \sim KS에 대해 (KS는 콜모고로프 스미르노프 분포), \epsilon_{i} \sim \mathcal{N}(0, \lambda_{i})를 샘플링하는 것이다. 더 간단한 방법은 로지스틱 분포를 스튜던트 분포로 근사하는 것이다. \epsilon_{i} \sim \mathcal{T}(0, 1, \nu), \nu \simeq 8로 근사를 한 뒤 스튜던트 분포의 가우시안 혼합 표현으로 추론을 간략화하는 것이다. 즉, 다음이 성립한다.

\lambda_{i} \sim \mathrm{Ga}(\frac{\nu}{2}, \frac{\nu}{2})

\epsilon_{i} \sim \mathcal{N}(0, \frac{1}{\lambda_{i}})

z_{i} = \mathbf{w}^{T} \mathbf{x}_{i} + \epsilon_{i}

y_{i} = 1 | z_{i} = \mathbf{1}_{z_{i} \geq 0}

이 때 \nu = 1로 두면 z_{i} \sim \mathcal{N}(\mathbf{w}^{T}\mathbf{x}_{i}, 1)이 되어 프로빗 회귀와 동치가 된다. 프로빗 회귀냐 로짓 회귀냐를 선택하는 대신에 \nu를 추정할 수 있다. 적합한 켤레사전분포는 없지만 유한한 범위의 가능한 값들을 고려한 뒤 사후분포를 p(\nu | \mathbf{\lambda}) \propto p(\nu) \prod_{i=1}^{N} \frac{1}{\Gamma(\frac{\nu}{2}) (\frac{\nu}{2})^{\frac{\nu}{2}}} \lambda_{i}^{\frac{\nu}{2} - 1} e^{-\frac{\nu \lambda_{i}}{2}}로 구하면 된다. \mathbf{V}_{0} = v_{0} \mathbf{I}로 놓으면 v_{0}도 샘플링할 수 있다.

24.5.2. Slice sampling

일변수 다최빈값 분포 \tilde{p}(x)를 생각해 보자. 보조 변수 u를 더해 이동의 크기를 크게 할 수 있다. 결합분포를 \hat{p}(x, u) = \frac{1}{Z_{p}} if 0 \leq u \leq \tilde{p}(x), (Z_{p} = \int \tilde{p}(x) dx) 이외 0으로 정의하면 x에서의 주변분포는 p(x) = \frac{\tilde{p}(x)}{Z_{p}}이 된다. 그러므로 \hat{p}(x, u)로부터 샘플링한 뒤 u를 무시하면 p(x)로부터 샘플링할 수 있다. 전체 조건부분포는 다음과 같다.

p(u | x) = U_{[0, \tilde{p}(x)]}(u)

p(x | u) = U_{A}(x) (A = \{x : \tilde{p}(x) \geq u \})

이를 절단 샘플링이라 한다. 실제로는 집합 A를 구하기 어렵다. 그래서 적당한 점 x_{s} 근처의 구간 x_{\min} \leq x \leq x_{\max}를 잡은 뒤 각 끝점이 절단부에 위치하는지 본다. 위치할 경우에는 절단면 밖으로 빠져나올 때까지 그 끝점을 확장시킨다 (스텝 아웃). 이후 후보 값 x^{\prime}이 그 지역으로부터 균일하게 선택된다. 그 변수가 절단면 내에 있으면 x_{s + 1} = x^{\prime}으로 유지시킨다. 그렇지 않으면 영역을 축소해 x^{\prime}이 끝점 중 하나가 되면서 x_{s}을 포함하도록 한다. 이후 표본을 하나 더 추출한다. 이를 표본이 채택될 때까지 반복한다.

이를 다변수 분포로 확장하려면 각 변수마다 보조 변수를 써야 한다. 이 접근법이 깁스 샘플링에 비해 갖는 이점은 전체 조건부를 알아야 할 필요가 없으며 비표준화 결합분포만 알면 충분하다는 점이다. 메트로폴리스 헤이스팅스 알고리즘에 비해 갖는 이점은 제안 분포를 특정해 줘야 할 필요가 없다는 점이다.

1D 절단 샘플링.
2D 절단 샘플링.

24.5.3. Swendsen Wang

p(\mathbf{x}) = \frac{1}{Z} \prod_{e} f_{e}(\mathbf{x}_{e}), \mathbf{x}_{e} = (x_{i}, x_{j}) (e는 간선) 꼴의 아이싱 모델을 생각해 보자. 간선 세기 J가 클 때 인접한 노드끼리의 상호 연관성이 클 수 있기 때문에 깁스 샘플링은 느릴 수 있다. 스웬젠 왕 알고리즘은 이에 보조 변수를 도입해서 훨씬 속도를 빠르게 하는 마르코프 연쇄 몬테 카를로 법이다.

각 간선마다 이진 보조 변수 (결합 변수) \mathbf{z}를 도입한 뒤 확장 모델 p(\mathbf{x}, \mathbf{z}) = \frac{1}{Z^{\prime}} \prod_{e} g_{e}(\mathbf{x}_{e}, z_{e}), \sum_{z_{e} = 0}^{1} g_{e}(\mathbf{x}_{e}, z_{e}) = f_{e}(\mathbf{x}_{e})가 되도록 잡으면 \sum_{\mathbf{z}} p(\mathbf{x}, \mathbf{z}) = p(\mathbf{x})이 된다. 즉 이 확장 모델에서 샘플링을 한 뒤 \mathbf{z}를 버리면 된다. 이 모델에 대해서는 깁스 샘플링을 적용하기 쉽다. 각각의 결합 변수가 정점 변수에 대해 조건부독립이므로 전체 조건부 함수가 간선에 대해 인수분해되기 때문이다. 전체 조건부함수 p(z_{e} | \mathbf{x}_{e})의 계산도 쉽다. 간선의 양 끝 정점의 상태가 같다면 z_{e} = 1p = 1 - e^{-2J}의 확률로 놓는다. 아니면 z_{e} = 0으로 놓는다.

p(\mathbf{x} | \mathbf{z})을 샘플링하려면 가동된 결합변수들로부터 정의되는 그래프의 연결요소를 찾은 뒤 그 중 하나를 고르고 그 연결 요소 전체를 무작위 상태로 바꿔버린다.

직관적으로 스웬젠 왕 알고리즘은 깁스 샘플링보다 상태 공간에서 갖는 이동의 크기가 크다. 간선 세기가 온도 T에 대해 \frac{J}{T}일 때 온도가 크면 깁스 샘플링이나 스웬젠 왕 알고리즘이나 똑같이 잘 동작하지만, 온도가 임계 온도 T_{c}에 도달할 때에는 상호 연관 길이가 대단히 길어지고 깁스 샘플링이 독립적인 표본을 추출하는 데 매우 오래 걸리게 된다. 반면, 스웬젠 왕 알고리즘은 모든 온도에서 상태를 빠르게 혼합한다.

단, 간선 세기 중 단 하나라도 음수이면 이 계는 좌절계가 되어 스웬젠 왕 알고리즘이 잘 작동하지 못한다. 인접한 변수들끼리 같은 변수를 갖도록 강요하는데, 이것이 자연스럽지 않기 때문이다.

24.5.4. Hybrid/Hamiltonian MCMC

연속된 상태 공간에서 로그 사후분포의 경사를 구함으로써 마르코프 연쇄 몬테 카를로 샘플링을 하는 방법이 있다. 기본 발상은 매개변수를 공간 내 입자로 보고 보조 변수를 이 입자의 모멘텀으로 정의하는 것이다. 그리고 이 매개변수-모멘텀 쌍을 특정한 규칙에 따라 업데이트하는 것이다. 이를 혼성 마르코프 연쇄 몬테 카를로 또는 해밀토니언 마르코프 연쇄 몬테 카를로법이라 한다. 이 때 사용자는 업데이트를 할 널뛰기 횟수와 그 크기를 지정해 줘야 한다. 이는 추계적 경사 하강법과 혼합되어 큰 데이터셋을 다룰 수도 있다.

24.6. Annealing methods

많은 분포는 최빈값이 여러 개여서 샘플링하기 쉽지 않다. 하지만, 계산적 온도 변수를 통해 분포를 다듬질하는 방법이 있다.

24.6.1. Simulated annealing

모의 담금질법은 모르는 함수 f(\mathbf{x})의 전역 최적해를 찾는 추계적 알고리즘이다. 이는 메트로폴리스-헤이스팅스 알고리즘과 연관이 있으며, 이산적/연속적 최적화에 모두 쓰일 수 있다.

이는 통계물리학에서 비롯되었다. 핵심 통계량은 상태 \mathbf{x}에 존재할 확률을 특정하는 볼츠만 분포 p(\mathbf{x}) \propto e^{-\frac{f(\mathbf{x})}{T}} 인데 (f(\mathbf{x})은 계의 에너지, T는 온도) 온도가 0에 접근할 수록 계는 최소 에너지 상태에 머무르게 된다. 높은 온도에서는 표면이 평면에 가까워지므로 이동을 하기 쉬워진다. 이후 온도를 식힘으로써 작은 고점들을 사라지게 하고 큰 고점들을 돌출시킨다. 이를 통해 가장 큰 고점을 추적함으로써 전역 최적해를 찾을 수 있다. 이를 연속화법이라 한다.

이를 서술하면 다음과 같다. 각 단계마다, 적절한 제안 분포 \mathbf{x}^{\prime} \sim q(\dot | \mathbf{x}_{k})에서 새 상태를 샘플링한다. 이후, \alpha = e^{\frac{f(\mathbf{x}) - f(\mathbf{x}^{\prime})}{T}}을 정의한 뒤, 새 상태를 \min(1, \alpha) 확률로 채택하고 이외 확률로 기각한다 (현상태 유지)

시간에 따라 온도를 변화시키는 냉각 계획을 정확히 정하기는 쉽지 않은데 (이것이 담금질 법의 중요한 단점 중 하나다), 대개 지수적 냉각 계획 T_{k} = T_{0} C^{k}을 사용한다. 대개 T_{0} \sim 1, C \sim 0.8이다.

모의 담금질법 예제.

24.6.2. Annealed importance sampling

모의 담금질 법과 중요도 샘플링을 혼합해 어려운 분포로부터 독립적 표본들을 추출하는 법을 담금질된 중요도 샘플링이라 한다. p_{0}(\mathbf{x}) \propto f_{0}(\mathbf{x})이 샘플링하기 어려운 함수이고 p_{n}(\mathbf{x}) \propto f_{n}(\mathbf{x})은 샘플링하기 쉬운 함수라고 가정하자. 이 때, 온도의 역수 \beta_{j}에 대해 중간 분포 f_{j}(\mathbf{x}) = f_{0}(\mathbf{x})^{\beta_{j}}f_{n}(\mathbf{x})^{1 - \beta_{j}}, 1 > \beta_{0} > \cdots > \beta_{n} = 0를 만들어 p_{n}으로부터 p_{0}에 접근하는 수열들을 만드는 것이다. 거기다 p_{j} 각각에 불변한 마르코프 연쇄 T_{j}(\mathbf{x}, \mathbf{x}^{\prime})들을 가진다고 하자. 이 때, p_{0}으로부터 \mathbf{x}을 샘플링하는 것은 다음과 같이 수행할 수 있다.

먼저 \mathbf{z} = (\mathbf{z}_{n-1}, \cdots, \mathbf{z}_{0})을 다음과 같이 샘플링한다. \mathbf{z}_{n-1} \sim p_{n}을 샘플링하고, \mathbf{z}_{n-2} \sim T_{n-1}(\mathbf{z}_{n-1}, \cdot)을 샘플링하고, 이를 반복해 \mathbf{z}_{0} \sim T_{1}(\mathbf{z}_{1}, \cdot)을 샘플링한다. 마지막으로 \mathbf{x} = \mathbf{z}_{0}으로 두고 가중치 w = \frac{f_{n-1}(\mathbf{z}_{n-1}) }{f_{n}(\mathbf{z}_{n-1})} \cdots \frac{f_{0}(\mathbf{z}_{0})}{f_{1}(\mathbf{z}_{0})}을 부여하면 된다.

24.6.3. Parallel tempering

마르코프 연쇄 몬테 카를로와 담금질을 결합하는 다른 방법은 여러 연쇄를 서로 다른 온도에서 병렬로 수행한 뒤 한 연쇄에서 인접한 연쇄로부터 샘플링을 하는 것이다. 이 방법에서는 고온을 가진 연쇄는 저온을 가진 연쇄에 영향을 줄 수 있게 된다. 이를 병렬 템퍼링이라 한다.

24.7. Approximating the marginal likelihood

주변가능도 p(\mathcal{D} | M) = \int p(\mathcal{D} | \mathbf{\theta}, M) p(\mathbf{\theta} | M) d \mathbf{\theta}를 어떻게 근사할까?

24.7.1. The candidate method

간단한 방법은 후보 방법으로 p(\mathcal{D} | M) = \frac{p(\mathcal{D} | \mathbf{\theta}, M) p(\mathbf{\theta} | M)}{p(\mathbf{\theta} | \mathcal{D}, M)}의 특성을 이용하는 것이다. 이는 임의의 \mathbf{\theta}에 대해 성립하기 때문에 적절한 값을 선택해 p(\mathcal{D} | \mathbf{\theta}, M) p(\mathbf{\theta} | M)을 쉽게 계산할 수 있다. \mathbf{\theta} 근처에서 사후분포 추정을 한다면 분모도 계산할 수 있다. 이는 대개 마르코프 연쇄 몬테 카를로로 근사된다.

이 방법의 단점은 p(\mathbf{\theta} | \mathcal{D}, M)이 사후분포의 모든 최빈값에 대해 주변화된다는 불확실한 가정을 하고 있다는 점인데, 그래서 실 성능은 그다지 좋지 못하다.

24.7.2. Harmonic mean estimate

다음과 같이 근사하는 방법이 있다.

\frac{1}{p(\mathcal{D})} \simeq \frac{1}{S} \sum_{s=1}^{S} \frac{1}{p(\mathcal{D} | \mathbf{\theta}_{s})}

그러나 이 방법도 사후분포로부터 추출된 표본에만 의존한다는 단점이 있다. 사후분포는 사전분포에 민감하지 않지만 가능도는 사전분포에 민감하기 때문에 이도 실 성능은 매우 떨어진다.

24.7.3. Annealed importance sampling

중요도 샘플링을 통해 분할 함수의 비율을 구할 수 있다.

\frac{Z_{0}}{Z_{n}} = \frac{\int f(\mathbf{z}) d \mathbf{z}}{\int g(\mathbf{z}) d \mathbf{z}} \simeq \frac{1}{S} \sum_{s=1}^{S} w_{s} (w_{s} = \frac{f(\mathbf{z}_{s})}{g(\mathbf{z}_{s})})

f_{n}을 사후분포, f_{0}을 사전분포로 넣으면 사전분포의 표준화 상수 Z_{0}을 알 때 Z_{n} = p(\mathcal{D})을 추정할 수 있다.

요점 정리

  • 마르코프 연쇄 몬테 카를로(MCMC) : 고차원 분포에서 샘플링하는 가장 널리 쓰이는 방법으로, 목표하는 확률분포를 정지 분포로 갖는 마르코프 연쇄를 만든 뒤 상태 공간에서 무작위 걸음을 해서 각 상태에서 소비하는 시간이 목표 분포에 비례하도록 하고 이 연쇄로부터 표본을 추출한다.
  • 깁스 샘플링 : 각 단계마다 하나의 변수를 다른 변수들은 분포에 존재한다는 조건하에 샘플링을 하는 것이다. 관측가능한 변수면 샘플링하지 않는다.
  • 메트로폴리스 헤이스팅스 알고리즘 : 각 단계에서 제안 분포의 확률로 현 상태에서 다음 상태 로 이동하는 것
  • 마르코프 연쇄 몬테 카를로 법의 속도와 정확성을 측정하는 여러 지표가 있다.
  • 더미 변수를 추가해 샘플링하는 보조 변수 마르코프 연쇄 몬테 카를로 법이 존재한다.
  • 계를 고온으로 만들어 표면을 평평하게 한 뒤 저온으로 식혀 전역 해를 찾는 담금질 법이 존재한다.
  • 마르코프 연쇄 몬테 카를로법에서 여러 방법으로 주변가능도를 근사할 수 있다.

답글 남기기

아래 항목을 채우거나 오른쪽 아이콘 중 하나를 클릭하여 로그 인 하세요:

WordPress.com 로고

WordPress.com의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Google photo

Google의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Twitter 사진

Twitter의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Facebook 사진

Facebook의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

%s에 연결하는 중