从极大似然说起
问题描述
假设我们需要调查我们学校学生的身高分布。我们先假设学校所有学生的身高服从正态分布 $N\left(\mu, \sigma^{2}\right)$ 。(注意：极大似然估计的前提一定是要假设数据总体的分布，如果不知道数据分布，是无法使用极大似然估计的)，这个分布的均值 $\mu$ 和方差 $\sigma^{2}$ 未知，如果我们估计出这两个参数，那我们就得到了结果。那么怎样估计这两个参数呢？
学校的学生这么多，我们不可能挨个统计吧？这时候我们需要用到概率统计的思想，也就是抽样，根据样本估算总体。我们可以先对学生进行抽样。假设我们随机抽到了 200 个人（也就是 200 个身高的样本数据，为了方便表示，下面，“人”的意思就是对应的身高）。然后统计抽样这 200 个人的身高。根据这 200 个人的身高估计均值 $\mu$ 和方差 $\sigma^{2}$ 。
那么问题来了怎样估算参数 $\theta(\theta=\mu, \sigma)$ 呢？
估计参数
我们先回答几个小问题：
问题一：抽到这 200 个人的概率是多少呢？
由于每个样本都是独立地从 $p(x | \theta)$ 中抽取的，换句话说这 200 个学生随便捉的，他们之间是没有关系的，即他们之间是相互独立的。假如抽到学生 A（的身高）的概率是 $p\left(x_{A} | \theta\right)$，抽到学生B的概率是 $p\left(x_{B} | \theta\right)$，那么同时抽到男生 A 和男生 B 的概率是 $p\left(x_{A} | \theta\right) \times p\left(x_{B} | \theta\right)$，同理，我同时抽到这 200 个学生的概率就是他们各自概率的乘积了，即为他们的联合概率用下式表示：
$$ L(\theta)=L\left(x_{1}, x_{2}, \cdots, x_{n} ; \theta\right)=\prod_{i=1}^{n} p\left(x_{i} | \theta\right), \quad \theta \in \Theta $$
这个概率反映了，在概率密度函数的参数是 $\theta$ 时，得到 X 这组样本的概率。等式右侧只有 $\theta$ 是未知数，所以它是 $\theta$ 的函数。
这个函数反映的是在不同的参数 $\theta$ 取值下，取得当前这个样本集的可能性，因此称为参数 $\theta$ 相对于样本集 X 的似然函数（likelihood function），记为 $ L(\theta) $ 。
为了便于分析，我们定义一个对数似然函数，将其变成连加的，称为对数似然函数：
$$ H(\theta)=\ln L(\theta)=\ln \prod_{i=1}^{n} p\left(x_{i} ; \theta\right)=\sum_{i=1}^{n} \ln p\left(x_{i} ; \theta\right) $$
问题二：学校那么多学生，为什么就恰好抽到了这 200 个人 ( 身高) 呢？
在学校那么学生中，我一抽就抽到这 200 个学生（身高），而不是其他人，那是不是表示在整个学校中，这 200 个人（的身高）出现的概率极大啊，也就是其对应的似然函数 $ L(\theta) $ 极大，即
$$ \hat{\theta}=\arg \max L(\theta) $$
$ \hat{\theta}= $ 这个叫做 $ \theta $ 的极大似然估计量，即为我们所求的值。
问题三：那么怎么极大似然函数？
求 $ L(\theta) $ 对所有参数的偏导数，然后让这些偏导数为 0，假设有 $n$ 个参数，就有 $n$ 个方程组成的方程组，那么方程组的解就是似然函数的极值点了，从而得到对应的 $ \theta $ 了。
极大似然估计
极大似然估计可以看作是一个反推。多数情况下我们是根据已知条件来推算结果，而极大似然估计是已经知道了结果，然后寻求使该结果出现的可能性极大的条件，以此作为估计值。
比如说，
- 假如一个学校的学生男女比例为 9:1 (条件)，那么你可以推出，你在这个学校里更大可能性遇到的是男生 (结果)；
- 假如你不知道男女比例，你走在路上，碰到100个人，发现男生就有90个 (结果)，这时候你可以推断这个学校的男女比例更有可能为 9:1 (条件)，这就是极大似然估计。
极大似然估计，只是一种概率论在统计学的应用，它是参数估计的方法之一。说的是已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，通过若干次试验，观察其结果，利用结果推出参数的大概值。
极大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率极大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。
求极大似然函数的一般步骤
（1）写出似然函数；
（2）对似然函数取对数，并整理；
（3）求导数，令导数为 0，得到似然方程；
（4）解似然方程，得到的参数。
极大似然函数的应用
应用一 ：回归问题中的极小化平方和（极小化代价函数）
假设线性回归模型具有如下形式: $h(x)=\sum_{j=1}^{d} \theta_{j} x_{ij}+\epsilon=\theta^{T} x_i+\epsilon$ ，其中 $x_i \in R^{1 \times d}, \theta \in R^{1 \times d}$ ， 误 差 $\epsilon \in R$ ， 当前已知 $X=\left(x_{1}, \cdots, x_{m}\right)^{T} \in R^{m \times d}$ ， 如何求 $ \theta $ 呢？
- 最小二乘估计：最合理的参数应该使得模型能最好地拟合样本数据，也就是估计值和观测值之差的平方和最小，其推导过程如下所示：
$$ J(\theta)=\sum_{i=1}^{n}\left(h_{\theta}\left(x_{i}\right)-y_{i}\right)^{2} $$
求解方法是通过梯度下降算法，通过训练数据不断迭代得到最终的值。 - 极大似然法：最合理的参数应该使得从模型中抽取该 $m$ 组样本观测值的概率极大，也就是似然函数极大。假设误差项 $\epsilon \in N\left(0, \sigma^{2}\right)$，则 $y_{i} \in N\left(\theta x_{i}, \sigma^{2}\right)$
$$p\left(y_{i} | x_{i} ; \theta\right)=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y_{i}-\theta^{T} x_{i}\right)^{2}}{2 \sigma^{2}}\right)$$
$$\begin{aligned} L(\theta) &= \prod_{i=1}^{m} p\left(y_{i} | x_{i} ; \theta\right) \\
&=\prod_{i=1}^{m} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y_{i}-\theta^{T} x_{i}\right)^{2}}{2 \sigma^{2}}\right) \\ H(\theta) &=\log (L(\theta)) \\
&=\log \prod_{i=1}^{m} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y_{i}-\theta^{T} x_{i}\right)^{2}}{2 \sigma^{2}}\right) \\
&=\sum_{i=1}^{m} \log \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y_{i}-\theta^{T} x_{i}\right)^{2}}{2 \sigma^{2}}\right) \\
&=-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{m}\left(y_{i}-\theta^{T} x_{i}\right)^{2}-m \ln \sigma \sqrt{2 \pi} \end{aligned}$$
令 $J(\theta)=\frac{1}{2} \sum_{i=1}^{m}\left(y_{i}-\theta^{T} x_{i}\right)^{2}$ 则 $ \arg \max _ {\theta} H(\theta) \Leftrightarrow \arg \min _ {\theta} J(\theta)$，即将极大似然函数等价于极小化平方和。
这时可以发现，此时的极大化似然函数和最初的最小二乘损失函数的估计结果是等价的。但是要注意这两者只是恰好有着相同的表达结果，原理和出发点完全不同。
应用二：分类问题中极小化交叉熵 （极小化代价函数）
在分类问题中，交叉熵的本质就是似然函数的极大化，逻辑回归的假设函数为：
$$h(x)=\hat{y}=\frac{1}{1+e^{-\theta^{T} x+b}}$$
根据之前学过的内容我们知道 $\hat{y}=p(y=1 | x, \theta)$，
当 y = 1 时，$p_{1}=p(y=1 | x, \theta)=\hat{y}$
当 y = 0 时，$p_{0}=p(y=0 | x, \theta)=1-\hat{y}$
合并上面两式子，可以得到 $p(y | x, \theta)=\hat{y}^{y}(1-\hat{y})^{1-y}$
$$\begin{aligned} L(\theta) &=\prod_{i=1}^{m} p\left(y_{i} | x_{i} ; \theta\right) \\ &=\prod_{i=1}^{m} \hat{y}_ {i}^{y_{i}}\left(1-\hat{y}_ {i}\right)^{1-y_{i}} \\
H(\theta) &=\log (L(\theta)) \\
&=\log \prod_{i=1}^{m} \hat{y}_ {i}^{y_{i}}\left(1-\hat{y}_ {i}\right)^{1-y_{i}} \\
&=\sum_{i=1}^{m} \log \hat{y}_ {i}^{y_{i}}\left(1-\hat{y}_ {i}\right)^{1-y_{i}} \\
&=\sum_{i=1}^{m} y_{i} \log \hat{y}_ {i}+\left(1-y_{i}\right) \log \left(1-\hat{y}_ {i}\right) \end{aligned}$$
令 $J(\theta)=-H(\theta)=-\sum_{i=1}^{m} y_{i} \log \hat{y}_ {i}+\left(1-y_{i}\right) \log \left(1-\hat{y}_ {i}\right)$ 则 $\arg \max_{\theta} H(\theta) \Leftrightarrow \arg \min_{\theta} J(\theta)$， 即将极大似然函数等价于极小化交叉熵。
EM算法
问题描述
上面我们先假设学校所有学生的身高服从正态分布 $N\left(\mu, \sigma^{2}\right)$。实际情况并不是这样的，男生和女生分别服从两种不同的正态分布，即男生 $\in N\left(\mu_{1}, \sigma_{1}^{2}\right)$，女生 $\in N\left(\mu_{2}, \sigma_{2}^{2}\right)$，(注意：EM算法和极大似然估计的前提是一样的，都要假设数据总体的分布，如果不知道数据分布，是无法使用EM算法的)那么该怎样评估学生的身高分布呢？
简单啊，我们可以随便抽 100 个男生和 100 个女生，将男生和女生分开，对他们单独进行极大似然估计。分别求出男生和女生的分布。
假如某些男生和某些女生好上了，纠缠起来了。咱们也不想那么残忍，硬把他们拉扯开。这时候，你从这 200 个人（的身高）里面随便给我指一个人（的身高），我都无法确定这个人（的身高）是男生（的身高）还是女生（的身高）。用数学的语言就是，抽取得到的每个样本都不知道是从哪个分布来的。那怎么办呢？
EM 算法
这个时候，对于每一个样本或者你抽取到的人，就有两个问题需要估计了，一是这个人是男的还是女的，二是男生和女生对应的身高的正态分布的参数是多少。这两个问题是相互依赖的：
- 当我们知道了每个人是男生还是女生，我们可以很容易利用极大似然对男女各自的身高的分布进行估计。
- 反过来，当我们知道了男女身高的分布参数我们才能知道每一个人更有可能是男生还是女生。例如我们已知男生的身高分布为 $N\left(\mu_{1}=172, \sigma_{1}^{2}=5^{2}\right)$，女生的身高分布为 $N\left(\mu_{2}=162, \sigma_{2}^{2}=5^{2}\right)$，一个学生的身高为 180，我们可以推断出这个学生为男生的可能性更大。
但是现在我们既不知道每个学生是男生还是女生，也不知道男生和女生的身高分布。这就成了一个鸡生蛋，蛋生鸡的问题了。鸡说，没有我，谁把你生出来的啊。蛋不服，说，没有我，你从哪蹦出来啊。为了解决这个你依赖我，我依赖你的循环依赖问题，总得有一方要先打破僵局，我先随便假设一个值出来，看你怎么变，然后我再根据你的变化调整我的变化，然后如此迭代着不断互相推导，最终就会收敛到一个解。这就是 EM 算法的基本思想了。
EM 的具体做法为：
- 先设定男生和女生的身高分布参数(初始值)，例如男生的身高分布为 44， 女生的身高分布 $$，当然了，刚开始肯定没那么准；
- 然后计算出每个人更可能属于第一个还是第二个正态分布中的（例如，这个人的身高是180，那很明显，他极大可能属于男生的那个分布），这个是属于
Expectation一步； - 我们已经大概地按上面的方法将这 200 个人分为男生和女生两部分，我们就可以根据之前说的极大似然估计分别对男生和女生的身高分布参数进行估计。这个是
Maximization； - 然后，当我们更新这两个分布的时候，每一个学生属于女生还是男生的概率又变了，那么我们就再需要调整 E 步；
- …… 如此往复，直到参数基本不再发生变化或满足结束条件为止。
总结
上面的学生属于男生还是女生我们称之为隐含参数，女生和男生的身高分布参数称为模型参数。
EM 算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含参数（EM 算法的 E 步），接着基于观察数据和猜测的隐含参数一起来极大化对数似然，求解我们的模型参数（EM 算法的 M 步)。由于我们之前的隐含参数是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。我们基于当前得到的模型参数，继续猜测隐含参数（EM 算法的 E 步），然后继续极大化对数似然，求解我们的模型参数（EM 算法的 M 步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。
EM算法的详细推导
基础知识
凸函数
设是定义在实数域上的函数，如果对于任意的实数，都有：
$$f^{\prime \prime} \geq 0$$
那么是凸函数。若不是单个实数，而是由实数组成的向量，此时，如果函数的 Hesse 矩阵是半正定的，即
$$H^{\prime \prime} \geq 0$$
那么是凸函数。特别地，如果 $f^{\prime \prime}>0$或者 $H^{\prime \prime}>0$，那么称为严格凸函数。
Jensen不等式
如下图，如果 $f$ 是凸函数，$x$ 是随机变量，那么有 0.5 的概率是 a，有 0.5 的概率是 b，$x$ 的期望值就是 a 和 b 的中值了那么：
$$E(f(x))\ge f(E(x))$$
其中，$E[f(x)]=0.5 f(a)+0.5 f(b), f(E(x))=f(0.5 a+0.5 b)$，这里 $a$ 和 $b$ 的权值为 0.5, $f(a)$ 与 $a$ 的权值相等，$f(b)$ 与 $b$ 的权值相等。
特别地，如果 $f$ 是严格凸函数，当且仅当 $p(x=E(x))=1$，(即随机变量是常量) 时等号成立。
注：若函数 $f$ 是凹函数，Jensen 不等式符号相反。
期望
对于离散型随机变量 $X$ 的概率分布为 $p_{i}=p{X=x_{i}}$，数学期望 $E(X)$ 为：
$$E(X)=\sum_{i} x_{i} p_{i}$$
$p_i$ 是权值，满足两个条件 $1 \geq p_{i} \geq 0, \sum_{i} p_{i}=1$ 。
若连续型随机变量 $X$ 的概率密度函数为 $f(x)$，则数学期望 $E(X)$ 为：
$$E(X)=\int_{-\infty}^{+\infty} x f(x) d x$$
设 $Y=g(X)$， 若 $X$ 是离散型随机变量，则:
$$E(Y)=\sum_{i} g\left(x_{i}\right) p_{i}$$
若 $f(x)$ 是连续型随机变量，则：
$$E(X)=\int_{-\infty}^{+\infty} g(x) f(x) d x$$
EM 算法的推导
对于 $m$ 个相互独立的样本 $x=\left(x^{(1)}, x^{(2)}, \ldots x^{(m)}\right)$，对应的隐含数据 $z=\left(z^{(1)}, z^{(2)}, \ldots z^{(m)}\right)$，此时 $(x, z)$ 即为完全数据，样本的模型参数为 $\theta$, 则观察数据 $x^{(i)}$ 的概率为 $P\left(x^{(i)} | \theta\right)$，完全数据 $\left(x^{(i)}, z^{(i)}\right)$ 的似然函数为 $P\left(x^{(i)}, z^{(i)} | \theta\right)$ 。
假如没有隐含变量 $z$，我们仅需要找到合适的 $\theta$ 极大化对数似然函数即可：
$$\theta=\arg \max _ {\theta} L(\theta)= \arg \max _ {\theta} \sum_{i=1}^{m} \log P\left(x^{(i)} | \theta\right)$$
增加隐含变量 $z$ 之后，我们的目标变成了找到合适的 $\theta$ 和 $z$ 让对数似然函数极大：
$$\theta, z=\arg \max _ {\theta, z} L(\theta, z)=\arg \max _ {\theta, z} \sum_{i=1}^{m} \log \sum_{z^{(i)}} P\left(x^{(i)}, z^{(i)} | \theta\right)$$
不就是多了一个隐变量 $z$ 吗？那我们自然而然会想到分别对未知的 $\theta$ 和 $z$ 分别求偏导，这样做可行吗？
理论上是可行的，然而如果对分别对未知的 $\theta$ 和 $z$ 分别求偏导，由于 $\log P\left(x^{(i)} | \theta\right)$ 是 $P\left(x^{(i)}, z^{(i)} | \theta\right)$ 边缘概率，转化为 $\log P\left(x^{(i)} | \theta\right)$ 求导后形式会非常复杂（可以想象下 $\log \left(f_{1}(x)+f_{2}(x)+\ldots\right.)$) 复合函数的求导，所以很难求解得到 $z$ 和 $\theta$ 。那么我们想一下可不可以将加号从 $log$ 中提取出来呢？我们可以对这个式子进行缩放如下：
$$ \begin{aligned} \sum_{i=1}^{m} \log \sum_{z^{(i)}} P\left(x^{(i)}, z^{(i)} | \theta\right) &=\sum_{i=1}^{m} \log \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{P\left(x^{(i)}, z^{(i)} | \theta\right)}{Q_{i}\left(z^{(i)}\right)} \\
& \geq \sum_{i=1}^{m} \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)} \end{aligned} $$
上面第 (1) 式引入了一个未知的新的分布 $Q_{i}\left(z^{(i)}\right)$，满足：
$$\sum_{z} Q_{i}(z)=1,0 \leq Q_{i}(z) \leq 1$$
上面第 (2) 式用到了 Jensen 不等式 (对数函数是凹函数)：
$$\log (E(y)) \geq E(\log (y))$$
其中：
$$ \begin{aligned} E(y) &= \sum_{i} \lambda_{i} y_{i}, \lambda_{i} \geq 0, \sum_{i} \lambda_{i}=1 \\
y_{i} &= \frac{P\left(x^{(i)} , z^{(i)} | \theta\right)}{Q_{i}\left(z^{(i)}\right)} \\
\lambda_{i} &= Q_{i}\left(z^{(i)}\right) \end{aligned} $$
也就是说 $\frac{P\left(x^{(i)}, z^{(i)} | \theta\right)}{Q_{i}\left(z^{(i)}\right)}$ 为第 $i$ 个样本，$Q_{i}\left(z^{(i)}\right)$ 为第 $i$ 个样本对应的权重，那么：
$$E\left(\frac{P\left(x^{(i)} , z^{(i)} | \theta\right)}{Q_{i}\left(z^{(i)}\right)}\right)=\sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{P\left(x^{(i)} , z^{(i)} | \theta\right)}{Q_{i}\left(z^{(i)}\right)}$$
上式实际上是我们构建了 $L(\theta, z)$ 的下界 (E步)，下一步要做的就是寻找一个合适的 $Q_{i}(z)$ 最优化这个下界(M步)。
假设 $\theta$ 已经给定，那么 $\log L(\theta)$ 的值就取决于 $Q_{i}\left(z^{(i)}\right)$ 和 $p\left(x^{(i)}, z^{(i)}\right)$了。我们可以通过调整这两个概率使下界逼近 $\log L(\theta)$ 的真实值，当不等式变成等式时，说明我们调整后的下界能够等价于 $\log L(\theta)$ 了。由 Jensen 不等式可知，等式成立的条件是随机变量是常数，则有：
$$\frac{P\left(x^{(i)}, z^{(i)} | \theta\right)}{Q_{i}\left(z^{(i)}\right)}=c$$
其中 $c$ 为常数，对于任意 $i$，我们得到：
$$P\left(x^{(i)}, z^{(i)} | \theta\right)=c Q_{i}\left(z^{(i)}\right)$$
方程两边同时累加和：
$$\sum_{z} P\left(x^{(i)}, z^{(i)} | \theta\right)=c \sum_{z} Q_{i}\left(z^{(i)}\right)$$
由于 $\sum_{z} Q_{i}\left(z^{(i)}\right)=1$ 。从上面两式，我们可以得到：
$$\sum_{z} P\left(x^{(i)}, z^{(i)} | \theta\right)=c$$
$$Q_{i}\left(z^{(i)}\right)=\frac{P\left(x^{(i)}, z^{(i)} | \theta\right)}{c}=\frac{P\left(x^{(i)}, z^{(i)} | \theta\right)}{\sum_{z} P\left(x^{(i)}, z^{(i)} | \theta\right)}=\frac{P\left(x^{(i)}, z^{(i)} | \theta\right)}{P\left(x^{(i)} | \theta\right)}=P\left(z^{(i)} | x^{(i)}, \theta\right)$$
其中：
边缘概率公式：$P\left(x^{(i)} | \theta\right)=\sum_{z} P\left(x^{(i)}, z^{(i)} | \theta\right)$
条件概率公式：$\frac{P\left(x^{(i)} , z^{(i)} | \theta\right)}{P\left(x^{(i)} | \theta\right)}=P\left(z^{(i)} | x^{(i)}, \theta\right)$
从上式可以发现 $Q(z)$ 是已知样本和模型参数下的隐变量分布。
如果 $Q_{i}\left(z^{(i)}\right)=P\left(z^{(i)} | x^{(i)}, \theta\right)$， 则第 (2) 式是我们的包含隐藏数据的对数似然的一个下界。如果我们能极大化这个下界，则也在尝试极大化我们的对数似然。即我们需要极大化下式：
$$\arg \max _ {\theta} \sum_{i=1}^{m} \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)}$$
$$=\arg \max _ {\theta} \sum_{i=1}^{m} \sum_{z^{(i)}}\left(Q_{i}\left(z^{(i)}\right) \log P\left(x^{(i)}, z^{(i)} | \theta\right)-Q_{i}\left(z^{(i)}\right) \log Q_{i}\left(z^{(i)}\right)\right)$$
去掉上式中为常数的部分，则我们需要极大化的对数似然下界为：
$$\arg \max _ {\theta} \sum_{i=1}^{m} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log P\left(x^{(i)}, z^{(i)} | \theta\right)$$
至此，我们推出了在固定参数 $\theta$ 后分布 $Q_{i}\left(z^{(i)}\right)$ 的选择问题，从而建立了 $\log L(\theta)$ 的下界。 这是 E 步，接下来的 M 步骤就是固定 $Q_{i}\left(z^{(i)}\right)$ 后，调整 $\theta$，去极大化 $\log L(\theta)$ 的下界。
EM算法流程
现在我们总结下 EM 算法的流程。
输入：观察数据 $x=\left(x^{(1)}, x^{(2)}, \ldots x^{(m)}\right)$，联合分布 $p(x, z | \theta)$，条件分布 $p(z | x, \theta)$，极大迭代次数 $J$。
- 随机初始化模型参数 $\theta$ 的初值 $\theta^0$
- for $i$ from $1$ to $J$
a) E步：计算联合分布的条件概率期望：
$Q_{i}\left(z^{(i)}\right) :=P\left(z^{(i)} | x^{(i)}, \theta\right)$
b) M步：极大化 $L(\theta)$ ,得到 $\theta$:
$\theta :=\arg \max _ {\theta} \sum_{i=1}^{m} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log P\left(x^{(i)}, z^{(i)} | \theta\right)$
c) 重复E、M步骤直到 $\theta$ 收敛 - 输出：模型参数 $\theta$
EM算法另一种理解
- 梯度下降 (或上升)法 ：沿函数值增幅最大 (即梯度方向) 的反方向走，就能使函数值减小 (或增大) 幅度越大。主要应用在逻辑回归算法，深度学习等算法中。如下图所示：
- 坐标上升 (或下降法) 法 (EM算法) ：每次只选择一个维度，将原始的问题变成一个一元函数，然后针对这个一元函数求极值，如此反复轮换不同的维度进行迭代。
什么时候用梯度下降，什么时候用EM算法呢？
这里举个简单的例子，在 x-y 坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步：固定 θ，优化 Q；M步：固定 Q，优化 θ；交替将极值推向极大。
EM算法应用
如果我们从算法思想的角度来思考EM算法，我们可以发现我们的算法里已知的是观察数据，未知的是隐含数据和模型参数，在E步，我们所做的事情是固定模型参数的值，优化隐含数据的分布，而在M步，我们所做的事情是固定隐含数据分布，优化模型参数的值。EM的应用包括：
- 支持向量机的 SMO 算法
- 混合高斯模型（GMM）
- K-means
在 K-Means 聚类时，每个聚类簇的质心是隐含数据。我们会- 假设 K 个初始化质心；
- 然后计算得到每个样本最近的质心，并把样本聚类到最近的这个质心，即 EM 算法的 M 步；
- 根据聚类后的样本，求出每个簇的质心，即 EM 算法的 E 步；
- 重复这个 E 步和 M 步，直到质心不再变化为止，这样就完成了 K-Means 聚类。
- 隐马尔可夫模型
EM算法案例-两硬币模型
假设有两枚硬币 A、B，以相同的概率随机选择一个硬币，进行如下的掷硬币实验：共做 5 次实验，每次实验独立的掷十次，结果如图中 a 所示，例如某次实验产生了 H、T、T、T、H、H、T、H、T、H (H代表正面朝上)。a 是在知道每次选择的是 A 还是 B 的情况下进行，b 是在不知道选择的是 A 还是 B 的情况下进行，问如何估计两个硬币正面出现的概率？
CASE a
已知每个实验选择的是硬币 A 还是硬币 B，重点是如何计算输出的概率分布，这其实也是极大似然求导所得。
$$\begin{aligned} \arg\max\log P(Y | \theta) &=\log \left(\left(\theta_{B}^{5}\left(1-\theta_{B}\right)^{5}\right)\left(\theta_{A}^{9}\left(1-\theta_{A}\right)\right)\left(\theta_{A}^{8}\left(1-\theta_{A}\right)^{2}\right)\left(\theta_{B}^{4}\left(1-\theta_{B}\right)^{6}\right)\left(\theta_{A}^{7}\left(1-\theta_{A}\right)^{3}\right)\right) \\ &=\log \left[\left(\theta_{A}^{24}\left(1-\theta_{A}\right)^{6}\right)\left(\theta_{B}^{9}\left(1-\theta_{B}\right)^{11}\right)\right] \end{aligned}$$
上面这个式子求导之后发现，5 次实验中 A 正面向上的次数再除以总次数即为 $\hat{\theta}_ A $，5 次实验中 B 正面向上的次数再除以总次数即为 $\hat{\theta}_ B $，即:
$$\hat{\theta}_ {A}=\frac{24}{24+6}=0.80$$
$$\hat{\theta}_ {B}=\frac{9}{9+11}=0.45$$
CASE b
由于并不知道选择的是硬币 A 还是硬币 B，因此采用 EM 算法。
E步：初始化 $\hat {\theta}_ {A}^{(0)} = 0.60$ 和 $\hat {\theta}_ {B}^{(0)} = 0.50$ ，计算每个实验中选择的硬币是 A 和 B 的概率，例如第一个实验中选择 A 的概率为：
$$P(z=A | y_{1}, \theta)=\frac{P(z=A, y_{1} | \theta)}{P(z=A, y_{1} | \theta)+P(z=B, y_{1} | \theta)}=\frac{(0.6)^{5} \ast (0.4)^{5}}{(0.6)^{5} \ast (0.4)^{5}+(0.5)^{10}}=0.45$$
$$P(z=B | y_{1}, \theta)=1-P(z=A | y_{1}, \theta)=0.55$$
计算出每个实验为硬币 A 和硬币 B 的概率，然后进行加权求和。
M步：求出似然函数下界 $Q(\theta, \theta^i)$，$y_j$ 代表第 $j$ 次实验正面朝上的个数，$\mu_j$ 代表第 $j$ 次实验选择硬币 A 的概率，$1-\mu_j$ 代表第 $j$ 次实验选择硬币 B 的概率 。
$$\begin{aligned} Q\left(\theta, \theta^{i}\right) &=\sum_{j=1}^{5} \sum_{z} P\left(z | y_{j}, \theta^{i}\right) \log P\left(y_{j}, z | \theta\right) \\ &=\sum_{j=1}^{5} \mu_{j} \log \left(\theta_{A}^{y_{j}}\left(1-\theta_{A}\right)^{10-y_{j}}\right)+\left(1-\mu_{j}\right) \log \left(\theta_{B}^{y_{j}}\left(1-\theta_{B}\right)^{10-y_{j}}\right) \end{aligned}$$
针对 $L$ 函数求导来对参数求导，例如对 $\theta_A$求导：
$$\begin{aligned} \frac{\partial Q}{\partial \theta_{A}} &=\mu_{1}\left(\frac{y_{1}}{\theta_{A}}-\frac{10-y_{1}}{1-\theta_{A}}\right)+\cdots+\mu_{5}\left(\frac{y_{5}}{\theta_{A}}-\frac{10-y_{5}}{1-\theta_{A}}\right)=\mu_{1}\left(\frac{y_{1}-10 \theta_{A}}{\theta_{A}\left(1-\theta_{A}\right)}\right)+\cdots+\mu_{5}\left(\frac{y_{5}-10 \theta_{A}}{\theta_{A}\left(1-\theta_{A}\right)}\right) \\ &=\frac{\sum_{j=1}^{5} \mu_{j} y_{j}-\sum_{j=1}^{5} 10 \mu_{j} \theta_{A}}{\theta_{A}\left(1-\theta_{A}\right)} \end{aligned}$$
求导等于 0 之后就可得到图中的第一次迭代之后的参数值:
$$\hat{\theta}_ {A}^{(1)}=0.71$$
$$\hat{\theta}_ {B}^{(1)}=0.58$$
当然，基于Case a 我们也可以用一种更简单的方法求得：
$$\hat{\theta}_ {A}^{(1)}=\frac{21.3}{21.3+8.6}=0.71$$
$$\hat{\theta}_ {B}^{(1)}=\frac{11.7}{11.7+8.4}=0.58$$
第二轮迭代：基于第一轮 EM 计算好的 $\hat{\theta}_ A^{(1)}, \hat{\theta}_ B^{(1)}$, 进行第二轮 EM，计算每个实验中选择的硬币是 A 和 B 的概率（E步），然后在计算M步，如此继续迭代……迭代十步之后 $\hat{\theta}_ A^{(10)} = 0.8, \hat{\theta}_ B^{(10)} = 0.52$
