[貝氏] Paper digest (4) 煩死人的Variational Bayes續集 -- Blackbox Variational Inference

這篇接著很久很久寫的煩死人的 Variational Bayes 初探往下寫的。如果用食記部落客的語法，應該改成「再訪！煩死人的Variational Bayes」(超無聊)。

好，會寫這篇的原因是因為最近拜讀了一篇由大神David Blei團隊所寫的這篇文章，順便邊讀邊寫，吸收快又好。題目叫做Blackbox Variational Inference，你沒看錯Variaional Bayes也有黑箱版的，這篇文章刊出的時間是2014年，那年台灣高中生也剛好在反黑箱課綱，真的蠻潮的。

這個方法的賣點是不需要手動計算複雜的積分也可以做Variational Inference。

先複習一下，Variational Bayes的精神就是找一個替代的替身 $q(\theta)$ 並且讓它和posterior $p(\theta, x)$ 越像越好，也就是盡全力把 $KL[q(\theta), p(\theta, x)]$ 變得很小就對了！接下來，

$KL \left[q \left( \theta \right),p \left( \theta \right|x) \right] $

$=\int q \left( \theta \right) \log \frac{q \left( \theta \right)}{p \left( \theta,x \right)} d \theta- \int q \left( \theta \right) \log p \left(x\right) d\theta$

$=\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta-\log p\left(x\right)$

$=E_{q}[q(\theta)] - E_{q}[\log p(\theta, x)]-C$

上式的$ E_{q}[q(\theta)] - E_{q}[\log p(\theta, x)]$取負號以後叫做 ELBO (Evidence lower bound)，也就是

$ ELBO = E_{q}[\log p(\theta, x)]-E_{q}[q(\theta)]$

只要不斷讓ELBO越來越大，  $KL[q(\theta), p(\theta, x)]$ 就會越來越小。

好現在假設我們的其中一個替身叫做$q(\theta|\lambda)$，最佳化 $\lambda$ 的梯度(gradient)就是

$\nabla_{\lambda}\int (\log p(\theta, x)-\log q(\theta|\lambda))q(\theta|\lambda)d\theta$

$=\int \nabla_{\lambda}(\log p(\theta, x)-\log q(\theta|\lambda))q(\theta|\lambda)d\theta$

$+\int \nabla_{\lambda}q(\theta|\lambda)(\log p(\theta, x)-\log q(\theta|\lambda))d\theta$

其中，

$\int \nabla_{\lambda}(\log p(\theta, x)-\log q(\theta|\lambda))q(\theta|\lambda)d\theta$

$=-E_{q}[\nabla_{\lambda}\log q(\theta|\lambda)]=0$

剩下另外一個，因為

$\nabla_{\lambda}q(\theta|\lambda)=q(\theta|\lambda) \nabla_{\lambda}\log q(\theta|\lambda)$

就可以寫成

$\int \nabla_{\lambda}\log q(\theta|\lambda)(\log p(\theta, x)-\log q(\theta|\lambda))q(\theta|\lambda)d\theta$

$=E_{q}[ \nabla_{\lambda}\log q(\theta|\lambda)(\log p(\theta, x)-\log q(\theta|\lambda))]$

這個東西是可以用Monte Carlo integration估計的，只要從 $q$ 中抽樣就可以了。

因此說了一堆，黑箱版本的好處其實就是不需要手動計算複雜的積分(再次強調)也可以找出梯度，然後一次做最佳化就可以了。

[貝氏] 煩死人的 Variational Bayes 初探

先說，搞懂這個東西的過程令人崩潰的程度大概是這樣

羽生結弦，日本天才花式滑冰選手
圖片來源http://abcnews.go.com/Sports/photos/tough-straight-face-midair-22463006/image-22525718

貝氏統計推論中，有些Posterior實在太複雜太難算，用sampling (如MCMC)抽一大堆貌似合理的樣本來進行參數推論是經常使用的替代方式。然而sampling也有既存的問題，比如說常常讓運算速度不夠快，或是完全來自隨機以致收斂性等表現難以掌握等等。總之，學界開始尋找其他逼近(approximation)的方法在某些場合下取代MCMC。

Variational Bayes算是最紅的，許多貝氏建模的文章都會提供相對應的演算細節。

大致上的概念是這樣的，我們感興趣的東西是$p\left ( \theta \right |x)$，並且

$p\left ( \theta \right |x) \propto p\left ( \theta \right )p\left ( x |\theta \right)$

找一個替代品$q\left(\theta\right)$，並且讓這個替代品和posterior很像很像，換句話說就是努力把

$KL\left [ q\left(\theta\right) , p\left( \theta \right | x)\right]$

變得很小很小(沒錯，Kullback-leibler divergence這個小壞壞又出現了)

然後，因為我們有好多個parameter要推論，mean field assumption會是一個把情況變得簡單的假設，也就是所有的參數都是獨立的。也就是

$q\left(\theta \right )=\prod_iq\left(\theta_i \right )$

把這個該變小的Kullback-leibler divergence寫開

$KL \left[q \left( \theta \right),p \left( \theta \right|x) \right] =\int q \left( \theta \right) \log \frac{q\left( \theta \right)}{p \left( \theta|x \right)}d\theta$

$=\int q \left( \theta \right) \log \frac{q\left( \theta \right)p \left(x \right)}{p \left( \theta,x \right)} d\theta$

$=\int q \left( \theta \right) \log \frac{q \left( \theta \right)}{p \left( \theta,x \right)} d \theta- \int q \left( \theta \right) \log p \left(x\right) d\theta$

$=\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta-\log p\left(x\right)$

由於對$\log p\left(x\right)$已經使不上力(和$\theta$無關)，所以接下來專心解決$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta$，並解我們知道$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta$和$KL\left [ q\left(\theta\right) , p\left( \theta \right | x)\right]$成正比，所以要盡量把他變小

接下來引入mean field assumption

$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta =\int \prod_i q\left(\theta_i \right) \log \frac{\prod_i q\left(\theta_i \right)}{p\left(\theta,x\right)} d\theta$

$=\int \prod_i q\left(\theta_i \right) \left( \sum_i \log q\left(\theta_i \right) - \log p\left(\theta, x \right ) \right)$

$=\int q\left(\theta_j \right) \prod_{i\neq j} q\left(\theta_i \right) \left( \sum_i \log q\left(\theta_i \right) - \log p\left(\theta, x \right ) \right) d\theta$

$=\int q\left(\theta_j \right) \prod_{i\neq j} q\left(\theta_i \right) \left( \log q\left(\theta_j \right) - \log p\left(\theta, x \right ) \right) d\theta$ $- \int q\left(\theta_j \right) \prod_{i\neq j} q\left(\theta_i \right) \left( \sum_{i\neq j} \log q\left(\theta_i \right) \right) d\theta$

$= \int q\left(\theta_j \right) \left( \log q\left(\theta_j \right) - \int \prod_{i\neq j} q\left(\theta_i \right) \log p\left(\theta, x \right ) d\theta_{i\neq j} \right ) d\theta_j$$+c$

$= \int q\left(\theta_j \right) \left( \log q\left(\theta_j \right) - E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right) d\theta_j +c$

$= \int q\left(\theta_j \right) \log \frac{\log q\left(\theta_j \right)}{\exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right )} d\theta_j +c$

$=KL\left[ q\left(\theta_j \right), \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right ) \right ] + c$


$KL\left[ q\left(\theta_j \right), \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right ) \right ]$最小就是讓$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta$最小，也就是讓$KL\left [ q\left(\theta\right) , p\left( \theta \right | x)\right]$最小，目標就完成啦!!!!!!! 哈哈哈哈~~~~~(已瘋)

而可以把$KL\left[ q\left(\theta_j \right), \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right ) \right ]$變得最小的$q\left(\theta_j \right)$不就是$\exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right )$本人嗎?

真是可喜可賀!!!

因此Variational Bayes演算法就大功告成了

$q^{*}\left(\theta_j \right) = \frac{1}{Z} \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right )$

或者寫成

$\log q^{*}=E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] - \log Z$

逐個迭代$q \left( \theta_{j} \right )$直到收斂為止

[無關程式] AIC和統計建模

這篇是我最近在讀這本書的一些些筆記，主題是AIC，除非你這個人夠無聊又極富耐性，不然請慎入(是說這個網誌除了我以外大概也沒有其他人會來吧吧吧吧)

沒耐性的就早早洗洗睡了吧

AIC就是Akaike Information Criterion的縮寫，在使用各種常見的統計方法進行建模時，常常都會用它來作為評估模型好壞的指標

顧名思義，這個指標是由一個叫做Akaike的人提出的，基本上不要試著用英文去唸這個字，因為他是日本人。這個人的名字應該要唸成『阿-咖-咿-擠』

噢扯遠了

言歸正傳，AIC基本上長這樣：

$-2log(L)+2p $

其中$log(L)$就是的log-likelihood然後$p$是參數的數目

先來看看看最大概似估計(maximum likelihood estimator)吧：

假設我們現在有一筆資料$x_n$，我們假設有一個模型$f(x|\theta)$是這筆資料最好的描述，我們必須解出$\theta$！

maximum likelihood estimator幫我們找到的是這樣的解 $\widehat{\theta}_n$：

$\widehat{\theta}_n = \operatorname*{arg\,max}_{\theta  \in  \Theta } log\big(f(x_n|\theta)\big)$

然而，夠大的log-likelihood就行了嗎?其實不然，在maximum likelihood estimate的架構下，一味追求把$log\big(f(x_n|\theta)\big)$變大其實會帶來非常不合理的結果

舉個非常簡單的例子，線性迴歸模型(linear regression model)：

$Y=\bf{X^T} \bf{\beta}+\epsilon$       $\epsilon\sim N(0,\sigma^2)$

likelihood function長這樣：

$f(y_i|\bf{x_i};\beta) = \frac{1}{2\pi\sigma^2} exp\big\{ -\frac{1}{2\sigma^2}(y_i-\bf{x^T_i} \bf{\beta})^2\big\}$

log-likelihood如下：

$l(\widehat{\beta}) = \sum_{i=1}^n logf(y_i|x_i;\beta) $

                                        $= - \frac{n}{2} 2log(\pi \sigma^2)  -\frac{1}{2\sigma^2} \sum_{i=1}^n(y_i-\bf{x^T_i}\beta)^2 $ 

如果我們的目的只是要讓$l(\widehat{\beta})$變小，我們就只要無止盡的增加$\bf{x^T_i}\beta $的維度，就可以讓每個$y_i-(x^T_{i1} \beta_1+x^T_{i2}\beta_2+x^T_{i3}\beta_3+...)$變得超小。但是這樣反而會使model裡有太多沒用的預測子，造成over-fitting的問題

因此，單靠夠大的log-likelihood是不夠的!!!

接下來插播Kullback-Leibler Divergence ($KL_{D}$)：

$KL_{D}$是一個表示兩個pdf距離的量，定義是：

$KL_{D} =\int g(z)\frac{g(z)}{f(z)} dx$

其中$g(z)$是true probability 

$f(z)$則是model

其實也就G(Z)是真實pdf的情況下，一個log-likehood的概念

$E_{G}  \left[log \frac{G(Z)}{F(Z)} \right]$

啊這到底和統計建模有筍磨關係？

$KL_{D}$、統計建模和log-likelihood

其實，統計建模的精神就是，我們永遠不知道真實的分配$G(Z)$長怎樣，所以我們就去創造一個叫做$F(Z)$的model來模擬$G(Z)$。換句話說，$F(Z)$和$G(Z)$越接近，我們對$F(Z)$就越有把握

所以我們希望$KL_{D}\big(G,F\big)$很小，最好小到不能再小，這樣就表示我們的$F(Z)$非常完美

$KL_{D}\big(G,F\big)= E_{G}\left[\textrm{log} \frac{G(Z)}{F(Z)} \right]$

$=E_{G}\left[\textrm{log} g(Z)\right]-E_{G}\left[\textrm{log}f(Z)\right]$

$E_{G}\left[\textrm{log} g(Z)\right]$是一個常數，基本上我們對他無能為力。我們只能盡量讓$E_{G}\left[\textrm{log}f(Z)\right]$變大，只要$E_{G}\left[\textrm{log}f(Z)\right]$夠大，就可以讓$KL_{D}\big(G,F\big)$就夠小，我們的model就越完美！

所以，請記得，越好的model的可以使得

$E_{G}\left[\textrm{log}f(Z)\right]$夠大！

log_likelihood的誤差估計

我們回到熟悉的log-likelihood：

$l(\widehat{\theta}) = \sum_i^n \textrm{log}f(x_i|\widehat{\theta})$

$\frac{1}{n} \sum_i^n \textrm{log}f(x_i|\widehat{\theta}) \rightarrow \int \textrm{log}f(z|\widehat{\theta})d \widehat{G}(z) = E_\widehat{G}\left[{log}f(Z|\widehat{\theta})\right]$

所以照理來說，$l(\widehat{\theta})$應該是$nE_\widehat{G}[{log}f(Z|\widehat{\theta})]$)的估計值，但是因為使用了同一筆資料兩次(先估計$\theta$再用來估計$nE_\widehat{G}[{log}f(Z|\widehat{\theta})]$)使得這樣的估計出現了偏誤(bias)

我們把這個bias寫成

$b(G) = E_G( \bf{x_n})\left[\textrm{log}f(\bf{X_n}|\widehat{\theta}(\bf{X_n}))-nE_G(z) \left[\textrm{log} f(Z|\widehat{\theta}(\bf{X_n}))\right]\right]$

假設當$n\rightarrow\infty$時$\theta_n$會逼近真實的值$\theta_0$

我們把剛剛的$b(G)$改寫一下，改成$D1$、$D2$和$D3$的相加

$b(G) = E_G( \bf{x_n} )\left[\textrm{log}f(\bf{X_n}|\widehat{\theta}(\bf{X_n}))-nE_G(z) \left[\textrm{log}f(Z|\widehat{\theta}(\bf{X_n}))\right]\right]$

$=E_G( \bf{x_n})\left[\textrm{log}f(\bf{X_n}| \widehat{\theta}( \bf{X_n} ))-\textrm{log}f(\bf{X_n}|\theta_0)\right]...D1$

$+E_G( \bf{x_n})\left[\textrm{log}f(\bf{X_n}|\theta_0)-nE_G(z) [\textrm{log}f(Z|\theta_0)]\right]...D2$

$E_G( \bf{x_n})\left[nE_G(z) \left[\textrm{log}f(Z|\theta_0)\right]-nE_G(z) [\textrm{log}f(Z|\widehat{\theta}(\bf{X_n}))]\right]...D3$

計算$D1$：

先把$logf(\bf{X_n}|\theta_0)$寫成$l(\theta_0)$

然後搬出泰勒展開式

$l(\theta_0) \sim l(\widehat{\theta}) + (\theta_0-\widehat{\theta})^T\frac{\partial l(\widehat{\theta})}{\partial \theta} + \frac{1}{2}(\theta_0-\widehat{\theta})^T\frac{\partial^2 l(\widehat{\theta})}{\partial \theta \partial \theta^T} (\theta_0-\widehat{\theta})$

因為maximum likelihood estimator基本上就是在找可以使$\frac{\partial l(\widehat{\theta})}{\partial \theta}$為0的$\theta$，所以

$(\theta_0-\widehat{\theta})^T\frac{\partial l(\widehat{\theta})}{\partial \theta}=0$

maximum likelihood estimator在大數法則下：

$\frac{\partial^2 l(\widehat{\theta})}{\partial \theta\partial \theta^T} \rightarrow nJ(\theta_0)$

然後$D1$就可以寫成

$E_G(x_n)\left[l(\widehat{\theta}) - l(\theta)\right]$

$=\frac{n}{2}E_G(x_n)\left[(\theta_0-\widehat{\theta})J(\theta_0)(\theta_0-\widehat{\theta})^T\right]$

$=\frac{1}{2}tr\big\{I(\theta_0)J(\theta_0)^{-1}\big\}$

其中：

$I(\theta_0)=\int g(x)\frac{\partial \textrm{log}f(x|\theta)}{\partial \theta}\frac{\partial \textrm{log}f(x|\theta)}{\partial \theta^T}dx$

$J(\theta_0)=\int g(x)\frac{\partial^2 \textrm{log}f(x|\theta)}{\partial \theta \partial \theta^T}dx$

計算$D2$：

$E_G( \bf{x_n})\left[\textrm{log}f(\bf{X_n}|\theta_0)-nE_G(z) [\textrm{log}f(Z|\theta_0)]\right]$

$=E_G( \bf{x_n})\left[\sum_{i=1}^n\textrm{log}f(\bf{X_i}|\theta_0)\right]-nE_G(z)\left[ [\textrm{log}f(Z|\theta_0)]\right]$

$=0$

計算$D3$：

先把$nE_G(z) \left[\textrm{log}f(Z|wide_theta{\theta})\right]$寫成$\eta(\widehat{\theta})$

然後再度搬出泰勒展開式：

$\eta(\widehat{\theta}) \sim \eta(\theta_0) + \sum_{i=1}^n(\widehat{\theta}_i-\theta_i^{(0)})\frac{\partial \eta(\theta_0)}{\partial \theta_i} +\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^n(\widehat{\theta}_i-\theta_i^{(0)})(\widehat{\theta}_j-\theta_j^{(0)})\frac{\partial^2 \eta(\theta_0)}{\partial \theta_i \partial \theta_j}$

和在計算$D1$時一樣，

$\frac{\partial \eta(\theta_0)}{\partial \theta_i}=E_G(z)\left[\frac{\partial}{\partial theta_i} \textrm{log} f(Z|\theta)\right|\theta_0] = 0,$   $i=1,2,...n$

因此

$\eta(\widehat{\theta}) \sim \eta(\theta_0) - \frac{1}{2}(\theta_0-\widehat{\theta})J(\theta_0)(\theta_0-\widehat{\theta})^T$

$D3$就可以寫成

$nE_G(z) \left[\eta(\theta_0)-\eta(\widehat{\theta})\right]$

$=nE_G(z) \left[\frac{1}{2}(\theta_0-\widehat{\theta})J(\theta_0)(\theta_0-\widehat{\theta})^T\right]$

$=\frac{n}{2}tr \big\{J(\theta_0)E_G(z) \left[(\theta_0-\widehat{\theta})(\theta_0-\widehat{\theta})^T\right]\big\}$

$=\frac{1}{2}tr\big\{I(\theta_0)J(\theta_0)^{-1}\big\}$

現在我們把$D1$、$D2$和$D3$加起來

$b(G)=\frac{1}{2}tr\big\{I(\theta_0)J(\theta_0)^{-1}\big\}+0+\frac{1}{2}tr\big\{I(\theta_0)J(\theta_0)^{-1}\big\}$

$=tr\big\{I(\theta_0)J(\theta_0)^{-1}\big\}$

現在我們回過頭來看看$I(\theta_0)$和$J(\theta_0)$

$-J(\theta_0)=\int g(x)\frac{\partial^2 \textrm{log}f(x|\theta)}{\partial \theta \partial \theta^T}dx$

$=E_G\left[\frac{1}{f(x|\theta)}\frac{\partial^2 }{\partial \theta \partial \theta^T}f(x|\theta)\right]-E_G\left[\frac{\partial}{\partial \theta_i} \textrm{log}f(x|\theta) \frac{\partial}{\partial \theta_j }\textrm{log} f(x|\theta)\right]$

當$g(x)=f(x|\theta_0)$時

$E_G\left[\frac{1}{f(x|\theta)}\frac{\partial^2 }{\partial \theta \partial \theta^T}f(x|\theta)\right]$

$=\int \frac{\partial^2}{\partial \theta_i \partial \theta_j} f(x|\theta_0)dx$

$=\frac{\partial^2}{\partial \theta_i \partial \theta_j}\int  f(x|\theta_0)dx=0$

這個時候$I(\theta_0)$就會等於$J(\theta_0)$

$b(G) =tr\big\{I(\theta_0)J(\theta_0)^{-1}\big\}$

$=tr(I_p) = p$

AIC、log-likelihood和$b(G)$

AIC其實就是

$-2(\textrm{log-likelihood})+2(\textrm{estimator of }b(G))$

因此AIC越小表示model越好

log-likelihood和$b(G)$代入就完成了！

$AIC = -2log(L)+2p$

參考資料：

Konishi, Sadanori, Kitagawa, Genshiro (2007). Information criteria and statistical modeling. Springer