算法思想：含有隐变量的极大似然估计

我们经常会从样本观察数据中，找出样本的模型参数。 最常用的方法就是极大化模型分布的对数似然函数。

但是在一些情况下，我们得到的观察数据有未观察到的隐含数据，此时我们未知的有隐含数据和模型参数，因而无法直接用极大化对数似然函数得到模型分布的参数。怎么办呢？这就是EM算法可以派上用场的地方了。那么先复习一下极大似然估计。

极大似然估计（MLE）

直接举个例子：

某位同学与一位猎人一起外出打猎，一只野兔从前方窜过。只听一声枪响，野兔应声到下，如果要你推测，这一发命中的子弹是谁打的？你就会想，只发一枪便打中，由于猎人命中的概率一般大于这位同学命中的概率，看来这一枪是猎人射中的。

这个例子所作的推断就体现了极大似然法的基本思想。

MLE就是利用已知的样本结果，反推最有可能(最大概率)导致这样结果的参数值 的计算过程。直白来讲，就是给定了一定的数据，假定知道数据是从某种分布中 随机抽取出来的，但是不知道这个分布具体的参数值，即“模型已定，参数未 知”，MLE就可以用来估计模型的参数。MLE的目标是找出一组参数(模型中的 参数)，使得模型产出观察数据的概率最大。

求最大似然函数估计值的一般步骤：

（1）写出似然函数；

（2）对似然函数取对数，并整理；

（3）求导数，令导数为0，得到似然方程；

（4）解似然方程，得到的参数即为所求；

EM算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含数据（EM算法的E步），接着基于观察数据和猜测的隐含数据一起来极大化对数似然，求解我们的模型参数（EM算法的M步)。由于我们之前的隐藏数据是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。不过没关系，我们基于当前得到的模型参数，继续猜测隐含数据（EM算法的E步），然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。

从上面的描述可以看出，EM算法是迭代求解最大值的算法，同时算法在每一次迭代时分为两步，E步和M步。一轮轮迭代更新隐含数据和模型分布参数，直到收敛，即得到我们需要的模型参数。

K-Means算法

下面回顾一下K-Means算法的基本思想。

K-means算法，也称为k-均值聚类算法，是一种非常广泛使用的聚类算法之一。

假定输入样本为S=x1,x2,x3,...,xm,则算法步骤为：

• 假设输入样本为T=X1,X2…,Xm则算法步骤为(使用欧几里得距离公式)

• 选择初始化的k个类别中心a1,...ak

• 对于每个样本X,将其标记位距离类别中心aj最近的类别

• 更新每个类别的中心点aj为隶属该类别的所有样本的均值

更新后的中心点为：

表示第j个簇中的样本点的均值为新的中心点。

• 重复上面两步操作,直到达到某个中止条件

中止条件：

迭代次数、最小平方误差MSE、簇中心点变化率。

一个最直观了解EM算法思路的是K-Means算法，见之前写的K-Means聚类算法原理。在K-Means聚类时，每个聚类簇的质心是隐含数据。我们会假设KK个初始化质心，即EM算法的E步；然后计算得到每个样本最近的质心，并把样本聚类到最近的这个质心，即EM算法的M步。重复这个E步和M步，直到质心不再变化为止，这样就完成了K-Means聚类。

当然，K-Means算法是比较简单的，实际中的问题往往没有这么简单。上面对EM算法的描述还很粗糙，我们需要用数学的语言精准描述。

EM算法
先从一个简单的例子开始：
现在一个班里有50个男生，50个女生，且男生站左，女生站右。我们假定男生的身高服从正态分布μ1，σ1，女生的身高则服从另一个正态分布μ2，σ2 。这时候我们可以用极大似然法（MLE），分别通过这50个男生和50个女生的样本来估计这两个正态分布的参数。
但现在我们让情况复杂一点，就是这50个男生和50个女生混在一起了。我们拥有100个人的身高数据，却不知道这100个人每一个是男生还是女生。
这时候情况就有点尴尬，因为通常来说，我们只有知道了精确的男女身高的正态分布参数我们才能知道每一个人更有可能是男生还是女生。但从另一方面去考量，我们只有知道了每个人是男生还是女生才能尽可能准确地估计男女各自身高的正态分布的参数。
下面以一个知乎上的例子来计算说明EM算法，地址：
www.zhihu.com/question/27…
背景：公司有很多领导=[A总，刘总，C总]，同时有很多漂亮的女职员=[小甲，小章，小乙]。（请勿对号入座）你迫切的怀疑这些老总跟这些女职员有问题。为了科学的验证你的猜想，你进行了细致的观察。于是，
观察数据：
1）A总，小甲，小乙一起出门了；
2）刘总，小甲，小章一起出门了；
3）刘总，小章，小乙一起出门了；
4）C总，小乙一起出门了；
收集到了数据，你开始了神秘的EM计算：
初始化，你觉得三个老总一样帅，一样有钱，三个美女一样漂亮，每个人都可能跟每个人有关系。所以，每个老总跟每个女职员“有问题”的概率都是1/3;
这样，（E step）
1） A总跟小甲出去过了 1/2 * 1/3 = 1/6 次，跟小乙也出去了1/6次；（所谓的fractional count）
2）刘总跟小甲，小章也都出去了1/6次
3）刘总跟小乙，小章又出去了1/6次
4）C总跟小乙出去了1/3次
总计，A总跟小甲出去了1/6次，跟小乙也出去了1/6次 ; 刘总跟小甲，小乙出去了1/6次，跟小章出去了1/3次；C总跟小章出去了1/3次；
你开始跟新你的八卦了(M step),
A总跟小甲，小乙有问题的概率都是1/6 / (1/6 + 1/6) = 1/2；
刘总跟小甲，小乙有问题的概率是1/6 / (1/6+1/6+1/6+1/6) = 1/4; 跟小章有问题的概率是(1/6+1/6)/(1/6 * 4) = 1/2;
C总跟小乙有问题的概率是 1。
然后，你有开始根据最新的概率计算了；（E-step）
1）A总跟小甲出去了 1/2 * 1/2 = 1/4 次，跟小乙也出去 1/4 次；
2）刘总跟小甲出去了1/2 * 1/4 = 1/12 次， 跟小章出去了 1/2 * 1/2 = 1/4 次；
3）刘总跟小乙出去了1/2 * 1/4 = 1/12 次， 跟小章又出去了 1/2 * 1/2 = 1/4 次；
4）C总跟小乙出去了1次；
重新反思你的八卦（M-step）:
A总跟小甲，小乙有问题的概率都是1/4/ (1/4 + 1/4) = 1/2；
B总跟小甲，小乙是 1/12 / (1/12 + 1/4 + 1/4 + 1/12) = 1/8 ; 跟小章是 3/4 ;
C总跟小乙的概率是1。
你继续计算，反思，总之，最后，你得到了真相！
EM算法过程
通过上面的计算我们可以得知，EM算法实际上是一个不停迭代计算的过程， 根据我们事先估计的先验概率A，得出一个结果B，再根据结果B，再计算得到结果 A，然后反复直到这个过程收敛。
可以想象饭店的后方大厨，炒了两盘一样的菜，现在，菜炒好后从锅中倒入盘， 不可能一下子就分配均匀，所以先往两盘中倒入，然后发现B盘菜少了，就从A中 匀出一些，A少了，从B匀.......不停迭代。
例如，小时候，老妈给一大袋糖果给你，叫你和你姐姐等分，然后你懒得去点糖果的个数，所以你也就不知道每个人到底该分多少个。咱们一般怎么做呢？先把一袋糖果目测的分为两袋，然后把两袋糖果拿在左右手，看哪个重，如果右手重，那很明显右手这代糖果多了，然后你再在右手这袋糖果中抓一把放到左手这袋，然后再感受下哪个重，然后再从重的那袋抓一小把放进轻的那一袋，继续下去，直到你感觉两袋糖果差不多相等了为止。呵呵，然后为了体现公平，你还让你姐姐先选了。
初始化分布参数 重复下列两个操作直到收敛：
E步骤：估计隐藏变量的概率分布期望函数；
M步骤：根据期望函数重新估计分布参数。
参数求解过程
假设我们有一个样本集{x(1),…,x(m)}，包含m个独立的样本。但每个样本i对应的类别z(i)是未知的（相当于聚类），也即隐含变量。故我们需要估计概率模型p(x,z)的参数θ，但是由于里面包含隐含变量z，所以很难用最大似然求解，但如果z知道了，那我们就很容易求解了。
对于参数估计，我们本质上还是想获得一个使似然函数最大化的那个参数θ，现在与最大似然不同的只是似然函数式中多了一个未知的变量z，见下式（1）。也就是说我们的目标是找到适合的θ和z让L(θ)最大。那我们也许会想，你就是多了一个未知的变量而已啊，我也可以分别对未知的θ和z分别求偏导，再令其等于0，求解出来不也一样吗？

但是可以看到里面有“和的对数”，求导后形式会非常复杂（自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)复合函数的求导），所以很难求解得到未知参数z和θ。那OK，我们可否对（1）式做一些改变呢？我们看（2）式，（2）式只是分子分母同乘以一个相等的函数，还是有“和的对数”啊，还是求解不了，那为什么要这么做呢？咱们先不管，看（3）式，发现（3）式变成了“对数的和”，那这样求导就容易了。我们注意点，还发现等号变成了不等号，为什么能这么变呢？这就是Jensen不等式的大显神威的地方。
步骤三用到了Jenson不等式，不等式的原理如下：
设f是定义域为实数的函数，如果对于所有的实数x。如果对于所有的实数x，f(x)的二次导数大于等于0，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的，那么f是凸函数。如果只大于0，不等于0，那么称f是严格凸函数。
Jensen不等式表述如下：
如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])
特别地，如果f是严格凸函数，当且仅当X是常量时，上式取等号。
Jensen不等式如下图所示：

OK，到这里，现在式（3）就容易地求导了，但是式（2）和式（3）是不等号啊，式（2）的最大值不是式（3）的最大值啊，而我们想得到式（2）的最大值，那怎么办呢？
上面的式（2）和式（3）不等式可以写成：似然函数L(θ)>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到它的最大值。
我们固定θ，调整Q(z)使下界J(z,Q)上升至与L(θ)在此点θ处相等，然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（θt到θt+1），然后再固定θ，调整Q(z)……直到收敛到似然函数L(θ)的最大值处的θ*。
在Jensen不等式中说到，当自变量X是常数的时候，等式成立。而在这里，即：

再推导下，由于

（因为Q是随机变量z⁽ⁱ⁾的概率密度函数），则可以得到：分子的和等于c（分子分母都对所有z⁽ⁱ⁾求和：多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），则：

至此，我们推出了在固定参数θ后，使下界拉升的Q(z)的计算公式就是后验概率，解决了Q(z)如何选择的问题。这一步就是E步，建立L(θ)的下界。接下来的M步，就是在给定Q(z)后，调整θ，去极大化L(θ)的下界J（在固定Q(z)后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：

EM算法（Expectation-maximization）： 期望最大算法是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。

EM的算法流程：

初始化分布参数θ；

重复以下步骤直到收敛：

E步骤：根据参数初始值或上一次迭代的模型参数来计算出隐性变量的后验概率，其实就是隐性变量的期望。作为隐藏变量的现估计值：

M步骤：将似然函数最大化以获得新的参数值：

这个不断的迭代，就可以得到使似然函数L(θ)最大化的参数θ了。

示例

简单给出一个EM算法的示例，代码如下：

#-*- conding:utf-8 -*-
# --encoding:utf-8 --
"""
实现GMM高斯混合聚类
根据EM算法流程实现这个流程
http://scipy.github.io/devdocs/generated/scipy.stats.multivariate_normal.html#scipy.stats.multivariate_normal
"""
import numpy as np
from scipy.stats import multivariate_normal


def train(x, max_iter=100):
    """
    进行GMM模型训练，并返回对应的μ和σ的值(假定x数据中的簇类别数目为2)
    :param x: 输入的特征矩阵x
    :param max_iter:  最大的迭代次数
    :return:  返回一个五元组(pi, μ1， μ2，σ1，σ2)
    """
    # 1. 获取样本的数量m以及特征维度n
    m, n = np.shape(x)

    # 2. 初始化相关变量
    # 以每一列中的最小值作为mu1，mu1中的元素数目就是列的数目(n)个
    mu1 = x.min(axis=0)
    mu2 = x.max(axis=0)
    sigma1 = np.identity(n)
    sigma2 = np.identity(n)
    pi = 0.5

    # 3. 实现EM算法
    for i in range(max_iter):
        # a. 初始化多元高斯分布（初始化两个多元高斯混合概率密度函数）
        norm1 = multivariate_normal(mu1, sigma1)
        norm2 = multivariate_normal(mu2, sigma2)

        # E step
        # 计算所有样本数据在norm1和norm2中的概率
        tau1 = pi * norm1.pdf(x)
        tau2 = (1 - pi) * norm2.pdf(x)
        # 概率做一个归一化操作
        w = tau1 / (tau1 + tau2)

        # M step
        mu1 = np.dot(w, x) / np.sum(w)
        mu2 = np.dot(1 - w, x) / np.sum(1 - w)
        sigma1 = np.dot(w * (x - mu1).T, (x - mu1)) / np.sum(w)
        sigma2 = np.dot((1 - w) * (x - mu2).T, (x - mu2)) / np.sum(1 - w)
        pi = np.sum(w) / m

    # 返回最终解
    return (pi, mu1, mu2, sigma1, sigma2)


if __name__ == '__main__':
    np.random.seed(28)

    # 产生一个服从多元高斯分布的数据（标准正态分布的多元高斯数据）
    mean1 = (0, 0, 0)  # x1\x2\x3的数据分布都是服从正态分布的，同时均值均为0
    cov1 = np.diag((1, 1, 1))
    data1 = np.random.multivariate_normal(mean=mean1, cov=cov1, size=500)

    # 产生一个数据分布不均衡
    mean2 = (2, 2, 3)
    cov2 = np.array([[1, 1, 3], [1, 2, 1], [0, 0, 1]])
    data2 = np.random.multivariate_normal(mean=mean2, cov=cov2, size=200)

    # 合并两个数据
    data = np.vstack((data1, data2))

    pi, mu1, mu2, sigma1, sigma2 = train(data, 100)
    print("第一个类别的相关参数:")
    print(mu1)
    print(sigma1)
    print("第二个类别的相关参数:")
    print(mu2)
    print(sigma2)

    print("预测样本属于那个类别(概率越大就是那个类别)：")
    norm1 = multivariate_normal(mu1, sigma1)
    norm2 = multivariate_normal(mu2, sigma2)
    x = np.array([0, 1, 0])
    print(pi * norm1.pdf(x)) # 属于类别1的概率为:0.0275  => 0.989
    print((1 - pi) * norm2.pdf(x))# 属于类别1的概率为:0.0003 => 0.011

运行结果如下：

第一个类别的相关参数:
[-0.03813579 -0.0265211   0.06435217]
[[ 0.89477712 -0.01408823 -0.0630486 ]
 [-0.01408823  1.04926061  0.05625028]
 [-0.0630486   0.05625028  1.01136971]]
第二个类别的相关参数:
[1.89859865 1.9345104  2.77821259]
[[0.84837261 1.11138647 1.12554675]
 [1.11138647 2.19189624 1.37848616]
 [1.12554675 1.37848616 3.43711088]]
预测样本属于那个类别(概率越大就是那个类别)：
0.02758258224417439
0.00039921371592994276

Process finished with exit code 0