【嵌入式拉普拉斯近似】Integrated Nested Laplace Approximation

Bayesian Inference with INLA

积分嵌套拉普拉斯近似（Integrated Nested Laplace approximation）是一种近似贝叶斯推理方法。

Introduction to Bayesian Inference

近年来，由于计算统计学的重要进步，贝叶斯推理得到了发展。本文章将重点介绍用于近似贝叶斯推理的集成嵌套拉普拉斯近似。INLA是贝叶斯统计学中最近的几个计算突破之一，它允许快速准确的模型拟合。

贝叶斯定理介绍

在贝叶斯范式中，模型中的所有未知量都被视为随机变量，其目的是计算（或估计）联合后验分布。这就是参数θ的分布，以观测数据y为条件。获得后验分布的方法依赖于贝叶斯定理：
$$P(\theta |x)=\frac{P(x|\theta )\cdot P(\theta )}{{\int }_{\theta }P(x|\theta )\cdot P(\theta )d\theta }$$

• 等式左边项$P(\theta |X)$表示后验概率，其中$P(\theta |X)\propto P(\theta )\cdot P(X|\theta )$
• $P(x|\theta )$表示似然
• $P(\theta )$表示先验

贝叶斯定理是统计的基础，提到贝叶斯统计，一个自然地问题就是先验分布该如何选择（问不同的人会给出不同的答案）。贝叶斯统计的一个缺点就是贝叶斯公式比频率派的参数估计更加复杂，那么如何选择合适的先验分布让贝叶斯统计变得更加容易计算就成为了一个非常现实的问题。
为什么说贝叶斯统计是难以计算的呢？一个最直接的困难就是我们难以保证后验分布具有解析解，这里的后验分布具有解析指的是后验分布的密度函数具有解析解。为什么后验分布有时会不具备解析解？因为在贝叶斯公式中分母是一个积分，那么显然就不一定具有解析解，一旦分母没有解析解后验分布就不一定具有解析解。若想让这个分母的积分具有解析解，由于分母的被积函数是模型和先验分布相乘，因此在这里你可能会简单的认为只要选取一个合适的模型和先验分布就能让后验分布具有解析解。但是，因为贝叶斯统计是不断收集数据和并更新参数分布的迭代过程。利用贝叶斯计算后验分布在更新的过程中，每次都会把上次迭代的后延分布设为新一轮迭代的先验分布。所以为了确保每次迭代的后验分布都能够得到一个解析解，为达此目的，须让先验分布和后延分布同属于一个参数化分布族，也就是让先验分布和后验分布的表达式形式相同。
注意，后验分布的形式是有先验分布和模型（似然函数）共同决定的。因此，需要同时选取特定的先验分布和模型才能保证后验分布和先验分布相同的形式。自然而然的我们引出共轭先验的定义：
共轭先验指的是一个分布族为模型$X\sim f_{x|\theta }(x|\theta ),\theta \in \Omega$的共轭先验，若只要先验分布$f_{\theta }(\theta )$是从该分布族中选取的，最终得到的后验分布$f_{\theta |x}(\theta |x)$也就属于该分布族。

共轭先验——Conjugate priors

如上所述，后验分布仅以闭合形式用于少数模型。具有共轭先验的模型是那些先验与似然具有相同形式的模型。例如，如果似然是具有已知精度的高斯分布，则均值上的共轭先验是高斯分布。这将确保平均值的后验分布也是高斯分布。
例如，我们有一系列的观测值 { y i } i = 1 n \{y_i\}_{i=1}^n {yi​}i=1n​服从高斯分布
$$y_i|\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)$其中
$$\begin{gathered} {\mu }_1={\mu }_0\frac{{\tau }_0}{{\tau }_0+\tau n}+\bar{y}\frac{\tau n}{{\tau }_0+\tau n} \\ {\tau }_1={\tau }_0+n\tau \end{gathered}$$
这很好地说明了贝叶斯定理。首先，后验平均值 是先验均值之间的折衷${\mu }_0$和观察到的平均值数据$\bar{y}$。当观测值数量较大时，则 后验平均值接近数据。如果 n 小，权重大鉴于我们先前对均值的信念。同样，后验精度为先验精度和似然精度的函数，它倾向于 无穷大与 n ，这意味着$\mu$的方差趋于零，因为 数据数量增加。

引用别人的博客的解释
在数学上，若先验概率与似然概率共轭，那么后验分布概率$p(w|t,X)$
与先验概率服从相同的分布。比如说：先验概率服从Beta分布，似然概率服从二项分布，这时先验概率分布与似然概率分布共轭了，那么后验概率也服从Beta分布。

计算方法

当后验分布不能求得一个解析解时，有必要诉诸其他方法来估计它（或者叫近似他），或者，从中抽取样本。给定来自后验分布的样本并遍历样本。通常，计算方法旨在估计出现的积分 在贝叶斯推理中。例如，参数的后验平均值${\theta }_i$（参数空间中的值${\theta }_i$） 的计算公式为:
$${\int }_{{\theta }_i}{\theta }_i\pi ({\theta }_i|y)d{\theta }_i$$
$\pi ({\theta }_i|y)$是参数${\theta }_i$的后验边际分布，类似地，也可以计算后验方差或者任何其他时刻。可以使用数值积分的方法和拉普拉斯近似法方便地近似这种类型的积分。

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

可以通过最大化似然与先验的乘积来获得模型参数的点估计。 这由最大后验 (MAP) 估计表示。 后验分布的最大化通常在对数缩放中执行，并且可以使用不同的方法有效地实现，例如 Newton-Raphson 算法或期望最大化算法 (Gelman et al. 2013)。

The integrated nested Laplace approximation拉普拉斯近似

在贝叶斯定理中我们介绍到由于无法直接求得后验概率的解析解，但是我们知道后验概率正比于似然函数乘以先验也就是$P(\theta |X)\propto P(\theta )\cdot P(X|\theta )$。于是乎我们就可以转而求解$P(\theta )\cdot P(X|\theta )$，在这里我们假设待求解的函数为$g(\theta ;X,t,{\sigma }^2)$。在这里，拉普拉斯近似就是首先假设函数$g$服从高斯分布，然后通过泰勒展开公式，将$log(g(\omega ;X,t,{\sigma }^2))$在${\omega }^{\ast }$处展开。

让我们回想一下基本的泰勒级数展开，其中一个函数关于一个点可以展开为项的和(有时是无限的)，使用有限个这些项可以作为近似值(通常使用前三项，直到二阶导数)。
$$f(x)=f(x_0)+\frac{f^{\prime }(x_0)}{1!}(x-x_0)+\frac{f^{\prime \prime }(x_0)}{2!}(x-x_0{)}^2+\frac{f^{\prime \prime \prime }(x-x_0}{3!}(x-x_0{)}^3+\dots \text{\hspace{1em}(2)}$$

根据公式（2），对函数$g$进行二阶泰勒展开如下：
$$\log g(\omega ;X,t,{\sigma }^2)\approx \log g(\hat{\omega };X,t,{\sigma }^2)+\frac{\partial \log g(\omega ;X,t,{\sigma }^2)}{\partial {\omega }^2}{\mid }_{\hat{\omega }}\frac{(\omega -\hat{\omega })}{1!}+\frac{{\partial }^2\log g(\omega ;X,t,{\sigma }^2)}{\partial {\omega }^2}{\mid }_{\hat{\omega }}\frac{(\omega -\hat{\omega }{)}^2}{2!}+\dots \text{\hspace{1em}(3)}$$
该模式的一阶导数为0，将公式（3）化简为：
$$\log g(\omega ;X,t,{\sigma }^2)\approx \log g(\hat{\omega };X,t,{\sigma }^2)-\frac{v}{2}(\omega -\hat{\omega }{)}^2$$
接着我们把表达式重写为：
$$\log g(x)\approx \log g(\hat{x})-\frac{1}{2{\sigma }^2}(x-\hat{x}{)}^2$$
使用拉普拉斯近似法，分布f(x)近似为均值$\hat{x}$的正态分布，可以通过求解 $f^{\prime }(x)=0$，方差 ${\sigma }^2=-1/f^{\prime \prime }(x)$。

参考链接

1.http://www.manongjc.com/detail/64-qmnljjsuobazkgc.html
2.https://ourcodingclub.github.io/tutorials/mcmcglmm/index.html
3.https://www.precision-analytics.ca/articles/a-gentle-inla-tutorial/
4.https://becarioprecario.bitbucket.io/inla-gitbook/index.html
5.https://www.cnblogs.com/hapjin/p/8848480.html