{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第十三讲：最大期望算法及其应用、因子分析模型\n",
    "\n",
    "上一讲我们介绍了无监督学习算法，现在使用的训练集都是没有类标记的数据集$\\left\\{x^{(1)},\\cdots,x^{(m)}\\right\\}$。我们想要对样本进行分类并得到每个类型的样本的的密度分布。于是，我们引入潜在变量$z$，用它来表示样本所属的分布，并将数据建模为$\\displaystyle P(x)=\\sum_zP(x\\mid z)P(z)$。在前半节课的介绍中，我们假设数据服从不同的高斯分布，并使用高斯混合模型来对样本建模，也就是令$z\\sim\\mathrm{Multinomial}(\\phi)$，令$x\\mid z=j\\sim\\mathcal N\\left(\\mu_j,\\varSigma_j\\right)$，使用EM算法的步骤迭代求出参数。\n",
    "\n",
    "后半节课介绍了泛化的EM算法，该算法的目标函数为$\\displaystyle\\operatorname*{max}_\\theta\\sum_i\\log p\\left(x^{(i)};\\theta\\right)$，因为潜在变量的存在，目标函数又可以表示为$\\displaystyle\\max_\\theta\\sum_i\\log\\sum_z^{(i)}\\log p\\left(x^{(i)},z^{(i)};\\theta\\right)$。接下来就是使用E步骤和M步骤迭代求参数了：\n",
    "\n",
    "* E步骤：$Q_i\\left(z^{(i)}\\right):=p\\left(z^{(i)}\\mid x^{(i)};\\theta\\right)$，构建目标函数$\\mathscr l$的下界；\n",
    "* M步骤：$\\displaystyle\\theta:=\\arg\\max_\\theta\\underbrace{\\sum_i\\sum_{z^{(i)}}Q_i\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta\\right)}{Q_i\\left(z^{(i)}\\right)}}_\\textrm{lower-bound function}$，最大化这个下界函数，求得新的参数。（如果把右边这个下界看做关于$\\theta$的函数，那么这个下界就是$\\mathscr l$的下界函数。）\n",
    "\n",
    "再介绍一下我们为什么要使用EM算法：因为直接最大化$\\mathscr l$，对参数求偏导、再将偏导置为零求参数的计算过程，会由于潜在变量的存在而难以求得解析解；而EM算法M步骤中的最大化的式子通常都是很容易由偏导解出解析解的形式。\n",
    "\n",
    "接着上一讲继续介绍EM算法。\n",
    "\n",
    "怎么能确定算法一定收敛呢？设$\\theta^{(t)}$和$\\theta^{(t+1)}$是EM算法两次成功迭代得到的参数，只要证明$\\mathscr l\\left(\\theta^{(t)}\\right)\\leq\\mathscr l\\left(\\theta^{(t+1)}\\right)$，就可以说明EM算法会随着迭代使目标函数值单调增加。而证明的关键在于$Q_i$的选择。当EM算法从参数$\\theta^{(t)}$开始迭代时，我们会选择$Q_i^{(t)}\\left(z^{(i)}\\right):=p\\left(z^{(i)}\\mid x^{(i)};\\theta^{(t)}\\right)$，在上一讲最后，我们看到这个选择保证了延森不等式成立，并且使[上一讲](chapter12.ipynb)中的$(3)$式取等号，而因为$\\displaystyle\\mathscr l\\left(\\theta^{(t)}\\right)=\\sum_i\\sum_{z^{(i)}}Q_i^{(t)}\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta^{(t)}\\right)}{Q_i^{(t)}\\left(z^{(i)}\\right)}$，则最大化等号右边的部分就可以得到参数$\\theta^{(t+1)}$。因此：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathscr l\\left(\\theta^{(t+1)}\\right)&\\geq\\sum_i\\sum_{z^{(i)}}Q_i^{(t)}\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta^{(t+1)}\\right)}{Q_i^{(t)}\\left(z^{(i)}\\right)}\\tag{4}\\\\\n",
    "&\\geq\\sum_i\\sum_{z^{(i)}}Q_i^{(t)}\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta^{(t)}\\right)}{Q_i^{(t)}\\left(z^{(i)}\\right)}\\tag{5}\\\\\n",
    "&=\\mathscr l\\left(\\theta^{(t)}\\right)\\tag{6}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "* $(4)$式来自$\\displaystyle\\mathscr l(\\theta)\\geq\\sum_i\\sum_{z^{(i)}}Q_i\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta\\right)}{Q_i\\left(z^{(i)}\\right)}$，不等式对任意$Q_i$和$\\theta$均成立，而对于$Q_i=Q_i^{(t)},\\theta=\\theta^{(t+1)}$时，则可以取到等号。\n",
    "* 而要得到$(5)$式，需要利用$\\theta^{(t+1)}$的定义式：$\\displaystyle\\mathrm{arg}\\operatorname*{max}_\\theta\\sum_i\\sum_{z^{(i)}}Q_i\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta\\right)}{Q_i\\left(z^{(i)}\\right)}$，正是这个式子使以$\\theta^{(t+1)}$为参数的目标函数总是等于或大于以$\\theta^{(t)}$为参数的目标函数。\n",
    "* $(6)$式已经在前面出现过，为了在参数为$\\theta^{(t)}$时，能够使延森不等式取到等号而选择的$Q_i^{(t)}$。\n",
    "\n",
    "因此，EM算法能够使似然函数单调收敛。在EM算法的描述中，我们提到“运行直到收敛”，依照上面的推导，一个合理的收敛检测可以是：检查两次成功迭代后$\\mathscr l(\\theta)$的增量，如果增量小于某个阈值，我们就可以说“因为EM算法增长已然很慢了，所以算法已经收敛”。\n",
    "\n",
    "**注意：**如果我们定义$\\displaystyle J(Q,\\theta)=\\sum_i\\sum_{z^{(i)}}Q_i\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta\\right)}{Q_i\\left(z^{(i)}\\right)}$，那么从前面的推导中可以得出$\\mathscr l(\\theta)\\leq J(Q,\\theta)$，则EM算法可以看做是函数J上的坐标上升过程，而E步骤就求$J$关于$Q$的函数的最大值，而M步骤就是求$J$关于$\\theta$的函数的最大值。\n",
    "\n",
    "## 3. 再看高斯混合模型\n",
    "\n",
    "有了一般化定义的EM算法，我们重新审视上一讲中的高斯混合模型，拟合模型中的参数$\\phi,\\mu,\\varSigma$。为了简洁一些，我们只看M步骤中关于$\\phi,\\mu_j$的参数更新过程，将$\\varSigma_j$的更新留给读者考虑。\n",
    "\n",
    "E步骤很简单，根据上面的推导，有：\n",
    "\n",
    "$$\n",
    "w_j^{(i)}=Q_i\\left(z^{(i)}=j\\right)=P\\left(z^{(i)}=j\\mid x^{(i)};\\phi,\\mu,\\varSigma\\right)=\\frac{p\\left(x^{(i)}\\mid z^{(i)}=j;\\mu,\\varSigma\\right)p\\left(z^{(i)}=j;\\phi\\right)}{\\sum_{i=1}^kp\\left(x^{(i)}\\mid z^{(i)}=l;\\mu,\\varSigma\\right)p\\left(z^{(i)}=l;\\phi\\right)}\n",
    "$$\n",
    "\n",
    "这里的$Q_i\\left(z^{(i)}=j\\right)$表示在$Q_i$分布中$z^{(i)}$取$j$的概率。\n",
    "\n",
    "接下来是M步骤，我们需要最大化关于参数$\\phi,\\mu,\\varSigma$的下界函数：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sum_{i=1}^m\\sum_{z^{(i)}}Q_i\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\phi,\\mu,\\varSigma\\right)}{Q_i\\left(z^{(i)}\\right)}&=\\sum_{i=1}^m\\sum_{j=1}^kQ_i\\left(z^{(i)}=j\\right)\\log\\frac{p\\left(x^{(i)}\\mid z^{(i)}=j;\\mu,\\varSigma\\right)p\\left(z^{(i)}=j;\\phi\\right)}{Q_i\\left(z^{(i)}=j\\right)}\\\\\n",
    "&=\\sum_{i=1}^m\\sum_{j=1}^kw_j^{(i)}\\log\\frac{\\frac{1}{\\sqrt{(2\\pi)^n\\left\\lvert\\varSigma_j\\right\\rvert}}\\exp\\left(-\\frac{1}{2}\\left(x^{(i)}-\\mu_j\\right)^T\\varSigma_j^{-1}\\left(x^{(i)}-\\mu_j\\right)\\right)\\cdot\\phi_j}{w_j^{(i)}}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "首先以$\\mu_l$为变量做最大化，如果我们对$\\mu_l$求偏导，则有：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\nabla_{\\mu_l}\\sum_{i=1}^m\\sum_{j=1}^kw_j^{(i)}\\log\\frac{\\frac{1}{\\sqrt{(2\\pi)^n\\left\\lvert\\varSigma_j\\right\\rvert}}\\exp\\left(-\\frac{1}{2}\\left(x^{(i)}-\\mu_j\\right)^T\\varSigma_j^{-1}\\left(x^{(i)}-\\mu_j\\right)\\right)\\cdot\\phi_j}{w_j^{(i)}}&=-\\nabla_{\\mu_l}\\sum_{i=1}^m\\sum_{j=1}^kw_j^{(i)}\\frac{1}{2}\\left(x^{(i)}-\\mu_j\\right)^T\\varSigma_j^{-1}\\left(x^{(i)}-\\mu_j\\right)\\\\\n",
    "&=\\frac{1}{2}\\sum_{i=1}^mw_l^{(i)}\\nabla_{\\mu_l}2\\mu_l^T\\varSigma_l^{-1}x^{(i)}-\\mu_l^T\\varSigma_l^{-1}\\mu_l\\\\\n",
    "&=\\sum_{i=1}^mw_l^{(i)}\\left(\\varSigma_l^{-1}x^{(i)}-\\varSigma_l^{-1}\\mu_l\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "将式子置为零再解出$\\mu_l$，就得到了参数$\\mu$的更新规则：\n",
    "\n",
    "$$\n",
    "\\mu_l:=\\frac{\\sum_{i=1}^mw_l^{(i)}x^{(i)}}{\\sum_{i=1}^mw^{(i)}}\n",
    "$$\n",
    "\n",
    "这就是上一讲中参数$\\mu$更新规则的由来。\n",
    "\n",
    "再来推导一个M步骤中的参数$\\phi$的更新规则。要对$\\phi$求偏导，则将原式中含有$\\phi$的项提出来，现在需要最大化：$\\displaystyle\\sum_{i=1}^m\\sum_{j=1}^kw_j^{(i)}\\log\\phi_j$。但是要注意，因为$\\phi_j=p\\left(z^{(i)}=j;\\phi\\right)$是一个概率表达式，所以参数还有一个限制条件——$\\phi_j$之和为$1$。为了处理限制条件$\\displaystyle\\sum_{j=1}^k\\phi_j=1$，我们构造拉格朗日算子：\n",
    "\n",
    "$$\n",
    "\\mathcal L(\\phi)=\\sum_{i=1}^m\\sum_{j=1}^kw_j^{(i)}\\log\\phi_j+\\beta\\left(\\sum_{j=1}^k\\phi_j-1\\right)\n",
    "$$\n",
    "\n",
    "式中的$\\beta$是拉格朗日乘数（我们不需要担心$\\phi_j\\geq0$的限制条件，在后面会看到，从这里推导出的结论将恰好能够保证该条件成立）。求偏导得：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial\\phi_j}\\mathcal L(\\phi)=\\sum_{i=1}^m\\frac{w_j^{(i)}}{\\phi_j}+\\beta\n",
    "$$\n",
    "\n",
    "将偏导置为零得到$\\displaystyle\\phi_j=\\frac{\\sum_{i=1}^mw_j^{(i)}}{-\\beta}$，也就是说$\\phi_j\\propto\\sum_{i=1}^mw_j^{(i)}$。由于$\\displaystyle\\sum_j\\phi_j=1=\\sum_{j=1}^k\\frac{\\sum_{i=1}^mw_j^{(i)}}{-\\beta}=\\frac{\\sum_{i=1}^m\\sum_{j=1}^kw_j^{(i)}}{-\\beta}=\\frac{m}{-\\beta}\\implies-\\beta=m$（因为$w_j^{(i)}=Q_i\\left(z^{(i)}=j\\right)$是一个概率密度，所以有$\\sum_jw_j^{(i)}=1$）。于是我们就得到了参数$\\phi$的更新规则：\n",
    "\n",
    "$$\n",
    "\\phi_j:=\\frac{1}{m}\\sum_{i=1}^mw_j^{(i)}\n",
    "$$\n",
    "\n",
    "## 4. 再看混合朴素贝叶斯模型\n",
    "\n",
    "简要介绍一下前面的（[第五讲](chapter5.ipynb)）朴素贝叶斯模型在EM算法下的应用。\n",
    "\n",
    "* 给定的训练集为$\\left\\{x^{(1)},\\cdots,x^{(m)}\\right\\},x^{(i)}\\in\\{0,1\\}^n$，即每个样本都是一个$n$维位向量，比如$x_j^{(i)}=1\\{\\textrm{单词}j\\textrm{是否出现在文档中}\\}$；\n",
    "* 模型中有潜在变量$z^{(i)}\\in\\{0,1\\}$（也就是说在本例中我们希望找到两个簇，而如果希望有更多分类的话可以继续扩展），$z^{(i)}\\sim\\mathrm{Bernoulli}(\\phi)$；\n",
    "* 同前面讲过的朴素贝叶斯模型一样，有$\\displaystyle p\\left(x^{(i)}\\mid z^{(i)}\\right)=\\prod_{j=1}^np\\left(x_j^{(i)}\\mid z^{(i)}\\right)$，特别的有$p\\left(x_j^{(i)}=1\\mid z^{(i)}=0\\right)=\\phi_{j\\mid z=0}$；\n",
    "\n",
    "如果把上面条件中的$z$换成$y$，就和前面介绍的朴素贝叶斯模型一样了。有了这些条件，应用EM算法：\n",
    "\n",
    "* E步骤：$w^{(i)}:=p\\left(z^{(i)}=1\\mid x^{(i)};\\phi_{j\\mid z},\\phi\\right)$；\n",
    "* M步骤：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\phi_{j\\mid z=1}&=\\frac{\\sum_{i=1}^mw^{(i)}1\\left\\{x_j^{(i)}=1\\right\\}}{\\sum_{i=1}^mw^{(i)}}\\\\\n",
    "\\phi_{j\\mid z=0}&=\\frac{\\sum_{i=1}^m\\left(1-w^{(i)}\\right)1\\left\\{x_j^{(i)}=1\\right\\}}{\\sum_{i=1}^mw^{(i)}}\\\\\n",
    "\\phi_z&=\\frac{\\sum_{i=1}^mw^{(i)}}{m}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "# 第十部分：因子分析（Factor analysis）\n",
    "\n",
    "当我们使用高斯混合模型对数据集$x^{(i)}\\in\\mathbb R^n$建模时，可以使用EM算法拟合模型求出参数。在这种情形下，我们通常都有足够的样本，可以用多个高斯分布将它们分开。这也就是训练样本数量$m$远大于数据维数$n$的情形。\n",
    "\n",
    "现在，来考虑另一种情形——当$n\\gg m$时。在这种情形下，因为样本太少，即使只用一个高斯分布也难以对数据进行建模，更别说用多个高斯分布将数据分开了。特别是因为样本短缺，仅$m$个样本在$\\mathbb R^n$空间中只能张成一个维度很低的子空间，如果此时使用高斯分布建模，并使用最大似然估计对期望和协方差做出估计时：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mu&=\\frac{1}{m}\\sum_{i=1}^mx^{(i)}\\\\\n",
    "\\varSigma&=\\frac{1}{m}\\sum_{i=1}^m\\left(x^{(i)}-\\mu\\right)\\left(x^{(i)}-\\mu\\right)^T\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "就会发现，$\\varSigma$是一个奇异矩阵，也就是说其逆$\\varSigma^{-1}$不存在，则$\\displaystyle\\frac{1}{\\sqrt{\\left\\lvert\\varSigma\\right\\rvert}}=\\frac{1}{0}$，而这几个量在我们计算多元高斯分布的概率密度时都需要使用。换一种方式描述这种困境，我们需要用某个高斯分布拟合数据集，而这个高斯分布的参数的最大似然估计的所有可能性将会落在一个由样本数据张成的仿射空间中（对某些满足$\\sum_{i=1}^m\\alpha_i=1$的$\\alpha_i$，能够使数据集中的样本$x$满足$x=\\sum_{i=1}^m\\alpha_ix^{(i)}$），这也意味着其协方差矩阵是奇异的。\n",
    "\n",
    "更普遍的说，如果样本数量$m$没有比特征数量$n$大出合理的倍数，则得到的期望和协方差参数都会很差。尽管如此，我们依然希望使用一个合理的高斯模型拟合数据集，也许能够得到一些有趣的关于数据的协方差矩阵，那么该怎么做呢？\n",
    "\n",
    "在下一节，我们会介绍两种加在$\\varSigma$上的约束，以保证它可逆。其中一种允许我们用很小的数据集拟合$\\varSigma$，不过这两种方法都不能很好的解决我们的问题；然后讨论一下某些会在在后面用到的高斯分布的性质，特别是如何求高斯分布的边缘分布及条件分布；最后将介绍因子分析模型，并对它应用EM算法。\n",
    "\n",
    "## 1. $\\varSigma$的约束条件\n",
    "\n",
    "如果没有足够的样本来拟合一个完整的协方差矩阵，我们就可以考虑在其矩阵空间上加一些约束。比如说，我们可以选择拟合一个对角矩阵形式的$\\varSigma$。在这种情形下，容易验证对角矩阵$\\varSigma$的最大似然估计满足$\\displaystyle\\varSigma_{jj}=\\frac{1}{m}\\sum_{i=1}^m\\left(x_j^{(i)}-\\mu_j\\right)^2$，其中$\\varSigma_{jj}$就是样本数据第$j$个分量的方差。\n",
    "\n",
    "回忆前面（[第五讲](chapter05.ipynb)）介绍的高斯分布概率密度的等高线图，普通的$\\varSigma$对应的等高线图的主轴大都是倾斜的，而当$\\varSigma$是一个对角矩阵时，其等高线图的主轴是平行于坐标轴的。\n",
    "\n",
    "有时，我们会加上更严格的约束条件，除了要令$\\varSigma$是对角矩阵外，还要令其对角线元素是相等的。在这种情形下，我们令$\\varSigma=\\sigma^2I$，其中的$\\sigma^2$是需要估计的参数：$\\displaystyle\\sigma^2=\\frac{1}{mn}\\sum_{j=1}^n\\sum_{i=1}^m\\left(x_j^{(i)}-\\mu_j\\right)^2$。在这种情形下的高斯分布概率密度的等高线图将是正圆形（在二维情况下，而非普通$\\varSigma$下的椭圆；在高维情况下则是球面或超球面）。\n",
    "\n",
    "如果我们要用数据集拟合一个完整的、不加约束条件的协方差矩阵$\\varSigma$，则至少要满足$m\\geq n+1$，以保证$\\varSigma$的最大似然估计是非奇异矩阵。在上面的两种情形下，只要$m\\geq2$，我们得到的$\\varSigma$就是非奇异的。\n",
    "\n",
    "不过，将$\\varSigma$限制为对角矩阵也就意味着样本的任意两个分量$x_i,x_j$之间是无关且独立的。通常，在不加约束条件时，通过数据拟合得到的$\\varSigma$将捕获到一些关于数据的有趣的潜在联系，而我们使用上面任意一种约束后，就相当于强行抹掉了样本各分量之间的联系。在后面对因子分析模型的介绍中，我们将在对角矩阵约束之外添加更多参数，用以捕获数据间的联系，但同样不会拟合一个完整的协方差矩阵。\n",
    "\n",
    "## 2. 高斯分布的边缘分布及条件分布\n",
    "\n",
    "在介绍因子分析模型之前，我们先了解一下计算多元高斯分布下随机变量的边缘分布及条件分布的计算。\n",
    "\n",
    "* 设有向量形式的随机变量$x=\\begin{bmatrix}x_1\\\\x_2\\end{bmatrix},\\ x_1\\in\\mathbb R^r,x_2\\in\\mathbb R^s,x\\in\\mathbb R^{r+s}$；\n",
    "* 设$x\\sim\\mathcal N\\left(\\mu,\\varSigma\\right)$，其中$\\mu=\\begin{bmatrix}\\mu_1\\\\\\mu_2\\end{bmatrix},\\ \\varSigma=\\begin{bmatrix}\\varSigma_{11}&\\varSigma_{12}\\\\\\varSigma_{21}&\\varSigma_{22}\\end{bmatrix},\\ \\mu_1\\in\\mathbb R^r,\\mu_2\\in\\mathbb R^s,\\varSigma_{11}\\in\\mathbb R^{r\\times r},\\varSigma_{12}\\in\\mathbb R^{r\\times s},\\varSigma_{21}\\in\\mathbb R^{s\\times r},\\varSigma_{22}\\in\\mathbb R^{s\\times s}$，而且协方差矩阵是对称的，所以有$\\varSigma_{21}=\\varSigma_{12}^T$。\n",
    "\n",
    "在上面的假设下，$x_1,x_2$就是多元高斯分布。如何计算$x_1$的边缘分布？不难发现$\\mathrm E[x_1]=\\mu_1$，$\\mathrm{Cov}\\left(x_1\\right)=\\mathrm E\\left[\\left(x_1-\\mu_1\\right)\\left(x_1-\\mu_1\\right)^T\\right]=\\varSigma_{11}$。可以亲自验证$x_1$的协方差，使用$x_1,x_2$间的协方差定义式：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathrm{Cov}&=\\varSigma\\\\\n",
    "&=\\begin{bmatrix}\\varSigma_{11}&\\varSigma_{12}\\\\\\varSigma_{21}&\\varSigma_{22}\\end{bmatrix}\\\\\n",
    "&=\\mathrm E\\left[\\left(x-\\mu\\right)\\left(x-\\mu\\right)^T\\right]\\\\\n",
    "&=\\mathrm E\\begin{bmatrix}\\begin{pmatrix}x_1-\\mu_1\\\\x_2-\\mu_2\\end{pmatrix}\\begin{pmatrix}x_1-\\mu_1\\\\x_2-\\mu_2\\end{pmatrix}^T\\end{bmatrix}\\\\\n",
    "&=\\mathrm E\\begin{bmatrix}\\left(x_1-\\mu_1\\right)\\left(x_1-\\mu_1\\right)^T&\\left(x_1-\\mu_1\\right)\\left(x_2-\\mu_2\\right)^T\\\\\\left(x_2-\\mu_2\\right)\\left(x_1-\\mu_1\\right)^T&\\left(x_2-\\mu_2\\right)\\left(x_2-\\mu_2\\right)^T\\end{bmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "多元高斯分布的边缘分布也是高斯分布，则随机变量$x_1$的边缘分布为$x_1\\sim\\mathcal N\\left(\\mu_1,\\varSigma_{11}\\right)$。\n",
    "\n",
    "那么，在给定$x_2$的条件下$x_1$的条件分布又该怎么计算？根据多元高斯分布的定义有$x_1\\mid x_2\\sim\\mathcal N\\left(\\mu_{1\\mid2},\\varSigma_{1\\mid2}\\right)$，其中：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mu_{1\\mid2}&=\\mu_1+\\varSigma_{12}\\varSigma_{22}^{-2}\\left(x_2-\\mu_2\\right)\\tag{1}\\\\\n",
    "\\varSigma_{1\\mid2}&=\\varSigma_{11}-\\varSigma_{12}\\varSigma_{22}^{-1}\\varSigma_{21}\\tag{2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "在求解因子分析模型时，需要经常使用这些式子计算边缘分布和条件分布（参见[多元正态分布的条件分布](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions)）。\n",
    "\n",
    "## 3. 因子分析模型\n",
    "\n",
    "在因子分析模型中，我们假设随机变量$(x,z)$上的联合分布如下（其中$z\\in\\mathbb R^k$是潜在随机变量）：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "z&\\sim\\mathcal N(0,I)\\\\\n",
    "x\\mid z&\\sim\\mathcal N(\\mu+\\varLambda z,\\varPsi)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "模型的参数为：\n",
    "* $\\mu\\in\\mathbb R^n$；\n",
    "* 矩阵$\\varLambda\\in\\mathbb R^{n\\times k}$；\n",
    "* $\\varPsi\\in\\mathbb R^{n\\times n}$，并且是一个对角矩阵。\n",
    "\n",
    "而通常情况下，我们选择会选择比$n$小的数作为$k$的值。\n",
    "\n",
    "因此，我们设想数据$x^{(i)}$来自源于$k$维多元高斯分布$z^{(i)}$的一次抽样；然后，通过$\\mu+\\varLambda z^{(i)}$式将上一步的抽样映射到$\\mathbb R^n$的一个仿射空间中；最后，通过向$\\mu+\\varLambda z^{(i)}$添加协方差噪音$\\varPsi$，最终得到$x^{(i)}$。\n",
    "\n",
    "等价的，我们也可以如下定义因子分析模型模型：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "z&\\sim\\mathcal N(0,I)\\\\\n",
    "\\epsilon&\\sim\\mathcal N(0,\\varPsi)\\\\\n",
    "x&=\\mu+\\varLambda z+\\epsilon\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "其中$z,\\epsilon$是相互独立的，$\\epsilon$同样是噪音项。\n",
    "\n",
    "现在我们来看看我们的模型究竟定义了一个怎样分布。我们的随机变量$z$和$x$服从联合高斯分布$\\begin{bmatrix}z\\\\x\\end{bmatrix}\\sim\\mathcal N\\left(\\mu_{zx},\\varSigma\\right)$，现在需要知道$\\mu_{zx}$和$\\varSigma$。\n",
    "\n",
    "我们从$z\\sim\\mathrm N(0,I)$知道$\\mathrm E[z]=0$，于是可以算出：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathrm E[x]&=\\mathrm E[\\mu+\\varLambda z+\\epsilon]\\\\\n",
    "&=\\mu+\\varLambda\\mathrm E[z]+\\mathrm E[\\epsilon]\\\\\n",
    "&=\\mu\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "综合上面的两个期望能够得到$\\mu_{zx}=\\begin{bmatrix}\\vec0\\\\\\mu\\end{bmatrix}$。\n",
    "\n",
    "接下来要找到$\\varSigma$需要计算（协方差矩阵是对称矩阵）：\n",
    "* $\\varSigma_{zz}=\\mathrm E\\left[(z-\\mathrm E[z])(z-\\mathrm E[z])^T\\right]$（参考上一节协方差推导结果$\\varSigma$左上角元素）；\n",
    "* $\\varSigma_{zx}=\\mathrm E\\left[(z-\\mathrm E[z])(x-\\mathrm E[x])^T\\right]$（参考上一节协方差推导结果$\\varSigma$右上角元素）；\n",
    "* $\\varSigma_{xx}=\\mathrm E\\left[(x-\\mathrm E[x])(x-\\mathrm E[x])^T\\right]$（参考上一节协方差推导结果$\\varSigma$右下角元素）。\n",
    "\n",
    "不过由于$z\\sim\\mathcal N(0,I)$，我们可以化简$\\varSigma=\\mathrm{Cov}(z)=I$，也可以化简$\\varSigma_{zx}$：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathrm E\\left[(z-\\mathrm E[z])(x-\\mathrm E[x])^T\\right]&=\\mathrm E\\left[z(\\mu+\\varLambda z+\\epsilon-\\mu)^T\\right]\\\\\n",
    "&=\\mathrm E\\left[zz^T\\right]\\varLambda^T+\\mathrm E\\left[z\\epsilon^T\\right]\\\\\n",
    "&=\\varLambda^T\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "关于上面的最后一步推导，一是根据[协方差定义](https://en.wikipedia.org/wiki/Covariance#Definition)有$\\mathrm{Cov}(zz^T)=\\mathrm E\\left[zz^T\\right]-\\mathrm E[z]\\mathrm E\\left[z^T\\right]=\\mathrm E\\left[zz^T\\right]$（因为$z$的期望为零）；二是根据$z$与$\\epsilon$相互独立有$\\mathrm E\\left[z\\epsilon^T\\right]=\\mathrm E[z]\\mathrm E\\left[\\epsilon^T\\right]$（仍是从协方差定义能够得到，独立事件协方差为零，所以独立事件之积的期望等于各自期望之积）。同样也可以化简$\\varSigma_{xx}$：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathrm E\\left[(x-\\mathrm E[x])(z-\\mathrm E[x])^T\\right]&=\\mathrm E\\left[(\\mu+\\varLambda z+\\epsilon-\\mu)(\\mu+\\varLambda z+\\epsilon-\\mu)^T\\right]\\\\\n",
    "&=\\mathrm E\\left[\\varLambda zz^T\\varLambda^T+\\epsilon z^T\\varLambda^T+\\varLambda z\\epsilon^T+\\epsilon\\epsilon^T\\right]\\\\\n",
    "&=\\varLambda\\mathrm E\\left[zz^T\\right]\\varLambda^T+\\mathrm E\\left[\\epsilon\\epsilon^T\\right]\\\\\n",
    "&=\\varLambda\\varLambda^T+\\varPsi\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "综合上面所求的期望和协方差：\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}z\\\\x\\end{bmatrix}\\sim\\mathcal N\\left(\\begin{bmatrix}\\vec0\\\\\\mu\\end{bmatrix},\\begin{bmatrix}I&\\varLambda^T\\\\\\varLambda&\\varLambda\\varLambda^T+\\varPsi\\end{bmatrix}\\right)\\tag{3}\n",
    "$$\n",
    "\n",
    "从这个式子也可以看出$x$的边缘分布为$x\\sim\\mathcal N\\left(\\mu,\\varLambda\\varLambda^T+\\varPsi\\right)$。因此，对于给定的训练集$\\left\\{x^{(i)};i=1,\\cdots,m\\right\\}$，我们可以得到参数的对数似然函数：\n",
    "\n",
    "$$\n",
    "\\mathscr l(\\mu,\\varLambda,\\varPsi)=\\log\\prod_{i=1}^m\\frac{1}{\\sqrt{(2\\pi)^n\\left\\lvert\\varLambda\\varLambda^T+\\varPsi\\right\\rvert}}\\exp\\left(-\\frac{1}{2}\\left(x^{(i)}-\\mu\\right)^T\\left(\\varLambda\\varLambda^T+\\varPsi\\right)^{-1}\\left(x^{(i)}-\\mu\\right)\\right)\n",
    "$$\n",
    "\n",
    "下一步就是计算最大似然估计，求似然函数关于每个参数的最大值。但是最大化这个式子非常难，我们无法直接计算其最大值的解析解。于是，我们转而使用EM算法：\n",
    "\n",
    "* E步骤：$Q_i\\left(z^{(i)}\\right):=p\\left(z^{(i)}\\mid x^{(i)};\\theta\\right)$\n",
    "* M步骤：$\\displaystyle\\theta:=\\mathrm{arg}\\operatorname*{max}_\\theta\\sum_i\\int_{z^{(i)}}Q_i\\left(z^{(i)}\\right)\\log\\frac{p\\left(x^{(i)},z^{(i)};\\theta\\right)}{Q_i\\left(z^{(i)}\\right)}\\mathrm dz^{(i)}$（因为这里的$z^{(i)}$是一个连续变量，所以使用积分）。\n",
    "\n",
    "下一节将详细介绍EM算法在因子分析模型上的应用。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}