变分贝叶斯方法

变分贝叶斯方法是近似贝叶斯推理与机器学习中较难积分的一系列技术，通常用于由观测变量及未知参数、潜变量组成的复杂统计模型，之间存在的各种关系可用图模式描述。贝叶斯推理中，参数和潜变量一般归为“未观察变量”。变分贝叶斯方法主要用于两个目的：

1. 提供未观测变量后验概率的分析近似，以便对变量统计推断
2. 得出观测数据的边缘似然的下界（即给定模型的数据的边缘概率，边缘化未观测变量）。通常用于模型选择，一般认为给定模型的边缘似然越高，拟合效果越好，产生数据的可能性越大（另见贝叶斯因子）。

就目的1（近似后验概率），变分贝叶斯是蒙特卡洛法（特别是马尔可夫链蒙特卡洛，如吉布斯采样）的替代，可采用完全贝叶斯的方法，对难以取值或抽样的复杂分布进行统计推断。蒙特卡洛技术利用一组样本对精确后验值进行数值近似，而变分贝叶斯可求近似后验的局部最优的精确解析解。

变分贝叶斯可视作最大期望算法（EM）的推广，从估计参数的单一最似然值的最大后验概率（MAP估计），推广到计算参数及潜变量的整个近似后验分布的完全贝叶斯估计。贝叶斯估计同EM一样，也能找到一组最优参数值，且与EM有相同的交替结构——基于不能求解析解的相依方程组。

很多应用中，变分贝叶斯能更快得到与吉布斯采样精度相当的解，但推导迭代方程组需要更多工作。即便是很多概念上非常简单的模型也如此，下面以只有2个参数、没有潜变量的基本不分层模型为例说明。

数学推导

问题

变分推断中，数据${\displaystyle \mathbf {X} }$的一组未观测变量${\displaystyle \mathbf {Z} =\{Z_{1}\dots Z_{n}\}}$的后验分布近似于所谓变分分布${\displaystyle Q(\mathbf {Z} )}$

${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )\approx Q(\mathbf {Z} ).}$

分布${\displaystyle Q(\mathbf {Z} )}$属于在形式上比${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$简单的分布族（如高斯分布族），这是为使${\displaystyle Q(\mathbf {Z} )}$与真实后验${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$更相似。

相似性用不相似函数${\displaystyle d(Q;P)}$衡量，推断（inference）通过求分布${\displaystyle Q(\mathbf {Z} )}$使${\displaystyle d(Q;P)}$最小化进行。

KL散度

最常见的变分贝叶斯法以Q与P的KL散度作为不相似函数，很适合最小化。KL散度的定义是

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)\triangleq \sum _{\mathbf {Z} }Q(\mathbf {Z} )\log {\frac {Q(\mathbf {Z} )}{P(\mathbf {Z} \mid \mathbf {X} )}}.}$

注意Q和P是相反的。反向KL散度在概念上类似于期望最大化算法（KL散度则产生期望传播算法）。

难解性

变分技术通常用于近似以下参数：

${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )={\frac {P(\mathbf {X} \mid \mathbf {Z} )P(\mathbf {Z} )}{P(\mathbf {X} )}}={\frac {P(\mathbf {X} \mid \mathbf {Z} )P(\mathbf {Z} )}{\int _{\mathbf {Z} }P(\mathbf {X} ,\mathbf {Z} ')\,d\mathbf {Z} '}}}$

对${\displaystyle \mathbf {Z} }$边缘化（以计算分母中的${\displaystyle P(\mathbf {X} )}$）通常是难解的，因为${\displaystyle \mathbf {Z} }$的搜索空间在组合上可能很大。因此，我们用${\displaystyle Q(\mathbf {Z} )\approx P(\mathbf {Z} \mid \mathbf {X} )}$求近似。

证据下界

主条目：证据下界

考虑到${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )={\frac {P(\mathbf {X} ,\mathbf {Z} )}{P(\mathbf {X} )}}}$，KL散度也可写作

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log {\frac {Q(\mathbf {Z} )}{P(\mathbf {Z} ,\mathbf {X} )}}+\log P(\mathbf {X} )\right]=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log P(\mathbf {X} )\right]}$

因为${\displaystyle Q(\mathbf {Z} )}$是分布，所以${\displaystyle \sum _{\mathbf {Z} }Q(\mathbf {Z} )=1}$；又由于${\displaystyle P(\mathbf {X} )}$对于${\displaystyle \mathbf {Z} }$是常数，所以

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\log P(\mathbf {X} )}$

根据（离散随机变量）期望的定义，上式可以写作

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)=\mathbb {E} _{\mathbf {Q} }\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\log P(\mathbf {X} )}$

重排为

${\displaystyle \log P(\mathbf {X} )=D_{\mathrm {KL} }(Q\parallel P)-\mathbb {E} _{\mathbf {Q} }\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]=D_{\mathrm {KL} }(Q\parallel P)+{\mathcal {L}}(Q)}$

由于对数证据${\displaystyle \log P(\mathbf {X} )}$对于Q是定值，最大化末项${\displaystyle {\mathcal {L}}(Q)}$会使Q对P的KL散度最小化。适当选择Q可以让${\displaystyle {\mathcal {L}}(Q)}$更容易计算，更容易最大化，因此既有后验${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$的解析近似Q，也有对数证据${\displaystyle \log P(\mathbf {X} )}$的下界${\displaystyle {\mathcal {L}}(Q)}$（KL散度非负）。

下界${\displaystyle {\mathcal {L}}(Q)}$与热力学自由能类似，也称作（负）变分自由能，因为它也可以表为负能量${\displaystyle \operatorname {E} _{Q}[\log P(\mathbf {Z} ,\mathbf {X} )]}$加Q的熵。${\displaystyle {\mathcal {L}}(Q)}$项也称作证据下界（ELBO），是数据的对数证据的下界。

证明

根据布雷格曼散度（KL散度是其特例）的广义勾股定理，可以证明^[1][2]

布雷格曼散度的广义勾股定理^[2]

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)\geq D_{\mathrm {KL} }(Q\parallel Q^{*})+D_{\mathrm {KL} }(Q^{*}\parallel P),\forall Q^{*}\in {\mathcal {C}}}$

其中${\displaystyle {\mathcal {C}}}$是凸集，且

${\displaystyle Q=Q^{*}\triangleq \arg \min _{Q\in {\mathcal {C}}}D_{\mathrm {KL} }(Q\parallel P).}$

时等号成立。这时，全局最小值${\displaystyle Q^{*}(\mathbf {Z} )=q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})q^{*}(\mathbf {Z} _{2})=q^{*}(\mathbf {Z} _{2}\mid \mathbf {Z} _{1})q^{*}(\mathbf {Z} _{1})}$，其中${\displaystyle \mathbf {Z} =\{\mathbf {Z_{1}} ,\mathbf {Z_{2}} \}}$可以如下求得：^[1]

${\displaystyle q^{*}(\mathbf {Z} _{2})={\frac {P(\mathbf {X} )}{\zeta (\mathbf {X} )}}{\frac {P(\mathbf {Z} _{2}\mid \mathbf {X} )}{\exp(D_{\mathrm {KL} }(q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})\parallel P(\mathbf {Z} _{1}\mid \mathbf {Z} _{2},\mathbf {X} )))}}={\frac {1}{\zeta (\mathbf {X} )}}\exp \mathbb {E} _{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}\left(\log {\frac {P(\mathbf {Z} ,\mathbf {X} )}{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}}\right),}$

其中归一化常数为

${\displaystyle \zeta (\mathbf {X} )=P(\mathbf {X} )\int _{\mathbf {Z} _{2}}{\frac {P(\mathbf {Z} _{2}\mid \mathbf {X} )}{\exp(D_{\mathrm {KL} }(q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})\parallel P(\mathbf {Z} _{1}\mid \mathbf {Z} _{2},\mathbf {X} )))}}=\int _{\mathbf {Z} _{2}}\exp \mathbb {E} _{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}\left(\log {\frac {P(\mathbf {Z} ,\mathbf {X} )}{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}}\right).}$

${\displaystyle \zeta (\mathbf {X} )}$项在实践中常称作证据下界（ELBO），因为${\displaystyle P(\mathbf {X} )\geq \zeta (\mathbf {X} )=\exp({\mathcal {L}}(Q^{*}))}$，^[1]如上所示。

交换${\displaystyle \mathbf {Z} _{1},\ \mathbf {Z} _{2}}$的角色，可以迭代计算真实模型边际${\displaystyle P(\mathbf {Z} _{1}\mid \mathbf {X} ),\ P(\mathbf {Z} _{2}\mid \mathbf {X} )}$的近似${\displaystyle q^{*}(\mathbf {Z} _{1}),\ q^{*}(\mathbf {Z} _{2})}$。这种迭代保证单调收敛，^[1]但收敛到的${\displaystyle Q^{*}}$只是${\displaystyle D_{\mathrm {KL} }(Q\parallel P)}$的局部极小值。

若约束空间${\displaystyle {\mathcal {C}}}$被限制在独立空间内，即${\displaystyle q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})=q^{*}(\mathbf {Z_{1}} )}$，则上述迭代将成为所谓平均场近似（mean field approximation）${\displaystyle Q^{*}(\mathbf {Z} )=q^{*}(\mathbf {Z} _{1})q^{*}(\mathbf {Z} _{2}),}$。

平均场近似

变分分布${\displaystyle Q(\mathbf {Z} )}$常被假定为潜变量的某划分的因式分解，即对潜变量${\displaystyle \mathbf {Z} }$的划分${\displaystyle \mathbf {Z} _{1}\dots \mathbf {Z} _{M}}$，

${\displaystyle Q(\mathbf {Z} )=\prod _{i=1}^{M}q_{i}(\mathbf {Z} _{i}\mid \mathbf {X} )}$

由变分法可以证明，对每个因子${\displaystyle q_{j}}$的最佳分布${\displaystyle q_{j}^{*}}$（就KL散度最小的分布而言，如上所述）满足^[3]

${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )={\frac {e^{\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}}{\int e^{\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}\,d\mathbf {Z} _{j}}}}$

其中${\displaystyle \operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}$是数据与潜变量的对数联合概率的期望值，是相对于不属于划分的所有变量的${\displaystyle q^{*}}$而取。关于${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )}$的推导见^[4]的引理 4.1。

实践中，通常用对数表示：

${\displaystyle \ln q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )=\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]+{\text{constant}}}$

上式中的常数与归一化常数（${\displaystyle q_{j}^{*}}$上述表达式中的分母）有关，通常通过检查来恢复，因为表达式的其余部分通常是已知的分布（如正态分布、伽马分布等）。

利用期望的性质，式${\displaystyle \operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}$通常可简为潜变量先验分布的固定超参数及不属于当前划分的潜变量（即不属于${\displaystyle \mathbf {Z} _{j}}$的潜变量）的期望值的函数（有时还包括方差等高阶矩）。这就在某一划分的变量分布参数同其他划分的变量的期望值之间产生了循环依赖，自然就需要类似期望最大化算法的迭代法，以某种方式（也许随机）初始化潜变量期望（可能还有高阶矩），再利用当前期望依次计算每个分布的参数、适当设置新算出的分布的期望。这种算法保证收敛。^[5]

换一种说法，对每个变量划分，简化其中变量分布的表达式、研究分布与相关变量的函数依赖关系，通常就能确定分布的族（反过来确定了常数）。分布的参数公式可用先验分布的超参数（已知常数）表示，也可用其他划分中变量函数的期望表示。通常这些期望可简为变量本身的期望值函数（即平均值）；有时也会出现变量平方（与变量方差有关）或更高次幂（）高阶矩）的期望。其他变量的分布一般来自已知族，相关期望的公式可查到；但公式取决于分布的参数，参数又取决于其他变量的期望，所以每个变量的分布参数公式都可表为一群变量间非线性相互依赖的方程。通常不能直接求解，不过可以用简单的迭代算法，大多数时候都能保证收敛。

变分推理对偶公式

用对偶公式图解坐标上升变分推理算法^[4]

下面的定理称作变分推理对偶公式，^[4]解释了变分贝叶斯方法中变分分布的一些重要性质。

Theorem 考虑两概率空间${\displaystyle (\Theta ,{\mathcal {F}},P),\ (\Theta ,{\mathcal {F}},Q)}$，其中${\displaystyle Q\ll P}$。假设存在共同的主导概率测度${\displaystyle \lambda }$使${\displaystyle P\ll \lambda ,\ Q\ll \lambda }$。令h是${\displaystyle (\Theta ,{\mathcal {F}},P)}$上的任意实值随机变量，满足${\displaystyle h\in L_{1}(P)}$，则下列等式成立

${\displaystyle \log E_{P}[\exp h]={\text{sup}}_{Q\ll P}\{E_{Q}[h]-D_{\text{KL}}(Q\parallel P)\}.}$

此外，当且仅当（supremum）

${\displaystyle {\frac {q(\theta )}{p(\theta )}}={\frac {\exp h(\theta )}{E_{P}[\exp h]}},}$

关于概率测度Q几乎是确定的，其中${\displaystyle p(\theta )=dP/d\lambda ,\ q(\theta )=dQ/d\lambda }$分别表示概率测度P、Q对${\displaystyle \lambda }$的Radon–Nikodym导数，右式才取上界。

基本例子

考虑简单的不分层贝叶斯模型，由来自正态分布（均值、方差未知）的独立同分布观测量组成。^[6]下面我们将详细介绍这模型，以说明变分贝叶斯方法的工作原理。

为了数学上的方便，下例中我们用精度（即方差倒数；多元正态分布时是协方差矩阵的逆。理论上精度和方差等价，因为两者间是双射）。

数学模型

将共轭先验分布置于未知均值${\displaystyle \mu }$、精度${\displaystyle \tau }$，即均值也遵循高斯分布，而精度遵循伽马分布。

${\displaystyle {\begin{aligned}\tau &\sim \operatorname {Gamma} (a_{0},b_{0})\\\mu |\tau &\sim {\mathcal {N}}(\mu _{0},(\lambda _{0}\tau )^{-1})\\\{x_{1},\dots ,x_{N}\}&\sim {\mathcal {N}}(\mu ,\tau ^{-1})\\N&={\text{number of data points}}\end{aligned}}}$

先验分布超参数${\displaystyle \mu _{0},\lambda _{0},a_{0}}$、${\displaystyle b_{0}}$是已知定值，可设为小正数，以得到较宽的先验分布，表明我们对${\displaystyle \mu }$、${\displaystyle \tau }$的先验分布一无所知。

已知N个数据点${\displaystyle \mathbf {X} =\{x_{1},\ldots ,x_{N}\}}$，而我们的目标是推断参数${\displaystyle \mu ,\ \tau }$的后验分布${\displaystyle q(\mu ,\tau )=p(\mu ,\tau \mid x_{1},\ldots ,x_{N})}$。

联合概率

所有变量的联合概率可重写为

${\displaystyle p(\mathbf {X} ,\mu ,\tau )=p(\mathbf {X} \mid \mu ,\tau )p(\mu \mid \tau )p(\tau )}$

个体因子是

${\displaystyle {\begin{aligned}p(\mathbf {X} \mid \mu ,\tau )&=\prod _{n=1}^{N}{\mathcal {N}}(x_{n}\mid \mu ,\tau ^{-1})\\p(\mu \mid \tau )&={\mathcal {N}}\left(\mu \mid \mu _{0},(\lambda _{0}\tau )^{-1}\right)\\p(\tau )&=\operatorname {Gamma} (\tau \mid a_{0},b_{0})\end{aligned}}}$

其中

${\displaystyle {\begin{aligned}{\mathcal {N}}(x\mid \mu ,\sigma ^{2})&={\frac {1}{\sqrt {2\pi \sigma ^{2}}}}e^{\frac {-(x-\mu )^{2}}{2\sigma ^{2}}}\\\operatorname {Gamma} (\tau \mid a,b)&={\frac {1}{\Gamma (a)}}b^{a}\tau ^{a-1}e^{-b\tau }\end{aligned}}}$

因式分解近似

假设${\displaystyle q(\mu ,\tau )=q(\mu )q(\tau )}$，即后验分布因式分解为${\displaystyle \mu ,\ \tau }$的独立因子。这种假设是变分贝叶斯方法的基础。事实上，真正的后验分布并非如此（这种简单情形我们知道它是正态-伽马分布），因此所得结果将是近似值。

导出.mw-parser-output .serif{font-family:Times,serif}q(μ)

则

${\displaystyle {\begin{aligned}\ln q_{\mu }^{*}(\mu )&=\operatorname {E} _{\tau }\left[\ln p(\mathbf {X} \mid \mu ,\tau )+\ln p(\mu \mid \tau )+\ln p(\tau )\right]+C\\&=\operatorname {E} _{\tau }\left[\ln p(\mathbf {X} \mid \mu ,\tau )\right]+\operatorname {E} _{\tau }\left[\ln p(\mu \mid \tau )\right]+\operatorname {E} _{\tau }\left[\ln p(\tau )\right]+C\\&=\operatorname {E} _{\tau }\left[\ln \prod _{n=1}^{N}{\mathcal {N}}\left(x_{n}\mid \mu ,\tau ^{-1}\right)\right]+\operatorname {E} _{\tau }\left[\ln {\mathcal {N}}\left(\mu \mid \mu _{0},(\lambda _{0}\tau )^{-1}\right)\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\ln \prod _{n=1}^{N}{\sqrt {\frac {\tau }{2\pi }}}e^{-{\frac {(x_{n}-\mu )^{2}\tau }{2}}}\right]+\operatorname {E} _{\tau }\left[\ln {\sqrt {\frac {\lambda _{0}\tau }{2\pi }}}e^{-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}}\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}\left({\frac {1}{2}}(\ln \tau -\ln 2\pi )-{\frac {(x_{n}-\mu )^{2}\tau }{2}}\right)\right]+\operatorname {E} _{\tau }\left[{\frac {1}{2}}(\ln \lambda _{0}+\ln \tau -\ln 2\pi )-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}-{\frac {(x_{n}-\mu )^{2}\tau }{2}}\right]+\operatorname {E} _{\tau }\left[-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}\right]+\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}{\frac {1}{2}}(\ln \tau -\ln 2\pi )\right]+\operatorname {E} _{\tau }\left[{\frac {1}{2}}(\ln \lambda _{0}+\ln \tau -\ln 2\pi )\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}-{\frac {(x_{n}-\mu )^{2}\tau }{2}}\right]+\operatorname {E} _{\tau }\left[-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}\right]+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right\}+C_{3}\end{aligned}}}$

上面的推导中，${\displaystyle C,\ C_{2},\ C_{3}}$指对${\displaystyle \mu }$为常数的值。注意${\displaystyle \operatorname {E} _{\tau }[\ln p(\tau )]}$项不是${\displaystyle \mu }$的函数，无论${\displaystyle \mu }$的值是多少，它都不变。因此，第三行中可以把它吸收到末尾的常数项，第七行也如此。

最后一行是${\displaystyle \mu }$的二次多项式。由于这是${\displaystyle q_{\mu }^{*}(\mu )}$的对数，我们可以看到${\displaystyle q_{\mu }^{*}(\mu )}$本身是正态分布。

通过一些繁琐的计算（展开括号里的平方、分离出涉及${\displaystyle \mu }$、${\displaystyle \mu ^{2}}$的项、对${\displaystyle \mu }$应用配方法），可得到正态分布的参数：

${\displaystyle {\begin{aligned}\ln q_{\mu }^{*}(\mu )&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\sum _{n=1}^{N}(x_{n}^{2}-2x_{n}\mu +\mu ^{2})+\lambda _{0}(\mu ^{2}-2\mu _{0}\mu +\mu _{0}^{2})\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\left(\sum _{n=1}^{N}x_{n}^{2}\right)-2\left(\sum _{n=1}^{N}x_{n}\right)\mu +\left(\sum _{n=1}^{N}\mu ^{2}\right)+\lambda _{0}\mu ^{2}-2\lambda _{0}\mu _{0}\mu +\lambda _{0}\mu _{0}^{2}\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\mu ^{2}-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu +\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\mu ^{2}-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu \right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)(\lambda _{0}+N)\mu \right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)\mu \right)\right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)\mu +\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}-\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}\right)\right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)\mu +\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}\right)\right\}+C_{5}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu -{\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}\right\}+C_{5}\\&=-{\frac {1}{2}}(\lambda _{0}+N)\operatorname {E} _{\tau }[\tau ]\left(\mu -{\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}+C_{5}\end{aligned}}}$

注意上面所有步骤用二次和公式都能化简。

即

${\displaystyle {\begin{aligned}q_{\mu }^{*}(\mu )&\sim {\mathcal {N}}(\mu \mid \mu _{N},\lambda _{N}^{-1})\\\mu _{N}&={\frac {\lambda _{0}\mu _{0}+N{\bar {x}}}{\lambda _{0}+N}}\\\lambda _{N}&=(\lambda _{0}+N)\operatorname {E} _{\tau }[\tau ]\\{\bar {x}}&={\frac {1}{N}}\sum _{n=1}^{N}x_{n}\end{aligned}}}$

推导q(τ)

${\displaystyle q_{\tau }^{*}(\tau )}$的推导大致相同，简洁起见略去细节。

${\displaystyle {\begin{aligned}\ln q_{\tau }^{*}(\tau )&=\operatorname {E} _{\mu }[\ln p(\mathbf {X} \mid \mu ,\tau )+\ln p(\mu \mid \tau )]+\ln p(\tau )+{\text{constant}}\\&=(a_{0}-1)\ln \tau -b_{0}\tau +{\frac {1}{2}}\ln \tau +{\frac {N}{2}}\ln \tau -{\frac {\tau }{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]+{\text{constant}}\end{aligned}}}$

两侧取指数，可见${\displaystyle q_{\tau }^{*}(\tau )}$服从伽马分布。

${\displaystyle {\begin{aligned}q_{\tau }^{*}(\tau )&\sim \operatorname {Gamma} (\tau \mid a_{N},b_{N})\\a_{N}&=a_{0}+{\frac {N+1}{2}}\\b_{N}&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]\end{aligned}}}$

计算参数

回忆前几节的结论：

${\displaystyle {\begin{aligned}q_{\mu }^{*}(\mu )&\sim {\mathcal {N}}(\mu \mid \mu _{N},\lambda _{N}^{-1})\\\mu _{N}&={\frac {\lambda _{0}\mu _{0}+N{\bar {x}}}{\lambda _{0}+N}}\\\lambda _{N}&=(\lambda _{0}+N)\operatorname {E} _{\tau }[\tau ]\\{\bar {x}}&={\frac {1}{N}}\sum _{n=1}^{N}x_{n}\end{aligned}}}$

${\displaystyle {\begin{aligned}q_{\tau }^{*}(\tau )&\sim \operatorname {Gamma} (\tau \mid a_{N},b_{N})\\a_{N}&=a_{0}+{\frac {N+1}{2}}\\b_{N}&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]\end{aligned}}}$

每种情形下，某变量的分布参数取决于对另一变量的期望。可用正态分布和伽马分布矩期望的标准公式推广期望：

${\displaystyle {\begin{aligned}\operatorname {E} [\tau \mid a_{N},b_{N}]&={\frac {a_{N}}{b_{N}}}\\\operatorname {E} \left[\mu \mid \mu _{N},\lambda _{N}^{-1}\right]&=\mu _{N}\\\operatorname {E} \left[X^{2}\right]&=\operatorname {Var} (X)+(\operatorname {E} [X])^{2}\\\operatorname {E} \left[\mu ^{2}\mid \mu _{N},\lambda _{N}^{-1}\right]&=\lambda _{N}^{-1}+\mu _{N}^{2}\end{aligned}}}$

大多数时候将这些公式应用到上述方程不困难，但${\displaystyle b_{N}}$需要更多工作：

${\displaystyle {\begin{aligned}b_{N}&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]\\&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[(\lambda _{0}+N)\mu ^{2}-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu +\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\\&=b_{0}+{\frac {1}{2}}\left[(\lambda _{0}+N)\operatorname {E} _{\mu }[\mu ^{2}]-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\operatorname {E} _{\mu }[\mu ]+\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\\&=b_{0}+{\frac {1}{2}}\left[(\lambda _{0}+N)\left(\lambda _{N}^{-1}+\mu _{N}^{2}\right)-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu _{N}+\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\\\end{aligned}}}$

这样就可以写出参数方程如下，不需要期望：

${\displaystyle {\begin{aligned}\mu _{N}&={\frac {\lambda _{0}\mu _{0}+N{\bar {x}}}{\lambda _{0}+N}}\\\lambda _{N}&=(\lambda _{0}+N){\frac {a_{N}}{b_{N}}}\\{\bar {x}}&={\frac {1}{N}}\sum _{n=1}^{N}x_{n}\\a_{N}&=a_{0}+{\frac {N+1}{2}}\\b_{N}&=b_{0}+{\frac {1}{2}}\left[(\lambda _{0}+N)\left(\lambda _{N}^{-1}+\mu _{N}^{2}\right)-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu _{N}+\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\end{aligned}}}$

注意${\displaystyle \lambda _{N}}$和${\displaystyle b_{N}}$存在相互依赖关系，自然产生了类似最大期望算法的算法：

1. 计算${\displaystyle \sum _{n=1}^{N}x_{n},\ \sum _{n=1}^{N}x_{n}^{2}.}$用所得值计算${\displaystyle \mu _{N},\ a_{N}.}$
2. 初始化${\displaystyle \lambda _{N}}$为任意值。
3. 用${\displaystyle \lambda _{N}}$的现值和其它参数的已知值计算${\displaystyle b_{N}}$。
4. 用${\displaystyle b_{N}}$的现值和其它参数的已知值计算${\displaystyle \lambda _{N}}$。
5. 重复最后两部，直到收敛（即直到两值的更新量不超过阈值）。

然后就有了后验参数近似分布的超参数值，可用于计算后验的任何属性，如均值、方差、95%最高密度区域（包含总概率95%的最小区间）等等。

可以证明这种算法保证收敛到局部极大值。

还要注意，后验分布与相应的先验分布有同样的形式。只需假设分布可以分解，分布形式便随之而来。事实证明，后验、先验形式相同并非巧合，而是先验分布为指数族成员时的一般结果，大多数标准分布都如此。

进一步讨论

分步法

上面例子展示了给定贝叶斯网络中推导后验概率密度的变分贝叶斯近似的方法：

1. 用图模式描述网络，确定观察变量（数据）${\displaystyle \mathbf {X} }$和未观察变量（参数${\displaystyle {\boldsymbol {\Theta }}}$、潜变量${\displaystyle \mathbf {Z} }$）及其条件概率分布。之后，变分贝叶斯构建后验概率${\displaystyle p(\mathbf {Z} ,{\boldsymbol {\Theta }}\mid \mathbf {X} )}$，这近似具有可分解分布的基本特性，即是多个独立分布在不相交的未观测变量子集上的积。
2. 将未观测变量划分为多个子集，在其上推导出独立因子。这种方法没有通用的程序，子集太多会使近似结果不佳，子集过少会使整个变分贝叶斯方法变得难以实现。第一种分割方法是将参数和潜变量分开，一般足以产生可操作的结果。划分记作${\displaystyle \mathbf {Z} _{1},\ldots ,\mathbf {Z} _{M}}$。
3. 对给定的划分${\displaystyle \mathbf {Z} _{j}}$，用基本方程${\displaystyle \ln q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )=\operatorname {E} _{i\neq j}[\ln p(\mathbf {Z} ,\mathbf {X} )]+{\text{constant}}}$写出最佳近似分布${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )}$。
4. 用图模式填写联合概率分布公式。任何不涉及${\displaystyle \mathbf {Z} _{j}}$中变量的分量条件分布都可忽略，它们将被折叠到常数项中。
5. 按上面例子，简化公式并应用期望算子。理想情况下这应该简化为不属于${\displaystyle \mathbf {Z} _{j}}$的变量的基本函数的期望（如第一或第二原始矩、对数期望等等）。为让变分贝叶斯方法顺利运行，这些期望值通常应可以解析地表为变量分布的参数和/或超参数的函数。这些期望项对当前划分的变量都是常数。
6. 式对当前划分中变量的函数形式表明了分布的类型。特别地，取公式指数，便得到分布的概率密度函数（PDF）（或至少是与之成正比的函数，带未知归一化常数）。要使方法可操作，应能识别出函数形式所属的已知分布；要将公式转为匹配已知分布的PDF的形式，可能需要大量计算。若能做到，就可根据定义恢复归一化常数，并提取公式中的适当部分得到已知分布的参数方程。
7. 期望都可用非当前划分变量的函数进行解析替换、并将PDF转为可与已知分布相对应的形式，此时结果应是一组方程，将最优参数值表为其他划分变量参数的函数。
8. 若这程序适用于所有划分，就会产生相互依赖的方程组，指明所有参数的最优值。
9. 最大期望算法型程序可为每个参数选取初值，再迭代。每一步中，我们都会在方程中循环，依次更新每个参数。可以确证收敛。

另见

• 变分消息传递：半分贝叶斯推理的模式化算法
• 变分自编码器：属于概率图模式和变分贝叶斯法的ANN
• 最大期望算法：与变分贝叶斯推理特例相对应的一种相关方法
• 广义滤波：非线性状态空间模型的变分滤波框架
• 变分法：数学分析中处理最值化函数的领域