[無關程式] 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