贝叶斯统计中的计算方法简介

来源1：Bayesian Data Analysis, Third Edition. By. Andrew Gelman. etl. Part Ⅲ Chapter10
来源2：[Rejection sampling - Wikipedia](https://en.wikipedia.org/wiki/Rejection_sampling)

贝叶斯计算主要涉及到两个步骤：计算后验分布$$p(\theta|y)$$，以及后验预测分布$$p(\tilde{y}|y)$$。针对标准的分布形式，如正态分布，gamma分布，beta分布，泊松分布等，我们可以直接使用公式计算。对于数值型的我们可以直接使用网格计算。然而，针对复杂的、不常见的或者是高维的分布，我们需要更好的计算方法，通常，最有效的方法都是联合不同的算法一起。首先，我们来简单介绍一下统计过程的积分的近似计算。为什么我们在统计中强调用仿真解决积分问题呢？因为**大多数积分都可以写成期望的形式**，所以实际上我们是用仿真解决估计量问题。

介绍之前我们先给一个例子，说明一下为何抽样的方法可以用来计算定积分的问题，该案例来自 [邓一硕《蒙特卡洛方法与定积分计算》](http://cos.name/2010/03/monte-carlo-method-to-compute-integration/ "邓一硕《蒙特卡洛方法与定积分计算》")。

首先，我们假设需要求解$$J=\int\_{0}^{1}f(x)dx$$，其中$$0 \le f(x) \le 1 $$。假设X、Y是服从正方形$$\\{ 0\le x \le 1, 0\le y \le 1 \\}$$的均匀分布，且X、Y相互独立。记事件 $$A=\\{Y\le f(x)\\}$$，那么A的概率是：

```math
P(A) = P( Y \le f(X)) = \int_{0}^{1}\int_{0}^{f(x)}dydx = \int _{0} ^{1} f(x) dx = J
```

也就是说事件A的概率就是定积分的结果。那么根据大数定理，我们重复足够多的试验，就可以把A事件出现的频率作为A的概率。事件A的频率如何求呢？这时我们就可以根据随机投影的结果结算了。我们可以随机选择n个x的点，然后计算$$f(x)$$的值，然后我们数$$y < f(x)$$的点的个数，把这个结果除以n就得到了事件A的概率了。举个例子，假设我们求$$\int \_{0} ^{1} e^{-x^2/2} / \sqrt{2\pi} dx$$的结果。我们模拟不同次数得到的抽样结果代码如下：

```r
n=10^4
x=runif(n)
y=runif(n)
plot(x,y)
f = function(x){ exp(-x^2/2)/sqrt(2*pi) }
mu_n = sum(y < f(x))
J = mu_n/n

```

我们模仿几次的结果如下：
  
![](http://www.datalearner.com/resources/blog_images/ba88b9c3-b6a7-4f17-9760-b20c606a12ae.png)

以上结果是[0,1]之间，对于非[0,1]区间可以采用映射方法。

[TOC]

###### *正规化和非正规化的密度函数*

###### *log密度*

为了避免计算溢出，如果可能的话，我们需要尽量计算后验密度的对数。对于求幂运算，我们尽量避免，且应该尽量放在后面运算。比如，在Metropolis算法中，我们需要计算两个密度的比值，这就应该作为log密度的计算后再指数化。

#### 1、数值积分

数值积分是指在有限的数据点上通过计算函数的值，来求该连续函数的值。通过增加点的数量，我们可以得到达到指定精确度的函数结果。数值积分的方法可以分成仿真类（如蒙特卡洛）和确定性计算（如利用各种积分规则）。

任意函数$$h(\theta)$$的后验期望可以定义成$$E(h(\theta)|y)=\int h(\theta)p(\theta|y)d \theta$$，其中积分的维度和$$\theta$$的维度一样多。相反的，我们可以把$$\theta$$空间上任意的积分都定义成一个合适的$$h(\theta)$$的后验期望。如果我们要从$$p(\theta|y)$$中抽取$$\theta^{s}$$，我们可以通过样本均值$$\frac{1}{S}\sum\_{s=1}^{S}h(\theta^{s})$$估计积分。对于任意有限个仿真样本，其估计精度可以通过$$h(\theta^{s})$$的标准差测量出来。如果从后验分布中抽取样本很困难，或者是$$h(\theta^{s})$$变化无常（也就是说样本均值太过于估计而没有用处），这时候就需要更多的抽样算法。

##### 1.1 仿真方法

仿真方法是从期望分布$$p(\theta)$$中获取随即样本$$\theta^{s}$$，然后估计任意函数$$h(\theta)$$的期望（公式1）：

```math
E(h(\theta)|y) = \int h(\theta)p(\theta|y) d\theta \approx = \frac{1}{S}\sum_{s=1}^{S}h(\theta^{s})
```

这个估计是随机的，它的结果取决于生成的随机数值，但是估计的精度可以通过获取更多的样本来提高。仿真的方法可以用在高维的分布中，也有通用的方法可以运用在广大的模型中，如果必要，也可以通过联合一般的想法，如仿真方法，确定性计算以及分布近似等。

##### 1.2 确定性方法

确定性数值积分方法是在给定的点$$\theta^{s}$$估计积分$$h(\theta)p(\theta|y)$$，是公式1的加权版本：

```math
E(h(\theta)|y) = \int h(\theta)p(\theta|y) d\theta \approx = \frac{1}{S}\sum_{s=1}^{S} w_{s}h(\theta^{s})p(\theta^{s}|y)
```

权重$$w\_{s}$$与点$$\theta^{s}$$所代表的空间大小有关。更加详细的方法如Simpson的，使用局部多项式来提高精确度。精确性数值积分规则比仿真方法有更好的可靠性，但是在高维情况下选择局部函数非常困难。

最简单的确定性计算是在等权重的情况下计算网格中的积分。网格方法可以适用于从后验模式中开始的网格公式。对于积分中包含的一些特殊的形式，如高斯形式，我们可以使用特殊的规则来计算结果。

##### 1.3、分布近似

分布近似是使用几个简单的参数分布来近似后验分布，从而使得积分可以被直接计算或者通过一些基本的仿真方法来近似起始点。

###### *通过忽略某些信息进行粗略估计*

在使用精心考虑的近似方法或者复杂的抽样方法求目标分布前，获取一个大概的目标分布的局部估计结果是非常有用的，也就是说模型中的参数进行点估计，这可以使用一些简单的非迭代的方法。首次估计的方法与问题很相关，但一般都是设计模型或者数据中的某部分，创建一个简单的问题，然后可以使用简单的估计方法。

在层次模型中，我们可以首先粗略估计超参数$$\phi$$，然后计算条件后验分布$$\gamma|\phi,y$$，即可以获取大致的主要参数$$\gamma$$的估计结果。

除了为更加精确的分析创造一个开始的点以外，粗略估计可以用来比较后验的结果。如果粗略估计的结果与完全分析的结果相差很大，那么后验的分析可能在模型或者程序上有问题。因为粗略估计可以用现有的计算机程序计算，所以它是比较方便可靠的。

#### 2、直接仿真和拒绝抽样

在简单的非层次贝叶斯模型中，直接从后验分布中抽取样本很简单，也就是在假设共轭先验的情况下。对于更加复杂的问题，将分布分解后再局部仿真是很有用的。我们可以首先抽取超参数的边缘后验分布的样本，然后在数据条件下抽取其他参数，以及仿真的超参数。有些情况下从大问题中分解几部分，直接仿真并分析其积分是可能的。

通常情况下，从标准分布或者低维非标准分布中抽取样本是必须的，不管是在简单问题中抽取估计量的后验分布样本，还是在复杂问题中某个步骤。

##### 2.1、通过计算坐标点直接近似

在最简单的离散近似中，我们可以首先在等距离的空间值$$\theta\_{1},\...,\theta\_{N}$$（它们可以覆盖参数空间$$\theta$$的范围）上计算目标密度，$$p(\theta|y)$$。然后通过在$$\theta\_{1},\...,\theta\_{N}$$点上的概率$$p(\theta\_{i}|y)/\sum\_{j=1}^{N}p(\theta\_{j}|y)$$，近似连续分布$$p(\theta|y)$$。因为近似的密度必须是正规化的，因此这个方法需要使用非正规化密度函数$$q(\theta|y)$$来代替$$p(\theta|y)$$。

一旦密度值的坐标计算出来，从[0,1]上的均匀分布中抽取一个随机样本U，通过这种方式可以获取$$p(\theta|y)$$的一个随机样本。当点$$\theta\_{i}$$间隔足够近的时候，没有缺少任何重要的界限，那么这个方法将会非常好用。离散近似在高维多元问题中的使用非常困难，因为计算多维坐标点的密度代价很大。

##### 2.2、从预测分布中近似

一旦我们有了后验分布的一个样本，我们就很容易从预测分布中抽取新的数据$$\tilde{y}$$。对于后验分布中任意一个$$\theta$$的抽取，我们都可以直接从预测分布$$p(\tilde{y}|\theta)$$中抽取一个$$\tilde{y}$$。从所有的$$\theta$$'s中抽取的近似$$\tilde{y}$$集合可以描述后验预测分布。

##### 2.3、拒绝抽样（Rejection Sampling）

拒绝抽样是从分布中产生观测值的一种基本方法，也称作**接受-拒绝抽样（acceptance-rejection sampling）**。它也是蒙特卡洛方法的一种。其主要思想是，目标分布的样本抽取比较困难，所以我们可以通过抽取另一个我们能抽取样本分布的分布来代替抽取目标分布的样本。这里需要解决的一个问题就是，如何构造另一个分布，使得我们对其抽样的结果也是来自于目标分布。我们可以假设替代分布为$$Y$$，其概率密度为$$g(\theta)$$。拒绝抽样是我们从$$Y$$抽取样本，接受的概率为$$p(\theta)/(Mg(\theta))$$。这里$$M$$是常量，且$$1 < M < \infty$$，它满足$$p(\theta) \le Mg(\theta)$$，也就是说$$Mg(\theta)$$的图像永远在$$f(\theta)$$的上方。拒绝抽样的有效性的验证可以通过**包络原则（envelop principle）**说明。我们抽取一个样本的时候，其实是会同时抽取两个数值$$(x,v=u \cdot Mg(x))$$。这里的$$x$$是从$$g(\theta)$$中抽取的，$$u$$是从[0,1]均匀分布中抽取的。只有当$$u < f(\theta)/Mg(\theta)$$的时候，我们才会接受这次抽取的结果，该样本就是来自$$p(\theta)$$的一个样本。这是什么原因呢？我们可以这样理解，每次抽取的时候我们都抽取了两个值，即$$\theta$$和来自[0,1]均匀分布的一个值$$u$$。当$$u < p(\theta)/Mg(\theta)$$时，我们认为该点落在$$p(\theta)$$下方，是来自于目标分布的，因此接受。当$$u > p(\theta)/Mg(\theta)$$时，我们认为该点落在$$p(\theta)$$和$$Mg(\theta)$$之间。我们拒绝它，因为它不是$$p(\theta)$$的样本。该过程如图10.1所示。外侧曲线是$$Mg(\theta)$$，其结果永远大于$$p(\theta)$$。该过程更加严谨的表述如下（下面所说的后验密度$$p(\theta|y)$$就是这里说的$$p(\theta)$$，这里我们把它当作任意的密度，下面是以贝叶斯推断中，估计后验概率为例，本质一样）：

假设我们希望从密度函数$$p(\theta|y)$$中（或者是一个非正规化的密度$$q(\theta|y)$$）获取一个随机抽取结果。我们使用$$p$$表示目标分布（但是我们也可以针对非正规化形势$$q$$做同样的操作）。为了构造拒绝抽样，我们需要一个正函数$$g(\theta)$$，其定义在所有的$$\theta$$上，其中$$p(\theta|y) > 0$$，其有如下特性：

1.我们可以按一定概率从$$g$$的概率密度中抽取样本（也就是上述中需要按照均匀分布的结果来考虑是否接受样本），我们不需要$$g(\theta)$$积分等于1，但是$$g(\theta)$$必须是有限积分。
　　2.重要率$$\frac{p(\theta|y)}{g(\theta)}$$必须有一个已知的界限，也就是说必须存在一个常量$$M$$，使得对所有的$$\theta$$存在$$\frac{p(\theta|y)}{g(\theta)} \le M$$。

拒绝抽样算法有两个主要的步骤：

1.按比例从$$g(\theta)$$的概率密度中抽取随机样本$$\theta$$。也就是产生样本$$\theta \sim g(\theta)$$，$$U \sim Uniform[0,1]$$
　　2.当$$U \le \frac{p(\theta|y)}{Mg(\theta)}$$，接收$$\theta$$是从p中抽取的一个样本。反之，$$\theta$$被拒绝，继续执行步骤1。

![](http://www.datalearner.com/resources/blog_images/5f4ef824-8546-472e-aacf-1a27d58522b1.png)

上图说明了拒绝抽样。接受的$$\theta$$有正确的分布$$p(\theta|y)$$。也就是说，在$$\theta$$被接受的条件下，$$\theta$$抽取的样本的分布是$$p(\theta|y)$$。有界的条件是必须的，才能保证步骤2不会产生大于1的结果。

对于拒绝抽样，一个好的近似密度$$g(\theta)$$应当大约与$$p(\theta|y)$$成比例。这个思想的条件是$$g \propto p$$。也就是在一个合适的值$$M$$下，我们可以以概率1接受所有的抽取结果。当$$g$$与$$p$$不成比例的时候，M的边界必须设置得很大以至于几乎所有的抽取都会被拒绝。拒绝抽样是自我监督的——如果方法不奏效的话，几乎很少有近似的样本会被接受。

函数$$g(\theta)$$是被选择用来近似$$p(\theta|y)$$。所以一般情况下它依赖于y。我们不使用$$g(\theta,y)$$或者是$$g(\theta|y)$$，然而，在实际中我们需要考虑一次对一个后验分布的近似，函数上$$g$$依赖于$$y$$我们不感兴趣。

拒绝抽样一般用于从标准一元分布中进行快速抽样。

我们用R语言来举个例子，假设我们希望抽取beta(4,4)的样本，根据图形我们可以构造一个M=7，$$g(\theta)=-1\*(x-0.5)^{2}+0.33$$的替代分布，产生的样本使用KS检验是否来自于beta分布，代码如下：

```r
n = 10^4

M < - 7
g = function(x){
	-1*(x-0.5)^2+0.33
}

p = function(x){
	dbeta(x,4,4)
}

x_t < - runif(n)

m < -  ifelse( x_t <= p(x)/M*g(x), x, NA)

plot(0,0,xlim=c(0,1),ylim=c(0,3))
lines(x,M*g(x),col='red')
lines(x,p(x),col='blue')

length(na.omit(m))
ks.test(na.omit(m),'pbeta',4,4)
```

其运行结果如下（p值大于0.05且D较小，证明来自于beta(4,4)的分布）：
![](http://www.datalearner.com/resources/blog_images/f0de7632-5acc-4609-a231-408fdf1dbf6a.png)
![](http://www.datalearner.com/resources/blog_images/bf738e57-0ddf-449a-b1e3-0982516d94c5.png)

我们可以通过如下方式从$$p(x)$$分布中抽取样本。这是一个比例常量，通过从另一个容易获取样本的分布$$q(x)$$中抽取，选择拒绝还是接受样本，只需要满足$$p(x)\leq Mq(x),M<\infty $$条件即可。拒绝或者接受的方式如下（可以参见图2）：

![](http://www.datalearner.com/resources/blog_images/ccdae751-0562-490e-b1f9-6eae15391a7c.jpg)
![](http://www.datalearner.com/resources/blog_images/c626d423-ed47-40b1-8c36-f122b0c2a4f3.jpg)

这个方法需要满足一个严格的条件。我们并不总能给$$p(x)/q(x)$$一个合理的常量M。当M很大的时候，接受的概率如下：
```math
Pr(x accepted) = Pr(u < \frac{p(x)}{Mq(x)}) = \frac{1}{M}
```
　　它的值将会很小使得高维场景中这个方法不现实。
  
#### 3、重要性抽样

重要性抽样与拒绝抽样有关，是Metropolis算法的前身，使用目标分布的近似样本的一个随机抽取计算期望。假设我们对$$E(h(\theta)|y)$$感兴趣，但是我们无法从$$p(\theta|y)$$中获取$$\theta$$的样本，进而无法通过近似值的简单平均来估算积分结果。

在很多应用中我们都希望计算$$\mu=E(f(x))$$，当$$P(X \in A)$$很小的时候，$$f(x)$$在区域A的外面几乎为0。集合A可能是很小的值，或者是在X分布的尾巴部分。一个简单的蒙特卡洛抽样可能连一个A区域的样本都抽不到。从感兴趣或者重要的区域抽取样本是我们必须要做的事情。我们可以通过抽取一个分布，其在我们重要的目标区域内具有较高的权重，因此成为重要性抽样。为了在这个区域抽样，我们需要调整估计对象。重要性抽样有很多好处，它使得原本无法解决的问题在蒙特卡洛里面变得可以解决。但它也可能使那些有限方差的简单的蒙特卡洛方法变成无限的方差。它是最难用好的方差缩减方法。重要性抽样不仅仅是一个方差缩减方法。它可以使我们从另外的分布中抽取样本来研究当前分布。我们可以使用重要性抽样替代接受-拒绝抽样，用来做敏感性分析或者作为计算概率密度正规化常数的基础。重要性抽样是蒙特卡洛序列的准备工作。

如果$$g(\theta)$$是一个我们能产生随机分布的一个概率密度，我们可以写作：

```math
E(h(\theta|y))=\frac{ \int h(\theta)q(\theta|y) d\theta }{ \int q(\theta|y) d\theta } = \frac { \int [h(\theta)q(\theta|y)/g(\theta)]g(\theta) d\theta }{ \int [q(\theta|y)/g(\theta)] g(\theta) d\theta }
```

我们可以通过如下表达式，根据从$$g(\theta)$$中抽取的S个样本$$\theta^{1},\...,\theta^{S}$$估计（公式3）：

```math

\frac { \frac{1}{S} \sum ^{S} _{s=1} h(\theta^{s})w(\theta^{s}) } { \frac{1}{S} \sum ^{S} _{s=1} w(\theta^{s}) }

```

其中因子：

```math

w(\theta^{s}) = \frac {q(\theta^{s}|y)} { g(\theta^{s}) }

```

称为重要性比率，或者是重要性权重。之前说过，$$q$$是非正规化密度的一般符号，也就是说$$q(\theta|y)$$等于$$p(\theta|y)$$乘以一些与$$\theta$$无关的因子。上式告诉我们，对原来的$$h(\theta)$$抽取样本可以转化成从$$g(\theta)$$中抽取$$h(\theta)\*w(\theta^s)$$的样本。当$$g(\theta_i)$$在某个点的值较大，那么这个位置的x抽取的概率就会增大。

一般情况下，使用相同的随机抽取计算公式3中的分子和分母是明智的，这样可以减少估计中的抽样错误。
如果$$g(\theta)$$可以选择使$$\frac{hq}{g}$$为一个大致的常量，那么我们就可以获得积分的相对较高的精度。如果重要性比率变化很大，那么重要性抽样并不实用。最坏的情况就是重要性比率很大概率较小，很小的概率较大，比如$$hq$$相比较$$g$$而言有很长的尾巴。

上式比较难以理解，我们可以看另外一个容易理解的说明。假设，我们的目标是计算均值$$\mu = E( f(x)) = \int\_{D}f(x)p(x)dx$$，其中p是概率密度函数，f是被积函数。那么我们有：
```math
\mu = \int_{D} f(x)p(x) dx=\int_{D} \frac { f(x)p(x) } { q(x) } q(x) dx = E_{q}( \frac{ f(X)p(X) } { q(X) } )
```
其中，$$E\_{q}(\cdot)$$表示$$X \sim q$$的期望。我们起初的目标是计算$$E\_{p}(f(x))$$。通过做乘法调整后，我们可以通过对$$q$$抽样来替代$$p$$。这里调整因子$$p(x)/q(x)$$是似然率（likelihood ratio）。分布$$q$$是重要性分布，$$p$$是名义分布（nominal distribution）。因此，对原始问题的估计结果就是：
```math
\hat{\mu_{q}} = \frac {1} {n} \sum_{i=1}^{n} \frac { f(X_i) p(X_i)} {q(X_{i})}, X_i \sim q
```

为了使用上式，我们需要计算$$fp/q$$。假设我们$$f$$是可计算的，那么我们只需要可以计算$$p(x)/q(x)$$。当p或者q有未知的正规化常数时，我们可以使用比例估计。目前我们假设它们是可估计的，我们来研究一下$$\hat{\mu\_{q}}$$的方差。我们这里直接给出该方差的计算结果，具体证明可参见[Importance Sampling By Art Owen](http://statweb.stanford.edu/~owen/mc/Ch-var-is.pdf "Importance Sampling By Art Owen")

```math
\sigma_{q}^{2} = \int_{D} \frac{ (f(x)p(x) - \mu q(x) )^{2}} { q(x)} dx
```

我们对其估计的结果是：

```math
\hat{\sigma}_{q}^{2} = \frac{1}{n} \sum_{i=1}^n (w_{i}f(x_i)-\hat{\mu}_{q})^2
```

从这里我们可以看出，好的q的选择应当是具有较小的$$\int\_{D}((fp)^2 /q)dx$$（均值）。同时上面第一个式子也说明了重要性抽样的能否成功的关键。当$$f(x)p(x)-\mu q(x)$$接近0的时候，被积函数的分子比较小，也就是说$$q(x)$$与$$f(x)p(x)$$是大约成比例的。
我们举个例子说明一下，假设对于所有的x，$$f(x) \ge 0$$。我们假设$$\mu >0 $$。那么$$q\_{opt} \equiv f(x)p(x) / \mu$$是一个概率密度，且$$\sigma\_{q\_{opt}}^2=0$$。这是最优的情况，因为$$\hat{\mu}\_{q\_{opt}}$$就是$$f(x\_{i})p(x\_{i})/q(x\_{i}) = \mu$$的均值，这意味着我们可以直接从$$f,p,q$$中计算出$$\mu$$而不需要任何抽样。尽管这是不可用的，但是我们可以据此看到，在同一个位置，q有spikes（这里的spikes原意是尖状物，可以理解为极大值点）比f或者p好，但最好的是fp有spikes。最好的q的选择应当是与fp成比例的。选择好的重要性抽样分布需要良好的猜测或者是可能的数值搜索。很多应用领域需要领域知识来发现spikes。在金融领域我们一般可以知道哪种股票的波动可以导致它出现最大值。对于查询系统我们可以知道什么样的请求会导致过载。

可参考 [重要性抽样案例](http://dept.stat.lsa.umich.edu/~jasoneg/Stat406/lab7.pdf "重要性抽样案例")

重要性抽样的精度和效率
一般情况下，当我们缺少精确或者近似密度的数学分析的形式，我们可以忽略一些不常见的且很大的重要性权重。然而，它却可以帮助我们检测抽取的重要性权重的分布以发现可能的问题。它可以用来检测最大的重要性比率的log直方图：如果最大的比率相比较均值过大的话估计结果通常不好。相反，我们不需要考虑较小的重要性比率的影响，因为它们对公式2的影响很小。如果权重的方差是有限的，那么有效的样本大小可以通过如下公式估计：

```math

S_{eff} = \frac {1} { \sum_{s=1}^{S} (\tilde{w}(\theta^{s}))^{2} }

```

其中，$$\tilde{w}(\theta^{s})$$是正规化权重，也就是说$$\tilde{w}(\theta^{s}) = w(\theta^{s})S/\sum^{S}\_{s,=1}w(\theta^{s,})$$。如果有很少的较大的权重会影响分布，那么有效样本的大小就会比较小。反之，这个估计就是自我噪音的。

##### 重要性重抽样
为了获取等权重的独立样本，使用重要性重抽样的方法是可以的。
一旦我们从近似分布$$g$$中抽取到了S个样本$$\theta^{1},\...,\theta^{S}$$，近似k（k< S）个抽取的方法可以按如下方式。
1、从集合$$\{\theta^{1},\...,\theta^{S}\}$$中抽取一个值$$\theta$$，每个$$\theta^{s}$$的抽取概率与权重成比例，$$w(\theta^{s}) = \frac { q(\theta^{s}|y) } { g(\theta^{s}) } $$
2、采用同样的方式抽取第二个样本值，但要除去第一个结果。
3、重复不替换抽样至(k-2)次。

为什么要不替换抽样？如果权重大小适中，那么是否替换抽样不影响结果。但是我们考虑一个坏情况，有一些很大的权重，还有一些较小的权重。替换抽样江湖重复选取相同的较少的值，想法不替换抽样会产生一个更加符合期望的近似结果。

#### 4、我们需要多少个仿真结果？

贝叶斯推断最方便的就是从模型参数的后验分布中随机抽样并汇总。一元估计量的分位数可以较好地传递分布的形状。比如某个估计量的样本分布的2.5%,25%,50%,75%和97.5%的点提供了50%和95%的区间，并描述了其边缘后验密度的偏度。仿真结果的散点图、轮廓图以及更高级的图形技术可以用来检验二维或者三维的后验分布。

我们也会使用后验仿真来推断预测量。给定每个抽取结果$$\theta^s$$，我们可以获取任意预测量，$$\hat{y}^s \sim p(\hat{y}|\theta^s)$$，或者对于回归模型，$$\hat{y}^s \sim p(\hat{y}|\hat{X},\theta^s)$$。

最后，给定一个仿真结果$$\theta^s$$，我们可以仿真出一个复制到数据集$$y^{rep-s}$$。我们可以比较数据集与我们复制的结果来评价模型的好坏。

在贝叶斯计算中，我们希望从后验分布中获取一组独立的抽取结果$$\theta^s$$，其中$$s=1,\...,S$$。当样本足够多的时候，我们就可以一定的精度来估计我们感兴趣的量。在大多数例子中，100个样本就足够用来计算后验汇总的量。我们可以考虑一个近似正态后验分布（均值为$$\mu\_{\theta}$$方差为$$\sigma\_{\theta}$$）的一个实量参数$$\theta$$。我们假设无法直接计算，需要通过仿真从均值$$\bar{\theta}$$和标准差$$s\_{\theta}$$上估计。那么后验均值的正确性近似$$s\_{\theta}/\sqrt{S}$$。标准差的为$$s\_{\theta}\sqrt{1+1/S}$$。当S=100时，因子$$\sqrt{1+1/S}=1.005$$意味着蒙特卡洛增加的误差几乎不增加来自实际后验方差的不确定性。当然，如果能很方便增加超过100个样本的画可以增加数值汇总的稳定性，尽管有时候这个稳定性增加没卵用。

在某些后验推断中，我们需要更多的仿真结果来获得期望中的准确率。

相关博客

最热博客