白家欣 2019302130057
变分贝叶斯是一种类似于EM的算法框架，用于逼近联合分布。EM通过交替评估对数密度的条件期望并使用这些来最大化一组超参数的函数（进而定义用于计算下一步期望的条件分布等等），收敛到超参数的点估计，从而逼近后验分布。在变分贝叶斯中，迭代导致封闭形式的近似值，该近似值与某些指定函数类中的后验分布（在某种意义上定义如下）最接近。
在变分贝叶斯中，参数逼近g(θ)是使用期望过程迭代构建的，正如我们将展示的，它具有最小化目标后验分布p(θ|y)的Kullback-Leibler 散度的效果，
$$
\mathrm{KL}(g | p)=-\mathrm{E}_{g}\left(\log \left(\frac{p(\theta \mid y)}{g(\theta)}\right)\right)=-\int \log \left(\frac{p(\theta \mid y)}{g(\theta)}\right) g(\theta) d \theta
$$(13.15)
这个散度的绝对最小值是0，这是在g≡p时得到的。困难在于，变分贝叶斯通常应用于我们不能直接总结p的情况下，也就是说，我们不能轻易地从p(θ|y)取后验，也不能轻易地计算兴趣期望Ep(h(θ))。在变分贝叶斯中，我们使用一些更简单的参数化分布g，它们更容易处理。
这里我们将用φ符号表示变分近似的超参数。因此，我们将逼近函数g写成g(θ|φ)，算法从φ的一些猜想开始，然后以数学上保证每一步减小Kullback-Leibler散度(13.15)的方式迭代更新它。与EM一样，在某一时刻φ不再产生任何可见的变化。
在这一点上，我们停止迭代，并使用g(θ|φ)给出最新更新的φ作为后验密度的近似。从不同的起点多次运行算法来检查结果是有意义的。
有多种方法可以定义变分近似g(θ|φ)的分布类别。标准方法是将θ的分量约束为独立的；因此，
$$
g(\theta \mid \phi)=\prod_{j=1}^{J} g_{j}\left(\theta_{j} \mid \phi_{j}\right)
$$
对于J维参数θ。在这种情况下，可以根据后验密度函数 p(θ|y) 的数学形式确定要优化的分布族g。
它是这样工作的:对于每个j，我们检查对数后验密度的期望，logp(θ|y)，把它看作θj的函数，对表示θ的其他j−1维的g−j的分布求平均值。
这与吉布斯抽样相似，只是我们感兴趣的是平均密度而不是条件密度。在建立变分贝叶斯算法的这一点上，我们还不需要评估期望Eg−j(logp(θ|y));我们只需要算出它的数学形式作为θj的函数。一旦我们计算出每个参数θj ，我们确定了近似分布的函数形式gj(θj|φj).
与EM一样，变分贝叶斯最适用于具有条件共轭先验分布的指数族模型，在这种情况下，近似分布族通常可以通过检验确定，必要的期望可以以封闭形式计算。在非共轭模型中，变分贝叶斯仍然可以通过在受限的函数形式(如正态分布)中工作来实现。
一旦近似分布gj(θj|φj)，计算从所有超参数φ的猜测开始。然后循环遍历分布gj，在每个步骤中更新超参数φj的向量，使
$$
\mathrm{E}{g-j}(\log p(\theta \mid y))=\int \log p(\theta \mid y) g{-j}\left(\theta_{-j} \mid \phi_{-j}\right) d_{\theta-j}
$$。
我们用符号表示Eg−j 表示除θj外所有参数的近似分布的平均值，条件为g−j的当前迭代。结果是一个J−1维积分，但对于许多模型，我们能够分析评估这些期望。
我们将简要证明变分贝叶斯的步骤降低KL(g||p)，从而逐渐使近似分布g(θ)更接近目标后验分布p(θ|y)。但首先，我们用一个简单但不平凡的例子来说明算法是如何工作的。
我们用5.5节中8所学校的层次模型来说明。与我们在12.5节中对HMC的演示一样，我们将八种学校效应（在第5章中定义为θj）标记为αj；参数θ的完整向量则有10个维度，对应于α1, ..., α8,µ,τ，对数后验密度为
$$
\log p(\theta \mid y)=-\frac{1}{2} \sum_{j=1}^{8} \frac{\left(y_{j}-\alpha_{j}\right)^{2}}{\sigma_{j}^{2}}-8 \log \tau-\frac{1}{2} \frac{1}{\tau^{2}} \sum_{j=1}^{8}\left(\alpha_{j}-\mu\right)^{2}+\text { const }
$$
(13.16)
我们将遵循变分贝叶斯的标准实践,并通过独立密度的乘积来近似 p(θ)；因此，
$$
g(\theta)=g\left(\alpha_{1}, \ldots, \alpha_{8}, \mu, \tau\right)=g\left(\alpha_{1}\right) \cdots g\left(\alpha_{8}\right) g(\mu) g(\tau)
$$
(13.17)
确定近似分布的形式。 我们开始变分贝叶斯,如下,对于模型中的十个变量中的每一个，我们检查对数后验密度，考虑其对其他九个变量的独立分布 g 的期望平均，并确定其参数形式：
•对于每个αj，我们查看Elogp，对其他七个α's、μ和τ求平均值；即，对(13.17)的除g(αj)之外的所有因子求平均 (13.16)。结果是αj的二次函数。对于这个检查，平均分布g的细节无关紧要；知道它们是独立的就足够了。我们需要做的就是查看（13.16）并考虑如果我们对除αj之外的所有参数求平均值会发生什么。结果是，
$$
\mathrm{E} \log p(\theta \mid y)=-\frac{1}{2} \frac{\left(y_{j}-\alpha_{j}\right)^{2}}{\sigma_{j}^{2}}-\frac{1}{2} \mathrm{E}\left(\frac{1}{\tau^{2}}\right) \mathrm{E}\left(\left(\alpha_{j}-\mu\right)^{2}\right)+\text { constant. }
$$
不同的期望不涉及αj 可以被带入常数项。更明确地说，我们可以把期望写成Eg−αj Log p，表示在(13.17)中除gαj外的所有g因子的平均值.在任何情况下，我们都知道 是α的二次函数j因此e是正比于法向密度当考虑作为g的函数j.
我们可以通过完成二次表达式的平方来确定这个正态分布的参数，或者从统计学的角度更直观地认为这个表达式等价于两条信息，一条以y为中心j与反变量σj−2和1个以E(E)为中心，方差为负的E(τ12)．我们通过加权均值和添加逆方差将这些组合起来，从而得到αj的变分贝叶斯分量的如下形式:
$$ g\left(\alpha_{j}\right)=\mathrm{N}\left(\alpha_{j} \mid \frac{\frac{1}{\sigma_{j}^{2}} y_{j}+\mathrm{E}\left(\frac{1}{\tau^{2}}\right) \mathrm{E}(\mu)}{\frac{1}{\sigma^{2}}+\mathrm{E}\left(\frac{1}{\tau^{2}}\right)}, \frac{1}{\frac{1}{\sigma_{3}^{2}}+\mathrm{E}\left(\frac{1}{\tau^{2}}\right)}\right) $$(13.18)
•对于µ，我们检测(13.16)。对除µ以外的所有参数求平均,表达式Elogp(θ|y)的形式为
$$
-\frac{1}{2} \mathrm{E}\left(\frac{1}{\tau^{2}}\right) \sum_{j=1}^{8}\left(\mathrm{E}\left(\alpha_{j}\right)-\mu\right)^{2}+\text { const. }
$$。如上文所述，这是一个正态密度函数的对数;这个分布的参数可以通过将其视为8条信息的组合来确定:
$$
g(\mu)=\mathrm{N}\left(\mu \mid \frac{1}{8} \sum_{j=1}^{8} \mathrm{E}\left(\alpha_{j}\right), \frac{1}{8} \frac{1}{\mathrm{E}\left(\frac{1}{\tau^{2}}\right)}\right)
$$(13.19)
•最后，对除τ以外的所有参数求平均给出了一个密度函数，该函数可以被识别为反伽马，或者在我们更喜欢的参数化中，在近似分布g上具有期望E.
$$
g\left(\tau^{2}\right)=\operatorname{Inv}-\chi^{2}\left(\tau^{2} \mid 7, \frac{1}{7} \sum_{j=1}^{8} \mathrm{E}\left(\left(\alpha_{j}-\mu\right)^{2}\right)\right)
$$(13.20)
上述表达式与第11.6节关于层次正态模型吉布斯采样器条件分布的推导和第13.6节EM算法的推导基本相同，唯一不同的是在8-schools的例子中，我们假设数据方差σj是已知的。
确定条件期望。 用通用符号重写上述因素，我们有:
$$
g\left(\alpha_{j}\right)=\mathrm{N}\left(\alpha_{j} \mid M_{\alpha_{j}}, S_{\alpha_{j}}^{2}\right), \text { for } j=1, \ldots, 8
$$(13.21)
$$
g(\mu)=N\left(\mu \mid M_{\mu}, S_{\mu}^{2}\right)
$$(13.22)
$$
g\left(\tau^{2}\right)=\operatorname{Inv}-\chi^{2}\left(\tau^{2} \mid 7, M_{\tau}^{2}\right)
$$(13.23)
我们需要这些来获得上述三个步骤的条件期望:
•为了确定(13.18)中αj的分布，我们需要Mµ来自(13.22)的E(µ)和1/Mτ2来自(13.23)的E(1/τ2)。
•为了确定(13.19)中µ的分布，我们需要Mαj来自(13.22)的E(αj)和1/Mτ2来自(13.23)的E(1/τ2)。
•为了确定(13.20)中τ的分布，我们需要E(αj −µ)2)，即(Mαj −Mµ)2+Sαj2+Sµ2 由(13.21)和(13.22)推导，并利用变分近似中密度g是独立的假设。
8所学校分层模型的变分 逼近参数的变分贝叶斯进展在一个随机的起点之后，参数需要大约50次迭代才能达到近似收敛。右下角的图显示了Kullback-Leibler散度KL(g||p)(计算到任意的加性常数);如果编制正确的变分算法，KL(g||p)能得到一致减小。
起始值 为了开始变分贝叶斯，我们需要初始化的不是变量α,µ,τ,而是它们分布中的参数g(α)，g(µ)，g(τ)。为简单起见，我们画出无界参数Mα1,...,Mα8,Mµ从独立的N(0,1)分布中得到有界参数Sα1,...,Sα8,Sµ从独立的U(0,1)分布。
运行算法 对α进行迭代α1,...,α8,µ,τ,在每次迭代时使用来自其他分布的当前值的期望更新分布(13.18)-(13.20)。也就是说，我们计算(13.18)-(13.20)中的参数，方法是使用M'm和S's的当前值，代入上面重点描述的期望。 然后我们将(13.18)-(13.20)中的新计算的均值和标准差标记为更新的M'm和S's。因此，该算法看起来很像EM，不同之处在于它是正在被更新的分布，而不是点估计。图13.5中的前5个图显示了分布(13.21)-(13.23)参数的进展。在这些特定的起点上，算法需要一段时间才能开始运行，但一旦它被解开，就会很快找到解决方案。
检查拟合是否在改善 如上所述，Kullback-Leibler散度(13.15)应该在变分Bayes的每一步中减小。在这个例子中，我们可以解析地计算这个表达式，我们这样做了。因为我们在比较这些值，而不关心它们的绝对水平，所以我们可以通过忽略不依赖于参数θ的常数来简化分析。
例如，Kullback-Leibler的分歧是,
图13.6对于100次变分贝叶斯迭代，A, B, C学校的效果推论进展线条和阴影区域显示了变异分布的中值、50%区间和90%区间。每个图的右侧显示的是通过模拟计算的完整贝叶斯推断的相应分位数。
图13.5右下角的图显示了随着算法的进行KL(g||p)的稳定下降。
比较变分贝叶斯解和全贝叶斯解。 图13.6给出了三个参数α的变分算法的进展α1,α2,α3，与数据中前三所学校的辅导计划的效果相对应。这里的函数形式是高斯函数，所以它必然无法捕捉后验分布的一些微妙之处，这可以从这个例子中与非对称全贝叶斯区间的比较中看出。此外，这种变分拟合不考虑α之间的依赖性j 因此在某些情况下是不合适的。也就是说，这种近似很好地拟合了边际分布，而变分贝叶斯代表了一种快速且可扩展的方法，用于对大数据集的这类问题进行推理。
不像MCMC，它最终收敛到后验密度，变分推理收敛到一个近似-在一个受限类内最接近的拟合。在这种情况下，我们也可以让MCMC运行足够长的时间以达到收敛，通过将其与实际后验密度进行比较来理解变分贝叶斯是有意义的。
Kullback-Leibler散度(13.15)是用后验密度p(θ|y)定义的。理解变分贝叶斯的第一步是用我们可以计算的非归一化密度p(θ，y) = p(θ)p(y|θ)重新表示。重新表达式是这样的:p(θ|y) = p(θ,y)/p(y),所以logp(θ|y)=logp(θ,y)−logp(y)，因此,
$$ \mathrm{KL}(g | p)=-\int \log \left(\frac{p(\theta \mid y)}{g(\theta)}\right) g(\theta) d \theta=-\int \log \left(\frac{p(\theta, y)}{g(\theta)}\right) g(\theta) d \theta+\int \log (p(y)) g(\theta) d \theta=-\mathrm{E}_{g}\left(\log \left(\frac{p(\theta, y)}{g(\theta)}\right)\right)+\log p(y) $$(13.24)
我们通常无法评估最后一项logp(y)，但为了变分推理的目的，我们需要意识到它不依赖于g。因此，对于任何给定的模型 p 和数据 y，减少 (13.24) 右边的第一项相当于减少 KL(g||p)。 这个表达式Eg(log(p( 被称为变分下界。
接下来，我们证明了变分算法的每一步都保证增加变分下界，或等效地减小全局Kullback-Leibler散度(13.24)。考虑分布gj(θj)正在更新。我们将用我们假设的近似分布的因子分解变分下界:g(θ)=gj(θj)g−j(θ−j)：
$$ \begin{aligned} -\mathrm{E}{g}\left(\log \left(\frac{p(\theta, y)}{g(\theta)}\right)\right)=& \iint g{j}\left(\theta_{j}\right) g_{-j}\left(\theta_{-j}\right)\left(-\log p(\theta, y)+\log g_{j}\left(\theta_{j}\right)+\log g_{-j}\left(\theta_{-j}\right)\right) d \theta_{j} d \theta_{-j} \ =&-\int g_{j}\left(\theta_{j}\right)\left(\int g_{-j}\left(\theta_{-j}\right) \log p(\theta, y) d \theta_{-j}\right) d \theta_{j}+\ & \quad+\int g_{j}\left(\theta_{j}\right) \log g_{j}\left(\theta_{j}\right) d \theta_{j}+\int g_{-j}\left(\theta_{-j}\right) \log g_{-j}\left(\theta_{-j}\right) d \theta_{-j} . \end{aligned} $$(13.25)
在上面最后一行中，我们能够将二重积分转化为单积分，因为我们假设gj和g-j是归一化概率密度，因此积分为1。
这里我们只考虑g的步长j Is被更新，所以我们可以忽略(13.25)中的最后一项，因为它只依赖于g−j ，我们可以认为第一项括号中的表达式(暂时)是一个常数，因为它不依赖于gj ：
$$ \mathrm{E}{g-j} \log p(\theta, y)=\int g{-j}\left(\theta_{-j}\right) \log p(\theta, y) d \theta_{-j} $$(13.26)
由于θ−j已被积分出来，式(13.26)仅是θj和y的函数，可以认为是θ非归一化概率密度的对数j ，我们称之为，
$$ \log \tilde{p}\left(\theta_{j}\right)=\mathrm{E}_{g-j} \log p(\theta, y)+\text { constant } $$(13.27)
即，logp是将logp(θ|y)考虑为θ的函数而得到的(非归一化)对数密度j把期望除以θ的所有其他分量，求g的当前迭代平均值−j (如本节前面8所学校的例子所详细说明的那样)。表达式(13.25)则变成,
$$ -\mathrm{E}{g}\left(\log \left(\frac{p(\theta, y)}{g(\theta)}\right)\right)=-\int g{j}\left(\theta_{j}\right) \log \left(\frac{\tilde{p}\left(\theta_{j}\right)}{g_{j}\left(\theta_{j}\right)}\right) d \theta_{j}+\text { const } $$(13.28)
最后一个表达式就是KL(gj||p~)， g的Kullback-Leibler散度j对于p~,当gj(θj)≡p˜(θj)，定义见(13.27)。
因此，当有可能评估(13.26)中的期望并由此确定g时j ，我们有我们的更新，并且变分算法保证在每一步减小全局Kullback-Leibler散度(13.15)。如果期望(13.26)不能以封闭的形式完成，则需要对g进行一些更新j 减少(13.28)，从而得到gj 在这一步中更接近p。
变分贝叶斯近似是一种生成模型，即参数θ的合适概率分布;因此，我们可以通过画一个样本θs来检验它的拟合性，s=1，...，S。然后，对于每个θs，绘制一个复制的数据集yrep,S。对于8个学校的例子，我们首先绘制1000个参数α的重复α1,…,α8 从我们变分贝叶斯计算的最终g，然后对8所学校的数据进行1000次重复采样。
或者，如果我们对新学校的预测感兴趣，我们将首先从g中画出1000个(µ，τ)的值，然后，对于每个模拟，从N(µ，τ2)画出8个新学校，然后从每个新的学派中得出一个新的观测数据(条件是假设的σ)。
一般来说，我们会期望，与完全贝叶斯相比，变分推理将提供更好的拟合观测数据。作为分布的点估计，变分贝叶斯可以过拟合数据。尽管如此，模型检查可能是一个好主意，因为它仍然可以揭示推论中的问题。
当基于模拟的全贝叶斯计算过于昂贵时，比如使用非常大的模型或数据集时，变分方法通常被用作一种近似方法。在这种情况下，使用变分估计作为随机算法的起点可能是有意义的，从而得到更好的目标分布的近似。
最简单的想法将重要性抽样:在各种问题的变分方法是驯良的,我们可以很容易地计算目标密度p(θ|y)和近似g(θ),我们也可以直接从g。我们可以快速模拟计算模拟了,从g和θ年代,对每个计算重要性权重,p(θs|y)/g(θs)。(通常，我们只需要这些权值达到一个任意的乘法常数;因此，使用非标准化密度来代替g或p是很好的)。然后我们可以使用加权平均计算期望。
不幸的是，从变分近似中直接使用重要性加权可能是灾难性的，因为变分贝叶斯拟合往往比目标分布的变量小，因此重要性比的分布可能有长尾，导致不稳定的平均值。因此，我们推荐重要性重采样，首先从g中提取，然后使用重要性权重进行不替换的重采样。(在不进行替换的情况下进行重新采样是至关重要的，这样任何具有极高权重的采样点都不会主导模拟。)在重要性重采样的一般情况下，在采样的两个阶段中，不总是清楚要抽取多少次。
更一般的方法是粒子滤波，同样使用变分贝叶斯的提取作为起点，然后使用Metropolis或哈密顿蒙特卡罗移动目标密度，并根据需要分裂和删除点。为任何特定的示例实现此方法可能需要大量的编程时间，但对于大问题，或者在已经设置了粒子过滤的计算环境中，这是有意义的。
变分推理以J步进行，每次更新一个条件分布gj EM有两个步骤(E-step和M-step)，交替估计参数φ和平均其他参数γ。EM可以被看作是一种特殊情况(a)的变分贝叶斯参数划分为两个部分,φ和γ(b)的近似分布φ需要一个质点(因此,更新g(φ)相当于更新φ的点估计值),和(c)γ的近似分布是不受约束的;因此g(γ)= p(γ|φ，y)，条件是φ的最新更新。
如上所述的变分推断是在某些近似分布g上近似目标分布p的一种方法，该方法使用迭代算法，每一步都减少Kullback-Leibler散度KL(g||p)。
如上所述，只有对于具有某些共轭性质的模型，这才能以封闭形式完成。这种模型包括许多重要的特殊情况(如正态分布和有限混合)，但更普遍的是，变分贝叶斯的思想可以通过用KL(g||p)的近似代替目标函数(要最小化的准则)来加以扩展。 对于许多问题，包括逻辑回归，都有良好的近似值，因此在这个新标准上进行优化的算法可以得到后验分布的良好近似值。通常，逼近本身随着每一步的变化而变化，它是根据逼近函数g的最新更新来定义的。