变分贝叶斯（Variational Bayes）

阅读量：5275 次

发布时间：2019-06-14

本文共 12916 字，大约阅读时间需要 43 分钟。

此文公式图片不全。详见博客：

【关键字】平均场理论，变分法，贝叶斯推断，EM算法，KL散度，变分估计，变分消息传递

引言

· 从贝叶斯推断说起

Question：如果我们有一组观测数据D，如何推断产生这些数据的模型m？

模型由1）模型的类别ξ（如高斯分布，伽马分布，多项式分布等）与2）模型的参数Θ共同决定，即 .

模型的选择

假设M为所有可能的模型集合（包括不同类别），那么选择

如何计算p(m | D)?
- 通常情况很难直接计算p(m | D)，根据贝叶斯公式有，p(m)表示模型的先验，p(D | m)表示证据；
- 先验：贝叶斯规则倾向于选择能解释数据的最简单模型：Occam剃刀原理。因为简单模型只在有限范围内做预测，复杂模型（如有更多自由参数）能对更宽范围做预测。

那么如何计算证据（evidence）？
- 参数θ的后验概率为
- 证据p(D | m)通常会在最可能的参数附近有一个很强的峰。
- 以一维参数为例：利用Laplace方法近似，即用被积函数乘以其宽度。即。
- 此处不在深究Occam因子。
- 从模型的选择可以看出参数的估计非常重要。

考虑同一个类别的模型。由于任何模型（函数）都可以由统一的数学形式给出，比如拉格朗日展开，傅里叶极数，高斯混合模型（GMM）等，因而通常我们更关心一个模型的参数Θ。换句话说，给出一组观测数据D，我们总是能够通过估计参数来推测模型，即。或者更简单的形式。

后验概率的估计

通常情况，取后验概率最大的参数值为估计值。根据贝叶斯公式，参数θ后验概率为

，

其中p(D)为归一化常数（normalizing constant）。

从经典的统计学角度看，概率是相对频率的，是真实世界的客观属性。因而每个模型被选择的概率是一样的，因而p(θ) =constant。此时问题转化为：，这便是极大似然法（ML, Maximum Likelihood）。

从贝叶斯学派的角度看，每一个模型都有一个先验概率p(θ)，但先验概率需事先给定。此时问题转化为：，这便是极大后验估计（MAP, Maximum A Posteriori）。

另一方面，许多科学问题的基本部分是计算一个目标函数的积分。Ω通常是高维空间中的一个区域，一般情况下f(x)稍微复杂一些，积分就难以计算。如果f(x)能被分解成一个函数g(x)与一个概率密度函数π(x)的乘积，那么上述积分可看做是g(x)在密度π(x)下的期望。

比如，

（1）计算后验概率：

（2）点估计：

（3）训练样本预测将来的数据的概率密度：假设D’与D条件独立， .

（4）新观测样本D’的隐藏变量(hidden variable) x’的后验分布：

上述积分最简单的近似方法就是通过估计参数θ来估计单点积分值，比如上述贝叶斯选择模型中极大后验估计（MAP）。由于ML、MAP只是估计概率密度而不是概率分布，因而省去了积分过程。ML, MAP估计最常用也最基本的方法是期望最大化算法（Expectation Maximization，EM）。

此外，可以通过蒙特卡洛方法（Monte Carlo），或马氏链蒙特卡洛法（Markov Chain Monte Carlo, MCMC）来模拟积分。此类方法具有较高的精度，但需要大量的计算。

本文介绍一种变分方法来近似积分。其主要思想是，对一个特定模型，构造一种简单（tractable）的数学形式来近似未观测变量的后验分布，同时给出观测数据的边缘似然（或者称证据，evidence）的下界（lower bound）。而积分过程转化为求下界的最优值问题。

· 传统方法的挑战

EM算法

EM算法主要是用于在不完全数据的情况下计算最大似然估计，是一个不断迭代优化的过程。以高斯混合模型为例。给定一些观察数据y,假设y符合如下高斯分布：

Laplace算法

隐马尔科夫模型HMM

再谈马尔科夫链蒙特卡洛(MCMC)

变分贝叶斯

· 前言

变分贝叶斯方法最早由Matthew J.Beal在他的博士论文《Variational Algorithms for Approximate Bayesian Inference》中提出，作者将其应用于隐马尔科夫模型，混合因子分析，非线性动力学，图模型等。变分贝叶斯是一类用于贝叶斯估计和机器学习领域中近似计算复杂（intractable）积分的技术。它主要应用于复杂的统计模型中，这种模型一般包括三类变量：观测变量(observed variables, data)，未知参数（parameters）和潜变量（latent variables）。在贝叶斯推断中，参数和潜变量统称为不可观测变量(unobserved variables)。变分贝叶斯方法主要是两个目的:

(1) 近似不可观测变量的后验概率，以便通过这些变量作出统计推断。

(2) 对一个特定的模型，给出观测变量的边缘似然函数（或称为证据，evidence）的下界。主要用于模型的选择，认为模型的边缘似然值越高，则模型对数据拟合程度越好，该模型产生Data的概率也越高。

对于第一个目的，蒙特卡洛模拟，特别是用Gibbs取样的MCMC方法，可以近似计算复杂的后验分布，能很好地应用到贝叶斯统计推断。此方法通过大量的样本估计真实的后验，因而近似结果带有一定的随机性。与此不同的是，变分贝叶斯方法提供一种局部最优，但具有确定解的近似后验方法。

从某种角度看，变分贝叶斯可以看做是EM算法的扩展，因为它也是采用极大后验估计(MAP)，即用单个最有可能的参数值来代替完全贝叶斯估计。另外，变分贝叶斯也通过一组相互依然（mutually dependent）的等式进行不断的迭代来获得最优解。

· 问题描述

现在重新考虑一个问题：1）有一组观测数据D，并且已知模型的形式，求参数与潜变量（或不可观测变量）的后验分布：P(Z|D)。

正如上文所描述的后验概率的形式通常是很复杂(Intractable)的,对于一种算法如果不能在多项式时间内求解，往往不是我们所考虑的。因而我们想能不能在误差允许的范围内，用更简单、容易理解(tractable)的数学形式Q(Z)来近似P(Z|D),即。

由此引出如下两个问题：

（1）假设存在这样的Q(Z),那么如何度量Q(Z)与P(Z|D)之间的差异性（dissimilarity）？

（2）如何得到简单的Q(Z)?

对于问题一，幸运的是，我们不需要重新定义一个度量指标。在信息论中，已经存在描述两个随机分布之间距离的度量，即相对熵，或者称为Kullback-Leibler散度。

对于问题二，显然我们可以自主决定Q(Z)的分布，只要它足够简单，且与P(Z|D)接近。然而不可能每次都手工给出一个与P(Z|D)接近且简单的Q(Z)，其方法本身已经不具备可操作性。所以需要一种通用的形式帮助简化问题。那么数学形式复杂的原因是什么？在“模型的选择”部分，曾提到Occam's razor，认为一个模型的参数个数越多，那么模型复杂的概率越大;此外，如果参数之间具有相互依赖关系(mutually dependent)，那么通常很难对参数的边缘概率精确求解。

幸运的是，统计物理学界很早就关注了高维概率函数与它的简单形式，并发展了平均场理论。简单讲就是：系统中个体的局部相互作用可以产生宏观层面较为稳定的行为。于是我们可以作出后验条件独立（posterior independence）的假设。即，。

· Kullback-Leibler散度

在统计学中，相对熵对应的是似然比的对数期望，相对熵D(p||q)度量当真实分布为p而假定分布为q时的无效性。

定义两个概率密度函数为p(x)和q(x)之间的相对熵定义为

KL散度有如下性质：

（1）；

（2），当且仅当p=q时为零；

（3）不满足三角不等式。

Q分布与P分布的KL散度为：

或者

由于对数证据logP(D)被相应的Q所固定，为了使KL散度最小，则只要极大化L(Q)。通过选择合适的Q，使L(Q)便于计算和求极值。这样就可以得到后验P(Z|D)的近似解析表达式和证据（log evidence）的下界L(Q)，又称为变分自由能（variational free energy）：

· 平均场理论（Mean Field Method）

数学上说，平均场的适用范围只能是完全图，或者说系统结构是well-mixed，在这种情况下，系统中的任何一个个体以等可能接触其他个体。反观物理，平均场与其说是一种方法，不如说是一种思想。其实统计物理的研究目的就是期望对宏观的热力学现象给予合理的微观理论。物理学家坚信，即便不满足完全图的假设，但既然这种“局部”到“整体”的作用得以实现，那么个体之间的局部作用相较于“全局”的作用是可以忽略不计的。

根据平均场理论，变分分布Q(Z)可以通过参数和潜在变量的划分（partition）因式分解，比如将Z划分为Z₁…Z_M,

注意这里并非一个不可观测变量一个划分，而应该根据实际情况做决定。当然你也可以这么做，但是有时候，将几个潜变量放在一起会更容易处理。

平均场方法的合理性

在量子多体问题中，用一个（单体）有效场来代替电子所受到的其他电子的库仑相互作用。这个有效场包含所有其他电受到的其他电子的库仑相互作用。这个有效场包含了所有其他电子对该电子的相互作用。利用有效场取代电子之间的库仑相互作用之后，每一个电子在一个有效场中运动，电子与电子之间的运动是独立的(除了需要考虑泡利不相容原理)，原来的多体问题转化为单体问题。

同样在变分分布Q(Z)这个系统中，我们也可以将每一个潜变量划分看成是一个单体，其他划分对其的影响都可以用一个看做是其自身的作用。采用的办法是迭代(Iterative VB(IVB) algorithm)。这是由于当变分自由能取得最大值的时候，划分Z_i与它的互斥集Z_-I(或者更进一步，马尔科夫毯(Markov blanket), mb(Z_i) )具有一个简单的关系：

（为保持文章的连贯性，此处先不证明，下文将详细说明）

于是，对于某个划分Z_i ,我们可以先保持其他划分Z_-i不变，然后用以上关系式更新Z_i。相同步骤应用于其他划分的更新，使得每个划分之间充分相互作用，最终达到稳定值。

具体更新边缘概率（VB-marginal）步骤如下：

（1）初始化Q⁽¹⁾(Z_i)，可随机取；

（2）在第k步，计算Z_-i的边缘密度

（3）计算Z_i的边缘密度

（4）理论上将会收敛，则反复执行(2), (3)直到Q(Z_i), Q(Z_-i)稳定，或稳定在某个小范围内。

（5）最后，得

平均场估计下边缘概率的无意义性（VB-marginals）

注意到Q(Z)估计的是联合概率密度，而对于每一个Q_i(Z_i)，其与真实的边缘概率密度P_i(Z_i)的差别可能是很大的。不应该用Q_i(Z_i)来估计真实的边缘密度，比如在一个贝叶斯网络中，你不应该用它来推测某个节点的状态。而这其实是很糟糕的，相比于其他能够使用节点状态信息来进行局部推测的算法，变分贝叶斯方法更不利于调试。

比如一个标准的高斯联合分布P(µ,x)和最优的平均场高斯估计Q(µ,x)。Q选择了在它自己作用域中的高斯分布，因而变得很窄。此时边缘密度Q_x(x)变得非常小，完全与P_x(x)不同。

· 边缘密度（VB-marginal）公式的证明

上文已经提到我们要找到一个更加简单的函数D(Z)来近似P(Z|D)，同时问题转化为求解证据logP(Z)的下界L(Q)，或者L(Q(Z))。应该注意到L(Q)并非普通的函数，而是以整个函数为自变量的函数，这便是泛函。我们先介绍一下什么是泛函，以及泛函取得极值的必要条件。

泛函的概念

【泛函】设对于(某一函数集合内的)任意一个函数y(x)，有另一个数J[y]与之对应，则称J [y]为y(x)的泛函。泛函可以看成是函数概念的推广。这里的函数集合，即泛函的定义域，通常要求y(x) 满足一定的边界条件，并且具有连续的二阶导数．这样的y(x)称为可取函数。

泛函不同于复合函数，例如g=g(f(x)); 对于后者，给定一个x值，仍然是有一个g值与之对应；对于前者，则必须给出某一区间上的函数y(x)，才能得到一个泛函值J[y]。(定义在同一区间上的)函数不同，泛函值当然不同，为了强调泛函值J[y]与函数y(x)之间的依赖关系，常常又把函数y(x)称为变量函数。

泛函的形式多种多样，通常可以积分形式：

泛函取极值的必要条件

泛函的极值

“当变量函数为y(x)时，泛函J [y]取极大值”的含义就是：对于极值函数y(x)及其“附近”的变量函数y(x) + δy(x)，恒有 ;

所谓函数y(x) + δy(x)在另一个函数y(x)的“附近”，指的是：

| δy(x) | < ε;

有时还要求| (δy)’(x) | < ε.

这里的δy(x)称为函数y(x)的变分。

Euler–Lagrange方程

可以仿造函数极值必要条件的导出办法，导出泛函取极值的必要条件，这里不做严格的证明，直接给出。泛函J[y]取到极大值的必要条件是一级变分δJ[y]为0，其微分形式一般为二阶常微分方程，即Euler-Largange方程：

泛函的条件极值

在约束条件下求函数J[y]的极值，可以引入Largange乘子λ，从而定义一个新的泛函，。仍将δy看成是独立的，则泛函在边界条件下取极值的必要条件就是，

问题求解

对于，将右式第一项定义为能量(Energy)，第二项看做是信息熵(Shannon entropy)。我们只考虑自然对数的形式，因为对于任何底数的对数总是可以通过换底公式将其写成自然对数与一个常量的乘积形式。另外根据平均场假设可以得到如下积分形式，

其中，且满足

考虑划分，其中，先考虑能量项(Energy)（第一项），

其中定义，C为的归一化常数。再考虑熵量(entropy)（第二项），

此时得到泛函，

注意到L(Q(Z))并非只有一个等式，如果不可观测变量有M个划分。那么将有M个方程。为了使得L(Q(Z))达到最大值，同时注意到约束条件，根据泛函求条件极值的必要条件，得，

直接求解将得到Gibbs分布，略显复杂;实际上，注意到KL散度，我们可以直接得到KL散度等于0的时候，L(D)达到最大值，最终得到

C为归一化常数，为联合概率函数在除Z_i本身外的其他划分下的对数期望。又可以写为

例子:高斯混合模型（GMM）

Question2 假设现在有独立同分布（i.d.d.）的训练样本X符合下列混合高斯分布

如何求解高斯混合分布的三组参数？

步骤一：选择无信息先验分布（non-informative prior）

先验分布的不是随便取的，一般可以根据共轭分布方法，Jefferys原则，最大熵原则等来确定。一般要求先验分布应取共轭分布（conjugate distribution）才合适，即先验分布h(θ)与后验分布h(θ|x)属于同一分布类型。本文不展开讨论，直接给出

k表示单高斯分布的个数，N表示样本个数，每个分布的解释，

SymDir() 表示K维对称 Dirichlet分布；它是卡方分布（categorical）或多项式分布（multinomial）的共轭先验分布。

W() 表示Wishart分布；对一个多元高斯分布（multivariate Gaussian distribution），它是Precision矩阵（逆协方差矩阵）的共轭先验。

Mult() 表示多项分布（此处也称卡方分布）；多项式分布是二项式分布的推广，表示在一个K维向量中只有一项为1，其它都为0.

N() 为高斯分布，在这里特别指多元高斯分布。

对变量的解释

是N个训练样本，每一项都是服从多元高斯分布的K维向量。

是一组潜在变量，每一项用于表示对应的样本属于哪个混合部分（mixture component）。

表示每个单高斯分布混合比例（mix proportion）

和分别表示每个单高斯分布参数的均值（mean）和精度（precision）

另外，为了区分联合分布的参数，以上分布的参数如又称为超参数（hyperparameter），并且都是已知量。

步骤二：写出联合概率密度函数

用“盘子表示法”（plate notation）表示贝叶斯多元高斯混合模型，如图所示。

小正方形表示不变的超参数，如β₀ ,ν₀ ,α₀ ,µ₀ ,W₀；圆圈表示随机变量，如；圆圈内的值为已知量。其中[K],[D]表示K、D维的向量，[D,D]表示的矩阵，单个K表示一个有K个值的categorical变量；波浪线和一个开关表示变量x_i通过一个K维向量z_i 来选择其他传入的变量(µ_k,Ʌ_k)。

假设各参数与潜在变量条件独立，则联合概率密度函数可以表示为

每个因子为：

其中，

D为各观测点的维度。

步骤三：计算边缘密度(VB- marginal)

（1）计算Z的边缘密度，根据平均场假设，，则

其中，

两边分别取对数可得，

归一化，即对于观测变量的属于某个单高斯分布的概率相加应等于1,则有

，其中

可见是多个单观测多项式分布(single-observation multinomial distribution)的乘积，可以因式分解成一个个以为参数的单观测多项式分布。更进一步，根据categorical分布，有 .

另外，需特别注意的是，在求期望的过程中，由于联合密度可以写成几个因子乘积的形式，而方程是关于Z的函数，因此对于不包含Z的密度函数可以当做常数处理。我们可以用马尔科夫毯（Markov blanket）描述，在一个贝叶斯网络中，它表示一个节点的父节点(parents)，子节点和子节点的父节点（co-parents）,如图所示。我们将在后文中详细说明。

（2）计算π的概率密度，

两边取对数，可见是Dirichlet分布，

其中， .

（3）最后同时考虑，对于每一个单高斯分布有，

经过一系列重组化解将得到Gaussian-Wishart分布，

其中定义，

步骤四：迭代收敛

最后，注意到对π,µ,Ʌ的边缘概率都需要且只需要r_nk；另一方面，r_nk的计算需要ρ_nk，而这又是基于 , , ,即需要知道π,µ,Ʌ的值。不难确定这三个期望的一般表达式为：

这些结果能够导出：。由于需要归一化使得，这样可以从线性相对值转化为绝对值。

再次分析各参数，

参数变量更新方程中的超参数都依赖与统计量，而这些统计量又依赖于 .

参数变量更新方程中的超参数都依赖于统计量，即 .

潜在变量的更新方程对超变量有直接的依赖关系，同时对通过有间接的依赖关系。

这样迭代过程便很清楚了，可以总结为如下两个迭代步骤：

1. 在VBE-Step，用参数和超参数的旧值计算潜在变量；

2. 在VBM-Step，用潜在变量计算参数和超参数的新值。

· 与EM算法的比较

以上迭代步骤与EM算法用ML或MAP解决高斯混合模型很相似。在E-step中，潜在变量r_nk对应于潜在变量关于数据样本的后验概率，如P(Z|X); 统计量对应于EM算法中“soft-count”统计量；用这些统计量去计算参数的新值与EM算法中用“soft-count”计算新参数值一致。

但经过如此，VBEM算法与EM算法还有有很多不同之处的。比如迭代中，逼近最优值的过程是不一样的，如下图所示。在有限制的情况下，EM算法极大似然值是动态变化的。刚开始与当前最优值相差一个KL.在E步骤，下界逼近最大似然值（或者由于条件限制相差一点）；然后在E步骤中，根据新参数重新确定新的似然值。如此往复，直到达到稳定。而在VBEM算法中，极大似然值是不变的。VBE与VBM步骤，都是逼近极大似然值的过程。

EM算法（with constrained）

VBEM算法

进一步讨论

· 再谈EM与VBEM

EM算法计算随机变量（或归类于参数）后验分布的点估计，但估计潜在变量的真实后验分布（至少在soft EM算法，并且经常只当潜在变量是离散化的情况）。这些参数的众数（modes）作为点估计，无任何其他信息。

而在VB算法，计算所有变量的真实后验分布的估计，包括参数和潜在变量。计算点估计的过程中，一般使用在贝叶斯推断中常用的均值（mean）而非众数。与此同时，应该注意的是计算参数在VB中与EM有不同的意义。EM算法计算贝叶斯网络本身的参数的最优值。而VB计算用于近似参数和潜在变量的贝叶斯网络的参数最佳值。正如之前的高斯混合模型例子，对于每一个混合部分都需要计算其参数。EM算法将会直接估计这些参数的最优值；而VB会先找一个合适的参数分布，通常是一个先验分布的形式，然后计算这个分布的参数值，更准确说是超参数，最后得到联合分布的各参数。

· 算法复杂度

· 信息耦合(Intractable Coupling)

改进：变分消息传递(Variational Message Passing, VMP)

从高斯混合模型的例子可以看出，传统的变分贝叶斯方法对模型的推导是繁琐而复杂的。而在推导边缘概率的时候，我们也提到对数联合概率的在一些参数下的期望可以简化，我们只需要关心所求参数的马尔科夫毯上的节点。另外，又认识到许多先验和条件概率属于指数分布族，而对数可以将乘积形式展开为加法。那么，我们是不是可以找到一些简单计算方法或者统一的形式呢？

这便是变分消息传递（Variational Message Passing; John Winn, Bishop 2005）他们考虑了贝叶斯网络中的共轭指数网络（conjugate-exponential networks）,这种方法使得充分统计量（sufficient statistics）与自然参数（natural parameter）都有一个标准形式（standard form），现在该方法已经取代了手工推导(derivation by hand),成为标准的变分贝叶斯推断方法。而对于非共轭指数网络（比如混合模型），也能通过进一步的近似转化为标准形式。

· 贝叶斯网络

变分信息传递方法是建立在贝叶斯网络上的，如图所示，对于一个节点H_j,它的父节点为pa_j，子节点为ch_j,子节点X_k的父节点为。所有节点可以统称为H_j的马尔科夫毯，对于变分贝叶斯推理，我们只需要关系这个模型，H为参数或潜在变量，其父节点为它的超参数，子节点为数据样本， co-parents为其他参数或潜在变量。

具体的，对于混合高斯模型，我们有如下图模型，

· 指数分布族

定义：设是可控参数统计结构，加入其密度函数可表示为如下形式：

并且它的支撑不依赖于θ，则称此结构为指数型的统计结构，简称指数结构，其中的分布族为指数分布族，这里的都与θ无关，且取有限值的B可测函数，k为正整数，h(x)>0，常见指数分布族，如二项分布，二元正态分布，伽马分布。

对于一个条件分布，如果它能写成如下形式，则称它属于指数分布族，

其中称为自然参数（natural parameter）向量，称为自然统计（natural statistic）向量。作为归一化函数使得对于任意的Y都能整合到统一的形式。指数分布族的好处是它的对数形式是可计算的（be tractable to compute）并且它的状态可以用自然参数向量所概括。

共轭指数模型（Conjugate-Exponential Model）当变量X关于父节点Y的条件概率分布P(X|Y)为指数分布族，且为父节点分布P(Y)的共轭先验，那么我们称这样的模型是共轭指数模型。

我们考虑共轭指数模型，这样后验的每个因子与它的先验都有相同的形式，因而我们只需要关心参数的变化，而无需整个函数。所谓相同的形式是指属于同样的分布，比如都属于正态分布，伽马分布，多项式分布等，后面我们将详细说明。

· 自然统计量的期望

如果我们知道自然参数向量，那么我们就能找到自然统计量的期望。重写指数分布族，用ϕ作为参数，g重新参数化( reparameterisation )为则,

对X积分有，

然后对ϕ微分，

得自然统计量的期望，

（1）如何计算？u(X)表示什么意思？

· 变分分布Q在共轭指数模型下的最优化

不失一般性，考虑变分分布的一个因子Q(Y),Y为马尔科夫毯上一个节点,子节点为X，如下图所示。

根据指数族条件分布的一般形式，则变量Y关于父节点的条件概率为，

（2）

的下标Y用于区分不同节点对数条件概率中的各成员。考虑Y的子节点，则X的关于其父节点的条件概率为，

（3）

可以将可以看出Y的先验，作为Y的似然函数。共轭的要求是这两个条件分布具有关于Y相同的函数形式，因而可以通过定义和函数将后者改写成

（4）

为了更新Q(Y),需要找到(2),(3)式关于除Y外其他因子的期望。对任何指数族的自然统计量u的期望都可以用自然参数向量ϕ带入（1）式得到。即对于任何变量Ａ，都可以找到。特别的，当A为被观测量时，我们能直接计算得 .

从(3)，(4)式可以看出与分布成线性关系。而共轭要求对数条件分布也会与所有的成线性，。因而看得出是一个关于u的多线性函数(multi-linear function)。

重新考虑Y的变分更新方程，

其中（5）

正如以上所解释的，和的期望都是关于自然统计量的多线性函数。因而有可能将以上期望重新参数化为

举例：如果X服从，那么

其中

可以重参数化为

· 下界L(Q)的计算

在贝叶斯网络中,由于Q可因式分解，则有

L(Q)被分解为每一个节点上的贡献值（contribution）{L_i}，如节点H_j的贡献值为

在变分消息传递算法中，和在找的后验分布时就已经计算了；在传出消息的时候也已经计算了。这样下界可以是分开计算的，大大降低了计算复杂性。

特别地，对于每个观测变量对下界的贡献值更简单，

· 变分消息传递算法

现在我们已经准确知道节点之间的消息应该由什么样的形式传递，那么定义变分消息传递算法：

来自父节点的消息（Message from parents）：父节点传递给子节点的消息只是自然统计量的期望：

（6）

消息传递给父节点（Message to parents）：依赖于X之前从Y的co-parents接收到的消息；对任何节点A，如果A是被观测量，那么 ,

（7）

用Y接收来自父节点与子节点的所有消息来计算，然后我们就能通过计算更新后的自然参数向量来找到Y的更新后的后验分布，的计算公式如下，

（8）

该式与（5）式一致。从（1）式可以看出自然统计量的期望是的单一函数，这样我们就可以用它来计算期望的新值。变分消息传递算法通过反复迭代的消息传递来最优化变分分布。

完整的算法描述如下，

Step1. 通过初始化相关的矩向量来初始化每个因子分布 .

Step2. 对于每一个节点，

从父节点和子节点接收 (6),(7)式所定义的消息。前提是子节点已经从的co-parents接收到消息。

使用（8）式更新自然参数向量

根据新的参数向量更新距向量

Step3. 计算新的下界（如果需要）

Step4. 如果经过数次迭代已经无法增加下界值，或者各边缘分布达到稳定值，则结束；否则回到第二步。

举例：对于单一高斯模型（univariate Gaussian Model）消息传递过程如下图(a)(b)(c)(d)，

· 混合模型

到目前为止，我们只考虑了来自指数族的分布。而通常来讲混合分布(Mixture Distribution)并非来自指数族，比如高斯混合模型。那么我们是否能将这些混合分布转化为指数族的分布呢？

考虑高斯混合模型，通常有如下形式，

我们可以引入一个离散型潜在变量λ,表示每个观测点是属于哪个单高斯分布。在变分贝叶斯方法中，我们已经举过该例子，故不加以描述。重写分布函数为：

加入该λ变量后该分布属于指数分布族，可写成

如果X有孩子Z，那么共轭条件要求每一个成分都有相同的自然统计向量，统一定义为。另外，我们可能要使模型的其他部分也有相同的形式，虽然不要求共轭，即。在这种情况下，混合模型的每个成分都有相同的形式，我们写成，

其中定义。这样对于每个成分来说条件分布都有了与指数分布族一样的形式。我们便可以应用变分消息传递算法。从某个节点X传递个子节点的消息为，而这是通过混合参数向量计算的。相似地，节点X到父亲节点的消息是那些以它为父节点的子节点发出的，而节点X中哪些属于是由指标变量的后验确定的。最后，从X到的消息是一个K维向量，其中第k个元素为 .