贝叶斯模型评价
赤池信息准则（AIC）
记 $M_j={f(x;\theta_j):\theta_j\in\Theta_j}$，$j=1,\dots,r$，$\hat{\theta}_j$ 为模型 $M_j$ 下参数的极大似然估计。简记 $\hat{f}_j(x)=f(x;\hat{\theta}_j)$，我们可以用KL距离衡量其与真实分布 $g(x)$ 之间的差距
$\displaystyle KL(g\|\hat{f}_j)=\int g(x)\log\dfrac{g(x)}{\hat{f}_j(x)}\text{d}x=\int g(x)\log g(x)\text{d}x-\int g(x)\log\hat{f}_j(x)\text{d}x,$
其中第一项与 $j$ 无关，因此对 $j=1,\dots,r$，最小化 $KL(g|\hat{f}_j)$ 等价于最大化
$\displaystyle K_j=\int g(x)\log\hat{f}_j(x)\text{d}x.$
它的一个直观估计为
$\displaystyle\overline{K}_j=\dfrac{1}{n}\sum_{i=1}^n\log f(x_i;\hat{\theta}_j)=\dfrac{\ell_j(\hat{\theta}_j)}{n}.$
但是这个估计的偏倚很大，因为样本被使用了两次：一次估计参数 $\theta_j$，一次估计积分。
Akaike（1974）注意到
$\begin{array}{rl}
K\approx&\!\!\!\!\displaystyle\int g(x)\left(\log f(x;\theta_0)+(\hat{\theta}-\theta_0)^{\text{T}}s(x,\theta_0)+\dfrac{1}{2}(\hat{\theta}-\theta_0)^{\text{T}}H(x,\theta_0)(\hat{\theta}-\theta_0)\right)\text{d}x\\
=&\!\!\!\! K_0-\dfrac{1}{2n}Z_n^{\text{T}}JZ_n,
\end{array}$
其中得分函数 $s(x,\theta)=\dfrac{\partial\log f(x;\theta)}{\partial\theta}$ 在真值 $\theta_0$ 处为零，$H(x,\theta)=\dfrac{\partial^2\log f(x;\theta)}{\partial\theta\partial\theta^\text{T}}$ 为Hessian矩阵，而
$K_0=\displaystyle\int g(x)\log f(x;\theta_0)\text{d}x,\quad Z_n=\sqrt{n}(\hat{\theta}-\theta_0),\quad J=-\mathbb E[H(X,\theta_0)].$
另一方面，
$\begin{array}{rl}
\overline{K}\approx&\!\!\!\!\displaystyle\dfrac{1}{n}\sum_{i=1}^n\left(\log f(x_i;\theta_0)+(\hat{\theta}-\theta_0)^{\text{T}}s(x_i,\theta_0)+\dfrac{1}{2}(\hat{\theta}-\theta_0)^{\text{T}}H(x_i,\theta_0)(\hat{\theta}-\theta_0)\right)\\
=&\!\!\!\! K_0+\Lambda_n+\dfrac{1}{n}Z_n^\text{T}S_n-\dfrac{1}{2n}Z_n^{\text{T}}J_nZ_n\\
\approx&\!\!\!\! K_0+\Lambda_n+\dfrac{1}{n}Z_n^\text{T}S_n-\dfrac{1}{2n}Z_n^{\text{T}}JZ_n,
\end{array}$
其中 $\displaystyle\Lambda_n=\dfrac{1}{n}\sum_{i=1}^n(\log f(x_i;\theta_0)-K_0)$，而
$\displaystyle J_n=-\dfrac{1}{n}\sum_{i=1}^n H(x_i,\theta_0)\overset{P}{\longrightarrow}J,\quad S_n=\dfrac{1}{\sqrt{n}}\sum_{i=1}^ns(x_i,\theta_0)\overset{d}{\longrightarrow}N(0,V).$
因此利用 $Z_n\approx J^{-1}S_n$，有
$\mathbb E(\overline{K}-K)\approx\overline{E}(\Lambda_n)+\mathbb E\left(\dfrac{1}{n}Z_n^\text{T}S_n\right)=0+\dfrac{1}{n}\textsf{tr}(J^{-1}V).$
如果模型是正确的，则 $J=V$，因而
$K\approx\overline{K}-\dfrac{d}{n},$
这里 $d$ 是参数 $\theta$ 的维数。
下面给出AIC的表达式，即
$AIC(j)=-2n\hat{K}_j=-2\ell_j(\hat{\theta_j})+2d_j,$
其中 $\ell(\theta)$ 是似然函数，$d$ 是参数的维数。我们希望其值越小越好。
贝叶斯信息准则（BIC）
基于贝叶斯的想法，由于
$\displaystyle f(\pmb x\mid M_j)=\int_{\Theta_j}f(\pmb x\mid M_j,\theta_j)\pi_j(\theta_j)\text{d}\theta_j=\int_{\Theta_j}L(\theta_j)\pi_j(\theta_j)\text{d}\theta_j,$
而 $\mathbb P(M_j\mid\pmb x)=f(\pmb x\mid M_j)\mathbb P(M_j)$，因此我们只要考虑选择最大化
$\displaystyle\log\int_{\Theta_j}L(\theta_j)\pi_j(\theta_j)\text{d}\theta_j+\log\mathbb P(M_j)$
的那个 $j$ 即可。使用Laplace近似有
$\displaystyle\int_{\Theta}L(\theta)\pi(\theta)\text{d}\theta\approx L(\hat{\theta})(2\pi)^dn^{-d/2}\vert\Delta(\hat{\theta})\vert^{-1/2}\pi(\hat{\theta}),$
忽略 $n\rightarrow\infty$ 时的有界项后有
$\displaystyle\log\int_{\Theta_j}L(\theta_j)\pi_j(\theta_j)\text{d}\theta_j+\log\mathbb P(M_j)\approx\ell_j(\hat{\theta}_j)-\dfrac{d_j}{2}\log n.$
下面给出BIC的表达式，即
$BIC(j)=-2\ell_j(\hat{\theta_j})+d_j\log n,$
其中 $\ell(\theta)$ 是似然函数，$d$ 是参数的维数。我们希望其值越小越好。
贝叶斯预测信息准则（BPIC）
Ando（2007）提出了最大化对数似然的后验均值
$\displaystyle\eta(G)=\int_{Z}\left\{\int_{\Theta}\log f(z\mid\theta)\pi(\theta\mid\pmb x)\text{d}\theta\right\}\text{d}G(z),$
其中 $\pmb x$ 为样本观测值，$G$ 为真实模型。利用经验分布函数，一个自然的估计是
$\displaystyle\eta(\hat{G})=\dfrac{1}{n}\int_\Theta\log f(\pmb x\mid\theta)\pi(\theta\mid\pmb x)\text{d}\theta.$
这里估计模型参数和计算期望对数似然的后验均值使用了相同的样本数据，导致其结果是有偏的，其偏差定义为
$b(G)=\displaystyle\int[\eta(\hat{G})-\eta(G)]\text{d}G(\pmb x),$
这里 $G(\pmb x)$ 为样本 $\pmb x$ 的密度。如果有偏差的估计 $\hat{b}(G)$，那么 $\eta(G)$ 的一个偏差修正的估计为
$\displaystyle\eta(G)=\dfrac{1}{n}\int_\Theta\log f(\pmb x\mid\theta)\pi(\theta\mid\pmb x)\text{d}\theta-\hat{b}(G).$
这个估计量常常表示为
$\displaystyle IC=-2\int_\Theta\log f(\pmb x\mid\theta)\pi(\theta\mid\pmb x)\text{d}\theta+2n\hat{b}(G).$
这里第一项度量了模型的拟合程度，第二项则是对模型复杂程度的一个惩罚。
假设参数模型 $f(\pmb x\mid\theta)$ 包含了真实的模型 $g(x)=f(\pmb x;\theta_0)$，以及 $\log\pi(\theta)=O_p(1)$，Ando（2007）证明了渐近偏差为 $\hat{b}=p/n$，因此提出了贝叶斯预测信息准则（BPIC）：
$\displaystyle BPIC=-2\int_\Theta\log f(\pmb x\mid\theta)\pi(\theta\mid\pmb x)\text{d}\theta+2p,$
其中 $p$ 为模型中的参数个数。我们希望其值越小越好。
在实际应用中，对数似然的后验均值往往没有解析表达，此时可以使用蒙特卡洛逼近：
$\displaystyle\int_\Theta\log f(\pmb x\mid\theta)\pi(\theta\mid\pmb x)\approx\dfrac{1}{L}\sum_{i=1}^L\log f(\pmb x\mid\theta^{(i)}),$
其中 $\theta^{(1)},\dots,\theta^{(L)}$ 是从后验分布 $\pi(\theta\mid\pmb x)$ 中抽取的样本。
偏差信息准则（DIC）
偏差（Deviance）的经典定义为
$D(\theta)=-2\log f(\pmb x\mid\theta)+2\log f_S(\pmb x),$
其中 $f_S(\pmb x)$ 是仅依赖于样本的标准化项，可以视为饱和模型的似然函数最大值，在模型选择中不起作用，常忽略之。
Spiegelhalter等（2002）指出对数似然的后验均值 $\overline{D}=\mathbb E[D(\theta)\mid\pmb x]$ 可以作为模型拟合程度的一个贝叶斯度量，其值越大表示模型拟合数据的程度越差。此外，还定义了有效参数个数来刻画模型的复杂程度：
$\displaystyle p_D=\overline{D}-D(\overline{\theta})=2\log f(\pmb x\mid\overline{\theta})-2\int_\Theta\log f(\pmb x\mid\theta)\pi(\theta\mid\pmb x)\text{d}\theta,$
其中 $\overline{\theta}$ 为后验均值。$p_D$ 越大，则模型拟合数据越容易。
Spiegelhalter等（2002）定义偏差信息准则（DIC）为
$DIC=p_D+\overline{D}=D(\overline{\theta})+2p_D.$
这实际上可以看作是AIC准则的一个推广：对非分层模型而言，$p\approx p_D$，$\hat{\theta}\approx\overline{\theta}$，从而 $DIC\approx AIC$。另外，DIC与BIC有如下比较：
- BIC试图来识别真实的模型，而DIC并没有假设“真实模型”，其目的在于考察短期的预测能力；
- BIC要求指定一些参数，而DIC估计有效参数个数；
- BIC提供了一种方法来进行模型平均，而DIC则没有。
模型平均
假设我们要基于数据 $D={Y_1,\dots,Y_n}$ 预测一个新观测 $Y$，有不同的候选模型 $M_j$，各模型先验概率相同，则贝叶斯模型平均方法使用下式进行预测：
$\displaystyle p(y\mid D)=\sum_j p(y\mid D,M_j)\mathbb P(M_j\mid D),$
其中
$\mathbb P(M_j\mid D)=\dfrac{\int_{\Theta_j}L(\theta_j)\pi_j(\theta_j)\text{d}\theta_j}{\sum_i\int_{\Theta_i}L(\theta_i)\pi_i(\theta_i)\text{d}\theta_i}\approx\dfrac{\exp(BIC(j))}{\sum_i\exp(BIC(i))}.$
记 $\hat{Y}_j=\mathbb E[Y\mid D,M_j]$，则后验均值和方差为
$\displaystyle\mathbb E[Y\mid D]=\sum_j\hat{Y}_j\mathbb P(M_j\mid D),\\
\displaystyle\text{Var}[Y\mid D]=\sum_j\left(\text{Var}[Y\mid D,M_j]+\hat{Y}_j^2\right)\mathbb P(M_j\mid D)-\mathbb E^2[Y\mid D].$
例如，考虑线性模型
$y=\alpha_\gamma+X_\gamma\beta_\gamma+\varepsilon,\quad\varepsilon\sim N_n(0,\sigma^2I),$
其中 $X_\gamma\in{X}$，假设 $X$ 包含了 $K$ 个可能有用的协变量，且做了标准化处理。
贝叶斯模型平均估计所有 $2^K$ 个变量组合，然后使用后验模型概率进行平均。这里，第 $\gamma$ 个模型的后验概率为
$\mathbb P(M_\gamma\mid y,X)=\dfrac{p(y\mid M_\gamma,X)\mathbb P(M_\gamma)}{\sum_{\gamma'=1}^{2^K}p(y\mid M_{\gamma'},X)\mathbb P(M_{\gamma'})},$
从而对任何量 $\theta$（比如回归系数 $\beta$），有
$\displaystyle p(\theta\mid y,X)=\sum_{\gamma=1}^{2^K}p(\theta\mid M_\gamma,y,X)\mathbb P(M_\gamma\mid y,X).$
在正态线性模型框架下，$\beta_\gamma$ 的先验分布常取为Zellner’s g：
$\beta_\gamma\mid\alpha_\gamma,\sigma^2\sim N_p(\tilde{\beta}_\gamma,g\sigma^2(X_\gamma^\text{T}X_\gamma)^{-1}).$
而 $(\alpha_\gamma,\sigma^2)$ 的一个无信息先验为 $\pi(\alpha_\gamma,\sigma^2)\propto\sigma^{-2}$，于是参数的联合先验密度为
$\pi(\alpha_\gamma,\beta_\gamma,\sigma^2)\propto(\sigma^2)^{-p/2-1}\exp\left\{-\dfrac{1}{2g\sigma^2}(\beta_\gamma-\tilde{\beta}_\gamma)^\text{T}X_\gamma^\text{T}X_\gamma(\beta_\gamma-\tilde{\beta}_\gamma)\right\}.$
而样本的联合密度函数为
$\begin{array}{rl}
f(y\mid\alpha_\gamma,\beta_\gamma,\sigma^2)\propto&\!\!\!\!(\sigma^2)^{-n/2}\exp\left\{-\dfrac{1}{2\sigma^2}(y-\alpha_\gamma1_n-X_\gamma\beta_\gamma)^\text{T}(y-\alpha_\gamma1_n-X_\gamma\beta_\gamma)\right\}\\
\propto&\!\!\!\!(\sigma^2)^{-n/2}\exp\left\{-\dfrac{1}{2\sigma^2}(y-\overline{y}1_n-X_\gamma\beta_\gamma)^\text{T}(y-\overline{y}1_n-X_\gamma\beta_\gamma)\right\}\\
&\quad\times\exp\left\{-\dfrac{n}{2\sigma^2}(\overline{y}-\alpha_\gamma)^2\right\}\\
\propto&\!\!\!\!(\sigma^2)^{-n/2}\exp\left\{-\dfrac{1}{2\sigma^2}(\beta_\gamma^\text{T}X_\gamma^\text{T}X_\gamma\beta_\gamma-2y^\text{T}X_\gamma\beta_\gamma)\right\}\\
&\quad\times\exp\left\{-\dfrac{1}{2\sigma^2}(y-\overline{y}1_n)^\text{T}(y-\overline{y}1_n)\right\}\times\exp\left\{-\dfrac{n}{2\sigma^2}(\overline{y}-\alpha_\gamma)^2\right\},
\end{array}$
所以参数的联合后验密度为
$\begin{array}{rl}
\pi(\alpha_\gamma,\beta_\gamma,\sigma^2\mid y)\propto&\!\!\!\!(\sigma^2)^{-(n+p)/2-1}\exp\left\{-\dfrac{1}{2\sigma^2}(\beta_\gamma^\text{T}X_\gamma^\text{T}X_\gamma\beta_\gamma-2y^\text{T}X_\gamma\beta_\gamma)\right\}\\
&\quad\times\exp\left\{-\dfrac{1}{2\sigma^2}(y-\overline{y}1_n)^\text{T}(y-\overline{y}1_n)\right\}\\
&\quad\times\exp\left\{-\dfrac{n}{2\sigma^2}(\overline{y}-\alpha_\gamma)^2\right\}\times\exp\left\{-\dfrac{1}{2g\sigma^2}\tilde{\beta}_\gamma^\text{T}X_\gamma^\text{T}X_\gamma\tilde{\beta}_\gamma\right\}\\
&\quad\times\exp\left\{-\dfrac{1}{2g\sigma^2}(\beta_\gamma^\text{T}X_\gamma^\text{T}X_\gamma\beta_\gamma-2\beta_\gamma^\text{T}X_\gamma^\text{T}X_\gamma\tilde{\beta}_\gamma)\right\},
\end{array}$
从中可以看到
$\alpha_\gamma\mid\sigma^2,y\sim N_1\left(\overline{y},\dfrac{\sigma^2}{n}\right),\\
\beta_\gamma\mid\sigma^2,y\sim N_p\left(\dfrac{g}{g+1}\hat{\beta}_\gamma+\dfrac{1}{g+1}\tilde{\beta}_\gamma,\dfrac{\sigma^2g}{g+1}(X_\gamma^\text{T}X_\gamma)^{-1}\right),$
其中 $\hat{\beta}\gamma=(X\gamma^\text{T}X_\gamma)^{-1}X_\gamma^\text{T}y$ 为 $\beta_\gamma$ 的极大似然估计，当然也可以看到 $\sigma^2$ 的后验分布为
$\sigma^2\mid y\sim\Gamma^{-1}\left(\dfrac{n-1}{2},s^2+\dfrac{(\tilde{\beta}_\gamma-\hat{\beta}_\gamma)^\text{T}X_\gamma^\text{T}X_\gamma(\tilde{\beta}_\gamma-\hat{\beta}_\gamma)}{g+1}\right),$
这里 $s^2=(y-y1n-X\gamma\hat{\beta}\gamma)^\text{T}(y-y1_n-X\gamma\hat{\beta}_\gamma)$。
进一步可以求得 $y$ 的边际密度为
$f(y\mid M_\gamma,X,g)=\dfrac{\Gamma(\frac{n-1}{2})}{\sqrt{n\pi^{n-1}}}(g+1)^{-p/2}\left(s^2+\dfrac{(\tilde{\beta}_\gamma-\hat{\beta}_\gamma)^\text{T}X_\gamma^\text{T}X_\gamma(\tilde{\beta}_\gamma-\hat{\beta}_\gamma)}{g+1}\right)^{-\frac{n-1}{2}}$
对超参数 $g$ 的选择常默认选单位信息先验（UIP），即 $g=n$，这可以使得先验和样本观测值占相同比重。
后记
本文内容参考自中国科学技术大学《贝叶斯分析》（2022Fall）课件。
如对本文内容有任何疑问，请联系 watthu@mail.ustc.edu.cn。
