9Bayesian Computation.md

贝叶斯计算

简介

贝叶斯计算具体围绕两个步骤展开，分别为：计算后验分布$p(θ| y)$和计算后验预测分布$p(\tilde{y} | y)$。

阅读相关材料可知，对于一些简单的案例，我们可以以闭合形式进行分析计算；或直接使用标准分布（正态分布、伽马分布、贝塔分布、泊松分布等）的预编程案例和相应的数值进行计算。然而，对于复杂或不寻常的模型或高维模型，往往需要更复杂的算法来近似后验分布，通常通过组合不同的算法才能实现最高效的计算。使用Bayes定理时，我们经常使用乘积$p(θ)p(y |θ)$，即“先验乘以似然”，它与后验密度成正比。本文简要概述这些近似计算积分的原理与过程。

知识补充：后验预测分布与后验分布的区别与联系

1. 后验预测分布$p(\tilde{y} | y)$，指的是在给定已知样本数据下，预测新数据所使用的分布。因此，这个分布是关于数据，而非参数的分布。通常运用于使用机器学习方法进行数据预测的场景。在使用机器学习方法进行预测的过程中，我们会使用训练集拟合一个模型，并利用模型去预测测试集中的数据。而在这里，后验预测分布实际上就是我们结合先验分布$p(θ)$与已知数据$y$得到的一个模型。
2. 后验分布$p(θ|y)$或$p_{posterior}(θ)$，指的是我们用已有数据对参数先验分布进行更新得到的分布。

虽然二者的含义完全不同，但他们之间存在着紧密的联系。可以说，后验预测分布正是通过后验分布得到的。我们发现，通过全概率公式，我们可以得到： $$p(\tilde{y} | y)=\int_{-∞}^{+∞}p(\tilde{y} | y,θ)p_{posterior}(θ)dθ$$ 即： $$p(\tilde{y} | y)=\int_{-∞}^{+∞}p(\tilde{y} | y,θ)p(θ|y)dθ$$ 在上式中，如果我们假设训练集的数据与测试集的数据在给定参数 $θ$ 下独立（实际上，在大多场景下这个假设都成立），那么我们可以整理上式得到: $$p(\tilde{y} | y)=\int_{-∞}^{+∞}p(\tilde{y} | θ)p_{posterior}(θ)dθ$$ 即： $$p(\tilde{y} | y)=\int_{-∞}^{+∞}p(\tilde{y} | θ)p(θ|y)dθ$$ 因此，在这样的前提下，我们成功将后验预测分布的形式转化成了似然函数与后验分布乘积的积分。

在这里，当我们仔细审视这个推导过程时可以发现，使用 $θ$ 的前验分布来推导关于 $y$ 的预测分布函数时实际上也能够得到一个预测分布。我们通常称这个预测分布为先验预测分布。

事实上，我们在贝叶斯统计中并不一定需要严格区分前验分布与后验分布，在对参数 $θ$ 的分布进行多次更新的过程中，这一轮更新的后验分布总会成为下一轮更新的先验分布。因此，先验分布与后验分布总是相对来说的。

---

下面介绍一些近似计算积分的原理和过程：

1.数值积分（Numerical integration)

数值积分，用于求定积分的近似值，也称为“求积”，是指通过计算有限个点的函数值来计算连续函数上的积分的方法。在数值分析中，数值积分是计算定积分数值的方法和理论。在数学分析中，给定函数的定积分的计算不总是可行的，许多定积分不能用已知的积分公式得到精确值。

数值积分是利用黎曼积分等数学定义，用数值逼近的方法近似计算给定的定积分值。借助于电子计算设备，数值积分可以快速而有效地计算复杂的积分。 通过增加评估函数的点数，可以获得所需的精度。数值积分方法可分为模拟（随机）方法（如蒙特卡罗）和确定性函数方法（如许多求积规则方法）。

（1）模拟方法（Simulation methods）：

模拟（随机）方法是指从期望分布 $p(θ)$ 中获得随机样本 $θ^s$ 并估计任何函数 $h(θ)$ 的期望。 $$E(h(θ)| y)=\int h(θ)p(θ|y)dθ≈ \frac{1}{S}\sum_{s=1}^{S} h(θ^s)$$ 估计是随机的，取决于生成的随机数，可以通过获取更多样本来提高模拟的准确性。产生独立样本的方法有基本蒙特卡罗方法，马尔可夫链蒙特卡罗方法等，马尔可夫链蒙特卡罗方法在使贝叶斯推理适用于通用分层模型方面非常重要。模拟方法可用于高维分布，并且有适用于各种模型的通用算法；在必要时，通过将这些一般性的思想与特定的模拟方法、确定性算法及分布近似相结合，可以获得更有效的计算。

（2）确定性函数方法（Deterministic methods）：

确定性数值积分方法基于评估选定点 $θ^s$ 处的被积函数 $h(θ)p(θ|y)$ ，即上文模拟方法公式的加权化处理： $$E(h(θ)| y)=\int h(θ)p(θ|y)dθ≈ \frac{1}{S}\sum_{s=1}^{S} w^sh(θ^s)p(θ^s|y)$$ 权重 $w_s$ 对应于点 $θ^s$ 表示的空间体积。更复杂的规则，例如 Simpson 规则，则使用局部多项式来提高准确性。确定性数值积分规则通常比模拟方法具有更低的方差，但在高维中位置选择变得困难。

最简单的确定性方法是在具有相等权重的网格中评估被积函数。网格方法（ Grid methods） 可以自适应地从后验模式开始网格形成。对于一个部分具有某种特定形式的高斯被积函数，存在特定的求积规则可以用更少的被积函数评估给出更准确的估计。其有界和无界区域都存在正交规则。

2.分布近似（Distributional approximations）

分布（解析）近似，其用一些更简单的参数分布近似后验，可以直接计算积分；也可以使用分布近似作为模拟方法的起点。可以称之为“忽略一些信息的粗略估计”。

在进行更精细的近似或采用其他复杂方法之前，使用一些简单的非迭代方法获得目标分布位置的粗略估计（即模型中参数的点估计）是很有用的。丢弃掉模型和数据的一部分，然后创建一个简单的问题，可以方便地找到参数估计值。例如在分层模型中，可以先粗略估计超参数 $φ$，然后使用 $γ|φ, y$ 的条件后验分布来粗略估计主要参数 $γ$。

当某些数据丢失时，分布近似是补齐数据的一种好的方法，即根据可用数据简单地估算缺失值。这种粗略估计通常既方便又可靠，因为它们可以直接通过使用计算机相应的程序进行计算。

3.直接模拟法和拒绝采样法（Direct simulation and rejection sampling）

在简单的非分层贝叶斯模型中，相应参数通常很容易直接从后验分布中得出，尤其是在假设共轭先验分布的情况下。而对于更复杂的问题，直接模拟法和拒绝采样法可以帮助我们分析、分解分布并对这些分布进行部分模拟。首先从超参数的边缘后验分布中采样，然后根据数据和模拟的超参数绘制其他参数。很多时候我们可以对“较大”问题的“部分”执行直接模拟和分析积分方法，

（1）通过在点网格处计算直接逼近（Direct approximation by calculating at a grid of points）

对于最简单的离散近似，首先计算目标密度函数 $p(θ|y)$，在一组均匀间隔中的值 $θ_1, \ldots , θ_N $, 其涵盖了 $θ$ 的广泛参数空间。然后通过$θ_1, \ldots , θ_N $,处的离散密度估计近似的连续密度$p(θ|y)$，其概率为 $p(θ_i|y)/ \sum_{j=1}^{N} p(θj |y)$。可以看出，无论如何都必须对近似密度进行归一化，所以使用未归一化的密度函数 $q(θ|y)$ 代替 $ p(θ|y)$ 也能正常计算。

(2)从预测分布进行模拟(Simulating from predictive distributions)

一旦我们从后验分布 $p(θ|y)$ 中获得样本，就很容易得到未观察到的或预测的未来数据的分布 $\tilde{y}$ 。 对于从后验分布中每次抽取 $θ$ ，只需从预测分布 $p(\tilde{y}|θ)$ 中抽取一个$\tilde{y}$。 来自所有 θ 的一组模拟 $\tilde{y}$ 表征了后验预测分布。

(3)拒绝采样(Rejection sampling)

假设我们想获得概率密度 $p(θ|y)$ 或一个未归一化的密度 $q(θ|y)$（其中 $p(θ|y) = q(θ|y)/ \int q(θ |y)dθ)$。 在下面的描述中，我们使用 $p$ 来表示目标分布，同样我们也可以使用非标准化形式 $q$ 来代替。 为了执行拒绝抽样，我们需要为所有 $θ$ 定义一个正函数(positive function) $g(θ)$，其中 $p(θ|y) > 0$。 这样的形式具有以下属性特点：

• 我们可以从与 $g$ 成比例的概率密度中得出。 不要求 $g(θ)$ 积分为 1，但 $g(θ)$ 必须有一个有限积分。
• 重要性比（importance ratio） $\frac{p(θ|y)}{ g(θ)}$ 必须有一个已知的界限； 也就是说，必须存在已知的常数 $M$，对于所有的 $θ$，$\frac{p(θ|y)}{ g(θ)} ≤ M$。

拒绝采样算法分两步进行：

1. 从与 $g(θ)$ 成比例的概率密度中随机采样 $θ$。
2. 以概率$\frac{p(θ|y)}{ Mg(θ)}$ ，接受 $θ$ 作为 $p$ 的一个抽样。 如果描述的 $θ$ 被拒绝，则返回步骤 1。

用于拒绝采样的近似密度 $g(θ)$ 应大致与 $p(θ|y)$ 成正比（可将其视为 $θ$ 的函数）。 较为理想的情况是 $g ∝ p$，在这种情况下，如果 $M$ 的值合适，我们可以认为在第 1 步中每一次抽取的概率为 1，其也将在第 2 步中被拒绝。拒绝抽样的一个优点是它是 自我监督（self-monitoring） 的——如果该方法不能有效地工作，则很少会接受模拟抽取。

使用符号选择函数 $g(θ)$ 来逼近 $p(θ|y)$时，通常取决于 $y$。 然而，当我们使用 $g(θ, y)$ 或 $g(θ|y)$时，在具体案例中往往需要每一次考虑一个后验分布的近似值，其中 $g$ 对 $y$ 的函数依赖性并不重要。

拒绝抽样一般用于从标准单变量分布中进行抽样，是一种极其方便与快截的方法。 如果截断部分中密度质量的比例不接近 1，它也可以用于常规的截断多元分布（ truncated multivariate distributions）。

4.重要性采样（Importance sampling）

重要性采样是一种与拒绝抽样相关的方法，也是 Metropolis 算法的前身，用于对从目标分布的近似值中抽取的随机样本计算期望值。

如果 $g(θ)$ 是我们随机抽取的概率密度，那么我们可以写出以下公式，

$$ E(h(θ|y)) =\frac{\int h(θ)q(θ|y)dθ} {\int q(θ|y) )dθ} = \frac{\int [h(θ)q(θ|y)/g(θ)] g(θ)dθ} {\int [q(θ|y)/g(θ)] g(θ)dθ} $$

其中可以使用估计值$S$绘制 $θ_1, \ldots, θ_S$ ，从 $g(θ)$ 表达式的表达式可得，

$$ \frac{\frac{1}{S} \sum_{s=1}^{S}h(θ^s)w(θ^s)} {\frac{1}{S} \sum_{s=1}^{S} w(θ^s)} $$

其中因子：

$$ w(θ^s ) = \frac{q(θ^s|y)}{ g (θ^s)}$$ 称为重要性比率或重要性权重。 回想一下，$q$是我们对非归一化密度的一般表示法； 也就是说，$q(θ|y)$ 等于 $p(θ|y)$ 乘以某个不依赖于 $θ$ 的因子。

通常建议对上面表达式中的分子和分母使用相同的随机抽取集，以减少估计中的抽样误差。如果可以选择 $g(θ)$ 使得$ \frac{hq} {g}$ 大致恒定，则可以获得相当精确的积分估计。如果重要性比率变化很大，那么重要性采用就不是一种有用的方法。最糟糕的情况发生在重要性比率很小的概率很高，而概率很低时又很大，比如：如果 $hq$ 与 $g$ 相比具有宽尾，作为$θ$的函数，就会发生这样的情况。

重要性抽样估计的准确性和效率:

一般来说，如果没有对精确和近似密度进行某种形式的数学分析，我们总是有可能遗漏一些非常大但很少见的重要性的权重。 但是，检查采样重要性权重的分布可能有助于发现可能的问题。 它可以帮助我们检查最大重要性比率的对数的直方图：如果最大比率相对于平均值太大，估计通常会很差。 相比之下，我们不必担心重要性比小的行为，因为它们对上面方程的影响很小。

---

计算环境与应用

在现存的科学研究中，存在许多用于常用模型（例如分层线性和逻辑回归）以及一些非参数模型的完全贝叶斯推理的程序。这些实现使用了其余同学们编写的贝叶斯计算算法的各种组合。

我们看到（至少）四个原因需要一个用于拟合贝叶斯模型的自动化通用程序。

• 首先，许多应用统计学家和主题研究人员想要拟合贝叶斯模型，但不具备数学、统计和计算机技能来对推理步骤进行编程。经济学家和政治学家可以通过一行代码或单击菜单进行回归，流行病学家可以进行逻辑回归，社会学家可以拟合结构方程模型，心理学家可以拟合方差分析，教育研究人员可以拟合分层线性模型。我们希望所有这些人都能够拟合贝叶斯模型（包括前面提到的所有特殊情况。但也允许泛化，例如鲁棒误差模型、混合分布和各种任意函数形式，更不用说框架了 包括先验信息）。

• 一般贝叶斯程序的第二个用途是用于教学。 学生可以先运用程序自动推理，专注于模型的结构而不是计算，之后再学习代数和计算。当然，学习代数及计算方面更深入的理解是有用的，最终了解推理和计算背后的原因也是有帮助的，因为我们理解模型的一种方式是通过与稍微简单或更复杂的类似模型进行比较，而我们理解模型拟合过程的一种方式是将其视为一个从数据和假设映射到推论。 即使在使用黑匣子或“推理引擎”（黑箱理论）时，我们也经常想回过头来看看我们的后验分布的实质性重要特征是从哪里来的。

• 编写自动化通用程序的第三个动机是作为一个用于实现新模型的编程环境，以及让用户对自己的模型进行编程，从而能够专注于更重要的统计问题。

• 最后，通用贝叶斯程序的第四个潜在好处是它可以比自定义代码更快。用经济学术语讲，即为存在规模经济。由于自动程序将被多次使用，我们便可以通过各种方式对其进行优化、并行实施，并包含需要更多编码工作且在困难问题上运行速度更快的算法。

总的来说，没有任何程序可以真正通用。任何这样的软件都应该是开放的和可访问的（开源的），然后存在一些地方让骨灰级用户（大佬）可以改变、更新程序或 “握住它的手”（hold its hand ，以确保程序能够做它应该做的事情。