[貝氏] 有關ABC的一些筆記

ABC是Approximate Bayesian Computation的縮寫，名字真的還蠻俏皮的。由於實習的關係，目前為止已經被迫和這個東西相處了三個月，寫一篇筆記 。

與本文無關，純粹覺得好笑

基本上，貝氏統計基本上可以濃縮下面成一行：

$\pi(\theta|x)\propto\pi(\theta)p(y|\theta)$

$\pi(\theta|x)$是所謂的後驗分配(posterior distribution)，基本上是得用手算的，但通常來說不太可能，尤其是模型非常複雜的時候，於是就出現了MCMC (Monte-Carlo Markov chain)還有Variational Bayes之類的方法，讓大家日子在電腦的方著之下可以好過一點(應該說很多)。

但隨著大家的胃口越來越大，像是是在處理基因或是疾病散播問題時，很多時候模型已經複雜到連概似函數(likelihood function) $p(y|\theta)$ 都寫不出來。ABC就是為了這樣的情境而出現的。

Rejection sampling ABC

這是一個最基本的ABC演算法，假設我們有一個simulator $S$，基本上可以把$S$想成一個由參數空間 $\Theta$ 投射到資料空間 $Y$ 的函數$S: \Theta \rightarrow Y$，如果Simulator $S$裡有潛在變數$v$的話simulator則寫成，$S: \Theta \times V \rightarrow Y$

若是手上所收集到的資料是 $y_{0}$ ，而先驗分配prior是 $\pi(.)$ ，演算法的步驟如下

for $i=1$ to $N$ do
    repeat
        從 $\pi(.)$ 抽出$\theta$
        由 $S(\theta)$ 產生 $y_{\theta}$
    until $d(y_{0}, y_{\theta})<\epsilon$
    $\theta^{(i)} \leftarrow \theta$
 end for

 $d: Y \times Y \rightarrow \mathbb{R}$ 是一個距離函數，用來定義生成的 $y_{\theta}$ 和觀察到的 $y_{0}$ 兩筆資料的異同。這個函數可以是Mean square error或是其他特定的函數，而$\epsilon$則是一個很小的容許值。值得注意的是，演算法中完全不需要定義所謂的likelihood函數。

由於likelihood函數的功能在上述的演算法中已經被$d$ 函數所取代，所以如何定義一個好的$d$在ABC中就是一個非常重要的課題了。

普遍來說，直接比較兩筆資料很沒效率，所以一般而言會比較統計量，summary statistics $s_{obs}$ 和$s_{simulated}$，而理論上，如果所選取的$s$是充分統計量(sufficient statistics)的話，$p(\theta|y_{obs}) = p(\theta|s_{obs})$。而統計量的選曲又是一個很大的題目了，這裡先跳過哈哈哈哈。

MCMC-ABC

Rejection sampling是一個很直覺的做法，但是收斂速度往往很慢，因此MCMC-ABC之類方法就發展起來彌補其不足。

Metropolis-Hastings sampler

Metropolis-Hastings sampler [1] (以下簡稱MH)對熟悉貝氏統計的人應該都不算陌生，在進行MH的時候，必須要選定一個proposal distribution $q(.|\theta)$ 用以決定參數在抽樣時移動的行為。

設定起始值$\theta^{(0)}$

for $t=1$ to $N$ do
    從 $q(.|\theta^{(t-1)})$ 抽出$\theta$
    計算 $a = \frac{\pi(\theta) p(y|\theta) q(\theta^{(t-1)}|\theta)}{\pi(\theta^{(t-1)}) p(y|\theta^{(t-1)}) q(\theta|\theta^{(t-1)})}$
    抽出 $u \sim Uniform(0, 1)$
    if $u < a$ then
        $\theta^{(t)} \leftarrow \theta$
    end if
    $\theta^{(t)} \leftarrow \theta^{(t-1)}$
end for

MCMC-ABC

MCMC-ABC [2]等於模仿了MH只是用rejection sampling中比較生成資料的概念取代了需要概似函數的部分

for $t=1$ to $N$ do
    從 $q(.|\theta^{(t-1)})$ 抽出$\theta$
    由 $S(\theta)$ 產生 $y_{\theta}$
    if $d(y_{0}, y_{\theta})<\epsilon$ then
        計算 $a = \frac{\pi(\theta) q(\theta^{(t-1)}|\theta)}{\pi(\theta^{(t-1)}) q(\theta|\theta^{(t-1)})}$
        抽出 $u \sim Uniform(0, 1)$
        if $u < a$ then
            $\theta^{(t)} \leftarrow \theta$
        end if
    end if
    $\theta^{(t)} \leftarrow \theta^{(t-1)}$
end for

搞了半天其實說穿了也沒什麼，呵呵。當然還有很多好事之徒發展了很多變形，這裡就不討論了。

提升效率

因為Simulator通常很複雜，也是最耗時的步驟，所以也有一些方法發展出來減少跑Simulator的次數以增進效率。這裡介紹Regression adjustment和BOLFI

Regression adjustment

一般來說，$\epsilon$是越小越好，最好可以接近零。但是如果真的設定$\epsilon = 0$的話又很花時間，因為會製造成堆沒路用的simulated data。

因此，與其這樣，不如直接統計建模，描述探討$\theta$和summery statistics $s$的關係。

$\theta^{i}=m(s^{i})+e^{i}$

其中$e^{i}$是error term。Beaumont [3]等人提出的方法是用一個簡單的線性回歸建模，

$m(s^{i})=\alpha + \beta^{T} s^{i}$

因此，在建模完成之後，就可以利用該模型回推回去若當 $s^{i} = s_{obs}$ 時， $\theta$ 會在哪裡，意即：

$\theta^{*i} = m(s_{obs})+(e^{i})$

而根據模型，可以估計$e^{i}$

$e^{i} = \theta^{i}-\hat{m}(s^{i})$

我們可以回推

$\theta^{*i}=\hat{m}(s_{obs}) + (\theta^{i} - \hat{m}(s^{i}))$

當然，市面上也有其他更複雜的模型用來描述 $\theta$ 和 $s$ 的關係，這裡就不詳述了。

BOLFI

這個東西的全稱是叫做Bayesian Optimization for Likelihood-Free Inference of Simulator-Based Statistical Method [4]。個人目前不是很喜歡這個東西，但似乎大家都在用，所以稍微提一下(?)

其實就是使用Gaussian Process Optimisation的技巧，去針對$\theta$和$d(y_{\theta}, y_{0})$的關係進行Gaussian process建模，並進行最佳化。

並且，在某些條件下，概似函數(likelihood)是可以由中央極限定理或是無母數的方法估計的，所以BOLFI也提供了一個直接由$\theta$取得likelihood value的架構。


[1] Chib, Siddhartha, and Edward Greenberg. "Understanding the metropolis-hastings algorithm." The american statistician 49.4 (1995): 327-335.
[2] Marjoram, Paul, et al. "Markov chain Monte Carlo without likelihoods."Proceedings of the National Academy of Sciences 100.26 (2003): 15324-15328.
[3] Beaumont, Mark A., Wenyang Zhang, and David J. Balding. "Approximate Bayesian computation in population genetics." Genetics 162.4 (2002): 2025-2035.
[4] Gutmann, Michael U., and Jukka Corander. "Bayesian optimization for likelihood-free inference of simulator-based statistical models." arXiv preprint arXiv:1501.03291 (2015).