PRML笔记（十）

10. Approximate Inference

在probabilistic models中的一个核心任务是，在给定observed（visible）data variables $\mathbf{X}$ 的时候去计算关于latent variables $\mathbf{Z}$的posterior distribution $p(\mathbf{Z}|\mathbf{X})$。并且去在该概率分布下计算一些expectations。这个model可能会包含一些deterministic parameters，我们此时将设这些parameters为implicit，或是说，它会是一个fully Bayesian model，其中任何的unknown parameters都会有相应的prior distribution，并且这些parameters将会容纳到由vector $\mathbf{Z}$所标记的latent variables中去。例如，在EM算法中，我们需要计算在latent variables的posterior distribution下的complete-data log likelihood的expectation。对于许多具有实际应用的models，我们无法计算得到相应的posterior distribution，或者去计算关于这些distribution的expectations。这可能是由于latent space的dimensionality过于高，导致我们无法直接处理它，或者是由于posterior distribution具有一种非常复杂的形式，从而使得这里的expectations无法求得解析解。在continuous variables的情况中，相应的integrations可能不具有closed-form analytical solutions，同时space的维数以及被积函数的复杂度可能会抑制numerical integration。对于discrete variables来说，marginalizations涉及到对关于hidden variables的所有可能取值的求和，从原则上来说，尽管它总是可能的，但是我们在实际中经常会发现，这里的hidden states是指数级别的，从而使得exact calculation的代价十分大。

在这样一种情况下，我们需要将目光转向一种approximation schemes，这样的schemes主要分为两类（根据具体的方法是依赖于stochastic或是deterministic approximations）。第11章中所叙述的stochastic 方法（例如Markov chain Monte Carlo）是一种在许多领域中都广泛使用的Bayesian methods。它们普遍具有一个性质，即在给定无限的计算资源的时候，该类方法可以得到exact results，而当我们使用a finite amount of processor time的时候，我们将得到一个approximation。在实际中，sampling methods非常消耗时间，从而使得该类方法通常被局限于small-scale problems中。同样的，我们难以知道，一个sampling scheme是否能够从required distribution中生成independent samples。

在本章中，我们引入一系列的deterministic approximation schemes，其中一些能够较好地扩展到大规模的应用中。这些方法是基于对posterior distribution的analytical approximations，例如通过假设它可能以一种特定的方式进行分解，或者它具有一种特定的形式（例如Gaussian）。同样的，它们不能生成exact results，因此它们的strengths和weaknesses与sampling methods是互补的。

在4.4节中，我们讨论了Laplace approximation，其中它是利用local Gaussian approximation来逼近一个distribution的mode而实现的。在这里，我们转而考虑另外一类近似方法，被称为是variational inference（或是variational Bayes），其中我们使用更为global的criteria，并且这类方法已经得到了广泛的使用。最后我们将介绍另外一种variational framework，被称为是expectation propagation。

10.1 Variational Inference

Variational methods源自18世纪Euler，Lagrange和其他人关于变分法的研究。标准的微积分关注于寻找functions的derivatives。我们可以视一个function为一个mapping，它将一个variable的value当作input，然后返回该function的value，作为output。因此function的derivative描述了当我们对input value进行了无限小的改变之后，output value是如何变化的。类似的，我们定义一个functional为，一个mapping，它以一个function为input，然后返回该functional的value作为output。一个关于functional的例子为entropy $H[p]$, 其中它以一个probability distribution $p(x)$ 为input，然后它返回如下的这个量：
$$H[p]=\int p(x)\text{ln}p(x)\text{d}x(10.1)$$
作为output。然后我们可以引入functional derivative的概念，其中它表达的是，当input function发生微小的改变的时候，functional的值是如何改变的。变分法的规则与标准微积分规则的表述在附录D中有所介绍。许多问题可以被表示成一个optimization problem，其中我们所要去优化的quantity为一个functional。我们可以通过探索所有可能的input function，并且寻找一个可能最大化或最小化functional的solution。Variational methods具有较为广泛的应用，并且在有限元方法和最大熵中都有所应用。

尽管变分法本质上没有什么近似的特性，但是该方法天然适合寻找近似解。这可以通过约束我们所需要优化的函数范围的方式来实现。例如，我们仅仅考虑二次函数，或是只考虑具有固定基函数的线性组合形式的函数，其中线性组合系数是可以变化的。在probabilistic inference的应用中，restriction可能需要具有factorization assumption的假设。

现在我们来从细节方面考虑variational optimization的概念可以以何种方式应用在inference problem中。假设我们具有一个fully Bayesian model，其中所有的parameters都具有相应的prior distributions。这个model可能具有latent variables以及parameters，并且我们将所有的observed variables的set表示为$\mathbf{X}$。例如，我们可能具有一个含有$N$个独立同分布data points的data set，其中有 X = { x 1 , . . . , x N } \mathbf{X}=\{\mathbf{x}_1,...,\mathbf{x}_N\} X={ x1​,...,xN​}, 并且有 Z = { z 1 , . . . , z N } \mathbf{Z}=\{\mathbf{z}_1,...,\mathbf{z}_N\} Z={ z1​,...,zN​}。我们的probabilistic model对joint distribution $p(\mathbf{X},\mathbf{Z})$进行了明确，并且我们的目标在于寻找到一个关于posterior distribution $p(\mathbf{Z}|\mathbf{X})$的approximation，以及model evidence $p(\mathbf{X})$。正如我们在EM算法中的讨论一样，我们可以使用如下的方式对log likelihood probability进行分解：
$$\text{ln}p(\mathbf{X})=\mathcal{L}(q)+\text{KL}(q||p)(10.2)$$
其中我们有定义：
$$\mathcal{L}(q)=\int q(\mathbf{Z})\text{ln}\left\{\frac{p(\mathbf{X},\mathbf{Z})}{q(\mathbf{Z})}\right\}\text{d}\mathbf{Z}(10.3)$$
$$\text{KL}(q||p)=-\int q(\mathbf{Z})\text{ln}\left\{\frac{p(\mathbf{Z}|\mathbf{X})}{q(\mathbf{Z})}\right\}\text{d}\mathbf{Z}(10.4)$$
这与我们对EM算法的讨论的不同之处在于，parameter vector $\mathbf{\theta }$不再出现了，因为此时parameters都是stochastic variables，因而被包含在$\mathbf{Z}$中。因为在本章中，我们将主要考虑continuous variables，所以我们使用integrations而不是summations来进行这里的分解。然后，如果这里的一些或所有的variables都是discrete，那么我们可以将integrations代替为summations的方式来使得这里的分析不发生改变。与之前相同的是，我们可以通过优化distribution $q(\mathbf{Z})$的方式，来对lower bound $\mathcal{L}(q)$进行maximize，其中它等价于最小化KL divergence。如果我们允许任何可能的关于$q(\mathbf{Z})$的选择，那么对lower bound的最大化将在KL divergence消失的时候发生，该情况在$q(\mathbf{Z})$与posterior distribution $p(\mathbf{Z}|\mathbf{X})$相等时发生。然而，我们假设，直接计算真实的posterior distribution是intractable的。

因此，我们考虑一类受限的distributions $q(\mathbf{Z})$, 然后寻找到这类distributions中那个能使得KL divergence最小的distribution。我们的目标在于对概率分布族进行约束，使得该族中仅包含tractable distributions，同时允许该概率分布族能够充分多样化，并且足够灵活，从而能够为真实的posterior distribution提供一个较好的approximation。我们需要强调的是，restriction的作用仅仅是为了能使distribution具有tractability，在受到这样的约束的前提下，我们应该尽可能使用一族含有概率分布形式数量丰富的approximating distributions。特别的，在这里，较高灵活度的distributions并不与‘over-fitting’产生关系。使用更为灵活的approximations仅仅使得我们可以更为接近真实的posterior distribution。

一种限制approximating distributions的方法是使用parametric distribution $q(\mathbf{Z}|\mathbf{w})$, 它受到着一系列parameters $\mathbf{w}$的控制。此时lower bound $\mathcal{L}(q)$将会成为是一个关于$\mathbf{w}$的function，并且我们可以利用标准的非线性优化方法来确定parameters的最优值。一个关于该方法的例子（其中variational distribution是Gaussian）如下图所示，其中我们对概率分布的的mean和variance进行了优化：

10.1.1 Factorized distributions

在这里，我们考虑另外一种方法去对概率分布族 $q(\mathbf{Z})$ 进行限制。假设我们将$\mathbf{Z}$中的elements划分成disjointed groups，并被标记为${\mathbf{Z}}_i,i=1,...,M$。此时我们假设distribution $q$可以依照这些groups而进行划分，从而有：
$$q(\mathbf{Z})=\prod _{i=1}^M{q}_i({\mathbf{Z}}_i)(10.5)$$
我们需要强调的是，我们并没有对distributions进行进一步的假设。特别的，我们并没有对每一个factor $q_i({\mathbf{Z}}_i)$进行限制。这里的factorized form of variational inference对应于一种approximation framework，被称为是mean field theory。

在所有的具有(10.5)形式的$q(\mathbf{Z})$ distributions中，我们现在需要寻找的是一个能够使得lower bound $\mathcal{L}(q)$最大的distribution (注意到，这里我们并没有对KL divergence直接进行最小化，因为我们无法求得该项的表达式，所以就没法使用变分法来求得$q(\mathbf{Z})$的值)。因此，我们希望能够构造一个关于所有distributions $q_i({\mathbf{Z}}_i)$ 的free form (variational) 的$\mathcal{L}(q)$的optimization，其中我们轮流对每一个factor进行优化。为了实现这件事，我们首先将(10.5)代入(10.3)，然后分析出关于factor $q_j({\mathbf{Z}}_j)$的依赖。我们将$q_j({\mathbf{Z}}_j)$简单标记为$q_j$, 从而简化符号，因而我们可以得到：
$$\begin{gathered} \mathcal{L}(q)=\int \prod _i{q}_i\left\{\text{ln}p(\mathbf{X},\mathbf{Z})-\sum _i\text{ln}{q}_i\right\}\text{d}\mathbf{Z} \\ =\int q_i\left\{\int \text{ln}p(\mathbf{X},\mathbf{Z})\prod _{i\ne j}{q}_i\text{d}{\mathbf{Z}}_i\right\}\text{d}{\mathbf{Z}}_j-\int q_j\text{ln}{q}_j\text{d}{\mathbf{Z}}_j+\text{const} \\ =\int q_j\text{ln}\tilde{p}({\mathbf{X},\mathbf{Z}}_j)\text{d}{\mathbf{Z}}_j-\int q_j\text{ln}{q}_j\text{d}{\mathbf{Z}}_j+\text{const}(10.6) \end{gathered}$$

其中我们通过如下的关系定义了一个新的distribution $\tilde{p}(\mathbf{X},{\mathbf{Z}}_j)$：
$$\text{ln}\tilde{p}(\mathbf{X},{\mathbf{Z}}_j)={\mathbb{E}}_{i\ne j}[\text{ln}p(\mathbf{X},\mathbf{Z})]+\text{const}(10.7)$$
在这里，符号${\mathbb{E}}_{i\ne j}[\dots ]$表示q distributions在所有variables ${\mathbf{z}}_i,i\ne j$上的expectation，从而有：
$${\mathbb{E}}_{i\ne j}[\text{ln}p(\mathbf{X},\mathbf{Z})]=\int \text{ln}p(\mathbf{X},\mathbf{Z})\prod _{i\ne j}{q}_i\text{d}{\mathbf{Z}}_i(10.8)$$
现在，假设我们保持 { q i ≠ j } \{q_{i\neq j}\} { qi​=j​}是固定的，并且关于所有可能的distribution $q_j({\mathbf{Z}}_j)$的形式，最大化(10.6)中的$\mathcal{L}(q)$。这一件事情其实是很容易的，因为我们发现，(10.6)其实是在$q_j({\mathbf{Z}}_j)$和$\tilde{p}(\mathbf{X},{\mathbf{Z}}_j)$之间的negative KL divergence。因此，对(10.6)的最大化等价于对KL divergence的最小化，因此，最小值出现在当$q_j({\mathbf{Z}}_j)=\tilde{p}({\mathbf{X},\mathbf{Z}}_j)$满足时。因此，我们得到了一个$q_j^{\ast }({\mathbf{Z}}_j)$最优解的一般性表达：
$$\text{ln}{q}_j^{\ast }({\mathbf{Z}}_j)={\mathbb{E}}_{i\ne j}[\text{ln}p(\mathbf{X},\mathbf{Z})]+\text{const}(10.9)$$
我们花费一些时间来研究这个表达式是十分值得的，因为它提供了variational method的一种应用。所以说，对于factor $q_j$来说，它的最优解的log形式可以简单地通过考虑在所有的hidden以及visible variables的joint distribution的log的形式，然后对其关于所有的其他factors { q i } , i ≠ j \{q_i\},i\neq j { qi​},i​=j求期望的方式得到。

在(10.9)中的addictive constant是通过对distribution $q_j^{\ast }({\mathbf{Z}}_j)$进行normalizing而得到的。因此，如果我们对上式两边同时取指数，然后进行归一化，则有：
$$q_j^{\ast }({\mathbf{Z}}_j)=\frac{\text{exp}({\mathbb{E}}_{i\ne j}[\text{ln}p(\mathbf{X},\mathbf{Z})])}{\int \text{exp}({\mathbb{E}}_{i\ne j}[\text{ln}p(\mathbf{X},\mathbf{Z})])\text{d}{\mathbf{Z}}_j}$$

在实际中，我们将考虑一种更为简单的方式，来处理(10.9), 然后（在需要的时候）通过观察的方式得到normalization constant。在接下来的例子中，这件事将变得更为清楚。

由(10.9)所给出的一系列的equations（$j=1,...,M$）表达着一系列的在受到着factorization constraint的时候的maximum of the lower bound的一致性条件。然后，它们并没有表达成一种显式解的形式，因为(10.9)右式关于$q_j^{\ast }({\mathbf{Z}}_j)$的最优解依赖于在其他factors $q_i({\mathbf{Z}}_i),i\ne j$上的期望。因此我们将探究一种consistent solution，其中我们首先对所有的factors $q_i({\mathbf{Z}}_i)$进行适当的初始化，然后在这些factors上进行循环，并利用由(10.9)右式中关于所有其他的factors相应的current estimates所计算出来的revised estimate，来替换原来的estimate。由于与每一个factor $q_i({\mathbf{Z}}_i)$相关的bound是convex的，所以convergence是可以得到保证的。

10.1.2 Properties of factorized approximations

我们关于variational inference的方法是基于对true posterior distribution的factorized approximation。我们首先考虑用一个factorized distribution来逼近一个general distribution的问题。首先，我们讨论使用一个factorized Gaussian来逼近一个Gaussian distribution的问题，其中我们将深入理解在使用factorized approximation时会带来的inaccuracy的类型。考虑一个关于两个correlated variables $\mathbf{z}=(z_1,z_2)$的Gaussian distribution $p(\mathbf{z})=\mathcal{N}(\mathbf{z}|\mathbf{\mu },{\mathbf{\Lambda }}^{-1})$，其中mean和precision为：
$$\mathbf{\mu }=\left(\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right),\mathbf{\Lambda }=\left(\begin{array}{cc}{\Lambda }_{11}& {\Lambda }_{12}\\ {\Lambda }_{21}& {\Lambda }_{22}\end{array}\right)(10.10)$$
其中由于precision matrix具有对称性，所以有${\Lambda }_{21}={\Lambda }_{12}$。现在，假设我们希望去利用一个具有$q(\mathbf{z})=q_1(z_1)q_2(z_2)$形式的factorized Gaussian去近似这个概率分布。我们首先利用(10.9)所述的general result，从而找到一个关于optimal factor $q_1^{\ast }(z_1)$的expression。为了实现这件事，我们注意到，在该等式的右侧，我们仅仅需要保留那些依赖于$z_1$的项，因为所有其他的项将被容纳在normalization constant中。因此，我们有：
$$\begin{gathered} \text{ln}{q}_1^{\ast }(z_1)={\mathbb{E}}_{z_2}[\text{ln}p(\mathbf{z})]+\text{const} \\ ={\mathbb{E}}_{z_2}\left[-\frac{1}{2}(z_1-{\mu }_1{)}^2{\Lambda }_{11}-(z_1-{\mu }_1){\Lambda }_{12}(z_2-{\mu }_2)\right]+\text{const} \\ =-\frac{1}{2}{z}_1^2{\Lambda }_{11}+z_1{\mu }_1{\Lambda }_{11}-z_1{\Lambda }_{12}(\mathbb{E}[z_2]-{\mu }_2)+\text{const}(10.11) \end{gathered}$$
然后，我们观察到，这个表达式的右侧是一个关于$z_1$的二次函数，因此我们可以将$q^{\ast }(z_i)$视为一个Gaussian。我们需要强调的是，我们并没有假设$q(z_i)$是Gaussian的，而是说我们通过对KL divergence关于所有可能的distribution $q(z_i)$进行variational optimization而得到了这个结果。我们同样注意到，我们其实并不需要显式考虑(10.9)中的addictive constant，因为它表示着我们可以在操作的末尾通过观察而得到的normalization constant。使用配方法的技巧，我们可以识别出这个Gaussian的mean和precision，从而有：
$$q^{\ast }(z_1)=\mathcal{N}(z_1|m_1,{\Lambda }_{11}^{-1})(10.12)$$
其中有：
$$m_1={\mu }_1-{\Lambda }_{11}^{-1}{\Lambda }_{12}(\mathbb{E}[z_2]-{\mu }_2)(10.13)$$
由对称性，$q_2^{\ast }(z_2)$同样是一个Gaussian，并且它可以写成是：
$$q_2^{\ast }(z_2)=\mathcal{N}(z_2|m_2,{\Lambda }_{22}^{-1})(10.14)$$
其中有：
$$m_2={\mu }_2-{\Lambda }_{22}^{-1}{\Lambda }_{21}(\mathbb{E}[z_1]-{\mu }_1)(10.15)$$
我们注意到这里的解都是相互耦合的，从而使得$q^{\ast }(z_1)$依赖于由$q^{\ast }(z_2)$所计算出来的期望值，并且反之亦然。通常来说我们通过将variational solutions看作re-estimation equations，然后循环这些variables并对它们进行更新，直到某些收敛准则得到满足。我们将在之后看到一个相关的例子。然后，在这里，我们注意到这个问题其实是十分简单的，因为相应的closed form solution可以找到。特别的，因为$\mathbb{E}[z_1]=m_1,\mathbb{E}[z_2]=m_2$, 我们可以知道的是，如果我们设$\mathbb{E}[z_1]={\mu }_1,\mathbb{E}[z_2]={\mu }_2$的话，这两个equations可以得到满足，并且在假定该distribution为nonsingular的时候，这其实就是唯一的solution。这样的结果如下图(a)所示：

我们知道的是，这里的mean可以被正确找出，但是$q(\mathbf{z})$的variance是由$p(\mathbf{z})$的最小的variance的direction所控制的，并且在与之正交的方向上，variance被极大地低估了。这其实是一个一般性的结果，即factorized variational approximation倾向于给posterior distribution提供一个过于compact的approximation。

为了进行比较，假设我们已经对reverse KL divergence $\text{KL}(p||q)$进行了最小化。我们将会看到的是，这种形式的KL divergence在另外一种被称为是expectation propagation的approximate inference framework中会被使用出来。因此我们考虑一个一般性的问题，即当$q(\mathbf{Z})$是通过(10.5)的形式来作为一个factorized approximation的时候，我们对$\text{KL}(p||q)$进行最小化。此时KL divergence可以被写成是如下的形式：
$$\text{KL}(p||q)=-\int p(\mathbf{Z})\left[\sum _{i=1}^M\text{ln}{q}_i({\mathbf{Z}}_i)\right]\text{d}\mathbf{Z}+\text{const}(10.16)$$
其中constant term仅仅是$p(\mathbf{Z})$的entropy，因为它并不依赖于$q(\mathbf{Z})$。此时我们可以关于每一个factors $q_j({\mathbf{Z}}_j)$来进行优化，其中它可以利用Lagrange multiplier而轻易地得到：
$$q_j^{\ast }({\mathbf{Z}}_j)=\int p(\mathbf{Z})\prod _{i\ne j}\text{d}{\mathbf{Z}}_i=p({\mathbf{Z}}_j)(10.17)$$
在这种情况中，我们可以知道的是，关于$q_j({\mathbf{Z}}_j)$的最优结果可以通过相应的关于$p(\mathbf{Z})$的marginal distribution而给出。注意到这其实是一个closed-form solution，因此它并不需要任何的iteration。

为了对将这个结果应用在关于vector $\mathbf{z}$的Gaussian distribution $p(\mathbf{z})$上的情况进行说明，我们可以利用(2.98), 其中这个结果在图10.2(b)中体现了出来。我们再一次发现了关于这个approximation的mean是正确的，但是它将大量的probability mass放置在了具有非常低概率值的variable space region中。这两个结果之间的区别为：我们注意到KL divergence中，就算是region of $\mathbf{Z}$ space中$p(\mathbf{z})$接近于0，但是如果$q(\mathbf{Z})$并不接近于0的话，那么相应的KL divergence中就会有一个large positive contribution：
$$\text{KL}(q||p)=-\int q(\mathbf{Z})\text{ln}\left\{\frac{p(\mathbf{Z})}{q(\mathbf{Z})}\right\}\text{d}\mathbf{Z}(10.18)$$
因此，对这样一种形式的KL divergence的最小化会导致$q(\mathbf{Z})$忽略掉那些$p(\mathbf{Z})$比较小的regions。反过来说，KL divergence $\text{KL}(p||q)$关于$q(\mathbf{Z})$的最小化在$q(\mathbf{Z})$非0且$p(\mathbf{Z})$非零的区域内。

如果我们考虑利用一个unimodal distribution来拟合一个multimodal distribution的话，我们就可以对这两个KL divergence之间的区别有一个深入的理解，如下图所示：

在实际的应用中，真实的posterior distribution往往是multimodal的，其中大部分的posterior mass集中在一些相对较小的regions of parameter space。这些multiple modes可能来源于latent space的不可识别性，或是来源于parameters的非线性依赖。这两种形式的multimodality都在第9章关于Gaussian mixture的部分遇到了，其中它们以likelihood function中的multiple maxima的形式展现了出来，并且基于minimizing $\text{KL}(q||p)$的一种variational treatment倾向于寻找出这其中的一个mode。反过来说，如果我们要minimize $\text{KL}(p||q)$, 最终结果的approximation将会对所有的modes取平均，并且，在mixture model的情况下，将会导致poor predictive distributions（因为两个较好的参数值的平均一般来说并非是一种较好的参数值）。利用$\text{KL}(p||q)$来定义一种有用的inference procedure是可能的，但是这样的方法与我们在这里所讨论的方法是相当不同的，我们将在讨论expectation propagation的时候详细说明这个问题。

我们所讨论的这两种形式的KL divergence是alpha family of divergences的两个成员，该family定义如下：
$${\text{D}}_{\alpha }(p||q)=\frac{4}{1-{\alpha }^2}\left(1-\int p(x{)}^{(1+\alpha )/2}q(x{)}^{(1-\alpha )/2}\text{d}x\right)(10.19)$$
其中$-\infty <\alpha <\infty$是一个continuous parameter。KL divergence $\text{KL}(p||q)$对应于极限$\alpha \to 1$, 然而$\text{KL}(q||p)$对应于极限$\alpha \to -1$。对于$\alpha$的所有值，我们有${\text{D}}_{\alpha }(p||q)\ge 0$, 其中等号成立当且仅当$p(x)=q(x)$。假设$p(x)$是一个fixed distribution，并且我们通过在一个set of distributions $q(x)$中寻找一个distribution，来最小化${\text{D}}_{\alpha }(p||q)$。那么，对于$\alpha \le -1$来说，divergence是zero forcing的，即任何能够使得$p(x)=0$成立的$x$的值，都会使得$q(x)=0$成立，此时有$q(x)$会underestimate $p(x)$的support，并且会倾向于去寻找那个具有largest mass的mode。反过来说，对于$\alpha \ge 1$来说，这里的divergence是zero-avoiding的，即满足$p(x)>0$的value也同样会使得$q(x)>0$，此时$q(x)$将会覆盖$p(x)$的所有区域，从而致使对support of $p(x)$的高估。当$\alpha =0$的时候，我们可以得到一个symmetric divergence，它与Hellinger distance之间是线性相关的：
$${\text{D}}_{\text{H}}(p||q)=\int \left(p(x{)}^{1/2}-q(x{)}^{1/2}\right)\text{d}x(10.20)$$
Hellinger distance的平方根是一个valid distance metric。

10.1.3 Example: The univariate Gaussian

此时我们利用一个在single variable $x$上的Gaussian distribution来对factorized Variational approximation进行说明。我们的目标在于，在给定一个关于observed values $x$的data set D = { x 1 , . . . , x N } \mathcal{D}=\{x_1,...,x_N\} D={ x1​,...,xN​}（我们假设这些data points都是从一个Gaussian distribution中独立同分布地采样出来的）的情况下，对mean $\mu$和precision $\tau$的posterior distribution进行推断。此时likelihood function为：
$$p(\mathcal{D}|\mu ,\tau )={\left(\frac{\tau }{2\pi }\right)}^{N/2}\text{exp}\left\{-\frac{\tau }{2}\sum _{n=1}^N(x_n-\mu {)}^2\right\}(10.21)$$
此时我们引入关于$\mu$和$\tau$的conjugate prior为：
$$p(\mu |\tau )=\mathcal{N}(\mu |{\mu }_0,({\lambda }_0\tau {)}^{-1})(10.22)$$
$$p(\tau )=\text{Gam}(\tau |a_0,b_0)(10.23)$$
其中$\text{Gam}(\tau |a_0,b_0)$是gamma distribution（由(2.146)定义）。这些distributions组成了一个Gaussian-Gamma conjugate prior distribution。

对于这个简单的问题来说，posterior distribution可以被精确地找出，同样的，它们也是Gaussian-gamma distribution。然而，为了教学的目的，我们将考虑一个关于该posterior distribution的factorized variational approximation：
$$q(\mu ,\tau )=q_{\mu }(\mu )q_{\tau }(\tau )(10.24)$$
我们注意到，真实的posterior distribution不是以这种方式进行分解的。而factors的最优解$q_{\mu }(\mu )$和$q_{\tau }(\tau )$可以从一般性的结果(10.9)中，以如下的方式得到。对于$q_{\mu }(\mu )$, 我们有：
$$\begin{gathered} \text{ln}{q}_{\mu }^{\ast }(\mu )={\mathbb{E}}_{\tau }[\text{ln}p(\mathcal{D}|\mu ,\tau )+\text{ln}p(\mu |\tau )]+\text{const} \\ =-\frac{\mathbb{E}[\tau ]}{2}\left\{{\lambda }_0(\mu -{\mu }_0{)}^2+\sum _{n=1}^N(x_n-\mu {)}^2\right\}+\text{const}(10.25) \end{gathered}$$
对关于$\mu$的项进行配方，我们有$q_{\mu }(\mu )$是一个Gaussian $\mathcal{N}(\mu |{\mu }_N,{\lambda }_N^{-1})$，其中mean和precision为：
$${\mu }_N=\frac{{\lambda }_0{\mu }_0+N\bar{x}}{{\lambda }_0+N}(10.26)$$
$${\lambda }_N=({\lambda }_0+N)\mathbb{E}[\tau ](10.27)$$
注意到对于$N\to \infty$，这提供了一个maximum likelihood result，其中${\mu }_N=\bar{x}$, 并且此时precision是infinite的。

类似的，对于factor $q_{\tau }(\tau )$来说，optimal solution为：
$$\begin{gathered} \text{ln}{q}_{\tau }^{\ast }(\tau )={\mathbb{E}}_{\mu }[\text{ln}p(\mathcal{D}|\mu ,\tau )+\text{ln}p(\mu |\tau )]+\text{ln}p(\tau )+\text{const} \\ =(a_0-1)\text{ln}\tau -b_0\tau +\frac{N}{2}\text{ln}\tau \\ -\frac{\tau }{2}{\mathbb{E}}_{\mu }\left[\sum _{n=1}^N(x_n-\mu {)}^2+{\lambda }_0(\mu -{\mu }_0{)}^2\right]+\text{const}(10.28) \end{gathered}$$
因此，$q_{\tau }(\tau )$是一个gamma distribution $\text{Gam}(\tau |a_N,b_N)$, 其中有参数：
$$a_N=a_0+\frac{N}{2}(10.29)$$
$$b_N=b_0+\frac{1}{2}{\mathbb{E}}_{\mu }\left[\sum _{n=1}^N(x_n-\mu {)}^2+{\lambda }_0(\mu -{\mu }_0{)}^2\right](10.30)$$
同样的，当$N\to \infty$的时候，这也表现出期望的性质。

我们需要强调的是，我们并没有为这里的optimal distributions $q_{\mu }(\mu )$和$q_{\tau }(\tau )$来假设特定的functional form。它们的形式是从structure of the likelihood function以及相应的conjugate priors中自然而然地产生出来的。

因此，我们就已经拥有了关于optimal distributions $q_{\mu }(\mu )$和$q_{\tau }(\tau )$的optimal distributions，其中它们每一个都依赖于另外一个distribution的矩的值。因此，一种寻找solution的方法是首先先进行猜测，例如对$\mathbb{E}[\tau ]$的值进行猜测，然后使用它来对distribution $q_{\mu }(\mu )$来进行re-compute。在给定这个revised distribution之后，我们此时可以计算得到相应的moments $\mathbb{E}[\mu ],\mathbb{E}[{\mu }^2]$，并利用这些结果来对distribution $q_{\tau }(\tau )$进行重新计算，以此类推。因为这个例子中的hidden variables的space是2维的，所以我们可以通过画出这里的true posterior以及factorized approximation的contours的方式来对该variational approximation与相应的posterior distribution的逼近进行说明，如下图所示：

一般来说，我们将需要一种例如上述所示的迭代方式来求解optimal factorized posterior distribution。然而，对于非常简单的例子（如我们在这里讨论的例子）来说，我们可以通过同时求解关于optimal factors $q_{\mu }(\mu )$和$q_{\tau }(\tau )$的方程，从而得到其显式解。在做这件事之前，我们可以通过考虑较为broad，noninformative priors，其中有${\mu }_0=a_0=b_0={\lambda }_0=0$。尽管这些parameter settings对应于improper priors，我们可以知道的是，此时posterior distribution仍然是welldefined。对一个gamma distribution的mean使用标准的结果$\mathbb{E}[\tau ]=a_N/b_N$，以及使用(10.29)和(10.30), 我们有：
$$\frac{1}{\mathbb{E}[\tau ]}=\mathbb{E}\left[\frac{1}{N}\sum _{n=1}^N(x_n-\mu {)}^2\right]=\overline{x^2}-2\bar{x}\mathbb{E}[\mu ]+\mathbb{E}[{\mu }^2](10.31)$$
然后，使用(10.26)以及(10.27), 我们可以得到$q_{\mu }[\mu ]$的一阶矩和二阶矩，如下所示：
$$\mathbb{E}[\mu ]=\bar{x},\mathbb{E}[{\mu }^2]={\bar{x}}^2+\frac{1}{N\mathbb{E}[\tau ]}(10.32)$$
此时我们可以将这些矩代入(10.31)中，然后求解出$\mathbb{E}[\tau ]$:
$$\frac{1}{\mathbb{E}[\tau ]}=\frac{1}{N-1}\sum _{n=1}^N(x_n-\bar{x}{)}^2(10.33)$$
我们可以看出来的是，该式的右侧与一个univariate Gaussian distribution的variance的unbiased estimator是相似的，因此我们可以发现的是，使用一个Bayesian approach可以避免maximum likelihood solution中的bias。

10.1.4 Model comparison

除了在hidden variables $\mathbf{Z}$上进行inference之外，我们同样也希望去对一系列candidate models进行比较，它们被标记为index $m$，并且具有prior probabilities $p(m)$。此时我们的目的在于对posterior probabilities $p(m|\mathbf{X})$进行近似，其中$\mathbf{X}$是observed data。这比我们之前讨论的情况要稍微复杂一些，因为不同的models可能具有不同的structure，以及不同的关于hidden variables $\mathbf{Z}$ 的dimensionality。因此我们不能仅仅考虑分解$q(\mathbf{Z})q(m)$, 而需要明确的是，关于$\mathbf{Z}$的posterior需要conditioned on $m$。因此我们必须要考虑$q(\mathbf{Z},m)=q(\mathbf{Z}|m)q(m)$。此时我们就已经准备好验证如下的基于该variational distribution的decomposition：
$$\text{ln}p(\mathbf{X})=\mathcal{L}-\sum _m\sum _{\mathbf{Z}}q(\mathbf{Z}|m)q(m)\text{ln}\left\{\frac{p(\mathbf{Z},m|\mathbf{X})}{q(\mathbf{Z}|m)q(m)}\right\}(10.34)$$
其中$\mathcal{L}$是$\text{ln}p(\mathbf{X})$的一个lower bound：
$$\mathcal{L}=\sum _m\sum _{\mathbf{Z}}q(\mathbf{Z}|m)q(m)\text{ln}\left\{\frac{p(\mathbf{Z},\mathbf{X},m)}{q(\mathbf{Z}|m)q(m)}\right\}(10.35)$$
其中我们假设$\mathbf{Z}$是discrete的，但是同样的分析应用在continuous latent variables之后会需要我们将summations代替为integrations。此时我们可以对${\mathcal{L}}_m$关于distribution $q(m)$进行最大化（利用Lagrange multiplier）。
为了实现这个目的，我们首先对lower bound $\mathcal{L}$的表达式进行变换：
L = ∑ m ∑ Z q ( Z ∣ m ) q ( m ) ln { p ( Z , X , m ) q ( Z ∣ m ) q ( m ) } = ∑ m ∑ Z q ( Z ∣ m ) q ( m ) { ln  p ( Z , X ∣ m ) + ln  p ( m ) − ln  q ( Z ∣ m ) − ln  q ( m ) } = ∑ m q ( m ) ( ln  p ( m ) − ln  q ( m ) + ∑ Z q ( Z ∣ m ) { ln  p ( Z , X ∣ m ) − ln  q ( Z ∣ m ) } ) = ∑ m q ( m ) { ln ( p ( m ) exp { L m } ) − ln  q ( m ) }      ( 264 ) \mathcal{L}=\sum_m\sum_{\mathbf{Z}}q(\mathbf{Z}|m)q(m)\text{ln}\left\{ \frac{p(\mathbf{Z,X},m)}{q(\mathbf{Z}|m)q(m)} \right\}\\ =\sum_m\sum_{\mathbf{Z}}q(\mathbf{Z}|m)q(m)\{\text{ln}\ p(\mathbf{Z,X}|m)+\text{ln}\ p(m)-\text{ln}\ q(\mathbf{Z}|m)-\text{ln}\ q(m)\}\\ =\sum_m q(m)\left( \text{ln}\ p(m)-\text{ln}\ q(m) \right.\\ \left. +\sum_{\mathbf{Z}}q(\mathbf{Z}|m)\{\text{ln}\ p(\mathbf{Z,X}|m)-\text{ln}\ q(\mathbf{Z}|m)\} \right)\\ =\sum_mq(m)\{\text{ln}(p(m)\text{exp}\{\mathcal{L}_m\})-\text{ln}\ q(m)\}\ \ \ \ (264) L=m∑​Z∑​q(Z∣m)q(m)ln{ q(Z∣m)q(m)p(Z,X,m)​}=m∑​Z∑​q(Z∣m)q(m){ ln p(Z,X∣m)+ln p(m)−ln q(Z∣m)−ln q(m)}=m∑​q(m)(ln p(m)−ln q(m)+Z∑​q(Z∣m){ ln p(Z,X∣m)−ln q(Z∣m)})=m∑​q(m){ ln(p(m)exp{ Lm​})−ln q(m)}    (264)
其中有:
$${\mathcal{L}}_m=\sum _{\mathbf{Z}}q(\mathbf{Z}|m)\text{ln}\left\{\frac{p(\mathbf{Z},\mathbf{X}|m)}{q(\mathbf{Z}|m)}\right\}$$
其中我们注意到，(264)式其实是在$q(m)$与（不保证归一化的） p ( m ) exp { L m } p(m)\text{exp}\{\mathcal{L}_m\} p(m)exp{ Lm​}之间的negative KL divergence。该值的最大化意味着相应KL divergence的最小化，因而为了使得我们所取得的那个$q(m)$能够让这个式子最大化，需要满足：
q ( m ) ∝ p ( m ) exp { L m }      ( 10.36 ) q(m)\propto p(m)\text{exp}\{\mathcal{L}_m\}\ \ \ \ (10.36) q(m)∝p(m)exp{ Lm​}    (10.36)
然而，如果我们对$\mathcal{L}$关于$q(\mathbf{Z}|m)$进行最大化, 我们将会看到的是，关于不同的$m$的solutions将会耦合在一起，这正如我们所预料的，因为它们都conditioned on $m$。因此，我们转而首先通过优化(10.35)的方式来对每一个$q(\mathbf{Z}|m)$单独进行优化，然后接下来使用(10.36)来确定$q(m)$的值。在normalization之后，关于$q(m)$的resulting values可以用在model selection或是model average中。

10.2 Illustration: Variational Mixture of Gaussians

我们现在回到我们关于Gaussian mixture model的讨论中，并且将我们在之前章节中得到的variational inference machinery应用进来。这将会为Variational methods的应用提供一个较好的解说，并且同样也会阐明Bayesian treatment是如何优雅地解决maximum likelihood approach中的一些困难点的。本书作者建议读者细致地研究这个例子，因为它提供了关于Variational methods的许多practical application。许多Bayesian models（对应于许多更为复杂的distributions）可以通过这里所叙述的分析的extension以及generalization而得到解决。

我们首先从Gaussian mixture model的likelihood function开始，如图9.6所示。对于每一个observation ${\mathbf{x}}_n$, 我们具有一个对应的latent variable ${\mathbf{z}}_n$，它包含着一个1-of-K binary vector, 其中有elements $z_{nk},k=1,...,K$。同之前一样，我们将observed data set标记为 X = { x 1 , . . . , x N } \mathbf{X}=\{\mathbf{x}_1,...,\mathbf{x}_N\} X={ x1​,...,xN​}, 类似的，我们将latent variables标记为$\mathbf{Z}={\mathbf{z}}_1,...,{\mathbf{z}}_N$。由(9.10), 我们可以在给定mixing coefficients $\mathbf{\pi }$的时候写出关于$\mathbf{Z}$的conditional distribution：
p ( Z ∣ π ) = ∏ n = 1 N ∏ k = 1