模型选择问题是要在给定的样本下，在一类候选模型中依照某个准则选择最优的模型。对模型的选择和评价依赖于抽样分布的结构、模型参数的先验分布指定等等。
贝叶斯因子
设总体 $X\sim f(x\mid\theta)$，其中 $\theta$ 为一未知参数且 $\theta\in\Theta$。比较两个模型
等价于假设检验问题 $H_0:\theta\in\Theta_0\leftrightarrow H_1:\theta\in\Theta_1$，其中 $\Theta_0=\Theta-\Theta_1$。
令 $g_i(\theta)$ 表示给定真实模型 $M_i$ 下 $\theta$ 的先验密度（$i=0,1$），则当有了样本 $\pmb x$ 后，可以使用贝叶斯因子来比较 $M_0$ 和 $M_1$：
其中 $\pi_0=\mathbb P^\pi(M_0)=\mathbb P^\pi(\Theta_0)$，$\mathbb P^\pi(M_1)=\mathbb P^\pi(\Theta_1)=1-\pi_0$，而
因此
如果
- 先验密度 $g_0,g_1$ 可以被指定，则可以仅仅使用贝叶斯因子 $\text{BF}_{01}(\pmb x)$ 进行模型选择；
- 进一步 $\pi_0$ 也可以被指定，则可以直接计算模型 $M_0$ 和 $M_1$ 的后验机会比进行模型选择，
但是大部分情况下，这些都不容易计算，甚至有时先验分布也难以指定时，其计算就更为复杂。
正常先验下的贝叶斯因子
假设候选模型有 $M_1,\dots,M_r$，在每个模型 $M_k$ 下抽样密度为 $f_k(x\mid\theta_k)$，其中 $\theta_k\in\Theta_k\subset\mathbb R^p$ 为未知的 $p$ 维参数向量。记 $\pi_k(\theta_k)$ 为模型 $M_k$ 下参数 $\theta_k$ 的先验密度，则在样本 $\pmb x$ 下，模型 $M_k$ 的后验密度为
其中 $f_k(\pmb x\mid\theta_k)$ 为样本 $\pmb x$ 在模型 $M_k$ 下的似然函数，$\mathbb P(M_k)$ 为模型 $M_k$ 的先验概率。
按照后验模型概率最大原则，最佳模型选择为
其中积分部分 $\displaystyle\mathbb P(\pmb x\mid M_k)=\int_{\Theta_k}f_k(\pmb x\mid\theta_k)\pi_k(\theta_k)\text{d}\theta_k$ 是样本 $\pmb x$ 在模型 $M_k$ 下的边际似然函数，表示指定先验分布对样本的拟合程度。
- 在对模型无任何先验信息参考时，常常认为
$\mathbb P(M_k)=\dfrac{1}{r},$ 此时后验密度与边际似然成正比；
- 有时希望对简单的模型赋予更多的先验概率，可取泊松先验
$\mathbb P(M_k)\propto\lambda^{p_k}\text{e}^{-\lambda},$ 其中 $p_k$ 为模型 $M_k$ 的复杂度（其值越大表示越复杂），而 $\lambda$ 衡量了模型中期望的复杂度。
- 对线性回归模型 $y=\pmb x^\textsf{T}\pmb\beta+\varepsilon$，上面的 $p_k$ 可理解为模型 $M_k$ 下的解释变量个数，另一种先验分布为
$\displaystyle\mathbb P(M_k)\propto\prod_{j=1}^p\pi_j^{r_j}(1-\pi_j)^{1-r_j},$ 其中 $\pi_j$ 表示解释变量 $x_j$ 被包含在模型 $M_k$ 中的概率，$r_j=1$ 当且仅当 $x_j$ 被包含进模型 $M_k$。
类似地，可以计算模型 $M_k$ 和 $M_j$ 的贝叶斯因子
从而模型 $M_k$ 的后验概率为
因此可以根据模型的后验概率或所有候选模型进行两两比较贝叶斯因子，从中选出最优模型。下表给出了Jeffreys对贝叶斯因子作为证据的程度的解释。
| 贝叶斯因子 | 解释 |
|---|---|
| $\text{BF}_{kj}<1$ | 否定模型 $M_k$ |
| $1<\text{BF}_{kj}<3$ | 对模型 $M_k$ 的支持证据微乎其微 |
| $3<\text{BF}_{kj}<10$ | 较强的证据支持模型 $M_k$ |
| $10<\text{BF}_{kj}<30$ | 强烈的证据支持模型 $M_k$ |
| $30<\text{BF}_{kj}<100$ | 非常强烈的证据支持模型 $M_k$ |
| $\text{BF}_{kj}>100$ | 肯定支持模型 $M_k$ |
非正常先验下的贝叶斯因子
由非正常先验的定义 $\pi(\theta)\propto h(\theta)$，其中 $\displaystyle\int_\Theta h(\theta)\text{d}\theta=\infty$，因此对任意正常数 $C$，可以使用 $q(\theta)=C\pi(\theta)$ 来作为先验，此时使用先验 $q_k(\theta_k)$ 和 $q_j(\theta_j)$ 的贝叶斯因子为
由于 $C_k,C_j$ 可以任意选取，因此贝叶斯因子不唯一。
潜在贝叶斯因子
潜在贝叶斯因子的想法是考虑使用部分样本作为训练样本估计先验分布（实际上是基于部分样本的后验分布），再使用剩余的样本来进行模型比较。
记样本 $\pmb X_n=(X_1,\dots,X_n)$，$\pmb X_{n(l)}$ 为 $\pmb X_n$ 的子集，再记 $\pmb X_{-n(l)}$ 为 $\pmb X_n$ 除去 $\pmb X_{n(l)}$ 外剩余的样本，则在先验 $\pi$ 下，
其中子集 $\pmb X_{n(l)}$ 的选择要使得此后验分布为正常分布。
以此上述分布 $\pi(\theta\mid\pmb X_{n(l)})$ 作为先验分布，结合剩余的样本 $\pmb X_{-n(l)}$，计算的贝叶斯因子称为在给定样本 $\pmb X_{n(l)}$ 时的潜在贝叶斯因子（IBF）或部分贝叶斯因子，其表达式为
此时该量是唯一的，不确定项 $C_k/C_j$ 在上式中被消去，但是依赖于部分样本 $\pmb X_{n(l)}$ 的选择。
- 算数潜在贝叶斯因子：将样本 $\pmb X_n$ 分为 $N$ 个最小子集 $\pmb X_{n(l)},l=1,\dots,N$，这些子集是使得 $\pi(\theta\mid\pmb X_{n(l)})$ 是正常密度的最小子集，定义算数潜在贝叶斯因子（AIBF）为
$\text{AIBF}_{kj}=\displaystyle\frac{1}{N}\sum_{l=1}^N\text{IBF}_{kj}(\pmb X_{n(l)}).y$ 它是 $N$ 个潜在贝叶斯因子的算术平均。
- 几何潜在贝叶斯因子：考虑这些潜在贝叶斯因子的几何平均即得几何潜在贝叶斯因子（GIBF）为
$\text{GIBF}_{kj}=\displaystyle\left[\prod_{l=1}^N\text{IBF}_{kj}(\pmb X_{n(l)})\right]^{1/N}.$
分数贝叶斯因子
为了避免潜在贝叶斯因子的数据子集选择问题，O’Hagan从另一角度提出了分数贝叶斯因子。
记观测样本为 $\pmb x_n=(\pmb z_{n-m},\pmb y_m)$，其中 $\pmb y$ 作为训练样本，则潜在贝叶斯因子为
注意到
而当 $n$ 和 $m$ 都很大时，在简单样本情况下似然 $f(\pmb y\mid\theta)$ 趋于 $[f(\pmb x\mid\theta)]^b$，其中 $b=m/n$，因此定义分数贝叶斯因子（FBF）为
其中 $q_i(\pmb x,b)=\dfrac{\int_{\Theta_i}\pi_i(\theta_i)f_i(\pmb x\mid\theta_i)\text{d}\theta_i}{\int_{\Theta_i}\pi_i(\theta_i)[f_i(\pmb x\mid\theta_i)]^b\text{d}\theta_i},~i=k,j$。可以看出
它也不受不确定项 $C_k/C_j$ 的影响。
后验贝叶斯因子
Aitkin提出后验贝叶斯因子（PBF）来克服不正常先验情形下贝叶斯因子定义的缺点，其定义为
其中 $\pi_i(\theta_i\mid\pmb x)$ 为 $\theta_i$ 在样本 $\pmb x$ 下的后验密度。显然，后验贝叶斯因子为不同模型下似然函数的后验均值之比，其使用与贝叶斯因子类似：$\text{PBF}_{kj}$ 的值小于 $1/20$、$1/100$ 和 $1/1000$ 分别表示有强、非常强、极强的证据来否定模型 $M_k$ 而支持模型 $M_j$。
基于交叉验证的拟贝叶斯因子
Gelfand等提出使用交叉验证预测密度（CVPD）
其中 $\pmb x_{-i}$ 表示除去 $x_i$ 后剩余的所有样本，从而定义拟贝叶斯因子（PSBF）为
显然 $\displaystyle f(x_i\mid\pmb x_{-i})=\int_\Theta f(x_i\mid\theta)\pi(\theta\mid\pmb x_{-i})\text{d}\theta$ 为预测密度。交叉验证方法的优点是适用于各种实际情形，而且当样本量较大时计算时间也很可观。
贝叶斯因子的Laplace近似
记 $\tilde{\theta}_i$ 为 $\theta_i$ 的后验众数，并且假定 $\tilde{\theta}_i$ 为 $\Theta_i$ 的内点，则
其中 $H_{\tilde{\theta}_i}$ 为相应的Hessian矩阵在 $\tilde{\theta}_i$ 处的值。根据Laplace近似公式知
这里 $p_i$ 是参数 $\theta_i$ 的维数。于是
（注：实践中常用 $2\ln\text{BF}_{kj}(\pmb x)$ 来表示样本 $\pmb x$ 对模型 $M_k$ 相对于 $M_j$ 的支持程度。）
类似地，可以用极大似然估计 $\hat{\theta}_i$ 代替 $\tilde{\theta}_i$，则得到
这里 $H_{\hat{\theta}_i}^{-1}$ 为观测的Fisher信息阵在 $\hat{\theta}_i$ 处的值。如果样本是独立同分布的，则
这里 $H_{1,\hat{\theta}_i}^{-1}$ 为期望的Fisher信息阵在 $\hat{\theta}_i$ 处的值。于是
因此
注意到Fisher信息阵对样本量 $n$ 的有界性，如果先验密度满足 $\ln g_i=O(1)$，则上式的一个近似为
这就是使用BIC逼近贝叶斯因子，其中 $(p_k-p_j)\ln n$ 可以视为模型复杂度的一个惩罚。
贝叶斯因子的模拟计算
重要性抽样方法
注意到在模型 $M_k$ 下，样本的边际密度
记被积函数（没有正则化的后验密度）为 $\tilde{f}_k(\theta_k\mid\pmb x)$，可以考虑对其作重要性抽样
其中 $q_k(\cdot)$ 为重要性抽样密度。由此可以计算贝叶斯因子
其中 $w_r(M_i)=\dfrac{\tilde{f}_i(\theta_i^{(r)}\mid\pmb x)}{q_i(\theta_i^{(r)})}$，$\theta_i^{(r)}\sim q_i$，$i=k,j$。
这种方法的计算精度显然依赖于重要性抽样密度的选取，尤其是参数为高维时选取就尤为困难。
MCMC方法
一般的MCMC方法将感兴趣的量表示为某个量的后验期望，从而抽取后验样本对其进行估计。Gelfand和Dey提供了一个基本的表达
其中 $q_k(\cdot)$ 为一个密度。从而贝叶斯因子可以表示为
利用MCMC抽取的后验样本，我们可以估计每个边际密度 $m_i(\pmb x)$：
其中 $\theta_i^{(r)}$ 为从后验密度 $f_i(\theta_i\mid\pmb x)$ 中抽取的样本。
可逆跳转MCMC
在一些问题下，MCMC过程应该能够在两个不同维数的参数空间之间移动。可逆跳转MCMC就实现了这一想法。假定我们有模型 $k~(k\in\mathcal K)$，$\cal K$ 为一可数集，且模型 $k$ 有连续的参数空间 $\Theta_k\in\mathbb R^{n_k}$，不同模型的参数维数可能是不相同的。假定模型的先验分布为
显然 $\sum\limits_{k\in\mathcal K}p_k=1$。对每个模型 $k$，参数 $\pmb\theta_k$ 的先验分布为 $\pi(\pmb\theta\mid k)$。记样本为 $\pmb y$，假设 $f(\pmb y\mid k,\pmb\theta_k)$ 为在模型 $k$ 下的似然函数，则 $(k,\pmb\theta_k)$ 的后验分布为
于是贝叶斯因子为
其中 $\displaystyle\mathbb P^\pi(k\mid\pmb y)=\int_{\Theta_k}\pi(k,\pmb\theta_k\mid\pmb y)\text{d}\pmb\theta_k$ 为模型 $k$ 的后验概率。
一个常规的想法是单独计算每个 $\mathbb P^\pi(k\mid\pmb y)$，然后选择模型后验概率最大的模型，或者使用后验概率 $\mathbb P^\pi(k\mid\pmb y)$ 为权进行模型平均。这种方法对存在大量的可选模型时是不可行的，且效率低下。
考虑可逆跳转MCMC方法。记 $\pmb x=(k,\pmb\theta)$，则 $\pmb x$ 的取值空间为
其分布为前述后验分布形式。按照M-H算法的要求，我们需要使用一个提议分布 $q(\cdot\mid\pmb x)$，在当前状态 $\pmb x=(k,\pmb\theta_k)$ 下，产生一个候选状态 $\pmb x’=(m,\pmb\theta_m)$，并以一定的接受概率接受其为下一步状态。由于 $q(m,\pmb\theta_m\mid\pmb x)=q(m\mid\pmb x)q(\pmb\theta_m\mid\pmb x,m)$，特别地 $q(m\mid\pmb x)=q(m\mid k)$，$q(\pmb\theta_m\mid\pmb x,m)=q(\pmb\theta_m\mid\pmb\theta_k)$，于是接受概率
但是由于 $\pmb\theta_m$ 和 $\pmb\theta_k$ 维数的差异，我们需要选取辅助变量 $\pmb u$ 和 $\pmb v$ 使得 $(\pmb\theta_k,\pmb u)$ 和 $(\pmb\theta_m,\pmb v)$ 维数一致。因此，假设我们从某个提议密度 $g(\pmb u\mid\pmb\theta_k,m)$ 中产生一个 $\pmb u$，然后通过一一映射 $\phi$ 建立
于是
其中 $\text{det}~\pmb J_{\phi}(\pmb\theta_k,\pmb u)$ 为一一映射 $\phi$ 在 $(\pmb\theta_i,\pmb u)$ 处的Jacobi行列式的值。由此可以简化接受概率的计算。
综上所述，可逆跳转MCMC算法为
迭代 $t$：若 $\pmb x_t=(k,\pmb\theta_k)$，则
- 以概率 $q(m\mid k)$ 选择模型 $m$；
- 从提议分布 $g(\pmb u\mid\pmb\theta_k,m)$ 中产生 $u$，令 $(\pmb\theta_m,\pmb v)=\phi(\pmb\theta_k,pmb u)$；
- 以概率
$\alpha(\pmb x_t,\pmb x_{t+1})=\min\left\{1,\dfrac{\pi(m,\pmb\theta_m\mid\pmb y)q(k\mid m)g(\pmb v\mid\pmb\theta_m)}{\pi(k,\pmb\theta_k\mid\pmb y)q(m\mid k)g(\pmb u\mid\pmb\theta_k)}\left\vert\text{det}~\pmb J_{\phi}(\pmb\theta_k,\pmb u)\right\vert\right\}$ 接受 $\pmb x_{t+1}=(m,\pmb\theta_m)$，否则 $\pmb x_{t+1}=\pmb x_t$。
例如，在对计数数据进行建模时，经常感兴趣的一个问题是数据对泊松分布来说是否过度分散。如果过度分散，则使用负二项分布来拟合可能更合适。对于给定的数据 $\pmb y$，使用可逆跳转MCMC方法进行模型选择。
给定容量为 $n$ 的i.i.d.样本 $\pmb y$ 时，在参数为 $\lambda>0$ 的泊松分布下似然函数为
而在参数 $\lambda>0$，$\kappa>0$ 的负二项分布下似然函数为
泊松分布和负二项分布的均值都为 $\lambda$，而负二项分布的方差为 $\lambda(1+\kappa\lambda)$，因而更适合过度分散的数据。
假设两个模型指示变量的先验分布为 $\mathbb P(k=1)=\mathbb P(k=2)=0.5$，参数 $\theta_1=\lambda$ 和 $\theta_2=(\theta_{21},\theta_{22})=(\lambda,\kappa)$ 的先验取为 $\theta_1,\theta_{21}\sim\Gamma(\alpha_\lambda,\beta_\lambda)$，而 $\theta_{22}\sim\Gamma(\alpha_\kappa,\beta_\kappa)$。于是后验分布为
其中 $p(\theta_1\mid k=1)=\gamma(\theta_1;\alpha_\lambda,\beta_\lambda)$，$p(\theta_{21},\theta_{22}\mid k=2)=\gamma(\theta_{21};\alpha_\lambda,\beta_\lambda)\gamma(\theta_{22};\alpha_\kappa,\beta_\kappa)$，这里 $\gamma(\cdot;\alpha,\beta)$ 为 $\Gamma(\alpha,\beta)$ 的密度函数。
下面构造从模型1到模型2的合适转移。记 $\pmb x=(1,\lambda)$ 为当前链的状态，由于模型1中没有和参数 $\kappa$ 等价的分量，考虑引入辅助变量 $u\sim N(0,\sigma^2)$（$\sigma^2$ 已知），然后令 $\pmb x’=(2,\theta)$，其中 $\theta=(\lambda,\kappa)=\phi(\lambda,u)=(\lambda,\mu\text{e}^u)$，其中 $\mu$ 是固定的。换言之，$\lambda$ 在变换中不变，而 $\kappa$ 是一个对数正态随机变量。因而变换的Jacobi行列式的值为
从而从模型1到模型2的接受概率为
再看模型2到模型1的转移。令 $(\lambda,u)=\phi’(\lambda,\kappa)=(\lambda,\ln(\kappa/\mu))$，因此模型2转移到模型1的接受概率为
后记
本文内容参考自中国科学技术大学《贝叶斯分析》（2022Fall）课件。
如对本文内容有任何疑问，请联系 watthu@mail.ustc.edu.cn。
