23. Monte Carlo inference

23.1. Introduction

이제까지는 사후분포 추정에 대한 결정론적 알고리즘을 알아보았다. 이들은 베이지안적 접근법의 장점을 누리면서 최적화 기반 점 근사만큼 빠르다. 문제점은 맞는 공식을 이론적으로 유도하기 어렵고 적용 범위가 협소하다는 점이다. 또한 정확도가 근사법의 선택에 크게 의존한다.

이 장에서는 몬테 카를로 근사에 기반한 대안을 다룬다. 발상은 사후분포로부터 표본 \mathbf{x}_{s} \sim p (\mathbf{x} | \mathcal{D})을 추출하고 이로부터 관심이 있는 통계량을 계산하는 것이다. 이런 통계량들은 적절한 함수 f에 대해 \mathbb{E}[f | \mathcal{D}] \simeq \frac{1}{S} \sum_{s=1}^{S} f(\mathbf{x}_{s})로 계산할 수 있다.

표본을 충분히 추출하면 정확도를 원하는 만큼 얻을 수 있다. 이 장에서는 표본을 추출하는 방법 자체도 다룬다.

23.2. Sampling from standard distributions

널리 쓰이는 분포에서 추출하는 법을 알아보자.

23.2.1. Using the cdf

일변수 분포로부터 추출하는 가장 간단한 방법은 역확률 변환이다. F를 원하는 분포의 누적밀도함수라 할 때 다음이 성립한다.

Theorem. U ~ U (0, 1)이 균등확률변수이면 F^{-1}(U) \sim F이다.

즉 일변수분포에 대해서는 누적밀도함수의 역함수를 계산할 수 있을 경우 u \sim U(0, 1)유사난수 생성기로부터 생성한 뒤 x = F^{-1}(u)를 계산하면 된다.

역확률 변환법.

23.2.2. Sampling from a Gaussian (Box-Muller method)

가우시안으로부터 샘플링할 때는 반경이 1인 원에서 무작위로 표본을 추출한 뒤 변수 변환을 써서 구형 2D 가우시안에서의 샘플링으로 변환하는 것이다. 이는 1D 가우시안에서 표본을 2개 추출한 것으로 볼 수 있다. 이를 박스-뮬러 법이라 한다. 다변수 가우시안에 대해서는 공분산 행렬에 대한 촐스키 분해 \mathbf{\Sigma} = \mathbf{L}\mathbf{L}^{T}를 계싼한 뒤 박스-뮬러 법으로 \mathbf{x} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})을 추출하고 \mathbf{y} = \mathbf{L}\mathbf{x} + \mathbf{\mu}을 적용하면 된다.

23.3. Rejection sampling

누적밀도함수의 역함수 법을 쓸 수 없을 경우에는 기각 추출법을 사용한다.

23.3.1. Basic idea

기각 추출법에서는 적당한 M에 대해 Mq(x) \geq \tilde{p}(x)을 만족하는 제안 분포 q(x)을 구성한다. 이후 x \sim q(x), u \sim U(0, 1)을 추출하고, u > \frac{\tilde{p}(x)}{Mq(x)}이면 이를 기각하고 아니면 이를 채택한다. 채택 확률은 \frac{1}{M} \int \tilde{p}(x) dx이므로 M은 부등식을 만족하는 한 최대한 작게 잡아도 상관없다.

23.3.2. Example

감마 분포 Ga(x | \alpha, \lambda)에 대해서는 제안 분포를 Ga(k, \lambda - 1)을 쓸 수 있다. 또는 코시 분포도 쓸 수 있다.

감마 분포에 대한 기각 추출법.

23.3.3. Application to Bayesian statistics

사후분포 p(\mathbf{\theta} | \mathcal{D}) = \frac{p(\mathcal{D} | \mathbf{\theta})p(\mathbf{\theta})}{p(\mathcal{D})}로부터 표본을 추출할 때 \tilde{p}(\mathbf{\theta}) = p(\mathcal{D} | \mathbf{\theta})p(\mathbf{\theta}), q(\mathbf{\theta}) = p(\mathbf{\theta}), M = p(\mathcal{D} | \hat{\mathbf{\theta}}_{\mathrm{MLE}})을 사용하여 기각 추출법을 적용할 수 있다. 사전분포와 사후분포가 크게 차이나면 이 방법은 비효율적이다.

23.3.4 Adaptive rejection sampling

p(x)가 로그 오목이면 타이트한 상한 q(x)을 잡을 수 있다. 로그 밀도함수에 대해 이에 접하는 구간적 선형 함수로 상한을 잡으면 된다. 이를 적응적 기각 추출법이라 한다. 이는 비표준화된 분포에도 적용할 수 있다.

적응적 기각 추출 상한의 예.
1D 가우시안에 대한 적응적 기각 추출의 예.

23.3.5. Rejection sampling in high dimensions

차원의 저주로 인해 기각 추출법은 고차원으로 확장하기 어렵다. 그래서 마르코프 연쇄 몬테 카를로 추출법을 쓰는데, 이와 혼합한 것을 적응적 기각 메트로폴리스 추출법이라 한다.

23.4. Importance sampling

I = \mathbb{E}[f] = \int f(\mathbf{x}) p(\mathbf{x}) d \mathbf{x}을 근사하는 중요도 추출법을 알아보자.

23.4.1. Basic idea

발상은 \mathbf{x}p(\mathbf{x})\lvert f(\mathbf{x}) \rvert와 모두 높은 범위에서 추출하는 것이다. 이는 초효과적이다. 즉, 표본 수가 적게 필요하다. 이는 희귀 사건을 추출할 때 특히 유용하다.

중요도 추출법은 제안 분포 q(\mathbf{x})로부터 추출한 뒤 이를 이용해 기대값 적분을 다음과 같이 추정한다.

\mathbb{E}[f] = \int f(\mathbf{x}) \frac{p(\mathbf{x})}{q(\mathbf{x})}d\mathbf{x} \simeq \frac{1}{S} \sum_{s=1}^{S} w_{s} f(\mathbf{x}_{s}) = \hat{I} (w_{s} = \frac{p(\mathbf{x}_{s})}{q(\mathbf{x}_{s})}중요도 가중치)

이 때 표본은 기각하지 않는다. 그렇다면 제안 분포는 어떻게 선택하나? 자연스러운 방법은 추정 \hat{I} = \sum_{s} \frac{p(\mathbf{x}_{s}) f(\mathbf{x}_{s})}{q(\mathbf{x}_{s})}을 최소화하는 분포를 선택하는 것이다. 이를 계산하면 다음의 최적 중요도 분포를 얻는다.

q^{\ast}(\mathbf{x}) = \frac{\lvert f(\mathbf{x}) \rvert p(\mathbf{x})}{\int \lvert f(\mathbf{x}^{\prime}) \rvert p(\mathbf{x}^{\prime}) d \mathbf{x}^{\prime}}

적절한 f(\mathbf{x})이 생각나지 않는다면 그냥 q(\mathbf{x})p(\mathbf{x})에 가장 가깝게 만들면 된다. 고차원에서는 이는 어렵지만 적응적 방법이 존재한다. 이를 적응적 중요도 추출이라 한다.

23.4.2. Handling unnormalized distributions

비표준화 대상 분포 \tilde{p}(\mathbf{x})는 알 수 있지만 표준화 상수 Z_{p}는 계산할 수 없을 때가 있다. 제안 분포 \tilde{q}(\mathbf{x})도 비표준화시켜야 할 때가 있다. 이 때는 먼저 다음을 계산한다.

\mathbb{E}[f] = \frac{Z_{q}}{Z_{p}} \int f(\mathbf{x}) \frac{\tilde{p}(\mathbf{x})}{\tilde{q}(\mathbf{x})} q(\mathbf{x})d \mathbf{x} \simeq \frac{Z_{q}}{Z_{p} S} \sum_{s=1}^{S} \tilde{w}_{s} f(\mathbf{x}_{s}) (\tilde{w}_{s} = \frac{\tilde{p}(\mathbf{x}_{s})}{\tilde{q}(\mathbf{x}_{s})})

이후 \frac{Z_{p}}{Z_{q}} \simeq \frac{1}{S} \sum_{s=1}^{S} \tilde{w}_{s}로 근사한다. 따라서 w_{s} = \frac{\tilde{w}_{s}}{\sum_{s^{\prime}} \tilde{w}_{s^{\prime}}}이라 할 때 \hat{I} = \sum_{s=1}^{S} w_{s} f(\mathbf{x}_{s})을 얻는다.

이 추정은 비편향 추정이 아니지만 S \to \infty일 때 \hat{I} \to I가 된다.

23.4.3. Importance sampling for a DGM: likelihood weighting

방향그래프모델로부터 중요도 추출을 이용해 표본을 추출하는 법을 알아보자. 먼저 루트 노드를 추출한 뒤, 이들의 자식들을 추출하고, 또 이들의 자식들을 추출하는 법을 생각할 수 있다. 이를 조상 샘플링이라 한다. 방향 비순환그래프의 경우 위상정렬이 가능하기 때문에 이는 가능하다.

증거도가 있는 경우에는 변수가 이산적이라면, 추출한 값이 관측값들과 일관성이 없는 경우에는 이를 기각할 수 있다. 이를 논리적 추출이라 한다. 이는 매우 비효율적이지만, 관측된 변수에 대해서는 그냥 따로 추출을 하지 않고 관측된 값들을 사용함으로써 이를 개선할 수 있다. 이를 가능도 가중이라 한다.

23.4.4. Sampling importance resampling (SIR)

p(x)로부터 중요도 추출을 이용해 비가중 표본을 추출해 p(\mathbf{x}) \simeq \sum_{s} w_{s} \delta_{\mathbf{x}_{s}} (\mathbf{x}) 형태의 분포를 얻을 수 있다. 이 분포로부터 재추출을 한 뒤 먼저 추출했던 표본들을 대체할 수 있는데, 이를 추출 중요도 재추출이라 한다.

이는 저차원 베이지안 추론에 사용할 수 있지만, 사전분포와 사후분포간 차이가 크다면 역시 비효율적이다.

23.5. Particle filtering

입자 거르기(PF)는 몬테 카를로, 또는 시뮬레이션 기반 알고리즘으로 재귀적 베이지안 추론을 수행한다. 즉, 예측-업데이트 루프를 근사한다.

23.5.1. Sequential importance sampling

기본 발상은 믿음 상태를 입자들의 가중치 집합으로 근사하는 것이다.

p(\mathbf{z}_{1 : t} | \mathbf{y}_{1 : t}) \simeq \sum_{s=1}^{S} \hat{w}_{t, s} \delta_{\mathbf{z}_{1 : t, s}}(\mathbf{z}_{1 : t})

이로부터 가장 최근 상태의 주변분포는 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}) \simeq \sum_{s=1}^{S} \tilde{w}_{t, s} \delta_{\mathbf{z}_{t, s}}(\mathbf{z}_{t})로 근사할 수 있다.

이는 간단하지만 잘 작동하지는 못한다.

23.5.2. The degeneracy problem

대부분의 입자들의 가중치가 매우 작으므로 몇 단계를 거치고 나면 기본적인 순차적 중요도 샘플링 알고리즘은 잘 작동하지 못한다. 이를 퇴화 문제라 한다. 이는 고차원에서특히 더 심해진다. 퇴화 정도를 유효 표본수 S_{\mathrm{eff}} = \frac{S}{1 + \mathrm{var}[w_{t, s}^{\ast}]} \simeq \frac{1}{\sum_{s=1}^{S} (w_{t, s})^{2}} (w_{t, s}^{\ast} = \frac{p(\mathbf{z}_{t, s}  | \mathbf{y}_{1 : t})}{q(\mathbf{z}_{t, s} | \mathbf{z}_{t-1, s}, \mathbf{y}_{t})})로 수치화할 수 있다.

퇴화 문제에 대한 해결책은 크게 2가지가 있다. 재추출 단계를 추가하거나, 좋은 제안 분포를 쓰는 것이다.

23.5.3. The resampling step

순차적 중요도 추출법을 개선하는 방법은 유효 표본수가 일정 기준 이하로 내려가면 낮은 가중치의 입자들을 버리고 남은 입자들의 복제를 생성하는 것이다. (그래서 이를 적자 생존, 회춘이라고도 한다) 이러한 재추출 법은 여러가지인데, 가장 간단한 것은 다항 재추출이고, 그 외에 체계적 재추출 잔여 재추출, 계층 재추출 등이 있다. 모두 선형 시간을 가진다.

이는 퇴화 문제를 완화시키지만 고유한 문제가 또 존재한다. 먼저 같은 표본이 계속 반복되어 추출되므로 표본의 다양성이 줄어든다(표본 빈곤화). 이를 해결하기 위한 여러 방법들이 있는데, 첫째는 매번 재추출을 하는 붓스트랩 필터를 하지 않고 필요할 때만 재추출을 하는 것이고, 둘째는 이전 입자들을 복제한 뒤 이를 그대로 사용하지 않고 사후 분포 조건을 유지하는 마르코프 연쇄 몬테 카를로 법으로 추출해 사용하는 것이고 (재추출-이동), 셋째는 적당한 다듬질 커널 \kappa에 대해 커널 밀도 추정 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}) \simeq \sum_{s=1}^{S} w_{t, s} \kappa(\mathbf{z}_{t} - \mathbf{z}_{t, s})을 사용한 뒤 이 다듬질된 분포에서 추출하는 것 (정규화 입자 거르기)이다. 넷째는 추론을 할 때 인공적인 노이즈를 더하는 것이다.

23.5.4. The proposal distribution

가장 자연스러운 제안 분포는 사전분포이다. q(\mathbf{z}_{t} | \mathbf{z}_{t-1, s}, \mathbf{y}_{t}) = p(\mathbf{z}_{t} | \mathbf{z}_{t-1, s})

이 경우 가중치 업데이트는 w_{t, s} \propto w_{t-1, s}p(\mathbf{y}_{t} | \mathbf{z}_{t, s})로 간소화된다. 이는 응결 알고리즘에 유용하다.

다만 가능도가 사전분포에 비해 협소하면 이는 매우 비효율적이 된다. 많은 입자들이 낮은 가중치를 배정받기 때문이다. 이를 해결하기 위해 데이터 \mathbf{y}_{t}를 고려해

q(\mathbf{z}_{t} | \mathbf{z}_{t-1, s}, \mathbf{y}_{t}) = \frac{p(\mathbf{y}_{t} | \mathbf{z}_{t}) p(\mathbf{z}_{t} | \mathbf{z}_{t-1, s})}{p(\mathbf{y}_{t} | \mathbf{z}_{t-1, s})}을 고려할 수 있다. 이 경우 가중치 업데이트는 w_{t, s} \propto w_{t-1, s} p(\mathbf{y}_{t} | \mathbf{z}_{t-1, s}) = w_{t-1, s} \int p(\mathbf{y}_{t} | \mathbf{z}_{t}^{\prime}) p(\mathbf{z}_{t}^{\prime} | \mathbf{z}_{t-1, s}) d \mathbf{z}_{t}^{\prime}로 이루어진다. 이 제안분포는 최적해이다.

일반적으로 이 최적해는 계산가능하지 않지만 \mathbf{z}_{t}가 이산적이거나 p(\mathbf{z}_{t} | \mathbf{z}_{t-1, s}, \mathbf{y}_{t})이 가우시안이면 계산가능하다.

모델이 선형-가우시안이 아니어도 p(\mathbf{z}_{t} | \mathbf{z}_{t-1, s}, \mathbf{y}_{t})에 대한 가우시안 근사를 수행한 뒤 이를 제안 분포로 쓸 수 있다. 이를 무취 입자 거르기라 한다. 더 일반적으로는 여러 데이터 기반 제안들을 할 수 있다.

23.5.5. Application: robot localization

점유 격자 위에서 움직이는 로봇에 대해 입자 거르기를 통해 위치를 근사할 수 있다. 이를 몬테 카를로 위치 추정이라 한다. 균일 사전분포로부터 시작하는 전역 위치 추정을 할 수도 있고 지각 에일리어싱을 수행할 수도 있다.

23.5.6. Application: visual object tracking

물체를 추적할 때에는 가능도 모델에 대한 색 히스토그램에 대해 바타차르야 거리를 이용해 이들을 비교해 모델링할 수 있다. 대안으로는 감지에 의한 추출이 있다.

23.5.7. Application: time series forecasting

추계적 휘발성 모델 같은 시계열 모델에도 적용할 수 있다.

23.6. Rao-Blackwellised particle filtering (RBPF)

은닉 변수를 \mathbf{q}_{t}, \mathbf{z}_{t}로 나누어 \mathbf{q}_{1 : t}를 알고 있을 때 해석적으로 \mathbf{z}_{t}를 적분해 빼낼 수 있다. 즉, \mathbf{q}_{1 : t}만 추출하면 된다는 것이다. 이로부터 유도되는 혼성 입자를 분산 입자 또는 붕괴 입자라 한다.

이 방식의 이점은 추출하는 차원수와 추정의 분산을 줄일 수 있다는 점이다. 그래서 이 방법은 라오-블랙웰화 입자 거르기(RBPF)라고 불린다.

23.6.1. RBPF for switching LG-SSMs

라오-블랙웰화 입자 거르기의 자연스러운 적용례는 전환 선형 동역학계 모델이다. 사전분포로부터 제안 q(q_{t} = k | q_{t-1, s})을 수행할 때 가중치 업데이트는 w_{t, s} \propto w_{t-1, s} p(\mathbf{y}_{t} | q_{t}= k, \mathbf{q}_{1 : t-1, s}, \mathbf{y}_{1 : t-1}) = w_{t-1, s} L_{t, k, s}이 된다. (L_{tk, s} = \int p(\mathbf{y}_{t} | q_{t} = k, \mathbf{z}_{t}) p(\mathbf{z}_{t} | q_{t} = k, \mathbf{y}_{1 : t-1} \mathbf{q}_{1 : t-1, s})d \mathbf{z}_{t})

이 때 L_{tk, s}는 새 관측값 \mathbf{y}_{t}을 목표와 q_{t}=k와 이전 기록 \mathbf{q}_{1 : t-1, s} 에 대해 조건을 건 예측밀도값이다. 이를 계산하는 법을 칼만 거르기 혼합이라 한다. K가 작을 때에는 최적의 제안 분포를 계산할 수 있는데, 이는 p(q_{t} = k | \mathbf{y}_{1 : t-1} \mathbf{q}_{1 : t-1, s}) = \frac{L_{tk, s}p(q_{t} = k | q_{t-1, s})}{\sum_{k^{\prime}} L_{tk^{\prime}, s} p(q_{t} = k^{\prime} | q_{t-1, s})}가 된다. 이로부터의 가중치 업데이트는 w_{t, s} \propto w_{t-1, s} p(\mathbf{y}_{t} | \mathbf{q}_{1 : t-1, s}, \mathbf{y}_{1 : t-1})이 된다.

이 때 가중치들은 q_{t}로부터 추출된 새 값과 독립적이기 때문에, 가중치를 먼저 계산하고 어떤 입자로 알고리즘을 진행할지 선택할 수 있다. 이를 예견 라오-블랙웰화 입자 거르기라 한다.

상태 공간이 이산적이라는 사실을 이용해 알고리즘을 더 최적화하는 것도 가능하다.

23.6.2. Application: tracking a maneuvering target

전환 선형 동역학계 모델은 궤도 조정 대상 추적 등에 쓰일 수 있다.

물체 위치 추적에 대한 입자 거르기 vs 라오-블랙웰화 입자 거르기.

23.6.3. Application: Fast SLAM

동시 위치 추정 및 매핑을 더 빠르게 수행하는 FastSLAM도 가능하다.

요점 정리

  • 사후분포로부터 표본을 추출하고 이로부터 관심이 있는 통계량을 계산하는 몬테 카를로 근사를 할 수 있다.
  • 가장 쉬운 방법은 역확률 변환이다.
  • 역확률 변환을 쓸 수 없으면 기각 추출법을 사용한다.
  • 함수값과 가능도가 모두 높은 영역에서 추출하는 중요도 추출법도 유용하다.
  • 입자 거르기는 재귀적 베이지안 추론을 수행하는 몬테 카를로 기반 알고리즘이다.
  • 입자 거르기에서 은닉 변수의 일부만 추출하고 나머지는 이로부터 추론하는 라오-블랙웰화 입자 거르기도 가능하다.

답글 남기기

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

WordPress.com 로고

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

Google photo

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

Twitter 사진

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

Facebook 사진

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

%s에 연결하는 중