从本章开始，关于贝叶斯的内容大多集中在贝叶斯推断的具体实施，这往往需要依赖于一些统计软件进行数值计算，但它们也都有良好的理论基础。
引入：贝叶斯推断例子
在3.1节贝叶斯统计推断（点估计与区间估计）中我们提到一个柯西分布总体在无信息先验下的最大后验HPD可信集的计算，如果我们关心参数的后验均值和后验方差，显然也没有显式表达。
又如，假设 $X_1,\dots,X_k$ 为独立的Poisson随机变量，且 $X_i\sim\mathcal P(\theta_i),i=1,\dots,k$。如果 $\theta_i$ 的先验分布为通过其对数的联合分布给出：
$\pmb\nu=(\ln\theta_1,\dots,\ln\theta_k)^\textsf{T}\sim N_k(\mu\pmb 1_k,\tau^2[(1-\rho)\pmb I_k+\rho\pmb J_k]),$
其中 $\pmb 1_k$ 为元素是1的 $k$ 维列向量，$\pmb I_k$ 为 $k$ 阶单位阵，$\pmb J_k$ 是元素全为1的 $k$ 阶方阵，$\mu,\tau^2,\rho$ 均为已知的常数。由
$\begin{array}{rl}
f(\pmb x\mid\pmb\nu)=&\!\!\!\!\displaystyle\exp\left\{-\sum_{i=1}^k(\text{e}^{\nu_i}-\nu_ix_i)\right\}\Big/\prod_{i=1}^kx_i!,\\
\pi(\pmb\nu)\propto&\!\!\!\!\displaystyle\exp\left\{-\dfrac{1}{2\tau^2}(\pmb\nu-\mu\pmb 1_k)^\textsf{T}[(1-\rho)\pmb I_k+\rho\pmb J_k]^{-1}(\pmb\nu-\mu\pmb 1_k)\right\},
\end{array}$
可得
$\begin{array}{rl}
\pi(\pmb\nu\mid\pmb x)\propto&\!\!\!\! g(\pmb\nu\mid\pmb x)\\
=&\!\!\!\!\displaystyle\exp\left\{-\sum_{i=1}^k(\text{e}^{\nu_i}-\nu_ix_i)-\dfrac{1}{2\tau^2}(\pmb\nu-\mu\pmb 1_k)^\textsf{T}[(1-\rho)\pmb I_k+\rho\pmb J_k]^{-1}(\pmb\nu-\mu\pmb 1_k)\right\}.
\end{array}$
于是后验均值为
${\mathcal E}^\pi(\theta_j\mid\pmb x)=\dfrac{\int_{\mathbb R^k}\text{e}^{\nu_j}g(\pmb\nu\mid\pmb x)\text{d}\pmb\nu}{\int_{\mathbb R^k}g(\pmb\nu\mid\pmb x)\text{d}\pmb\nu}.$
这是两个 $k$ 重积分的比值，当 $k$ 越大时就越难处理，而数值积分方法在这种场合下不再是有效的方法。近年来，又有许多新的方法提出，例如分析逼近方法、E-M方法、MCMC抽样方法（Gibbs抽样方法、M-H算法）等等。
分析逼近方法
在计算一些积分时，可以使用Laplace逼近方法来逼近所要计算的积分。在贝叶斯框架下，设 $\pmb\theta$ 的后验密度为 $\pi(\pmb\theta\mid\pmb x)$，我们关心函数 $g(\pmb\theta)$ 的后验期望，即计算积分
$\displaystyle{\mathbb E}^\pi[g(\pmb\theta)\mid\pmb x]=\int_{\mathbb R^k}g(\pmb\theta)\pi(\pmb\theta\mid\pmb x)\text{d}\pmb\theta=\dfrac{\int_{\mathbb R^k}g(\pmb\theta)f(\pmb x\mid\pmb\theta)\pi(\pmb\theta)\text{d}\pmb\theta}{\int_{\mathbb R^k}f(\pmb x\mid\pmb\theta)\pi(\pmb\theta)\text{d}\pmb\theta},$
其中 $g,f,\pi$ 均为 $\pmb\theta$ 的光滑函数。
参数维数 $k=1$ 的情形
当参数维数 $k=1$ 时，将分子分母中的积分统一表示为
$\displaystyle I=\int_{-\infty}^{\infty}q(\theta)\text{e}^{-nh(\theta)}\text{d}\theta,$
其中 $-nh(\theta)=\ln f(\pmb x\mid\theta)\pi(\theta)$。假定 $q$ 和 $h$ 为 $\theta$ 的光滑函数，且 $-h$ 有唯一极大值点 $\hat{\theta}$。由于 $\hat{\theta}$ 是 $\theta$ 的MLE或后验众数，因此积分的主要贡献来自 $\hat{\theta}$ 的一个邻域 $(\hat{\theta}-\delta,\hat{\theta}+\delta)$，且当 $n\rightarrow\infty$ 时，有
$\displaystyle I\sim I_1=\int_{\hat{\theta}-\delta}^{\hat{\theta}+\delta}q(\theta)\text{e}^{-nh(\theta)}\text{d}\theta,$
这里 $I\sim I_1$ 表示 $I/I_1\rightarrow1(n\rightarrow\infty)$。
Laplace逼近方法涉及 $q$ 和 $h$ 在 $\hat{\theta}$ 处的Taylor展开，有
$\begin{array}{rl}
q(\theta)=&\!\!\!\!q(\hat{\theta})+(\theta-\hat{\theta})q'(\hat{\theta})+\dfrac{1}{2}(\theta-\hat{\theta})^2q''(\hat{\theta})+\text{余项},\\
h(\theta)=&\!\!\!\!h(\hat{\theta})+\dfrac{1}{2}(\theta-\hat{\theta})h''(\hat{\theta})+\text{余项},
\end{array}$
所以
$q(\theta)\text{e}^{-nh(\theta)}\approx\text{e}^{-nh(\hat{\theta})}q(\hat{\theta})\left[1+(\theta-\hat{\theta})\dfrac{q'(\hat{\theta})}{q(\hat{\theta})}+\dfrac{1}{2}(\theta-\hat{\theta})^2\dfrac{q''(\hat{\theta})}{q(\hat{\theta})}\right]\exp\left\{-\dfrac{n}{2}(\theta-\hat{\theta})^2h''(\hat{\theta})\right\},$
令 $c=h’’(\hat{\theta})$，$t=\sqrt{nc}(\theta-\hat{\theta})$，则 $\text{d}\theta=\text{d}t/\sqrt{nc}$。利用
$\displaystyle\int_{-\sqrt{nc}\delta}^{\sqrt{nc}\delta}t\text{e}^{-\frac{t^2}{2}}=0,\quad\int_{-\sqrt{nc}\delta}^{\sqrt{nc}\delta}\dfrac{1}{\sqrt{2\pi}}t^2\text{e}^{-\frac{t^2}{2}}\approx1$（当 $n$ 充分大时），
则有
$\begin{array}{rl}
I_1\approx&\!\!\!\!\displaystyle\text{e}^{-nh(\hat{\theta})}q(\hat{\theta})\dfrac{\sqrt{2\pi}}{\sqrt{nc}}\int_{-\sqrt{nc}\delta}^{\sqrt{nc}\delta}\left[1+\dfrac{1}{\sqrt{nc}}\dfrac{q'(\hat{\theta})}{q(\hat{\theta})}+\dfrac{t^2}{2nc}\dfrac{q''(\hat{\theta})}{q(\hat{\theta})}\right]\dfrac{1}{\sqrt{2\pi}}\text{e}^{-\frac{t^2}{2}}\text{d}t\\
=&\!\!\!\!\text{e}^{-nh(\hat{\theta})}q(\hat{\theta})\dfrac{\sqrt{2\pi}}{\sqrt{nc}}\left[1+\dfrac{q''(\hat{\theta})}{2ncq(\hat{\theta})}\right]\\
=&\!\!\!\!\text{e}^{-nh(\hat{\theta})}q(\hat{\theta})\dfrac{\sqrt{2\pi}}{\sqrt{nc}}\left[1+O(n^{-1})\right].
\end{array}$
这就是积分 $I$ 的近似 $I\sim\text{e}^{-nh(\hat{\theta})}q(\hat{\theta})\dfrac{\sqrt{2\pi}}{\sqrt{nc}}\left[1+O(n^{-1})\right]$。
例如，利用Laplace逼近方法求 $n!$ 的Stirling逼近如下：
$\displaystyle n!=\Gamma(n+1)=\int_0^{\infty}x^n\text{e}^{-x}\text{d}x=\int_0^{\infty}\text{e}^{-n(x/n-\ln x)}\text{d}x,$
所以 $q(x)=1$，$h(x)=\dfrac{x}{n}-\ln x$，故 $\hat{x}=n$，$c=h’’(\hat{x})=\dfrac{1}{n^2}$，所以
$n!\approx\text{e}^{-n(1-\ln n)}\cdot\sqrt{2\pi n}\left[1+O(n^{-1})\right]\approx\sqrt{2\pi}n^{n+1/2}\text{e}^{-n},$
即 $n!\sim\sqrt{2\pi}n^{n+1/2}\text{e}^{-n}$。
参数维数 $k>1$ 的情形
此时 $\pmb\theta=(\theta_1,\dots,\theta_k)^\textsf{T}$ 为 $k\times1$ 向量，故积分变为如下形式的积分
$\displaystyle I=\int_{\mathbb R^k}q(\pmb\theta)\exp\{-nh(\pmb\theta)\}\text{d}\pmb\theta,$
其中 $h$ 为一光滑函数且 $-h$ 有唯一的最大值点 $\hat{\pmb\theta}$。记 $\pmb\Delta_q$ 和 $\pmb\Delta_h$ 分别为函数 $q$ 和 $h$ 的Hessian矩阵，则
$\begin{array}{rl}
q(\pmb\theta)=&\!\!\!\!q(\hat{\pmb\theta})+(\pmb\theta-\hat{\pmb\theta})^\textsf{T}q'(\hat{\pmb\theta})+\dfrac{1}{2}(\pmb\theta-\hat{\pmb\theta})^\textsf{T}\pmb\Delta_q(\hat{\pmb\theta})(\pmb\theta-\hat{\pmb\theta})+\text{余项},\\
h(\pmb\theta)=&\!\!\!\!h(\hat{\pmb\theta})+\dfrac{1}{2}(\pmb\theta-\hat{\pmb\theta})^\textsf{T}\pmb\Delta_h(\hat{\pmb\theta})(\pmb\theta-\hat{\pmb\theta})+\text{余项},
\end{array}$
从而类似地，有
$I\sim\text{e}^{-nh(\hat{\pmb\theta})}q(\hat{\pmb\theta})(2\pi)^{k/2}n^{-k/2}\vert\pmb\Delta_h(\hat{\pmb\theta})\vert^{-1/2}\left[1+O(n^{-1})\right].$
将这种方法应用到后验期望的计算中就有
$\mathbb E^\pi[g(\pmb\theta)\mid\pmb x]\sim g(\hat{\pmb\theta})\left[1+O(n^{-1})\right].$
如果在Taylor展开中有更多的项保留在积分逼近中，则二阶逼近也是容易得到的。
E-M方法
E-M算法是一种最优化方法，常用于不完全数据求极大似然估计的问题中。
假设 $\pmb Y\mid\pmb\theta$ 有密度 $f(\pmb y\mid\pmb\theta)$，且 $\pmb\theta$ 的先验密度为 $\pi(\pmb\theta)$，由此得到的后验分布记为 $\pi(\pmb\theta\mid\pmb y)$。大部分情况下，计算 $\pi(\pmb\theta\mid\pmb y)$ 的数字特征非常困难，可以用一些“数据扩张”（data augmentation）的方法来解决。数据扩张的基本思想是将观测到的数据 $\pmb y$ 和缺失数据或隐变量数据 $\pmb z$ 扩张为“完全”数据 $\pmb x=(\pmb y,\pmb z)$，使得扩张后的后验分布 $\pi(\pmb\theta\mid\pmb x)=\pi(\pmb\theta\mid\pmb y,\pmb z)$ 在计算上是容易处理的。
E-M算法是这类数据扩张方法中最简单的一种方法，它的目的是最大化后验分布，也就是只能计算后验众数（MAP）。不过，如果使用数据扩张方法能够得到一个易于计算的后验分布，则可以使用更有效的计算工具来计算后验分布的各种数字特征。
E-M算法
记 $p(\pmb z\mid\pmb y,\hat{\pmb\theta})$ 为在给定 $\pmb y$ 和 $\hat{\pmb\theta}$ 时 $\pmb Z$ 的预测分布，$\hat{\pmb\theta}^{(i)}$ 为在第 $i$ 次迭代时 $\pmb\theta$ 的估计值。E-M算法的基本步骤如下：
- E步：计算 $Q(\pmb\theta\mid\hat{\pmb\theta}^{(i)})=\mathbb E[\ln\pi(\pmb\theta\mid\pmb y,\pmb z)\mid\pmb y,\hat{\pmb\theta}^{(i)}]$；
- M步：对 $\pmb\theta$ 最大化 $Q(\pmb\theta\mid\hat{\pmb\theta}^{(i)})$ 以获得 $\hat{\pmb\theta}^{(i+1)}$。
这种先求期望，然后求最大值的步骤即为E-M算法名称的由来。
下面这个过程是随机化E-M算法，它通过从 $p(\pmb z\mid\pmb y,\hat{\pmb\theta})$ 中随机抽样来生成缺失变量 $\pmb z$ 的值，然后带入到完全后验分布中：
- 计算 $\pmb z^{(i)}=\mathbb E[\pmb Z\mid\pmb y,\hat{\pmb\theta}^{(i)}]$；
- 将观测数据 $\pmb y$ 扩张为 $(\pmb y,\pmb z^{(i)})$，最大化 $\pi(\pmb\theta\mid\pmb y,\pmb z^{(i)})$，记其最大值点为 $\hat{\pmb\theta}^{(i+1)}$；
- 重复上述两步直至达到收敛要求。
一般来说，从任何初始值出发，由E-M算法一般都能达到一个局部最大值。
例1. 考虑Rao（1973）关于特定基因组合率的数据如下：
| 计数 |
$y_1=125$ |
$y_2=18$ |
$y_3=20$ |
$y_4=34$ |
| 概率 |
$\dfrac{1}{2}+\dfrac{\theta}{4}$ |
$\dfrac{1}{4}(1-\theta)$ |
$\dfrac{1}{4}(1-\theta)$ |
$\dfrac{\theta}{4}$ |
在均匀先验 $U(0,1)$ 下，$\theta$ 的后验分布为
$\pi(\theta\mid\pmb y)\propto(2+\theta)^{y_1}(1-\theta)^{y_2+y_3}\theta^{y_4}.$
由于 $2+\theta$ 项的存在，此分布不是一个标准分布。考虑把第一类拆成两类，其概率分别为 $1/2$ 和 $\theta/4$，则完全数据为 $\pmb x=(x_1,x_2,x_3,x_4,x_5)$，其中 $x_1+x_2=y_1$，$x_k=y_{k-1}~(k=3,4,5)$。因此扩张的后验分布为
$\pi(\theta\mid\pmb x)\propto\theta^{x_2+x_5}(1-\theta)^{x_3+x_4},$
这是一个贝塔分布。
E-M算法中的E步为
$\begin{array}{rl}
Q(\theta\mid\hat{\theta}^{(i)})=&\!\!\!\!\mathbb E\left[(x_2+x_5)\ln\theta+(x_3+x_4)\ln(1-\theta)\mid\pmb y,\hat{\theta}^{(i)}\right]\\
=&\!\!\!\!\left[\mathbb E(x_2\mid y_1,\hat{\theta}^{(i)})\right]\ln\theta+(y_2+y_3)\ln(1-\theta).
\end{array}$
M步为最大化 $Q$，由
$\dfrac{\partial Q(\theta\mid\hat{\theta}^{(i)})}{\partial\theta}=0$ 可得 $\hat{\theta}^{(i+1)}=\dfrac{\mathbb E(x_2\mid y_1,\hat{\theta}^{(i)})+y_4}{\mathbb E(x_2\mid y_1,\hat{\theta}^{(i)})+y_2+y_3+y_4}.$
注意到
$x_2\mid[(x_1+x_2),\hat{\theta}^{(i)}]\sim B\left(x_1+x_2,\dfrac{\hat{\theta}^{(i)}/4}{1/2+\hat{\theta}^{(i)}/4}\right),$
所以
$\mathbb E(x_2\mid y_1,\hat{\theta}^{(i)})=\mathbb E(x_2\mid x_1+x_2=y_1,\hat{\theta}^{(i)})=y_1\cdot\dfrac{\hat{\theta}^{(i)}}{2+\hat{\theta}^{(i)}}.$
由此得到迭代公式
$\hat{\theta}^{(i+1)}=\dfrac{y_1\cdot\hat{\theta}^{(i)}+y_4(2+\hat{\theta}^{(i)})}{y_1\cdot\hat{\theta}^{(i)}+(y_2+y_3+y_4)(2+\hat{\theta}^{(i)})}.$
利用如下代码实现E-M算法的迭代过程。
EM <- function(y, max.it = 10000, eps = 1e-5) {
i <- 1
theta1 <- 1
theta2 <- 0.5
while (abs(theta1 - theta2) >= eps) {
theta1 <- theta2
x2 <- y[1] * theta1 / (2 + theta1)
theta2 <- (x2 + y[4]) / (x2 + y[2] + y[3] + y[4])
if (i == max.it) { break }
i <- i + 1
}
return(theta2)
}
y <- c(125, 18, 20, 34)
EM(y, max.it = 10000, eps = 1e-5)
代码输出0.6268207，这就是所得的估计 $\hat{\theta}=0.62682$。
蒙特卡洛抽样方法
蒙特卡洛抽样方法是一种概率化的技巧，它的思想是利用样本均值或样本分位数来估计相应的总体特征，而大数定律则保证了所得估计量的相合性。
具体来说，假设 $f$ 为一概率函数，我们感兴趣的量为如下形式的有限期望
$\displaystyle\mathbb E_f[h(X)]=\int_{\mathscr{X}}h(x)f(x)\text{d}x.$
如果可以从 $f$ 中产生i.i.d.观测值 $X_1,\dots,X_m$，则
$\displaystyle\overline{h}_m=\dfrac{1}{m}\sum_{i=1}^mh(X_i)$
依概率收敛（甚至几乎处处收敛）到感兴趣的量 $\mathbb E_f[h(X)]$。这种估计的误差可以用样本标准差来刻画：
$\displaystyle s_m=\sqrt{\frac{1}{m-1}\sum_{i=1}^m[h(X_i)-\overline{h}_m]^2}.$
进一步，还可以得到 $\mathbb E_f[h(X)]$ 的渐近水平为 $1-\alpha$ 的置信区间
$\left[\overline{h}_m-u_{\alpha/2}s_m/\sqrt{m},\overline{h}_m+u_{\alpha/2}s_m/\sqrt{m}\right].$
在贝叶斯框架下，我们往往关心参数的后验期望，则可以通过从后验分布中产生i.i.d.观测值，然后计算相应的样本均值作为估计。但是，大多数情况下后验分布不是标准分布，很难从中直接抽样。为此，上述抽样方法的一个变种——蒙特卡洛重要性抽样方法——被提出。
假设从 $f$ 中直接抽样很困难，而从与 $f$ 很靠近的一个分布 $g$ 中抽样比较容易。那么
$\displaystyle\mathbb E_f[h(X)]=\int_{\mathscr{X}}h(x)f(x)\text{d}x=\int_{\mathscr{X}}h(x)\dfrac{f(x)}{g(x)}g(x)\text{d}x=\mathbb E_g[h(X)w(X)],$
其中 $w(x)=f(x)/g(x)$。从 $g$ 中产生i.i.d.样本 $X_1,\dots,X_m$，则一个合适的估计为
$\displaystyle\widehat{\mathbb E_f[h(X)]}=\dfrac{1}{m}\sum_{i=1}^nh(X_i)w(X_i).$
抽样分布 $g$ 称为重要性函数。
例如，假设 $\pmb X=(X_1,\dots,X_n)$ 为从 $N(\theta,\sigma^2)$ 中抽取的i.i.d.样本，其中 $\theta,\sigma^2$ 均未知。取 $\theta$ 和 $\sigma^2$ 的先验为独立先验，其中 $\theta$ 服从双指数先验分布，密度为 $\text{e}^{-\vert\theta\vert}/2$，$\sigma^2$ 有先验密度 $(1+\sigma^2)^{-2}$。这两个先验分布都不是标准先验，但都是从稳健性角度考虑而选取的。由
$\begin{array}{rl}
f(\pmb x\mid\theta,\sigma^2)=&\!\!\!\!\displaystyle(2\pi\sigma^2)^{-\frac{n}{2}}\exp\left\{-\frac{1}{2\sigma^2}\bigg[\sum_{i=1}^n(x_i-\overline{x})^2+n(\overline{x}-\theta)^2\bigg]\right\}\\
\pi(\theta,\sigma^2)=&\!\!\!\!\pi_1(\theta)\pi_2(\sigma^2)\propto\text{e}^{-\vert\theta\vert}(1+\sigma^2)^{-2},
\end{array}$
可知
$\begin{array}{rl}
\pi(\theta,\sigma^2\mid\pmb x)\propto&\!\!\!\!(\sigma^2)^{-\frac{n}{2}}\exp\left\{-\dfrac{n}{2\sigma^2}[(\theta-\overline{x})^2+s_n^2]\right\}\text{e}^{-\vert\theta\vert}(1+\sigma^2)^{-2}\\
=&\!\!\!\![(\theta-\overline{x})^2+s_n^2]^{\frac{n}{2}+1}(\sigma^2)^{-(\frac{n}{2}+1)-1}\exp\left\{-\dfrac{n}{2\sigma^2}[(\theta-\overline{x})^2+s_n^2]\right\}\\
&\times[(\theta-\overline{x})^2+s_n^2]^{-\frac{n+2}{2}}\text{e}^{-\vert\theta\vert}\left(\dfrac{\sigma^2}{1+\sigma^2}\right)^2\\
\propto&\!\!\!\!g_1(\sigma^2\mid\theta,\pmb x)g_2(\theta\mid\pmb x)\text{e}^{-\vert\theta\vert}\left(\dfrac{\sigma^2}{1+\sigma^2}\right)^2,
\end{array}$
其中 $g_1$ 为形状参数为 $n/2+1$、尺度参数为 $n[(\theta-\overline{x})^2+s_n^2]/2$ 的逆伽马密度函数，$g_2$ 为自由度为 $n+1$、位置参数为 $\overline{x}$、尺度参数为 $s_n/\sqrt{n+1}$ 的 $t$ 分布密度函数，$s_n^2=\dfrac{1}{n}\sum\limits_{i=1}^n(x_i-\overline{x})^2$。
因此，可以选择 $g(\theta,\sigma^2\mid\pmb x)=g_1(\sigma^2\mid\theta,\pmb x)g_2(\theta\mid\pmb x)$ 作为重要性函数：在第 $i$ 步抽样时需要先从 $g_2$ 抽取一个 $\theta_i$，然后在给定 $\theta_i$ 的条件下，从 $g_1$ 抽取一个 $\sigma_i^2$，合在一起组成第 $i$ 步抽样 $(\theta_i,\sigma_i^2)$。如果 $\theta$ 的后验期望为感兴趣的量，即
$\displaystyle\mathbb E^\pi(\theta\mid\pmb x)=\int_{-\infty}^{\infty}\int_0^{\infty}\theta\pi(\theta,\sigma^2\mid\pmb x)\text{d}\sigma^2\text{d}\theta.$
在抽取了 $m$ 对 $(\theta_i,\sigma_i^2)~(i=1,\dots,m)$ 后得到后验期望 $\mathbb E(\theta\mid\pmb x)$ 的一个估计
$\displaystyle\widehat{\mathbb E^\pi(\theta\mid\pmb x)}=\sum_{i=1}^m\theta_iw(\theta_i,\sigma_i^2\mid\pmb x)\bigg/\sum_{i=1}^mw(\theta_i,\sigma_i^2\mid\pmb x),$
其中 $w(\theta,\sigma^2\mid\pmb x)=\dfrac{f(\pmb x\mid\theta,\sigma^2)\pi(\theta,\sigma^2)}{g(\theta,\sigma^2\mid\pmb x)}\propto\text{e}^{-\vert\theta\vert}\left(\dfrac{\sigma^2}{1+\sigma^2}\right)^2$。
后记
本文内容参考自中国科学技术大学《贝叶斯分析》（2022Fall）课件。
如对本文内容有任何疑问，请联系 watthu@mail.ustc.edu.cn。
