引言
过去讲过的这些模型可以分为两类：概率模型和非概率模型。
- 对于概率模型，我们假设数据满足某些概率分布，例如贝叶斯分类器、Logistic回归、最小二乘和岭回归（如果从ML和MAP的角度看）以及贝叶斯线性回归
- 非概率模型不牵涉有关概率分布的假设，例如感知机、SVM、决策树和k均值聚类
对于以上所有模型，其实我们都定义了各自相对应的目标函数来优化，尽管优化的方法有很多种。其中概率模型的目标函数都是最大似然。最基本假设为
- 我们有一组模型参数\(\theta\)
- 有一组数据\(\{x_1, \ldots, x_n\}\)
- 假设一个条件概率分布\(p(x|\theta)\)
- 假设样本是独立同分布的，即\(x_i \mathop{\sim}^{iid} p(x|\theta)\)
最大似然是要找到能使得似然函数最大化的\(\theta\)，即 \[ \begin{align*} \theta_{\rm ML} &= {\rm arg}\max_\theta p(x_1, \ldots, x_n|\theta) \\ &= {\rm arg}\max_\theta \prod_{i=1}^n p(x_i|\theta) \\ &= {\rm arg}\max_\theta \sum_{i=1}^n \ln p(x_i|\theta) \end{align*} \]
对于类似贝叶斯分类器这样的模型，\(\theta\)可以从目标函数中求得解析解。对于更复杂的模型，可以把参数分割成\(\theta_1\)、\(\theta_2\)两组。尽管不能同时求出最优的\(\theta_1\)和\(\theta_2\)，但是可以使用kmeans中提到的坐标下降方法，分别迭代求解。
第三种办法是当 \[ \theta_{1,{\rm ML}} = {\rm arg}\max_{\theta_1} \sum_{i=1}^n \ln p(x_i|\theta_1) \] 很难直接优化时，引入一个变量\(\theta_2\)，即去解决 \[ \theta_{1,{\rm ML}} = {\rm arg}\max_{\theta_1} \sum_{i=1}^n \ln p(x_i, \theta_2|\theta_1) \] 这个式子可能更好求解。这是EM算法的基础
EM算法（Expectation-Maximization, 期望最大化算法）
算法每次迭代大致可以分成两步：Expectation（期望步，类似于坐标下降对第一个变量进行优化）和Maximization（似然步，类似于坐标下降对第二个变量进行优化）
可以从一个例子入手：假设\(x_i \in \mathbb{R}^d\)是一个数据点，但是一些维度上的分量已经丢失了，那么可以把这个数据点分成两个部分
- \(x_i^o\)：数据点中已经观察到的部分
- \(x_i^m\)：数据点中缺失/未知的部分
注意不同的数据点缺失的维度可能不同
假设\(x_i \mathop{\sim}^{iid} N(\mu, \Sigma)\)，想要求解 \[ \mu_{\rm ML}, \Sigma_{\rm ML} = {\rm arg}\max_{\mu, \Sigma} \sum_{i=1}^n \ln p(x_i^o|\mu, \Sigma) \] 是比较困难的。但是如果我们找到了\(x_i^m\)，就好办得多，即 \[ \mu_{\rm ML}, \Sigma_{\rm ML} = {\rm arg}\max_{\mu, \Sigma} \sum_{i=1}^n \ln \underbrace{p(x_i^o, x_i^m|\mu, \Sigma)}_{=p(x_i|\mu,\Sigma)} \] 容易优化
EM算法在优化\(\sum_{i=1}^n \ln p(x_i^o|\mu, \Sigma)\)的同时还填入了缺失值\(\{x_1^m,\ldots, x_n^m\}\)。通常情况下假设我们有两个参数集\(\theta_1, \theta_2\)，其中 \[ p(x|\theta_1) = \int p(x, \theta_2|\theta_1)d\theta_2 \] 我们需要定义一个通用的目标函数来得到我们想要的结果：首先，需要可以以\(\theta_1\)为自变量优化边缘概率密度\(p(x|\theta_1)\)，其次，可以通过\(p(x,\theta_2|\theta_1)\)来简化计算。这个目标函数可以写为以下形式(*) \[ \ln p(x|\theta_1) = \int q(\theta_2)\ln\frac{p(x,\theta_2|\theta_1)}{q(\theta_2)}d\theta_2 + \int q(\theta_2) \ln \frac{q(\theta_2)}{p(\theta_2|x,\theta_1)}d\theta_2 \] 其中\(q(\theta_2)\)是任何概率分布，而且假设我们知道\(p(\theta_2|x,\theta_1)\)
首先证明上式的正确性：由积分相加的法则和对数法则，有 \[ \int q(\theta_2)\ln\frac{p(x,\theta_2|\theta_1)}{q(\theta_2)}d\theta_2 + \int q(\theta_2) \ln \frac{q(\theta_2)}{p(\theta_2|x,\theta_1)}d\theta_2 = \int q(\theta_2) \ln \frac{p(x,\theta_2|\theta_1)q(\theta_2)}{p(\theta_2|x,\theta_1)q(\theta_2)} d\theta_2 \] 又 \[ p(a,b|c) = p(a|b,c)p(b|c) \Rightarrow p(b|c) = \frac{p(a,b|c)}{p(a|b,c)} \] 因此有 \[ \int q(\theta_2) \ln \frac{p(x,\theta_2|\theta_1)q(\theta_2)}{p(\theta_2|x,\theta_1)q(\theta_2)} d\theta_2 = \int q(\theta_2) \ln p(x|\theta_1) d\theta_2 = \ln p(x|\theta_1) \int q(\theta_2)d\theta_2 = \ln p(x|\theta_1) \]
*式的两个积分项各自有各自的意义。第一个积分项是对\(\theta_2\)的积分，因此最后\(\theta_2\)是可以被积掉的，也就是说最后就剩下一个关于\(\theta_1\)的函数，定义这个函数为\(\mathcal{L}(x,\theta_1)\)；第二个积分项称为Kullback_Leibler差异，也称为相对熵，可以看作是衡量了分布\(q(\theta_2)\)和分布\(p(\theta_2|x,\theta_1)\)之间的重合程度。只有当\(p=q\)时，第二项才为0
总和项的意义在于，在迭代优化\(\ln p(x|\theta_1)\)的过程中，可以保证得到的\(\theta_1\)的序列能满足当前迭代的ln值大于前一个迭代的ln值，即\(\ln p(x|\theta_1^{(t)}) \ge \ln p(x|\theta_1^{(t-1)})\)，而且\(\theta_1^{(t)}\)最后会收敛于\(\ln p(x|\theta_1)\)的一个局部最大值
EM算法： 给定第\(t\)个迭代时的\(\theta_1^{(t)}\)，根据下列算法算出下一个迭代时该参数的值（如果是第一个迭代，可以随机生成一个初始值）： - E-Step：令\(q_t(\theta_2) = p(\theta_2|x,\theta_1^{(t)})\)，计算 \[ \mathcal{L}_t(x, \theta_1) = \int q_t(\theta_2)\ln p(x,\theta_2|\theta_1)d\theta_2 - \int q_t(\theta_2)\ln q_t(\theta_2)d\theta_2 \] 其中后一项由于是常数，可以忽略掉 - M-Step：令\(\theta_1^{(t+1)} = {\rm arg}\max_{\theta_1}\mathcal{L}_t(x,\theta_1)\)
关于EM算法可以提升\(\ln p(x|\theta_1)\)的证明
\[ \begin{align*} \ln p(x|\theta_1^{(t)}) &= \mathcal{L}(x, \theta_1^{(t)}) + \underbrace{KL\left(q(\theta_2) || p(\theta_2|x_1, \theta_1^{(t)})\right)}_{=0\ \because\ q = p} \\ &= \mathcal{L}_t(x, \theta_1^{(t)})\ \ \ \ \leftarrow {\rm E-Step} \\ &\le \mathcal{L}_t(x, \theta_1^{(t+1)})\ \ \ \leftarrow {\rm M-Step} \\ &\le \mathcal{L}(x, \theta_1^{(t+1)}) + \underbrace{KL\left(q_t(\theta_2) || p(\theta_2|x_1, \theta_1^{(t+1)})\right)}_{>0\ \because\ q \not= p} \\ &= \mathcal{L}(x, \theta_1^{(t+1)}) + {KL\left(q(\theta_2) || p(\theta_2|x_1, \theta_1^{(t+1)})\right)}\\ &= \ln p(x|\theta_1^{(t+1)}) \end{align*} \]
使用EM算法处理缺失数据
套用EM算法，缺失数据问题可以表示为 \[ \sum_{i=1}^n \ln p(x_i^o | \mu, \Sigma) = \sum_{i=1}^n \int q(x_i^m) \ln \frac{p(x_i^o, x_i^m|\mu, \Sigma)}{q(x_i^m)}dx_i^m + \sum_{i=1}^n \int q(x_i^m)\ln \frac{q(x_i^m)}{p(x_i^m|x_i^o, \mu, \Sigma)} dx_i^m \] 其中\(x_i^o\)是观察到的部分，\(x_i^m\)是缺失的部分。不失一般性地，可以假设所有观察到的数据维度都在该条数据的前部，缺失数据都在后部，因此有 \[ x_i = \left[\!\begin{array}{c} x_i^o \\ x_i^m \end{array}\!\right] \sim N\left(\left[\!\begin{array}{c} \mu_i^o \\ \mu_i^m \end{array}\!\right] , \left[\!\begin{array}{cc} \Sigma_i^{oo} & \Sigma_i^{om} \\ \Sigma_i^{mo} & \Sigma_i^{mm} \end{array}\!\right]\right) \]
略过证明，有结论 \[ \begin{align*} p(x_i^m | x_i^o, \mu, \Sigma) &= N(\hat{\mu}_i, \hat{\Sigma}_i) \\ \hat{\mu}_i &= \mu_i^m + \Sigma_i^{mo}(\Sigma_i^{oo})^{-1}(x_i^o - \mu_i^o) \\ \hat{\Sigma}_i &= \Sigma_i^{mm} - \Sigma_i^{mo}(\Sigma_i^{oo})^{-1}\Sigma_i^{om} \end{align*} \]
在E阶段，有 \[ \begin{align*} \mathbb{E}_{q(x_i^m)}[\ln p(x_i^o, x_i^m | \mu, \Sigma)] &= \mathbb{E}_q [(x_i-\mu)^\mathsf{T}\Sigma^{-1}(x_i-\mu)] \\ &= \mathbb{E}_q [{\rm trace}\{\Sigma^{-1}(x_i - \mu)(x_i - \mu)^\mathsf{T}\}] \\ &= {\rm trace}\{\Sigma^{-1}\mathbb{E}_q[(x_i - \mu)(x_i - \mu)^\mathsf{T}]\} \end{align*} \] 同时， \[ q(x_i^m) = p(x_i^m|x_i^o, \mu, \Sigma) = N(\hat{\mu}_i, \hat{\Sigma}_i) \] 记\(\hat{x}_i\)为将所有缺失的\(x_i\)替换为\(\hat{\mu}_i\)的向量，\(\hat{V}_i\)为全零矩阵加\(\hat{\Sigma}_i\)中缺失数据对应的部分，则 在M阶段，有 \[ \begin{align*} \mu_{\rm up},\Sigma_{\rm up} &= {\rm arg}\max_{\mu, \Sigma} \sum_{i=1}^n \mathbb{E}_q [\ln p(x_i^o, x_i^m|\mu, \Sigma)] \\ \mu_{\rm up} &= \frac{1}{n}\sum_{i=1}^n \hat{x}_i \\ \Sigma_{\rm up} &= \frac{1}{n}\sum_{i=1}^n \{(\hat{x}_i - \mu_{\rm up})(\hat{x}_i - \mu_{\rm up})^\mathsf{T} + \hat{V}_i\} \end{align*} \]
其中\(\mu\)和\(\Sigma\)的初始值可以设为0或随机值，然后逐步迭代，直到每步只能做出微小的改变
