Expectation Maximization Algorithm

EM算法和之前学的都不太一样，EM算法更多的是一种思想，所以后面用几个例子讲解，同时也会重点讲解GMM高斯混合模型。

①极大似然估计

极大似然估计这里面用的比较多。假设我们想要知道我们学生身高的分布，首先先假设这些学生都是符合高斯分布我们要做的就是要估计这两个参数到底是多少。学生这么多，挨个挨个来肯定是不切实际的，所以自然就是抽样了。
为了统计学生身高，我们抽样200个人组成样本
我们需要估计的参数首先估计一下抽到这两百人的概率一共是多少，抽到男生A的概率抽到学生B的概率所以同时抽到这两个学生的概率就是那么同时抽到这200个学生的G概率
最后再取一个对数就好了：

notation和log

上面有一条公式里面是同时存在了;和|，这两个符号差别其实有点大的。|一般我们是用来表示条件概率，比如就是表示x在θ的条件下发生的概率。也是一个意思。
分号;表示的就是表示后面的是待估计的参数，也就是说P(x;θ)意思就是后面的θ是需要估计的参数而不是条件，所以|也有另一层意思，如果不是表示条件概率，那么就是表示后面有待估计参数。
当然是在|不表示条件概率的情况下。
这两种表示法是源于两种学派的不同理解：
频率派认为参数为固定的值，是指真实世界中，参数值就是某个定值。
贝叶斯派认为参数是随机变量，是指取这个值是有一定概率的。当然，无论是;还是|，他们都是表示条件概率的意思，只不过这两个学派理解不一样而已。
notation的问题解决了之后就是log的问题了，为什么需要log化，讲道理，是不需要的。但是求log有这么几个好处：
1.计算简单，累乘是很难计算的，log之后可以变换成累加。
2.概率累乘是会出现数字非常小的情况，log之后就可以避免这种情况。
3.log之后函数的梯度方向是没有变化的，对于函数优化的方向影响很小。

似然函数的执行步骤：
1.得到似然函数
2.取对数整理
3.求导数，另导数为零
4.解方程得到解

②Jensen不等式

首先引出凸函数的概念那么就是凸函数，所以它的图像就是一个勾形的，看起来是一个凹函数，实际上是凸函数。

凸函数的性质很显而易见了

其实很明显的，看图就知道了，E(x)其实就是a到b的平均值，上面的结论很容易证实。那么如果想要取到等号，需要改变上面？取等号的意思就是相切，相切的意思就是a = b，既然a = b了，那么x自然就是常数了，所以当且仅当如果是凹函数，那么就是方向相反了。

③EM算法的推导

正常来看先是要引入一个最大似然函数：但这样其实是和难求的，P(x|θ)完全混在了一起，根本求不出来，所以我们要引入一个辅助变量z。

隐变量Z

隐变量是观测不到的，比如做一个抽样，从3个袋子里面抽取小球。而抽取这些小球的过程是不可见的，抽取的过程其实就是隐变量，而这些隐变量，也就是过程可以告诉我们这个x是从哪个袋子来的，由此来区分。这个隐变量和HMM里面的隐含序列不是一个东西，隐含序列是当前x的状态，而不是抽取x的过程。所以在这里，通俗点讲，这个z就是用来找到x的组类的，也就是说z来告诉这个x你是属于哪一组的。
另外需要注意的是，隐变量是不可以随便添加的，添加隐变量之后不能影响边缘概率。也就是说，原来是P(x)，添加之后就是P(x,z)，那么必须有：

所以我们引入隐变量的原因是为了转化成和这几个高斯模型相关的式子，否则无从下手。化简一下上式子：既然z可以指定x，那么我们只需要求解出z就好了。
注意上面凸函数所提到的一个期望性质，这里就可以使用了。因为虽然优化了上面的式子，还是不能求出来，因为z变量实在是太抽象了，找不到一个合适的公式来表示它。EM的一个方法就是用优化下界函数的方法来达到优化目标函数的目的。
既然z很抽象，那么我们就需要一个转变一下。对于每一个样例x都会对应一个z，那么假设一个分布Q(z)是满足了z的分布的，而Q(z)满足的条件是Qi意味着每一个x对应的z都会对应着一个Q了，这里有点复杂，再详细解释一下。一个x对应一组z，z是一个向量，但是每一个z又会分别对应一个一个分布Q。以为最后得到的z不会是一个数字，而是一个概率，也就是说Q(z)得到的是这个x样例属于这个类别的概率是多少。而z的数量，一个是当前有多少个分布混合在一起的数量。
再梳理一下：现在的样本是xi，那么每一个xi将会对应着一组的z，每一个xi同时也会对应着一个分布Qi，z其实就是反应了这个样本是来自于哪个分布的。比如这个x是A1分布做了3，A2分布做了5，那么z可能就是={3,5}。所以Qi(z)得到的是这个x属于这些个分布的概率，也就是说这些分布对x做了多少百分比的功，自然就是要等于1了。
还要注意的是，上面的这个并不能得到Qi(z)就是分布对x做了多少功的结论，得到这个结论是后面下界函数与目标函数相等得到的。这里只是知道了总和等于1，因为是分布的总和嘛。
现在就到了公式的化简：
仔细看一下这个式子这个式子其实就是求的期望，假设，那么可以利用上面。于是化简：
这个时候就得到了下界函数，上面也讲过了，想要相等，自然就是x要是常数，所以既然，而且z也是一样的，因为一个样本嘛。所以上下加和（如果是离散的，那就sum一下，连续的那就积分，这里是离散的，所以就是sum一下）。于是有
于是有：

所以Q(z)计算的是后验概率。计算在当前参数θ和已经抽样到x的条件下，这个x是从z来的概率。其实就是z对x做了多少贡献。
所以整个EM算法步骤就很清晰了：

EM算法计算步骤：

E-step：
对于每一个，求
M-step：

这时候就可以使用求导迭代的方法求解了。

这就是整一个EM算法的框架了，可以看到其实没有比较具体的算法，大致上就是一个框架。那么问题来了，怎么样证明这东西是一个收敛的？？

证明EM算法的收敛性

既然是极大似然估计，那么需要证明的自然就是，也就是说极大似然估计要单调递增。因为每一次都是极大，那么随着时间增加，自然就是要增大的了。
当给定一个的时候，相当于t时间的EM算法进行了E步了，所以有：
然后就走M步，M步之前还没有求最大，所以:这一步有点复杂，这里是M-step了，是L(θ)求极大值得到的，而最重要的是L和后面下界函数的关系。这两个东西是不相等的，下界函数等于目标函数的条件是x为constant，x == constant意味着Q得求出来，这里看这个Q貌似是出来了，但是这个Q是，是上一个时刻的，而不是这个时刻的，所以t+1时刻Q还没有被固定住，也就是没有满足x == constant的条件，所以下界函数小于等于目标函数。接下来就简单了，直接把换成就好了。

这样就证明了收敛性。

EM algorithm的优化方法：

之前讨论过，这种方法是迭代，使用极大似然估计和迭代的方法来进行优化，但实际上这更像是一种坐标上升的方法：

一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步：固定 θ，优化Q；M步：固定 Q，优化 θ；交替将极值推向极大。Kmeans也是这样更新的，固定中心点，优化Ein；优化完成，更新中心点。SMO也是，固定更新α1和α2。

④GMM的推导

可以直接把高斯混合模型代入EM框架里面。
存在多个高斯分布混合生成了一堆数据X，取各个高斯分布的概率是，第i个高斯分布的均值是，方差是，求法φ，μ，σ。
按照套路，第一个E-step求出Q，于是有：
意思就是求出第i个样本属于第j个分布的概率是多少。之后就是M-step了，就是化简了：

这里可能需要解释一下，根据至于条件，因为很明显，z是隐变量，只是指明了x是属于哪个类别，和μ，Σ没有什么关系，所以直接忽略那两个参数了，所以P(z)是没有那两个参数的，z是代表了分布，所以每一个分布的概率肯定是包括了，所以就只有一个概率的参数。P(x|z)是本身的概率，就是已经知道分布是那个了，求属于这个分布的概率是多少，既然已经选定了分布那么自然就不需要再看φ了，因为φ是各个分布的概率。

之后就是求导数了：

导数等于0，于是有：μ就出来了。
对于φ，那就简单了，实际上需要优化的式子：
求导等于0。
要注意的是φ还有一个条件：。有条件，那么问题来了，不能直接求导，所以只能用拉格朗日乘子法求解，于是构造式子：
接着求导。于是得到，再改变一下，两边加和，得到，这个就简单了，，所以-β = m。那么：

Σ也是一样：

这样整个推导就算结束了。

⑤硬币模型

现在有两个硬币AB，进行5次试验每一次投10次，并不知道是哪个硬币投的，求两种硬币的正面的概率。
首先E-step：
首先先初始化一下，
第一个试验选中A的概率：
同样求得
计算机出每一个试验的概率然后相加求均值。
之后就是M-step了：

表示选择了A的概率，表示选择了B的概率。之后就是求导更新了。

⑥代码实现

方差的求解就不玩了，主要就是迭代求解μ和φ的值了。
首先是生成数据，4个高斯分布，每一个高斯分布的sigma都是一样的，不一样的只有μ和α，也就是φ，习惯上把前面的一个参数叫做权值，所以用α来表示。

def generate_data(sigma, N, mu1, mu2, mu3, mu4, alpha):
    global X
    X = np.zeros((N, 2))
    X = np.matrix(X)
    global mu
    mu = np.random.random((4, 2))
    mu = np.matrix(mu)
    global expect
    expect = np.zeros((N, 4))
    global alphas
    alphas = [0.25, 0.25, 0.25, 0.25]
    for i in range(N):
        if np.random.random(1) < 0.1:
            X[i, :] = np.random.multivariate_normal(mu1, sigma, 1)
        elif 0.1 <= np.random.random(1) < 0.3:
            X[i, :] = np.random.multivariate_normal(mu2, sigma, 1)
        elif 0.3 <= np.random.random(1) < 0.6:
            X[i, :] = np.random.multivariate_normal(mu3, sigma, 1)
        else:
            X[i, :] = np.random.multivariate_normal(mu4, sigma, 1)
    plt.title('Generator')
    plt.scatter(X[:, 0].tolist(), X[:, 1].tolist(), c = 'b')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.show()

这四个模型的比例分别是1:2:3:4，使用EM来找到他们属于的类别。

其实如果用kmeans聚类的话更加快速，但是这里还是用EM。
E-step：

def e_step(sigma, k, N):
    global X
    global mu
    global expect
    global alphas
    for i in range(N):
        W = 0
        for j in range(k):
            W += alphas[j] * math.exp(-(X[i, :] - mu[j, :]) * sigma.I * np.transpose(X[i, :] - mu[j, :])) / np.sqrt(np.linalg.det(sigma))
        for j in range(k):
            w = math.exp(-(X[i, :] - mu[j, :]) * sigma.I * np.transpose(X[i, :] - mu[j, :])) / np.sqrt(np.linalg.det(sigma))
            expect[i, j] = alphas[j]*w/W
            pass

就是按照公式来求解w即可，求解每一个分布对样本点做了多少的功，之后求单个样本点求比例。
M-step：

def m_step(k, N):
    global expect
    global X
    global alphas
    for j in range(k):
        mathor = 0
        son = 0
        for i in range(N):
            son += expect[i, j]*X[i, :]
            mathor += expect[i, j]
        mu[j, :] = son / mathor
        alphas[j] = mathor / N

直接按照公式优化即可。

if __name__ == '__main__':
    iterNum = 1000
    N = 500
    k = 4
    probility = np.zeros(N)
    mu1 = [5, 35]
    mu2 = [30, 40]
    mu3 = [20, 20]
    mu4 = [45, 15]
    sigma = np.matrix([[30, 0], [0, 30]])
    alpha = [0.1, 0.2, 0.3, 0.4]
    generate_data(sigma, N, mu1, mu2, mu3, mu4, alpha)
    for i in range(iterNum):
        print('iterNum : ', i)
        err = 0
        err_alpha = 0
        Old_mu = copy.deepcopy(mu)
        Old_alpha = copy.deepcopy(alphas)
        e_step(sigma, k, N)
        m_step(k, N)
        for z in range(k):
            err += (abs(Old_mu[z, 0] - mu[z, 0]) + abs(Old_mu[z, 1] - mu[z, 1]))
            err_alpha += abs(Old_alpha[z] - alphas[z])
        if (err <= 0.001) and (err_alpha < 0.001):
            print(err, err_alpha)
            break
    color = ['blue', 'red', 'yellow', 'green']
    markers  = ['<', 'x', 'o', '>']
    order = np.zeros(N)
    for i in range(N):
        for j in range(k):
            if expect[i, j] == max(expect[i, ]):
                order[i] = j
        plt.scatter(X[i, 0], X[i, 1], c = color[int(order[i])], alpha=0.5, marker=markers[int(order[i])])
    plt.show()
    print('standedμ:',mu4, mu3, mu2, mu1)
    print('standedα:',alpha)
    print('new μ:', mu)
    print('new α:',alphas)

运行函数。看看结果：

结果其实还是相差不大。达到预期。

⑦EM的另一种理解方法

上面所讲的其实只是一种理解方法，在李航老师的统计学习方法里面是另一种比较厉害的解法：
1.E-step：求出Q函数。
2.M-step：利用Q函数求极大值。
其实这两种方法是完全一样的，Q函数就是下界函数，

这个公式和上面的： 是一样的，至于为什么有下面这个 这是李航老师书上的，既然θ是已经知道了，那么自然就可以去掉，因为在log下面做分母就相当于是常数了，拆开可以变减号嘛，在前面的就不可以。回来看看我们刚刚推导的，其实是一样的，下面的那个Q确实可以去掉，因为是事先知道的了。在使用Jensen不等式的时候，需要假设隐变量服从某种形式的概率分布，才可以将推导过程的一部分看成是期望的表达形式从而应用Jensen不等式。然而这个分布不是随便指定的。我们令Jensen不等式取等号的时候，可以计算出这个分布其实就是：已知观测数据的隐变量的后验概率分布。由于求Q函数需要先求出隐变量的后验概率的期望，因此，这就可以解释为什么EM算法的“通俗”理解角度的E步骤是求隐变量的期望了。有时候在用EM算法解决某个具体问题的时候，会发现M步骤极大化的居然是完全数据的对数似然函数。这是因为，Q函数虽然是完全数据的对数似然函数的某种期望，但是求这个期望的过程有时其实就是将隐变量的后验概率的期望代入就可以了。因此，本质上我们其实还是在求Q函数的极大。
事实上，李航老师的Estep求出Q函数，其实就是通俗理解里面的求出Q(z)并代入下界函数的过程。因为求出Q(z)就相当于这个Q(z)被固定了，可以去掉了，于是就和李航老师的式子一样了。有机会再补一些李航老师的推导，勃大精深。

⑧总结

EM和Kmeans算法其实很类似，事实上步骤基本可以用EM框架来替换，但是Kmeans算法是硬分类，说一不二，但是EM算法不太一样，是软分类，百分之几是那个，百分之几是这个。

缺点也还是有的：初值敏感，局部最优。因为存在了隐变量，所以导致了直接对x做极大似然是不可行的，log已经在sum的外面了。所以EM算法就转向了下界函数，而这种方法本来就不保证找到局部最优解。

如果将样本看作观察值，潜在类别看作是隐藏变量，那么聚类问题也就是参数估计问题。如果一个目标函数存在多个变量，那么梯度下降牛顿法这些逼近方法就用不了了。但我们可以使用坐标上升方法，固定一个变量，对另外一个求导数，然后替换最后逐步逼近极值点。对应到EM算法也是一样，E步求隐含的z变量，Mstep求解其他参数。

最后符上GitHub代码：

https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/ExpectationMaximization