当有一组来自多个高斯分布的混合数据 $x^{(i)} \in \mathbb{R}^n$ 时，EM算法可以用来拟合高斯混合模型。这个问题中，通常假设拥有足够的数据，可以描绘出多个高斯分布的结构，严格来说，需要训练集大小 $m$ 显著地大于 数据维度 $n$。
现在，考虑 $n \gg m$ 的情况。这时，即便要拟合单个高斯模型，都是十分困难的，更不必说高斯混合模型。具体来看，由于 $m$ 个数据点稀疏地散布在较低的 $n$ 维空间中，假如要用高斯建模，并用最大似然估计法，估计均值和协方差，可得 $$ \begin{split} \mu &= \frac{1}{m}\sum_{i=1}^m x^{(i)} \\ \Sigma &= \frac{1}{m}\sum_{i=1}^m (x^{(i)} - \mu)(x^{(i)} - \mu)^T \end{split} $$ 可以发现，求出的结果，$\Sigma$ 是奇异矩阵，这意味着 $\Sigma^{-1}$ 不存在，并且 $1/|\Sigma|^{1/2}$ = 1/0，而这两项都需要用来计算多元高斯分布的概率密度。另一个困难点在于，最大似然估计法估计出的参数，会使得高斯分布的概率分布在数据的远交空间内。
更广泛来看，除非 $m$ 比 $n$ 大出较多的数量，否则最大似然估计法估计出的参数会非常糟糕。然而，我们依然希望对数据给出高斯模型的比较好的拟合，甚至捕捉到数据中关于协方差的一些结构信息。这要怎么做？
接下来，我们会首先增加两种可能的针对协方差的约束，使得可以用更小的数据集拟合 $\Sigma$。但这两种方法都无法给出令人满意的解。于是，我们会引入一些关于高斯分析的特性——边缘分布和条件分布——这将用于构建因子分析模型，以及对应的EM算法。
本节包含以下内容：
当没有足够的数据来拟合完整的协方差矩阵，可以考虑给协方差增加一些约束，从而减少对数据量的要求。比如，约束协方差为对角矩阵。在这个约束下，最大似然估计法估计出的协方差矩阵 $\Sigma$ 满足 $$ \Sigma_{jj} = \frac{1}{m}\sum_{i=1}^m (x_j^{(i)} - \mu_j)^2 $$ 因此，$\Sigma_{jj}$ 是数据第 $j$ 坐标方差的经验估计。
回想高斯密度的轮廓图是椭圆形的。一个是对角矩阵的协方差矩阵，对应着所有对称轴都和坐标轴平行的椭圆。
有时，还可以进一步约束协方差矩阵不仅为对角矩阵，并且对角线上的所有值都相等，也就是说 $\Sigma = \sigma^2I$，这里 $\sigma^2$ 可以根据最大似然估计法获得 $$ \sigma^2 = \frac{1}{mn}\sum_{j=1}^n \sum_{i=1}^m (x_j^{(i)} - \mu_j)^2 $$
这个模型对应着轮廓图为圆形的高斯密度函数。（在超过二维的情形下，则是三维的球体，或者高维的超球面）
如果要用最大似然估计法拟合无约束的完整协方差矩阵 $\Sigma$，$\Sigma$ 可逆（不是奇异矩阵）的必要条件为 $m \geq n+1$。而引入上面两种约束的任意一种，则只需要 $m \geq 2$ 即可。
但是，约束 $\Sigma$ 为对角矩阵，也就意味着不等坐标的 $x_i$ 和 $x_j$ 相互独立不相关。更多的时候，我们还是希望尽可能地捕捉到数据维度的相关性信息。而因子分析模型，就是为了无法拟合整个协方差矩阵的前提下，来描述数据维度的相关性信息。
在介绍银子分析之前，首先回顾一下多元高斯联合分布的条件分布和边缘分布。
首先设使我们有一个值为向量的随机变量 $$ \begin{split} x &= \begin{bmatrix} x_1 \\ x_2 \\ \end{bmatrix} \end{split} $$ 其中，$x_1 \in \mathbb{R}^r, x_2 \in \mathbb{R}^s, x \in \mathbb{R}^{r+s}$，假设 $x \sim \mathcal{N}(\mu,\Sigma)$，其中 $$ \begin{split} \mu &= \begin{bmatrix} \mu_1 \\ \mu_2 \\ \end{bmatrix} \end{split} , \begin{split} \Sigma &= \begin{bmatrix} \Sigma_{11} \, \Sigma_{12} \\ \Sigma_{21} \, \Sigma_{22} \\ \end{bmatrix} \end{split} $$ 这里，$\mu_1 \in \mathbb{R}^r, \mu_2 \in \mathbb{R}^s, \Sigma_{11} \in \mathbb{R}^{r \times r}, \Sigma_{12} \in \mathbb{R}^{r \times s}$，以此类推。注意协方差矩阵是对称的，$\Sigma_{12}=\Sigma_{21}^T$
根据我们的假设，$x_1$ 和 $x_2$ 是联合多元高斯分布。那么 $x_1$ 的边缘分布是什么？不难看出，$E[x_1]=\mu_1$，而 $Cov(x_1)=E[(x_1-\mu_1)(x_1-\mu_1)]=\Sigma_{11}$，其中后者是按照一下步骤推导得知的： $$ \begin{split} Cov(x) &= \Sigma \\ &= \begin{bmatrix} \Sigma_{11} \, \Sigma_{12} \\ \Sigma_{21} \, \Sigma_{22} \\ \end{bmatrix} \\ &= E[(x-\mu)(x-\mu)^T] \\ &= E[\begin{bmatrix}x_1-\mu_1 \\ x_2-\mu_2 \end{bmatrix}\begin{bmatrix}x_1-\mu_1 \\ x_2-\mu_2 \end{bmatrix}^T] \\ &= E[\begin{bmatrix}(x_1-\mu_1)(x_1-\mu_1)^T (x_1-\mu_1)(x_2-\mu_2)^T \\ (x_2-\mu_2)(x_1-\mu_1)^T (x_2-\mu_2)(x_2-\mu_2)^T \end{bmatrix}] \end{split} $$ 匹配矩阵的左上角即得到上述结果。
既然多元高斯分布的边缘分布本身也是高斯分布，从而有 $x_1$ 的边缘分布为 $x_1 \sim \mathcal{N}(\mu_1,\Sigma_{11})$。
此外，当给定 $x_2$ 时，$x_1$ 的条件分布是什么？参考多元高斯分布的定义，可以得知 $x_1|x_2 \sim \mathcal{N}(\mu_{1|2}, \Sigma_{1|2})$，其中 $$ \begin{split} \mu_{1|2} &= \mu_1 + \Sigma_{12}\Sigma_{22}^{-1}(x_2-\mu_2) \\ \Sigma_{1|2} &= \Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} \end{split} $$
在推导因子分析模型的过程中，多元高斯分布的条件分布和边缘分布将会非常有用。
在因子分析模型中，假定 $(x,z)$ 的联合分布如下，其中 $z \in \mathbb{R}^k$ 是隐含变量： $$ \begin{split} z &\sim \mathcal{N}(0, I) \\ x|z &\sim \mathcal{N}(\mu+\Lambda z, \Psi) \end{split} $$
这里模型的参数为向量 $\mu \in \mathbb{R}^n$，矩阵 $\Lambda \in \mathbb{R}^{n \times k}$，对角矩阵 $\Psi \in \mathbb{R}^{n \times n}$。$k$ 的值通常选取得比 $n$ 小。
因此，我们想象每个数据点 $x^{(i)}$ 首先从 $k$ 维的多元高斯分布中取得 $z^{(i)}$，之后，通过 $\mu+\Lambda z^{(i)}$ 投射到高维空间 $\mathbb{R}^n$，最后，通过协方差矩阵 $\Psi$ 对 $\mu+\Lambda z^{(i)}$ 增加噪声，得到 $x^{(i)}$。
于是，等价地，也可以像下面这样定义因子分析模型 $$ \begin{split} z &\sim \mathcal{N}(0,I) \\ \epsilon &\sim \mathcal{N}(0, \Psi) \\ x &= \mu + \Lambda z + \epsilon \end{split} $$ 其中，$\epsilon$ 和 $z$ 相互独立。
接下来，我们推导一下因子分析模型定义的分布。随机变量 $z$ 和 $x$ 有联合高斯分布 $$ \begin{bmatrix} x \\ z \\ \end{bmatrix} \sim \mathcal{N}(\mu_{zx}, \Sigma) $$ 我们需要计算出 $\mu_{zx}, \Sigma$
由于 $z \sim \mathcal{N}(0, I)$，所以 $E[z]=0$。同时，有 $$ \begin{split} E[x] &= E[\mu + \Lambda z + \epsilon] \\ &= \mu + \Lambda E[x] + E[\epsilon] \\ &= \mu \end{split} $$ 结合在一起，得到 $$ \mu_{zx} = \begin{bmatrix} \vec{0} \\ \mu \end{bmatrix} $$
而要计算协方差矩阵 $\Sigma$，需要分别计算 $\Sigma_{zz}=E[(z-E[z])(z-E[z])^T]$（$\Sigma$的左上角），$\Sigma_{zx}=E[(z-E[z])(x-E[x])^T]$（$\Sigma$的右上角），$\Sigma_{xx}=E[(x-E[x])(x-E[x])^T]$（$\Sigma$的右下角）
由于 $z \sim \mathcal{N}(0, I)$，所以 $\Sigma_{zz}=Cov(z)=I$。同时，有 $$ \begin{split} E[(z-E[z])(x-E[x])^T] &= E[z(\mu+\Lambda z+\epsilon - \mu)^T] \\ &= E[zz^T]\Lambda ^ T + E[z\epsilon^T] \\ &= \Lambda^T \end{split} $$ 这里，最后一步用到了 $E[zz^T]=Cov(z)=I$ ($z$ 的均值为零)，$E[z\epsilon^T]=E[z]E[\epsilon^T]=0$（$z$ 和 $x$ 相互独立）。类似地，可以推导 $\Sigma_{xx}$ $$ \begin{split} E[(x-E[x])(x-E[x])^T] &= E[(\mu+\Lambda z+\epsilon - \mu)(\mu+\Lambda z+\epsilon - \mu)^T] \\ &= E[\Lambda zz^T\Lambda^T + \epsilon z^T\Lambda^T + \Lambda z \epsilon^T + \epsilon\epsilon^T] \\ &= \Lambda E[zz^T]\Lambda^T + E[\epsilon\epsilon^T] \\ &= \Lambda \Lambda^T + \Psi \end{split} $$
讲上述推导的所有结果结合起来，有 $$ \begin{bmatrix}z\\x\\ \end{bmatrix} \sim \mathcal{N}(\begin{bmatrix}\vec{0}\\ \mu \\ \end{bmatrix}, \begin{bmatrix}I &\, \Lambda^T\\ \Lambda &\, \Lambda \Lambda^T + \Psi \end{bmatrix}) $$
根据上面这个联合分布，可以知道 $x$ 的边缘分布为 $x \sim \mathcal{N}(\mu, \Lambda \Lambda^T+\Psi)$，给定训练集 $\{x^{(i)};i=1,\cdots,m\}$，可以将所有参数的对数最大似然估计为 $$ \ell(\mu, \Lambda, \Psi) = log\prod_{i=1}^m \frac{1}{(2\pi)^{n/2}|\Lambda \Lambda^T + \Psi|}exp(-\frac{1}{2}(x^{(i)}-\mu)^T(\Lambda \Lambda^T + \Psi)^{-1}(x^{(i)}-\mu)) $$
同样地，由于隐含变量的存在，尽管参数里不包含，这个对数最大似然函数依然很难直接求解，同样需要EM算法来求解。
E步骤相对简单，需要计算 $Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\mu,\Lambda,\Psi)$，根据章节2中高斯分布的条件分布的公式，替换到章节3里的最终公式，可知 $z^{(i)}|x^{(i)};\mu,\Lambda,\Psi \sim \mathcal{N}(\mu_{z^{(i)}|x^{(i)}}, \Sigma_{z^{(i)}|x^{(i)}})$，其中 $$ \begin{split} \mu_{z^{(i)}|x^{(i)}} &= \Lambda^T(\Lambda \Lambda^T + \Psi)^{-1}(x^{(i)}-\mu) \\ \Sigma_{z^{(i)}|x^{(i)}}) &= I - \Lambda^T(\Lambda \Lambda^T + \Psi)^{-1}\Lambda \\ \end{split} $$ 根据这两个值的定义，有 $$ Q_i(z^{(i)}) = \frac{1}{(2\pi)^{k/2}|\Sigma_{z^{(i)}|x^{(i)}}|^{1/2}}exp(-\frac{1}{2}(z^{(i)}-\mu_{z^{(i)}|x^{(i)}})^T \Sigma^{-1}_{z^{(i)}|x^{(i)}}(z^{(i)}-\mu_{z^{(i)}|x^{(i)}})) $$
In [ ]: