1.EM算法是什么

 EM算法的英文全称是Expectation Maximization Algorithm——期望极大化算法，它采用迭代的方式逼近带隐变量的似然函数。通过对似然函数的一个下界求偏导，得到每一步参数估计的过程。
 这个名称由于缺乏作用对象，让人一头雾水。这里的期望是什么？为什么我们要极大化这个期望，我们试图优化什么？
 这里的期望的含义其实是针对极大似然估计中的似然函数来说的，这里的期望就是似然函数的一个下界，我们的目的是求这样一个期望: 这个下界是根据詹森不等式(Jensen's inequality）放缩得到的，为什么要放缩呢？因为我们试图找出一个下界，极大化这个带参数的下界之后，可以无限近似于似然函数。你想，如果这个做法ok的话，意味着什么？意味着我们可以通过这个过程找出极大似然估计或最大后验估计的参数近似解。这也意味着我们可以搞一个迭代法来得到一个近似解。但是即便我说的天花乱坠，这个下界要是不收敛那也白搭。而下界要收敛必须满足两个条件：
 1.这个下界的取值要单调递增（因为每回迭代的结果要比上一次的取值更大）
 2.这个下界必须有上界(这个上界就是对数似然函数，且这一点可以由詹森不等式保证，这个也是EM的核心)
大白话就是单调有界必有极限。

EM算法收敛性证明：

我们来证明一下它确实是收敛的。
 首先，在极大似然估计中，我们的目的是根据手头上的个样本，最大化后，将参数估计出来；引入对数：;此时引入辅助变量;我们的对数似然函数就变成了：

设置变分函数：;那么：

根据琴生不等式,对数函数为凸函数(因为：等号在为常数时取到)：

上面的这个下界，就是用来逼近原对数似然函数的，这里我们已经证明了算法收敛的一个条件，有界性；但是在继续进行下一步的时候，我们还有一个问题没搞清楚，那就是变分函数的具体形式,实际上，我们可以通过琴生不等式等号成立的条件导出我们要的这个变分函数q。
令为常数：
接着我们代入变分函数的形式，定义这个下界的第一项：

定义下界的第二项：

对于第二项，我们看看随着迭代步数的增大，它是否是递增的，

我们将不同参数的与看作是两个分布，那么这个结果实际上是一个KL散度，它表征两个分布的相似度,其结果总是大于等于0的。
大于等于0的原因：

所以：

H为一个递增的序列，那么剩下的就是Q是否递增了，基于前面提到的这个下界是有上界的，上界就是我们的对数似然函数。在这个前提下，现在我们只要证明，Q递增是整个下界递增的充分必要条件即可。
必要性：

当整个下界递增，即：
那么：
所以 单调递增，必要性得证。
充分性：
因为：
前面已经证明：

又因为：

所以：

即，在递增的情况下，整个下界递增。
充分性得证。
证毕。

回答最初的问题：

 这个算法名称里提及的期望究竟是什么？
我们说了这么多，实际都是要做一件事，那就是：

由于前面证明过整个下界有界。且只要找到让第i次迭代的Q函数最大的作为下一次迭代的参数，我们就可以让Q递增，让算法收敛。
我们来看看Q的形式。

这就是为什么叫期望最大算法的原因。

2.概率PCA的EM推导与应用

 我们以概率PCA，来展开EM的应用：
 当然这里的应用只涉及变分函数已知情况下的应用，并不涉及广义EM的内容，日后看完文献在来唠唠广义EM，AVE，GAN等内容。
 我们先来算一下PPCA的EM期望的形式：

在概率PCA中，我们有提到：


所以：


所以期望里面是这个式子：

我们的目的是要估计出和;那么我们分别对它们求偏导:

所以：


因为：

代入偏导中

所以：

我们偏导得到的结果就是：

我们会发现我们还有两个估计量没解决，一个是一阶估计量,另一个是二阶估计量
在概率PCA中，我们提到过：

那么我们就有了一阶估计量：

二阶估计量可以通过简单的计算得到：

剩下的代入即可.

3.代码：

import numpy as np
class EM_Algorithm_For_PPCA:
    def __init__(self,data,q):
        self.X=data
        self.q=q
        self.n=np.shape(self.X)[0]
        self.p=np.shape(self.X)[1]
        self.H=np.identity(self.n)-(1/self.n)*np.ones([self.n,self.n])
        self.x_bar=np.average(self.X,axis=0)
        self.sigma_square=np.random.rand(1)
        self.W=np.random.random([self.p,self.q])
    def Moment_Calculation(self,x_i):
        self.inv_M=np.linalg.inv(self.W.T.dot(self.W)+self.sigma_square*np.identity(self.q))
        self.Ez=self.inv_M.dot(self.W.T).dot(np.array([x_i]).T-np.array([self.x_bar]).T)
        self.first_order_moment=self.Ez
        self.second_order_moment=self.sigma_square*self.inv_M+self.Ez.dot(self.Ez.T)
        return self.first_order_moment,self.second_order_moment
    def E_step(self,x_i): return self.Moment_Calculation(x_i)
    def M_step(self,x_i):
        self.x_i_minus_mu=(np.array([x_i]).T-np.array([self.x_bar]).T)
        self.Ez,self.Ez_EzT=self.E_step(x_i)
        self.W=(np.array([x_i]).T-np.array([self.x_bar]).T).dot(self.Ez.T).dot(np.linalg.inv(self.Ez_EzT))
        self.sigma_square=(1/self.p)*np.trace((1/self.n)*self.x_i_minus_mu.dot(self.x_i_minus_mu.T))-np.trace(self.Ez_EzT.dot(self.W.T.dot(self.W)))
        return self.W, self.sigma_square
    def iteration_step(self):
        for i,x_i in enumerate(self.X):
            self.W,self.sigma_square=self.M_step(x_i)
            print("The {} [info]: W = ".format(i),self.W,"  ","Sigma= ",self.sigma_square)
        return np.linalg.inv(self.W.T.dot(self.W)+self.sigma_square*np.identity(self.q)).dot(self.W.T).dot(self.X.T).dot(self.H)
if __name__ == '__main__':
    from sklearn.datasets import load_diabetes
    import seaborn as sns
    import matplotlib.pyplot as plt
    def plot(dfData):
        plt.subplots(figsize=(9, 9))
        sns.heatmap(dfData, annot=True, vmax=1, square=True, cmap="Blues")
        plt.show()
    data = load_diabetes()
    X = data["data"]
    EM=EM_Algorithm_For_PPCA(X,4)
    dimensional_reduction=EM.iteration_step()
    H=np.identity(np.shape(dimensional_reduction)[1])-(1/np.shape(dimensional_reduction)[1])*np.ones([np.shape(dimensional_reduction)[1],np.shape(dimensional_reduction)[1]])
    S=(1/np.shape(dimensional_reduction)[1])*dimensional_reduction.dot(H).dot(dimensional_reduction.T)

结果展示：

结果展示.png