贝叶斯模型比较

注意：本文主要来自于Joshua C.C.Chan 的notes on bayesian macroeconometrics，只为记录学习，请勿他用。

1 边际似然的计算

1.1 Savage-Dicky Density Ratio

使用这个此密度比，我们可以解决边际似然计算困难的问题。先验分布通常选择简单常用的分布，我们可以解析地计算$p\left(\mathbf{\omega }={\mathbf{\omega }}_0\mid M_U\right)$， 分母可以通过计算蒙特卡洛平均来估计：
$$\frac{1}{R}\sum _{r=1}^{R}p\left(\mathbf{\omega }={\mathbf{\omega }}_0\mid \mathbf{y},{\mathbf{\psi }}^{(r)},M_U\right)$$
其中，${\mathbf{\psi }}^{(1)},\dots ,{\mathbf{\psi }}^{(R)}$是从无约束模型$M_U$抽取的draws。

例如： 我们重新讨论一个带有MA(1)误差的线性回归模型。如果MA(1)模型的系数$\psi$为0，则MA(1)模型就退化为带有独立误差的标准的线性回归模型。因此，支持MA(1)模型的贝叶斯因子可以写为$p(\psi =0)/p(\psi =0\mid \mathbf{y})$。

1.2 Modified Harmonic Mean

估计边际似然的一种常用方法是Gelfand和Dey(1994)的修正调和平均。定义$f$是一个p.d.f，它的支撑包含在后验分布的支撑剂之中。修改的调和平均基于如下的等式：
$$\mathbb{E}\left[\frac{f(\mathbf{\theta })}{p(\mathbf{\theta })p(\mathbf{y}\mid \mathbf{\theta })}\mid \mathbf{y}\right]=\int \frac{f(\mathbf{\theta })}{p(\mathbf{\theta })p(\mathbf{y}\mid \mathbf{\theta })}\frac{p(\mathbf{\theta })p(\mathbf{y}\mid \mathbf{\theta })}{p(\mathbf{y})}\mathrm{d}\mathbf{\theta }=p(\mathbf{y}{)}^{-1}$$
这个期望是关于后验分布$p(\mathbf{\theta }\mid \mathbf{y})=p(\mathbf{\theta })p(\mathbf{y}\mid \mathbf{\theta })/p(\mathbf{y})$来求的。因此，我们可以使用如下的估计量来估计边际似然：
$${\widehat{p(\mathbf{y})}}_{\mathrm{G}\mathrm{D}}={\left[\frac{1}{R}\sum _{r=1}^R\frac{f\left({\mathbf{\theta }}^{(r)}\right)}{p\left({\mathbf{\theta }}^{(r)}\right)p\left(\mathbf{y}\mid {\mathbf{\theta }}^{(r)}\right)}\right]}^{-1}$$
其中，${\mathbf{\theta }}^{(1)},\dots ,{\mathbf{\theta }}^{(R)}$是后验draws。当R趋近于无穷时，边际似然$p(\mathbf{y})$依概率是相合的。但这个估计并不是无偏的。Geweke(1999)表明，如果函数f的尾部比后验分布的尾部薄，则修正调和均值具有有限方差。由于后验分布在一定的正则性条件下是渐近正态分布，Geweke(1999)建议使用尾截尾后验分布的正态近似。更准确地说， 定义$\hat{\mathbf{\theta }}$和${\mathbf{Q}}_{\mathbf{\theta }}$表示参数$\mathbf{\theta }$的后验均值和后验协方差矩阵。 考虑选择的$f$是截断的密度$\mathcal{N}\left(\hat{\theta },{\mathbf{Q}}_{\theta }\right)$，区域为：
$$\left\{\mathbf{\theta }\in {\mathbb{R}}^m:(\mathbf{\theta }-\hat{\mathbf{\theta }}{)}^{\prime }{\mathbf{Q}}_{\theta }^{-1}(\mathbf{\theta }-\hat{\mathbf{\theta }})<{\chi }_{\alpha ,m}^2\right\}$$
其中，$m$表示$\mathbf{\theta }$的维度，${\chi }_{\alpha ,m}^2$表示${\chi }_m^2$的$1-\alpha$分位数。

改进的谐波均值的主要优点是它易于实现-编程的工作量最小，只需要后验绘制。虽然这种估计器通常对低维问题(例如，小于12个参数)是准确的，但在高维问题中它可能在数值上不稳定。

1.3 Chib’s Method

另外一种计算边际似然的方法是Chib(1995)， 这个方法的基础是边际似然是后验分布的标准化常数。或者等价地，
$$p(\mathbf{y})=\frac{p(\mathbf{y}\mid \mathbf{\theta })p(\mathbf{\theta })}{p(\mathbf{\theta }\mid \mathbf{y})}$$
因此，$p(\mathbf{y})$的一个自然的估计如下表示（在对数尺度下）：
$$\log {\widehat{p(\mathbf{y})}}_{\text{Chib }}=\log p\left(\mathbf{y}\mid {\mathbf{\theta }}^{\ast }\right)+\log p\left({\mathbf{\theta }}^{\ast }\right)-\log p\widehat{\left({\mathbf{\theta }}^{\ast }\mid \mathbf{y}\right)}$$
其中，${\mathbf{\theta }}^{\ast }$是后验分布支撑集中的任意点。具体原因稍后解释，选择${\mathbf{\theta }}^{\ast }$是某个高密度点如后验均值或者后验众数是有用的。

通常情形下，我们都可以解析地计算似然函数和先验分布的值，唯一未知量为后验分布$p\left({\mathbf{\theta }}^{\ast }\mid \mathbf{y}\right)$。和之前一样，我们估计后验分布使用蒙特卡洛方法。特别地，如果所有的条件分布都是已知的，然后$p\left({\mathbf{\theta }}^{\ast }\mid \mathbf{y}\right)$可以通过一系列设计好的Gibbs samplers中抽样，这就是所谓的reduced runs。

给一个具体的例子： 假定我们将参数$\beta$划分为带有完全条件分布$p\left({\mathbf{\theta }}_1\mid \mathbf{y},{\mathbf{\theta }}_2,{\mathbf{\theta }}_3\right)$，$p\left({\mathbf{\theta }}_2\mid \mathbf{y},{\mathbf{\theta }}_1,{\mathbf{\theta }}_3\right)$，$p\left({\mathbf{\theta }}_3\mid \mathbf{y},{\mathbf{\theta }}_1,{\mathbf{\theta }}_2\right)$的三个块${\beta }_1,{\beta }_2,{\beta }_3$. 我们强调这些条件分布是完全已知的，并且他们可以准确地计算：

我们首先将$p\left({\mathbf{\theta }}^{\ast }\mid \mathbf{y}\right)$分解为：
$$\begin{array}{rl}\log p\left({\mathbf{\theta }}^{\ast }\mid \mathbf{y}\right)& =\log p\left({\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_2^{\ast },{\mathbf{\theta }}_3^{\ast }\mid \mathbf{y}\right)\\ & =\log p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y}\right)+\log p\left({\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)+\log p\left({\mathbf{\theta }}_3^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_2^{\ast }\right)\end{array}$$
我们分别获得上述公式中三项的每一项。最后一项$p\left({\mathbf{\theta }}_3^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_2^{\ast }\right)$ 是已知的， 因为我们假定所有的条件分布都是已知的。剩余两项可以通过蒙特卡洛估计得到。特别地，第一项$p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y}\right)$可以使用主要的Gibbs 抽样来获得估计，然而第二项$p\left({\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)$可以通过reduced runs 进行抽样估计。接下来，我们提供计算细节。

为了估计$p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y}\right)$，首先请注意
$$\begin{array}{rl}p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y}\right)& =\iint p\left({\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_2,{\mathbf{\theta }}_3\mid \mathbf{y}\right)\mathrm{d}{\mathbf{\theta }}_2\mathrm{d}{\mathbf{\theta }}_3\\ & =\iint p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_2,{\mathbf{\theta }}_3\right)p\left({\mathbf{\theta }}_2,{\mathbf{\theta }}_3\mid \mathbf{y}\right)\mathrm{d}{\mathbf{\theta }}_2\mathrm{d}{\mathbf{\theta }}_3.\end{array}$$
因此，$p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y}\right)$可以被估计通过MC平均：
$$\widehat{p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y}\right)}=\frac{1}{R}\sum _{r=1}^{R}p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_2^{(r)},{\mathbf{\theta }}_3^{(r)}\right)$$
其中，$\left({\mathbf{\theta }}_2^{(1)},{\mathbf{\theta }}_3^{(1)}\right),\dots ,\left({\mathbf{\theta }}_2^{(R)},{\mathbf{\theta }}_3^{(R)}\right)$是从边际后验分布$p\left({\mathbf{\theta }}_2,{\mathbf{\theta }}_3\mid \mathbf{y}\right)$抽取的。也就是说，这些样本室从 main Gibbs run 中获得的。

如上所述， 原则上${\mathbf{\theta }}_1^{\ast }$的可以被设置为任意值，只要它在后验分布的支撑集中即可，可是实际上，${\mathbf{\theta }}_1^{\ast }$选择的是某一个具有高密度的点。从上述的MC估计可以看出，原因是清楚的，对于大部分组$\left({\mathbf{\theta }}_2^{(r)},{\mathbf{\theta }}_3^{(r)}\right)$，$p\left({\mathbf{\theta }}_1^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_2^{(r)},{\mathbf{\theta }}_3^{(r)}\right)$是足够大，反过来，估计量在数值上将是非常稳定的。

相似地， 注意
$$\begin{array}{rl}p\left({\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)& =\int p\left({\mathbf{\theta }}_2^{\ast },{\mathbf{\theta }}_3\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)\mathrm{d}{\mathbf{\theta }}_3\\ & =\int p\left({\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_3\right)p\left({\mathbf{\theta }}_3\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)\mathrm{d}{\mathbf{\theta }}_3\end{array}$$
因此，$p\left({\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)$可以通过如下公式进行估计：
$$p\left(\widehat{{\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }}\right)=\frac{1}{R}\sum _{r=1}^{R}p\left({\mathbf{\theta }}_2^{\ast }\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_3^{(r)}\right)$$
其中，${\mathbf{\theta }}_3^{(1)},\dots ,{\mathbf{\theta }}_3^{(R)}$是从分布$p\left({\mathbf{\theta }}_3\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast }\right)$中抽样的。这些样本可以通过reduced run获得。

更具体地说， 我们初始化${\mathbf{\theta }}_2^{(0)}={\mathbf{a}}_0$ and ${\mathbf{\theta }}_3^{(0)}={\mathbf{b}}_0$。然后，我们重复如下步骤$r=1,\dots ,R$：

1. Draw ${\mathbf{\theta }}_2^{(r)}\sim p\left({\mathbf{\theta }}_2\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_3^{(r-1)}\right)$.
2. Draw ${\mathbf{\theta }}_3^{(r)}\sim p\left({\mathbf{\theta }}_3\mid \mathbf{y},{\mathbf{\theta }}_1^{\ast },{\mathbf{\theta }}_2^{(r)}\right)$.

如果在Gibbs samplers中有很多的blocks，则需要更多的reduced run。获得这些额外的样本需要更多的编程和计算时间，但是他们也降低了模拟误差，并且产生了更加准确的结果估计。

Chib方法的主要局限性是它需要完全了解吉布斯采样器中的所有条件分布。为了克服这一限制，Chib和Jeliazkov(2001)将基本方法扩展到处理某些条件分布是非标准的情况，并且需要Metropolis-Hastings步骤。然而，实现这些扩展要复杂得多。

为了阐述Chib‘s方法， 我们仍然使用带有MA(1)误差的AR(2)模型的边际似然来对PCE 膨胀建模。

1.4 Cross - Entropy Method

修正调和均值和Chib方法都是利用MCMC抽样来计算边际似然。通过构造，这些抽取的样本是自相关的，通常获取它们的成本很高。在这一节中，我们讨论了一个基于简单分布独立抽取样本的重要抽样估计。其基本思想是对抽样分布进行偏置，以便在模拟中生成更多的“重要值”。然后对样本进行加权以校正使用不同的分布，以给出无偏估计。在我们估计边际似然的背景下，考虑以下重要抽样估计量:。
$${\widehat{p(\mathbf{y})}}_{\mathrm{I}\mathrm{S}}=\frac{1}{R}\sum _{r=1}^R\frac{p\left(\mathbf{y}\mid {\mathbf{\theta }}^{(r)}\right)p\left({\mathbf{\theta }}^{(r)}\right)}{g\left({\mathbf{\theta }}^{(r)}\right)}$$
其中，${\mathbf{\theta }}^{(1)},\dots ,{\mathbf{\theta }}^{(R)}$从重要性密度$g(\cdot )$抽取的样本。 这个估计量是无偏并且相合的。但是它的性能主要取决于重要性密度的选择。

下面我们描述一种经典的交叉熵方法的变种来构造重要性密度。最初的交叉熵方法是由Rubinstein(1997,1999)开发的用于罕见事件模拟的方法，使用多层程序来获得最佳的重要性采样密度。Chan和Kroese(2012)的研究表明，MCMC方法可以在一步内更准确地获得最优重要性采样密度。我们跟随Chan和Eisenstat(2015)使用这个新变量进行边际似然估计。

交叉熵方法基于几个简单的观测值。首先，对于边际似然估计存在一个给出零方差估计的重要密度。 特别地，如果我们使用后验分布作为重要性密度，即$g(\mathbf{\theta })=p(\mathbf{\theta }\mid \mathbf{y})=p(\mathbf{y}\mid \mathbf{\theta })p(\mathbf{\theta })/p(\mathbf{y})$，这个对应的重要性抽样估计具有0方差。
$${\widehat{p(\mathbf{y})}}_{\mathrm{I}\mathrm{S}}=\frac{1}{R}\sum _{r=1}^R\frac{p\left(\mathbf{y}\mid {\mathbf{\theta }}^{(r)}\right)p\left({\mathbf{\theta }}^{(r)}\right)}{g\left({\mathbf{\theta }}^{(r)}\right)}=\frac{1}{R}\sum _{r=1}^R\frac{p\left({\mathbf{\theta }}^{(r)}\right)p\left(\mathbf{y}\mid {\mathbf{\theta }}^{(r)}\right)}{p\left({\mathbf{\theta }}^{(r)}\right)p\left(\mathbf{y}\mid {\mathbf{\theta }}_r\right)/p(\mathbf{y})}=p(\mathbf{y})$$
此时，我们们只需要生产R = 1的观测即可。我们把这个具有0方差的重要性密度记为$g^{\ast }$。

虽然原则上$g^{\ast }$给除了$p(y)$的最佳可能的估计量，但这在实际中是不能使用的，因为它的标准化常数也就是边际似然是我们想要估计的未知量。然而，这提出了一种获得最佳重要性密度的实用方法。直观地说，如果我们选择一个重要密度$g$，它“足够接近”$g^{\ast }$，以便两者的表现相似，那么由此产生的重要性抽样估计应该具有合理的准确性。因此，我们的目标是找到一个方便的密度，在定义良好的基础上与$g^{\ast }$很接近。

为此，我们考虑一个考虑由参数向量v索引的参数族 F = { f ( θ ; v ) } \mathcal{F}=\{f(\boldsymbol{\theta} ; \mathbf{v})\} F={f(θ;v)}，我们在该参数族中选择重要性密度$g$。一种方便的测量密度之间的距离的方法是K-L divergence或 Cross - Entropy disttance，测度公式如下：
$$\mathcal{D}\left(h_1,h_2\right)=\int h_1(\mathbf{x})\log \frac{h_1(\mathbf{x})}{h_2(\mathbf{x})}\mathrm{d}\mathbf{x}$$
给定此测度后，下一步我们寻找$f(\cdot ;\mathbf{v})\in \mathcal{F}$使得$\mathcal{D}\left(g^{\ast },f(\cdot ;\mathbf{v})\right)$最小。
$$\begin{array}{rl}{\mathbf{v}}_{\mathrm{c}\mathrm{e}}^{\ast }& =\underset{\mathbf{v}}{\mathrm{argmin}}\mathcal{D}\left(g^{\ast },f(\cdot ;\mathbf{v})\right)\\ & =\underset{\mathrm{v}}{\mathrm{argmin}}\left(\int g^{\ast }(\mathbf{\theta })\log g^{\ast }(\mathbf{\theta })\mathrm{d}\mathbf{\theta }-p(\mathbf{y}{)}^{-1}\int p(\mathbf{y}\mid \mathbf{\theta })p(\mathbf{\theta })\log f(\mathbf{\theta };\mathbf{v})\mathrm{d}\mathbf{\theta }\right)\end{array}$$
其中，我们使用这个事实，$g^{\ast }(\mathbf{\theta })=p(\mathbf{y}\mid \mathbf{\theta })p(\mathbf{\theta })/p(\mathbf{y})$。由于第一项不依赖于$v$，上述的最小化问题变为下述的最大化问题：
$${\mathbf{v}}_{\mathrm{c}\mathrm{e}}^{\ast }=\underset{\mathrm{v}}{\mathrm{argmax}}\int p(\mathbf{y}\mid \mathbf{\theta })p(\mathbf{\theta })\log f(\mathbf{\theta };\mathbf{v})\mathrm{d}\mathbf{\theta }$$
实际上，这个最优化问题很难得到解析解。相反我们考虑它的随机部分：
$${\hat{\mathbf{v}}}_{\mathrm{c}\mathrm{e}}^{\ast }=\underset{\mathbf{v}}{\mathrm{argmax}}\frac{1}{M}\sum _{m=1}^M\log f\left({\mathbf{\theta }}_m;\mathbf{v}\right)$$
其中，${\mathbf{\theta }}^{(1)},\dots ,{\mathbf{\theta }}^{(M)}$是后验样本。换句话说，如果我们把$f(\mathbf{\theta };\mathbf{v})$看做似然函数，则${\hat{\mathbf{v}}}_{\mathrm{c}\mathrm{e}}^{\ast }$是$v$的精确的极大似然估计。由于找到最大似然估计量是一个标准问题，求解(4.8)通常很容易。特别是，(4.8)的解析解可以明确地找到指数族(例如，Rubinstein和Kroese, 2004, p. 70)。

在实际中，选择的参数族$\mathcal{F}$通常要保证每一个成员$f(\mathbf{\theta };\mathbf{v})$是密度的乘积。例如：$f(\mathbf{\theta };\mathbf{v})=f\left({\mathbf{\theta }}_1;{\mathbf{v}}_1\right)\times \cdots \times f\left({\mathbf{\theta }}_B;{\mathbf{v}}_B\right)$，其中，$\mathbf{\theta }=\left({\mathbf{\theta }}_1,\dots ,{\mathbf{\theta }}_B\right)$ and $\mathbf{v}=\left({\mathbf{v}}_1,\dots ,{\mathbf{v}}_B\right)$。在这种情况下，可以将可能的高维最大化问题简化为B个低维问题，然后就可以很容易地解决了。

一旦最优密度找到，将其设置为$g(\cdot )=f\left(\cdot ;{\hat{\mathbf{v}}}_{\mathrm{c}\mathrm{e}}^{\ast }\right)$，并利用其构造重要性抽样估计。这种重要抽样方法的主要优点是它快速且易于实现。此外，由于它是基于独立的样本，估计器的数字标准误差很容易得到。