[cite]10.1007/978-1-4757-3437-9_1[/cite]

1 An Introduction to Sequential Monte Carlo Methods

1.1 Motivation

많은 데이터 분석 알고리즘은 실세계 어떤 수치를 관찰된 값으로부터 추정하는 일을 동반합니다. 그리고 어떠한 그와 관련된 사전 지식이 있다면 베이지안 모델을 적용하는 것이 가능해집니다. 베이지안 모델은 그 수치가 나타날 수 있는 사전 분포와 실제 관측된 수치로부터 계산된 likelihood를 이용하여 사후 확률(posterior)을 꾸민 뒤 문제를 풀어내게 되지요.

어떠한 분야에서는 관찰 값들이 시간이 흘러감에 따라 순차적으로 들어오는 경우가 있습니다. 따라서 새로운 관찰 값이 들어오면 그에 따라 사후 확률을 갱신해야 할 것입니다. 이러한 환경에서는 시간이 무한히 흘러가더라도 이전의 데이터를 모두 저장해놓을 필요없이, 계산 또한 간단히 할 수 있는 알고리즘이 필요할 것입니다.

만약 상태 모델이 linear Gaussian state-model으로 가정할 경우 Kalman filter를, 데이터가 일부만 관측되고 finite state-space Markov chain으로 가정하면 HMM filter를 사용할 수 있습니다. 하지만 여전히 실 세계의 데이터는 복잡하고 Guassian의 성질을 가지지 않으며 선형이지도 않습니다.

이제부터 이야기할 Sequential Monte Carlo (SMC) 방법은 시뮬레이션을 통하여 사후 확률 분포를 계산해내는 방법입니다. 문제에 대하여 유연하고 구현이 쉬우며 병렬처리에도 적합합니다. SMC는 bootstrap filters, condenstation, particle filters, Monte carlo filters, interacting particle approximations, survival of the fittest 등으로 불리우고 있습니다.

1.2 Problem statement

먼저 우리가 다룰 모델은 Markovian, nonlinear, non-Gaussian state-space model으로 제한하도록 합니다.

관측되지 않은 값, 숨겨진 상태 : $\{ x_t; t \in N \}, x_t \in \mathit{X}$, Markov process으로 초기 상태 $p(x_0)$로부터 $p(x_t | x_{t-1} )$로 변화합니다.

관측 값 : $\{ y_t ; t \in N^* \}, y_t \in \mathit{Y}$, $x_t$와는 종속적이지 않고, 오직 $p(y_t | x_t)$에만 관계가 있습니다. 관측 값은 초기 상태에는 없고 t=1부터 존재합니다.

그리고 시각 t까지의 순차적인 두 상태 값을 $x_{0:t} \triangleq \{x_0, \dots, x_t \}$, $y_{1:t} \triangleq \{y_1 \dots, y_t \}$ 으로 표현합니다.

우리가 풀고자 하는 최종적인 문제는 매 시각마다 반복적으로 사후 확률 $p(x_{0:t} | y_{1:t})$과 어떤 함수 $f_t$에 대하여 그 기대값(expectation) $I(f_t)$를 계산하는 것입니다.

베이지안 정리에 의하여 사후 확률은 다음과 같이 표현됩니다.

$p(x_{0:t} | y_{1:t}) = \frac{p(y_{1:t} | x_{0:t}) p(x_{0:t})}{\int p(y_{1:t}|x_{0:t})p(x_{0:t})\ dx_{0:t}}$

또한 Markov process이기 때문에 t+1 시각의 사후 확률은 이전 시각 t의 확률분포와의 곱으로 나타내어집니다.

$p(x_{0:t+1} | y_{1:t+1}) = p(x_{0:t} | y_{1:t}) \frac{p(y_{t+1} | x_{t+1}) p(x_{t+1} | x_t)}{p(y_{t+1} | y_{1:t})}$

그리고 Marginal distribution $p(x_t | y_{1:t})$은 다음과 같이 반복적으로 표현됩니다.

$\text{Prediction: } p(x_t | y_{1:t-1}) = \int p(x_t | x_{t-1}) p (x_{t-1} | y_{1:t-1}) dx_{t-1}$

$\text{Updating: } p(x_t | y_{1:t}) = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{\int p(y_t|x_t) p(x_t | y_{1:t-1}) dx_t}$

이 모델은 간단해보이지만, 실제로는 사후 확률 $p(x_{0:t} | y_{1:t})$의 marginal인 $p(y_{1:t})$을 계산할 수가 없습니다.

1.3 Monte Carlo methods

앞에서 언급한 문제를 풀기 위하여 Monte carlo (MC) 방법이 제안되었습니다. 이 방법은 linearity나 Gaussianity라는 제한에 걸려있지 않는 장점이 있습니다. 그 아이디어는 많은 수의 샘플을 우리가 원하는 posterior에서 추출하는 것으로 본래 분포를 추측할 수 있다는 것에서 시작합니다.

1.3.1 Perfect Monte carlo sampling

우리의 관심사인 posterior $p(x_{0:t} | y_{1:t})$으로부터 N개의 샘플 $\{ x_{0:t}^{(i)}, i = 1, \dots, N \}$ 을 이제 파티클이라 부르겠습니다. 이들 샘플은 서로 독립적이면서 같은 분포에서 추출된 것(i.i.d)들 입니다.

이들을 이용하여 추정한 분포는 다음과 같이 표현할 수 있습니다.

$P_N(dx_{0:t} | y_{0:t}) = \frac{1}{N} \sum^N_{t=1} \delta_{x_{0:t}^{(i)}} (dx_{0:t})$

그리고 이를 이용하여 $I(f_t)$를 추정할 수 있습니다.

$I_N(f_t) = \int f_t(x_{0:t}) P_N(dx_{0:t} | y_{0:t}) = \frac{1}{N} \sum^N_{t=1} f_t(x_{0:t}^{(i)})$

샘플의 수 N이 커질 수록 실제 값에 수렴하게 됩니다.

하지만 실제 posterior인 $p(x_{0:t} | y_{1:t})$로부터 샘플을 추출해 내기란 매우 어렵습니다. posterior는 변수가 여러 개일 뿐더러, non-standard이고, 정확한 값을 알기보다는 비례적인 관계만 알 수 있지요. 이를 구현하기 위한 방법으로 Markov chain Monte Carlo (MCMC) 가 있긴 하지만, 반복적으로 수행하기에는 알맞지 않습니다.

1.3.2 Importance sampling

Importance sampling

Posterior로부터 샘플링하지 힘든 점을 해결하기 위하여 Importance sampling 방법이 있습니다. 먼저 어떤 함수 하나를 importance sampling distribution (또는, proposal distribution, importance function) 으로 지정하고 $\pi(x_{0:t} | y_{1:t})$으로 표현합니다.

그렇다면 기대값 식을 아래와 같이 변화시킬 수 있습니다.

$I_N(f_t) = \frac{ \int f_t(x_{0:t}) w(x_{0:t}) \pi(x_{0:t} | y_{1:t}) dx_{0:t}}{\int w(x_{0:t}) \pi(x_{0:t} | y_{1:t}) dx_{0:t}}$

여기서 $w(x_{0:t})$ 를 importance weight라고 부릅니다.

$w(x_{0:t}) = \frac{ p(x_{0:t} | y_{1:t})} {\pi(x_{0:t}|y_{1:t})}$

마찬가지로 Monte Carlo에서의 기대값 $I(f_t)$ 또한 다음과 같이 추정할 수 있습니다.

$\hat{I}_n(f_t) = \frac{ \frac{1}{N} \sum^N_{i=1} f_t(x_{0:t}^{(i)}) w(x_{0:t}^{(i)}) } { \frac{1}{N} \sum^N_{i=1} w(x_{0:t}^{(i)}) } = \sum^n_{i=1} f_t(x_{0:t}^{(i)}) \tilde{w}_t^{(i)} \tilde{w}_t^{(i)} = \frac{w(x_{0:t}^{(i)})} {\sum^N_{j=1}w(x_{0:t}^{(j)})}$

역시 마찬가지로 N이 커질 수록 본래 기대값에 수렴하게 됩니다.

결국 posterior는 importance weight로 추정할 수 있게 됩니다.

$\hat{P}_N(dx_{0:t} | y_{1:t}) = \sum^N_{i=1} \tilde{w}_i^{(i)} \delta_{x_{0:t}^{(i)}} (dx_{0:t})$

Importance sampling은 posterior에서의 샘플링 문제를 해결하였지만, 새로운 관측 데이터가 들어올 때마다 매번 weight를 재계산해줘야 하므로 반복적 계산에는 적당하지 않습니다.

Sequential Importance Sampling

이제 Importance sampling을 이전 결과에 대해서는 수정하지 않도록 수정하여봅시다. 시각 t에서의 Importance function $\pi$는 이전 시각 t-1 의 importance function의 marginal distribution으로 생각할 수 있습니다.

$\pi(x_{0:t} | y_{1:t}) = \pi(x_0) \prod_{k=1}^t \pi(x_k | x_{0:k-1}, y_{1:k})$

Importance weight 또한 수정됩니다.

$\tilde{w}_t^{(i)} \propto \tilde{w}_{t-1}^{(i)} \frac{p(y_t | x_t^{(i)}) p(x_t^{(i)} | x_{t-1}^{(i)})}{ \pi(x_t^{(i)} | x_{i:t-1}^{(i)}, y_{1:t})}$

어떤 경우, prior를 $p(x_{0:t}) = \pi(x_{0:t} | y_{1:t})$로 간주할 수 있습니다. 이러한 경우, $\tilde{w}_t^{(i)} \propto \tilde{w}_{t-1}^{(i)} p(y_t | x_t^{(i)})$를 만족하게 됩니다. 이것은 특별한 경우로 꼭 이렇게 사용해야하는 것은 아니지만, 이제부터 설명할 내용들은 이러한 것을 가정하고 설명할 것이니다.

1.3.3 The Bootstrap filter

SIS 방법의 문제는 시각 t가 늘어날 수록 importance weight가 점점 왜곡된다는 점입니다. 몇번 반복을 하다보면 하나의 파티클만이 weight를 갖고 나머지는 0에 가깝게 되어버립니다. 이를 해결하기 위하여 새로운 방법을 추가합니다.

Notion

아이디어는 다음과 같습니다. 낮은 importance weight를 갖는 파티클은 제거해버리고, 높은 weight를 갖는 파티클을 그만큼 복제하도록 합니다. 식으로 표현하면,

$\hat{P}_N(dx_{0:t} | y_{1:t}) = \sum^N_{i=1} \tilde{w}_i^{(i)} \delta_{x_{0:t}^{(i)}} (dx_{0:t})$

이었던 posterior에서 weight 대신 정수 $N_t^{(i)}$를 사용한 식으로 교체합니다.

$\hat{P}_N(dx_{0:t} | y_{1:t}) = \frac{1}{N} \sum^N_{i=1} N_t^{(i)} \delta_{x_{0:t}^{(i)}} (dx_{0:t})$

여기서 $N_t^{(i)}$은 그에 해당하는 파티클이 다시 몇개로 복제될 것인지 결정하는 수 입니다. 정수이기 때문에 모든 $N_t^{(i)}$를 더하면 N이 되어야 합니다. 만약 파티클 j에 해당하는 $N_t^{(j)} = 0$ 이라면 그 파티클은 제거해야할 대상이 됩니다. 따라서 이러한 동작을 잘 하려면 $N_t^{(i)}$은 posterior에 근거하여 결정을 해야할 것입니다.

$N_t^{(i)}$의 결정단계를 거치고 나면 결과적으로 0보다 큰 값을 가지는 파티클은 그 파티클이 여러번 샘플링된 것과 같은 효과를 지니게 되어 추정하고자 하는 posterior에 계산이 될 것입니다.

Algorithm description

1. Initializsation t=0
   for i = 1, ... N
   $x_0^{(i)} \sim p(x_0)$ // N개의 파티클을 샘플링 t=1
2. Importance sampling step
   For i = 1, ... N
   $x_i^{(i)} \sim p(x_i | x_{i-1})$ // N개의 파티클을 샘플링
   $x_{0:t}^{(i)} = (x_{0:t-1}^{(i)}, x_{0:t}^{(i)})$
   For i = 1, ... N
   $\tilde{w}_t^{(i)} = p(y_t | x_t^{(i)})$ // importance weight를 계산

Normalize importance weights

3. Selection step
   Importance weight에 따라 파티클의 제거와 복제 작업을 수행
   t = t + 1

단계 2, 3을 반복한다.

---