MCMC diagnosis: convergence, correlations, CLT, effective sample size (ESS)
MCMC: Pros and Cons
- (+) high dimensional data에서 잘 동작한다.
- (+) Metropolis-Hastings 알고리즘과 같이 general-purpose sampler로 확장이 쉽다
- (+) 구현이 쉬운 편이다
- (-) sequential한 성질 때문에 대규보 데이터로 확장이 어렵다 (not really scalable)
- (-) 어떤 chain이 target distribution에 도달하는지 명확하지 않다.
- (-) 수렴 지표가 명확하지 않다.
그렇다면 무엇이 더 좋은 MCMC 알고리즘으로 만들까?
좋은 MCMC는, high-density 영역에 오래 머물면서(exploitation) i.i.d.를 위해 때때로 다른 영역도 탐색해야한다. (exploration)
예를 들어, chain이 local optima에 머물게 되면 exploration이 부족하여 기댓값의 근사를 잘 구하지 못할 것이다.
즉, markov chain을 어떻게 i.i.d.한지 평가할지 필요하다.
Variance and Correlation
먼저, 우리는 MCMC를 통해 다음의 기댓값을 계산하기 원한다.
\[ \mathbb{E}_{p(x)}[f(x)] \approx \cfrac{1}{n} \sum_{t=1}^{n}f(x^{(t)}) \]
간단히 $y_i := f(x^{(i)})$, $\bar{y} := \frac{1}{n}\sum_{t=1}^{n}y_t$로 정의하자.
그러면 $\bar{y}$의 분산은 다음과 같다.
\begin{align} Var(\bar{y}) &= \cfrac{1}{n^2}\sum_{i=1}^{n} \sum_{j=1}^{n} Cov(y_i, y_j) \\ &= \cfrac{1}{n^2}\sum_{i=1}^{n}Var(y_i) + \cfrac{1}{n^2}\sum_{i \neq j}Cov(y_i, y_j) \\ &= \cfrac{1}{n^2}\sum_{i=1}^{n}Var(y_i) + \cfrac{2}{n^2}\sum_{k=1}^{n-1} \sum_{i=1}^{n-k}Cov(y_i, y_{i+k}) \end{align}
chain이 stationary distribution으로 수렴하면, 즉 $y=f(x)$에 대하여 $(y_i)_{i=1}^{n} \overset{\text{i.i.d.}}{\sim} p(y)$이다.
이때 $\mu = \mathcal{E}[y]$, $\sigma^2 = Var(y)$이라 하자. 그러면 $\sigma^2 = \frac{1}{n} \sum_{i=1}^{n} Var(y_i)$이고 위의 식은 다음으로 근사된다. (chain이 stationary distribution으로 수렴하기 때문에 $=$가 아니라 $\approx$이다)
\begin{align} Var(\bar{y}) &\approx \cfrac{\sigma^2}{n} + \cfrac{2}{n^2} \sum_{k=1}^{n-1} \sum_{i=1}^{n-k}Cov(y_i, y_{i+k}) \\ &\approx \cfrac{\sigma^2}{n}\left( 1 + \cfrac{2}{n}\sum_{k=1}^{n-1} \sum_{i=1}^{n-k} \cfrac{Cov(y_i, y_{i+k})}{\sigma^2} \right) \\ &= \cfrac{\sigma^2}{n} \left( 1+2 \cfrac{n-k}{n}\sum_{k=1}^{n-1} \cfrac{Cov(y_1, y_{1+k})}{\sigma^2} \right) \end{align}
CLT for Markov chain
$(x^{(t)})_{t \ge 1}$가 (homogeneous) Markov chain이고 stationary distribution이 $\pi(x)$라 하자.
그리고 $y_i = f(x_i)$, $\bar{y}=\frac{1}{n}\sum y_i$라 하자.
만약 $p_1(x)=\pi(x)$라면, $\mathbb{E}[f(x)]=\mu < \infty$이고 $\text{Var}(f(x))=\sigma^2 < \infty$인 measurable function $f$에 대하여 다음이 성립한다.
\[ \cfrac{\bar{y} - \mu}{\tilde{\sigma} / \sqrt{n}} \overset{\mathrm{d}}{\to} \mathcal{N}(0, 1) \]
이때, $\tilde{\sigma^2} = \sigma^2 \left( 1 + 2\sum_{k=1}^{\infty}\cfrac{Cov(y_1, y_{1+k})}{\sigma^2} \right)$
Effective Sample Size (ESS)
autocorrelation time at lag $k$를 다음과 같이 정의하자.
\[ \rho_k := \cfrac{Cov(y_1, y_{1+k})}{\sigma^2} \]
그러면 $n \to \infty$이면 위의 $Var(\bar{y})$는 다음으로 수렴한다.
\[ Var(\bar{y}) \approx \cfrac{\sigma^2}{n} \left(1 + 2\sum_{k=1}^{\infty}\rho_k \right) = \cfrac{\sigma^2}{n_{\text{eff}}} \]
이때 $n_{\text{eff}}$를 Effective Sample Size(ESS)라 하고
\[ n_{\text{eff}} := \cfrac{n}{1 + 2\sum_{k=1}^{\infty}\rho_k}\]
이다. 분모를 특별이 Integrated Autocorrelation Time (IAT)라 부른다.
\[ \text{IAT} = 1 + 2\sum_{k=1}^{\infty}\rho_k \]
위 식의 의미는 다음과 같다.
chain이 수렴하여 sample이 independent하더라도, $Var(\bar{y})$는 $\sigma^2 / n$이 아니라 correlation이 존재하여 $\sigma^2 / n_{\text{eff}}$로 수렴한다.
실제로는(in practice) $\rho_k$를 직접 계산하지 않고 아래 식으로 근삿값을 이용한다.
\[ \hat{\rho_k} = \cfrac{\frac{1}{n-k}\sum_{i=1}^{n-k} (y_i-\bar{y})(y_{i+k} - \bar{y})} {\frac{1}{n}\sum_{i=1}^{n}(y_i - \bar{y})^2} \]
근데 $k$값이 커지면 $1/(n-k)$값이 매우 불안정해진다. 그래서 $1/n$으로 대체하여 다음의 식으로 사용한다.
\[ \hat{\rho_k} = \cfrac{ \sum_{i=1}^{n-k} (y_i-\bar{y})(y_{i+k} - \bar{y}) }{ \sum_{i=1}^{n}(y_i - \bar{y})^2 } \]
그러면 ESS의 근삿값은 다음과 같다.
\[ n_{\text{eff}}= \cfrac{n}{1+2\sum_{k=1}^{m} \hat{\rho_k} }, \quad m < n \]
※ 그럼에도, IAT와 ESS는 fast Fourier transform을 이용하여 효율적으로 계산할 수 있다고 한다.
Some terminologies
chain "mixes" well: chain이 stationary distribution으로 수렴해 보인다.
Burn-in: mixing 이전의 iteration. 일반적으로 burn-in 이전의 sample은 기댓값 계산에서 제외한다.
Thining: chain의 모든 $k$번째 sample을 keep하는 것. 이렇게 하는 이유는 더 적은 correlated sample을 얻기 위함이다.
References
