查看“︁变分贝叶斯方法”︁的源代码

{{贝叶斯统计}}
{{Proofreader needed}}
'''变分贝叶斯方法'''是近似[[贝叶斯推理]]与[[机器学习]]中较难[[积分]]的一系列技术，通常用于由观测变量及未知[[参数]]、[[潜变量]]组成的复杂[[统计模型]]，之间存在的各种关系可用[[图模式]]描述。贝叶斯推理中，参数和潜变量一般归为“未观察变量”。变分贝叶斯方法主要用于两个目的：
# 提供未观测变量[[后验概率]]的分析近似，以便对变量[[统计推断]]
# 得出观测数据的[[边缘似然]]的[[下界]]（即给定模型的数据的[[边缘分布|边缘概率]]，边缘化未观测变量）。通常用于[[模型选择]]，一般認為给定模型的边缘似然越高，拟合效果越好，产生数据的可能性越大（另见[[贝叶斯因子]]）。

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

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

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

==数学推导==
===问题===
[[变分法|变分]]推断中，数据<math>\mathbf{X}</math>的一组未观测变量<math>\mathbf{Z} = \{Z_1 \dots Z_n\}</math>的后验分布近似于所谓'''变分分布'''<math>Q(\mathbf{Z})</math>
: <math>P(\mathbf{Z}\mid \mathbf{X}) \approx Q(\mathbf{Z}).</math>

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

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

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

:<math>D_{\mathrm{KL}}(Q \parallel P) \triangleq \sum_\mathbf{Z}  Q(\mathbf{Z}) \log \frac{Q(\mathbf{Z})}{P(\mathbf{Z}\mid \mathbf{X})}.</math>

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

===难解性===
变分技术通常用于近似以下参数：

:<math>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'}</math>

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

===证据下界===
{{Main|证据下界}}
考虑到<math>P(\mathbf Z \mid \mathbf X) = \frac{P(\mathbf X, \mathbf Z)}{P(\mathbf X)}</math>，KL散度也可写作

:<math>
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]
</math>

因为<math>Q(\mathbf{Z})</math>是分布，所以<math>\sum_\mathbf{Z} Q(\mathbf{Z}) = 1</math>；又由于<math>P(\mathbf{X})</math>对于<math>\mathbf Z</math>是常数，所以
:<math>
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})
</math>

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

:<math>
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})
</math>

重排为
:<math>
\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)
</math>

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

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

===证明===
根据[[布雷格曼散度]]（KL散度是其特例）的广义勾股定理，可以证明<ref name=Tran2018>{{cite arXiv|title=Copula Variational Bayes inference via information geometry|first1=Viet Hung|last1=Tran|year=2018|eprint=1803.10998|class=cs.IT}}</ref><ref name="Martin2014"/>
[[File:Bregman_divergence_Pythagorean.png|right|300px|thumb|[[布雷格曼散度]]的广义勾股定理<ref name="Martin2014">{{cite journal |last1=Adamčík |first1=Martin |title=The Information Geometry of Bregman Divergences and Some Applications in Multi-Expert Reasoning |journal=Entropy |date=2014 |volume=16 |issue=12 |pages=6338–6381|bibcode=2014Entrp..16.6338A |doi=10.3390/e16126338 |doi-access=free }}</ref>]]
:<math>
D_{\mathrm{KL}}(Q\parallel P) \geq D_{\mathrm{KL}}(Q\parallel Q^{*}) + D_{\mathrm{KL}}(Q^{*}\parallel P), \forall Q^{*} \in\mathcal{C}
</math>

其中<math>\mathcal{C}</math>是凸集，且

:<math> Q = Q^{*} \triangleq \arg\min_{Q\in\mathcal{C}}D_{\mathrm{KL}}(Q\parallel P). </math>

时等号成立。这时，全局最小值<math>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)</math>，其中<math>\mathbf{Z}=\{\mathbf{Z_1},\mathbf{Z_2}\}</math>可以如下求得：<ref name=Tran2018/>

:<math> 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),</math>

其中归一化常数为

:<math>\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).</math>

<math>\zeta(\mathbf{X})</math>项在实践中常称作证据下界（ELBO），因为<math>P(\mathbf{X})\geq\zeta(\mathbf{X})=\exp(\mathcal{L}(Q^{*}))</math>，<ref name=Tran2018/>如上所示。

交换<math>\mathbf{Z}_1,\ \mathbf{Z}_2</math>的角色，可以迭代计算真实模型边际<math>P(\mathbf{Z}_1\mid\mathbf{X}),\ P(\mathbf{Z}_2\mid\mathbf{X})</math>的近似<math>q^{*}(\mathbf{Z}_1),\ q^{*}(\mathbf{Z}_2)</math>。这种迭代保证单调收敛，<ref name=Tran2018/>但收敛到的<math>Q^{*}</math>只是<math>D_{\mathrm{KL}}(Q\parallel P)</math>的局部极小值。

若约束空间<math>\mathcal{C}</math>被限制在独立空间内，即<math>q^{*}(\mathbf{Z}_1\mid\mathbf{Z}_2) = q^{*}(\mathbf{Z_1})</math>，则上述迭代将成为所谓平均场近似（mean field approximation）<math>Q^{*}(\mathbf{Z}) = q^{*}(\mathbf{Z}_1)q^{*}(\mathbf{Z}_2),</math>。

==平均场近似==
变分分布<math>Q(\mathbf{Z})</math>常被假定为潜变量的某[[集合划分|划分]]的因式分解，即对潜变量<math>\mathbf{Z}</math>的划分<math>\mathbf{Z}_1 \dots \mathbf{Z}_M</math>，

:<math>Q(\mathbf{Z}) = \prod_{i=1}^M q_i(\mathbf{Z}_i\mid \mathbf{X})</math>

由[[变分法]]可以证明，对每个因子<math>q_j</math>的最佳分布<math>q_j^{*}</math>（就KL散度最小的分布而言，如上所述）满足<ref name="Nguyen">{{cite web|last=Nguyen|first=Duy|title=AN IN DEPTH INTRODUCTION TO VARIATIONAL BAYES NOTE|date=2023-08-15|doi=10.2139/ssrn.4541076|ssrn=4541076|url=https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076|access-date=2023-08-15|archive-date=2024-05-17|archive-url=https://web.archive.org/web/20240517013410/https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076|dead-url=no}}</ref>

:<math>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}</math>

其中<math>\operatorname{E}_{q^*_{-j}} [\ln p(\mathbf{Z}, \mathbf{X})]</math>是数据与潜变量的对数[[联合分布|联合概率]]的[[期望值]]，是相对于不属于划分的所有变量的<math>q^*</math>而取。关于<math>q_j^{*}(\mathbf{Z}_j\mid \mathbf{X})</math>的推导见<ref name=Yoon2021>{{Cite journal |last=Lee|first=Se Yoon|  title = Gibbs sampler and coordinate ascent variational inference: A set-theoretical review|journal=Communications in Statistics - Theory and Methods|year=2021|volume=51 |issue=6 |pages=1–21|doi=10.1080/03610926.2021.1921214|arxiv=2008.01006|s2cid=220935477}}</ref>的引理 4.1。

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

:<math>\ln q_j^{*}(\mathbf{Z}_j\mid \mathbf{X}) = \operatorname{E}_{q^*_{-j}} [\ln p(\mathbf{Z}, \mathbf{X})] + \text{constant}</math>

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

利用期望的性质，式<math>\operatorname{E}_{q^*_{-j}} [\ln p(\mathbf{Z}, \mathbf{X})]</math>通常可简为潜变量[[先验分布]]的固定[[超参数]]及不属于当前划分的潜变量（即不属于<math>\mathbf{Z}_j</math>的潜变量）的期望值的函数（有时还包括[[方差]]等高阶[[矩 (数学)|矩]]）。这就在某一划分的变量分布参数同其他划分的变量的期望值之间产生了[[循环依赖]]，自然就需要类似[[期望最大化算法]]的[[迭代]]法，以某种方式（也许随机）初始化潜变量期望（可能还有高阶矩），再利用当前期望依次计算每个分布的参数、适当设置新算出的分布的期望。这种算法保证收敛。<ref>{{cite book|title=Convex Optimization|first1=Stephen P.|last1=Boyd|first2=Lieven|last2=Vandenberghe|year=2004|publisher=Cambridge University Press|isbn=978-0-521-83378-3|url=https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf|access-date=2011-10-15|archive-date=2021-05-09|archive-url=https://web.archive.org/web/20210509212522/https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf|dead-url=no}}</ref>

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

==变分推理对偶公式==
[[File:CAVI algorithm explain.jpg|600px|thumb|right|用对偶公式图解坐标上升变分推理算法<ref name=Yoon2021/>]]

下面的定理称作变分推理对偶公式，<ref name=Yoon2021/>解释了变分贝叶斯方法中变分分布的一些重要性质。

{{EquationRef|3|Theorem}} 考虑两[[概率空间]]<math>(\Theta,\mathcal{F},P),\ (\Theta,\mathcal{F},Q)</math>，其中<math>Q \ll P</math>。假设存在共同的主导[[概率测度]]<math>\lambda</math>使<math>P \ll \lambda,\ Q \ll \lambda</math>。令''h''是<math>(\Theta,\mathcal{F},P)</math>上的任意实值[[随机变量]]，满足<math>h \in L_1(P)</math>，则下列等式成立

:<math> \log E_P[\exp h] = \text{sup}_{Q \ll P}  \{ E_Q[h] - D_\text{KL}(Q \parallel P)\}.</math>

此外，当且仅当（supremum）

:<math> \frac{q(\theta)}{p(\theta)} = \frac{\exp h(\theta)}{E_P[\exp h]},</math>

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

==基本例子==
考虑简单的不分层贝叶斯模型，由来自[[正态分布]]（[[均值]]、[[方差]]未知）的[[独立同分布]]观测量组成。<ref name=bishop>{{Cite book |last=Bishop |first=Christopher M. |chapter=Chapter 10 |title=Pattern Recognition and Machine Learning |url=https://archive.org/details/patternrecogniti0000bish |year=2006 |publisher=Springer |isbn=978-0-387-31073-2 |author-link=Christopher M. Bishop}}</ref>下面我们将详细介绍这模型，以说明变分贝叶斯方法的工作原理。

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

===数学模型===
将[[共轭先验]]分布置于未知均值<math>\mu</math>、精度<math>\tau</math>，即均值也遵循高斯分布，而精度遵循[[伽马分布]]。

:<math>
\begin{align}
\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{align}
</math>

先验分布[[超参数]]<math>\mu_0, \lambda_0, a_0</math>、<math>b_0</math>是已知定值，可设为小正数，以得到较宽的先验分布，表明我们对<math>\mu</math>、<math>\tau</math>的先验分布一无所知。

已知''N''个数据点<math>\mathbf{X} = \{x_1, \ldots, x_N\}</math>，而我们的目标是推断参数<math>\mu,\ \tau</math>的[[后验分布]]<math>q(\mu, \tau)=p(\mu,\tau\mid x_1, \ldots, x_N)</math>。

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

:<math>p(\mathbf{X},\mu,\tau) = p(\mathbf{X}\mid \mu,\tau) p(\mu\mid \tau) p(\tau)</math>

个体因子是

:<math>
\begin{align}
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{align}
</math>

其中

:<math>
\begin{align}
\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{align}
</math>

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

===导出{{math|''q''(''&mu;'')}}===
则

:<math>
\begin{align}
\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{align}
</math>

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

最后一行是<math>\mu</math>的二次多项式。由于这是<math>q_\mu^*(\mu)</math>的对数，我们可以看到<math>q_\mu^*(\mu)</math>本身是[[正态分布]]。

通过一些繁琐的计算（展开括号里的平方、分离出涉及<math>\mu</math>、<math>\mu^2</math>的项、对<math>\mu</math>应用[[配方法]]），可得到正态分布的参数：
:<math>\begin{align}
\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{align}</math>

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

即

:<math>
\begin{align}
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{align}
</math>

===推导{{math|q(&tau;)}}===

<math>q_\tau^*(\tau)</math>的推导大致相同，简洁起见略去细节。

:<math>
\begin{align}
\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{align}
</math>

两侧取指数，可见<math>q_\tau^*(\tau)</math>服从[[伽马分布]]。

:<math>
\begin{align}
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{align}
</math>

===计算参数===
回忆前几节的结论：

:<math>
\begin{align}
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{align}
</math>
:<math>
\begin{align}
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{align}
</math>

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

:<math>
\begin{align}
\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{align}
</math>

大多数时候将这些公式应用到上述方程不困难，但<math>b_N</math>需要更多工作：

:<math>
\begin{align}
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{align}
</math>

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

:<math>\begin{align}
\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{align}</math>

注意<math>\lambda_N</math>和<math>b_N</math>存在相互依赖关系，自然产生了类似[[最大期望算法]]的算法：

# 计算<math>\sum_{n=1}^N x_n,\ \sum_{n=1}^N x_n^2.</math>用所得值计算<math>\mu_N,\ a_N.</math>
# 初始化<math>\lambda_N</math>为任意值。
# 用<math>\lambda_N</math>的现值和其它参数的已知值计算<math>b_N</math>。
# 用<math>b_N</math>的现值和其它参数的已知值计算<math>\lambda_N</math>。
# 重复最后两部，直到收敛（即直到两值的更新量不超过阈值）。

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

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

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

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

==参考文献==
{{reflist}}

==外部链接==
* [https://www.inference.phy.cam.ac.uk/mackay/itila/ The on-line textbook: Information Theory, Inference, and Learning Algorithms] {{Wayback|url=https://www.inference.phy.cam.ac.uk/mackay/itila/ |date=20170512025952 }}, by [[David J.C. MacKay]] provides an introduction to variational methods (p.&nbsp;422).
* [https://www.robots.ox.ac.uk/~sjrob/Pubs/fox_vbtut.pdf A Tutorial on Variational Bayes] {{Wayback|url=https://www.robots.ox.ac.uk/~sjrob/Pubs/fox_vbtut.pdf |date=20240527220106 }}. Fox, C. and Roberts, S. 2012. Artificial Intelligence Review, {{doi|10.1007/s10462-011-9236-8}}.
* [https://www.gatsby.ucl.ac.uk/vbayes/ Variational-Bayes Repository] A repository of research papers, software, and links related to the use of variational methods for approximate Bayesian learning up to 2003.
* [https://www.cse.buffalo.edu/faculty/mbeal/thesis/index.html Variational Algorithms for Approximate Bayesian Inference] {{Wayback|url=https://www.cse.buffalo.edu/faculty/mbeal/thesis/index.html |date=20240527220112 }}, by M. J. Beal includes comparisons of EM to Variational Bayesian EM and derivations of several models including Variational Bayesian HMMs.
* [https://www.cs.jhu.edu/~jason/tutorials/variational.html High-Level Explanation of Variational Inference] {{Wayback|url=https://www.cs.jhu.edu/~jason/tutorials/variational.html |date=20240603140130 }} by Jason Eisner may be worth reading before a more mathematically detailed treatment.
* [https://arxiv.org/abs/1803.10998 Copula Variational Bayes inference via information geometry (pdf)] {{Wayback|url=https://arxiv.org/abs/1803.10998 |date=20221022060746 }} by Tran, V.H. 2018. This paper is primarily written for students. Via [[Bregman divergence]], the paper shows that Variational Bayes is simply a generalized Pythagorean projection of true model onto an arbitrarily correlated (copula) distributional space, of which the independent space is merely a special case.
* [https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076 An in depth introduction to Variational Bayes note] {{Wayback|url=https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076 |date=20240612080933 }}. Nguyen, D.  2023

{{DEFAULTSORT:Variational Bayesian Methods}}
[[Category:贝叶斯统计]]