Introduction to Bayesian Inference

Introduction to Bayesian Inference

Introduction

近年来，由于计算统计学的重要进展，贝叶斯推理得到了蓬勃发展。本书将重点介绍用于近似贝叶斯推理的集成嵌套拉普拉斯近似（INLA、Havard Rue、Martino 和 Chopin 2009）。 INLA 是贝叶斯统计领域最近的几项计算突破之一，可实现快速、准确的模型拟合。

本简介的目的不是全面介绍贝叶斯推理，而是介绍本书其他章节的一些符号和上下文。对于那些可能需要的读者，最近的贝叶斯推理介绍性文本包括 Kruschke (2015) 和 McElreath (2016)。这两本书都涵盖了基础知识，并提供了大量使用多种编程语言进行贝叶斯计算的示例。更高级的文本包括 Gelman 等人的经典著作。 (2013) 以及 Carlin 和 Louis (2008)，它们提供了有关贝叶斯推理的广泛主题的更深入的细节。

Bayesian inference

在贝叶斯范式中，模型中的所有未知量都被视为随机变量，目的是计算（或估计）联合后验分布。就是说参数$\theta$的分布在给定观测结果数据$y$的情况下，后验分布的获得方式依赖于贝叶斯定理：
$$\pi (\mathbf{\theta }\mid \mathbf{y})=\frac{\pi (\mathbf{y}\mid \mathbf{\theta })\pi (\mathbf{\theta })}{\pi (\mathbf{y})}$$
这里，$\pi (y|\theta )$是数据$y$在给定参数$\theta$的情况下(也就是取参数空间$\Theta$中的值，$\pi (\theta )$是参数的先验分布而且$\pi (y)$是边际似然函数，它充当归一化常数。边际似然通常很难计算，因为它等于
$${\int }_{\Theta }\pi (\mathbf{y}\mid \mathbf{\theta })\pi (\mathbf{\theta })d\mathbf{\theta }$$
注意到后验分布$\pi (\theta |y)$通常是一个高度多元的分布，并且它仅适用于少数模型的封闭形式，因为边际似然函数$\pi (y)$通常很难估计。实际上，后验分布是在不计算边际似然的情况下估计的。这就是为什么贝叶斯定理经常被表述为
$$\pi (\theta \mid \mathbf{y})\propto \pi (\mathbf{y}\mid \mathbf{\theta })\pi (\mathbf{\theta })$$
这意味着可以通过重新缩放似然性和先验的乘积来估计后验，以便其积分达到 1。

模型的似然描述了给定参数的数据生成过程，而先验通常反映了有关模型参数的任何先前知识。当这种知识稀缺时，会假设模糊的先验，以便后验分布由观察到的数据驱动。

$\theta$参数的汇总统计可以从$\pi (\theta |y)$的联合后验分布中获得另外，每个元素的后验边缘分布
$\theta$将其余参数积分即可得到：
$$\pi ({\theta }_i\mid \mathbf{y})=\int \pi (\mathbf{\theta }\mid \mathbf{y})d{\theta }_{-i},i=1,\dots ,\dim (\mathbf{\theta })$$
这里，${\theta }_{-i}$代表$\theta$参数集中排除${\theta }_i$，后边缘分布对于总结各个参数和计算汇总统计量非常有用。重要的量是后验均值、标准差、众数和分位数。

Conjugate priors

如上所述，后验分布仅适用于少数模型的封闭形式。具有共轭先验的模型是先验与似然具有相同形式的模型。例如，如果似然度是具有已知精度的高斯分布，则均值上的共轭先验是高斯分布。这将确保均值的后验分布也是高斯分布。

例如，假设观测值 { y i } i = 1 n \{y_i\}^n_{i=1} {yi​}i=1n​遵循高斯分布
$$y_i\mid \mu ,\tau \sim N(\mu ,\tau ),i=1,\dots ,n$$
其中$\mu$是未知的均值和$\tau$是标准差。关于$\mu$的先验可以是关于${\mu }_0$和${\tau }_0$的标准分布.
$$\mu \sim N({\mu }_0,{\tau }_0)$$
关于$\mu$的后验分布，(例如，在对给定数据y的情况下$\mu$的分布)是$N({\mu }_1,{\tau }_1)$
$${\mu }_1={\mu }_0\frac{{\tau }_0}{{\tau }_0+\tau n}+\overline{y}\frac{\tau n}{{\tau }_0+\tau n}$$
这很好地说明了贝叶斯学习。首先，后验均值是先验均值${\mu }_0$和$\bar{y}$之间的折衷。当观测值数量较多时，后验均值接近数据的均值。如果n较小时，我们对平均值的先前信念会更加重视。类似地，后验精度是先验精度和似然精度的函数，并且它趋于无穷大n，这意味着方差μ随着数据量的增加趋于零。

其他众所周知的共轭模型包括 Beta-二项式模型和 Poisson-Gamma 模型（Gelman 等人，2013）。

Computational methods

当后验分布无法以封闭形式获得时，有必要采用其他方法来估计它，或者从中抽取样本。给定后验样本，遍历定理确保可以估计矩和其他感兴趣的量（Brooks et al. 2011）。

一般来说，计算方法旨在估计贝叶斯推理中出现的积分。例如，参数${\theta }_i$的后验均值(参数空间中的值${\Theta }_i$)按照以下方式计算：
$${\int }_{{\Theta }_i}{\theta }_i\pi ({\theta }_i\mid \mathbf{y})d{\theta }_i$$
$\pi ({\theta }_i|y)$的分布是单变量参数${\theta }_i$的边际后验分布。类似地，也可以计算后验方差或任何其他矩。这种类型的积分可以使用数值积分方法和拉普拉斯近似（Tierney 和 Kadane 1986）方便地近似。

此外，从已知常数的密度中进行采样的典型蒙特卡罗方法可用于从后验分布中获取样本。这些方法包括拒绝抽样和其他方法（Carlin 和 Louis 2008）。然而，大多数这些方法在高维空间中效果不佳。

模型参数的点估计可以通过最大化似然与先验的乘积来获得。这由最大后验 (MAP) 估计来表示。后验分布的最大化通常以对数尺度进行，并且可以使用不同的方法有效地实现，例如牛顿-拉夫森算法或期望最大化算法（Gelman et al. 2013）。

Markov chain Monte Carlo

马尔可夫链蒙特卡罗 (MCMC) 方法（Gilks 等人，1996 年；Brooks 等人，2011 年）是一类从联合后验分布中抽取样本的计算方法。这些方法基于构建具有平稳分布（后验分布）的马尔可夫链。因此，通过从该马尔可夫链中重复采样，在多次迭代后获得联合后验分布的样本。

已经提出了几种算法（请参阅 Brooks 等人 2011 年的最新摘要），但其中两种算法为贝叶斯推理提供了工具箱：Metropolis-Hastings 算法和吉布斯采样。

Metropolis-Hastings 算法使用提议分布来获取参数集合的新值。也就是说给定当前的参数$\theta$一个新的参数${\theta }^{\prime }$从概率密度函数$q(\cdot |\theta )$中获得。可以通过多种不同方式选择提案分配。然而，并非所有提案都会被接受。当提出新提案时，它有可能被接受
$$\min \left\{1,\frac{q(\mathbf{\theta }\mid {\mathbf{\theta }}^{\prime })\pi ({\mathbf{\theta }}^{\prime }\mid \mathbf{y})}{q({\mathbf{\theta }}^{\prime }\mid \mathbf{\theta })\pi (\mathbf{\theta }\mid \mathbf{y})}\right\}$$
这里，$q(\theta |{\theta }^{\prime })$表示建议的概率分布密度$q(\cdot \mid {\mathbf{\theta }}^{\prime })$所评估的$\theta$。另请注意，接受概率包含要估计的后验分布（且未知）。然而，先前的接受概率可以重写为
$$\min \left\{1,\frac{q(\mathbf{\theta }\mid {\mathbf{\theta }}^{\prime })\pi (\mathbf{y}\mid {\mathbf{\theta }}^{\prime })\pi ({\mathbf{\theta }}^{\prime })}{q({\mathbf{\theta }}^{\prime }\mid \mathbf{\theta })\pi (\mathbf{y}\mid \mathbf{\theta })\pi (\mathbf{\theta })}\right\}$$
通过使用贝叶斯规则。请注意，现在缩放常数$\pi (y)$是不需要计算比率，这就是 Metropolis-Hastings 算法如此流行的原因。

该比率本质上确保新提出的值以正确的比例被接受为后验分布的样本。此外，当很有可能回到当前值时，该比率也将有利于转向新值。这是有道理的，因为这将避免陷入低后验概率的区域。

Gibbs 采样（Geman 和 Geman 1984）可以被视为 Metropolis-Hastings 算法的一个特例，其中提议分布是模型参数的完整条件分布，即$\pi ({\theta }_i|y,{\theta }_{-i})$这确保了 Metropolis-Hastings 算法中的接受概率始终为 1，这意味着所有提案都会被自动接受。

在实践中，收敛到后验分布将需要多次迭代，并且第一批模拟（或老化）将被丢弃以进行推理。此外，考虑到样本通常是相关的，对剩余样本进行稀疏化，其中许多样本也将被丢弃（例如，参见 Link 和 Eaton 2012 以及其中的参考文献）。可以使用几个标准来评估模拟是否已收敛到静止状态（Gelman 和 Rubin 1992；Brooks 和 Gelman 1998）。

The integrated nested Laplace approximation

Havard Rue、Martino 和 Chopin (2009) 开发了一种新颖的贝叶斯推理计算方法。特别是，他们专注于估计模型参数的后验边缘。因此，不要估计高度多元的联合后验分布$\pi (\theta |y)$他们专注于获得单变量后验分布的近似值$\pi ({\theta }_i|y)$。他们将这种近似称为积分嵌套拉普拉斯近似（INLA）。

尽管 INLA 可以拟合的模型仅限于可以表示为潜在高斯马尔可夫随机场的模型类别（Rue 和 Held 2005），但这包括无数常用模型。此外，作者还开发了一个 R 包 INLA，用于轻松进行模型拟合。

An introductory example: U’s in Game of Thrones books

为了介绍上述贝叶斯推理的不同方法，这里开发了一个介绍性示例。请注意，目的是展示 INLA 方法在贝叶斯生态系统中的适用位置，并且需要一些贝叶斯推理知识才能遵循该示例。

我们将考虑《权力的游戏》系列中任何一本书的一页上“u”的平均数量。对于数据收集，仅考虑了 Bantam Books 于 2001 年出版的五本书中每本书的前言页面。不考虑带有标题或不完整（即缺少行）的页面。数据由 2018 年 7 月在巴伦西亚（西班牙）举办的第二届巴伦西亚国际贝叶斯暑期学校（VIBASS2，http://vabar.es）的参与者收集。

数据以 CSV 文件形式提供，可以按如下方式加载和汇总：

GoT <- read.csv2(file = "data/GoT.csv")
summary(GoT)

总共有 31 个观察结果。因为这是一个计数过程，所以泊松似然似乎是一个合理的选择：
$$U_i\sim Po(\lambda ),i=1,\dots ,31$$
参数$\lambda$是泊松分布的平均值和我们感兴趣的参数。请记住，泊松分布的概率分布是
$$P(X=x\mid \lambda )=\frac{\exp (-\lambda ){\lambda }^x}{x!}$$
注意到字母$u$不是英语中最常见的元音，这个事实可以用来对先验分布进行选择。每页大约有 20 行，假设我们确定每行大约有一到两个 u。该先验猜测可用于设置参数的先验$\lambda$。

Gamma 分布是泊松似然的共轭先验（Gelman et al. 2013）。这是由两个参数定义的，a和b，定义平均值$a/b$,和方差$a/b^2$。伽马分布的概率密度为
$$f_X(x\mid a,b)=b\exp (-b\cdot x)\frac{(b\cdot x{)}^{a-1}}{\Gamma (a)}$$
给定我们的先验假设，我们可以取一个均值为30方差为5的伽马值，这样就可以得到一个很高的先验概率$\lambda$介于20到40之间。因此先验伽马分布的参数为a=180和b=6.

​ 图 1.1：《权力的游戏》一页上 u 数量的可能性和先验。

可视化可能性和先验的常见方法是将它们绘制在一起。图 1.1 显示了先验和似然。请注意它们的比例不同，但模式很接近，这意味着我们之前的猜测毕竟没有那么糟糕。后验分布将通过使用贝叶斯规则结合先验和似然来获得。

最大似然估计就是数据的平均值：

mean(goT$Us)

## [1] 33.03

最大后验 (MAP) 估计可以通过最大化对数先验和对数似然之和来获得：

got.logposterior <- function(lambda) {
  dgamma(lambda, 180, 6, log = TRUE) + sum(dpois(GoT$Us, lambda,
    log = TRUE))
}

#Maximize log-posterior
got.MAP <- optim(30, got.logposterior, control = list(fnscale = -1))
got.MAP$par

鉴于伽玛分布与泊松似然共轭，后验以封闭形式提供，并且它是带有参数的伽玛$a+{\sum }_{i=1}^{31}{U}_i\mathrm{and}b+n$， 后验分布如图 1.2 所示。

​ 图 1.2：《权力的游戏》书中使用共轭先验（黑色实线）和 MAP 估计（红色虚线）的 u 数量的后验分布。

用于估计参数后验分布的 Metropolis-Hastings 算法$\lambda$在这种情况下可以实施。作为提案密度，我们将考虑以当前平均值为中心的对数正态分布$\lambda$小标准差等于0.05。

该算法会得到2100个$\lambda$的值，该值以20开始：

n.sim <- 2100
lambda <- rep(NA, n.sim)
lambda[1] <- 20

Metropolis-Hastings算法的实现如下：

set.seed(1)
for (i in 2:n.sim) {
  # Draw new value
  lambda.new <- rlnorm(1, log(lambda[i - 1]), 0.05)

  #Log-acceptance probability
  acc.prob <- dlnorm(lambda[i -1], log(lambda.new), 0.05, log = TRUE) + 
    dgamma(lambda.new, 180, 6, log = TRUE) +
    sum(dpois(GoT$Us, lambda.new, log = TRUE)) -
    dlnorm(lambda.new, log(lambda[i - 1]), 0.05, log = TRUE) -
    dgamma(lambda[i - 1], 180, 6, log = TRUE) -
    sum(dpois(GoT$Us, lambda[i - 1], log = TRUE))

  acc.prob <- min(1, exp(acc.prob))

  # Accept/reject new value  
  if(runif(1) < acc.prob) {
    lambda[i] <- lambda.new
  } else {
    lambda[i] <- lambda[i -1]
  }

}

图 1.3 显示了示例$\lambda$通过 Metropolis-Hastings 算法获得。第一次迭代提供的值太接近初始值，这些值通常被视为老化并被丢弃。其余值用于推理。然而，鉴于它们似乎有点相关，通常会对它们进行稀疏化。在这种情况下，100次迭代将被视为老化，并对剩余的观测值进行细化，以便仅保留五分之一。结果，将有 400 个值用于推理。

删除老化迭代和细化可以按如下方式完成：

lambda2 <- lambda[seq(101, n.sim, by = 5)]

后验分布的汇总统计λ通过对采样值进行汇总统计可以得到：

summary(lambda2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    30.0    32.0    32.7    32.6    33.2    36.2

核密度可用于获得后验边缘分布的估计λ从最终的一组值。这如图 1.3 底部的图表所示。值得一提的是，使用 Metropolis-Hastings 算法获得的估计值与迄今为止计算的其他方法获得的估计值非常接近，因为我们实际上拟合的是相同的模型。任何差异都是由于 MCMC 方法固有的蒙特卡罗误差造成的，但通常可以通过增加抽取的样本数量来减少该误差。无论如何，我们都会对本节末尾获得的不同估计进行比较和总结。

​ 图1.3：$\lambda$的样本通过Metropolis-Hasting算法(上)和使用样本的后延分布的核密度估计(下)获得

最后，我们展示如何使用 INLA 拟合这些模型。首先，模型拟合如下：
$$\begin{gathered} U_i\sim Po(\lambda ),i=1,\dots ,31 \\ \log (\lambda )=\alpha \\ \pi (\alpha )\propto 1 \end{gathered}$$
参数α是线性预测器的截距，并为其指定一个常数先验。这是一个不正确的先验，因为它不能整合到一个先验中，但有时这些对于提供模糊的先验很有用。注意λ不是直接建模，而是通过线性预测器的截距以对数尺度建模α。该模型可以使用以下代码通过 INLA 来实现（这里不会详细讨论）：

GoT.inla <- inla(Us ~ 1, data = GoT, family = "poisson",
  control.predictor = list(compute = TRUE)
) 

$\alpha$的后验分布的汇总统计是：

GoT.inla$summary.fixed

##              mean      sd 0.025quant 0.5quant 0.975quant  mode       kld
## (Intercept) 3.498 0.03127      3.436    3.498      3.558 3.498 4.908e-07

$\lambda$的后验分布可以通过将$\alpha$的后验分布转换得到$\lambda =exp(\alpha )$:

marg.lambda <- inla.tmarginal(exp, GoT.inla$marginals.fixed[[1]])
inla.est <- inla.zmarginal(marg.lambda)

## Mean            33.049 
## Stdev           1.02682 
## Quantile  0.025 31.059 
## Quantile  0.25  32.3464 
## Quantile  0.5   33.0376 
## Quantile  0.75  33.7378 
## Quantile  0.975 35.0904

为了比较不同的方法，我们制作了表 1.2，其中包括基于迄今为止描述的方法的估计。在这种情况下，λ参数的点估计非常相似。另请注意，最大似然估计非常接近使用所有贝叶斯方法获得的后验模式。后验标准差的点估计在计算它们的三种贝叶斯方法中也非常相似。计算出的分位数也非常接近。值得一提的是，在这种特殊情况下，模式很简单，我们有足够的数据来估计模型中的唯一参数。因此，预计所有结果都应该相似。

表 1.2： 使用不同贝叶斯估计方法的估计摘要。

Method	Mean	Mode	St. dev.	0.025quant.	0.975quant.
Max. lik.	NA	33.03	NA	NA	NA
MAP	NA	32.51	NA	NA	NA
Conjugate	32.45	32.51	0.9378	30.73	34.40
M-H	32.47	32.77	1.2625	30.55	34.33
INLA	33.05	33.02	1.0268	31.06	35.09

表：(#tab:GoTestimates) 使用不同贝叶斯估计方法的估计摘要。

图1.4展示了参数λ的后验边缘估计和MAP估计的不同估计使用此处介绍的不同方法获得。正如我们之前所看到的，所有方法都提供非常相似的后边缘估计。请注意，对于更复杂的模型，可能根本不可能获得封闭形式的边际，并且估计通常依赖于 MCMC 或 INLA 方法。

图 1.4：λ后边缘的不同估计的比较由本章描述的方法提供：MCMC（实线）、共轭先验（虚线）、INLA（虚线）和MAP（垂直红色虚线）。

1.8 final remarks

正如本章开头所述，目的是介绍贝叶斯推理。特别是，提供一些背景信息，以便读者可以将集成的嵌套拉普拉斯近似放置在贝叶斯生态系统中。出于多种原因，这一点很重要。首先，INLA 是进行贝叶斯推理的众多方法之一。 INLA 在许多情况下可能是正确的方法，但在其他一些情况下可能只是错误的方法。其次，这也意味着深入了解 INLA 有助于识别哪些模型表现良好，哪些模型表现不佳。

第 2 章介绍了集成嵌套拉普拉斯近似方法以及 INLA 包在贝叶斯模型拟合中的使用。第 3 章描述了模型构建的许多潜在影响。接下来，第 4 章讨论了多级模型。第 5 章总结了先验分布。第 6 章详细介绍了可用于拟合更复杂模型的其他特征。第 7 章描述了具有空间依赖性的模型类型。第 8 章介绍了通过时间相关模型和时空模型对时间序列进行分析。第 9 章介绍了样条曲线和其他平滑方法。时间序列的分析第 10 章检查了使用生存和联合模型的事件数据。第 11 章说明了拟合 INLA 包中未实现的新模型的不同方法。第 12 章描述了处理响应和协变量中的缺失值。最后，第 13 章解释了使用 INLA 进行混合模型的分析。