EM算法（期望最大算法、Expectation Maximization Algorithm)

引言

EM算法，全称为期望最大（Expectation Maximization）算法，是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。
我们采用Nature Biotech中EM tutorial中掷硬币的例子来对EM算法进行一个具体化的解释。

我们假设有两枚硬币A与B，他们随机抛掷的结果如图，当我们知道了每次抛的是A还是B就可以直接进行估计。

那么根据图中的计算，我们可以得到，两枚硬币抛出正面的概率分别为${\hat{\theta }}_A=0.80$,${\hat{\theta }}_B=0.45$
但如果此时不知道抛掷的硬币是A还是B，也就是隐藏了图中表格里的Coin A和Coin B，将所有结果存到一列表格中。此时相当于添加了隐变量。只能观测到5轮循环，每轮循环10次，共计50次投币。这个时候就无法直接估计A和B的正面概率，因为不知道哪一轮循环是哪个硬币投掷的。最后，图中的表格就变为下面的形式：

Coin	Statistics
Unknown	5H,5T
Unknown	9H,1T
Unknown	8H,2T
Unknown	4H,6T
Unknown	7H,3T

此时我们的目标没有变，仍旧是估计A和B正面的概率。虽然此时我们失去了A和B的标号，但是我们多了一个硬币种类的隐变量，设为Z,Z是一个5维的向量，可以表示为：
$Z=(z_1,z_2,z_3,z_4,z_5)$，代表每一轮所使用的硬币。但是，对于Z并不知晓，就无法去估计${\theta }_A$和${\theta }_B$，那么我们就必须先得出Z，然后才能进一步估计P(A)和P(B)，要得出Z又必须先知道${\theta }_A$和${\theta }_B$，然后采用最大似然去预测，那么就会陷入一个死循环。

为了解决上述问题，我们可以先随机初始化${\theta }_A$和${\theta }_B$，然后去估计Z，基于估计出来的Z按照最大似然去估计新的${\theta }_A$和${\theta }_B$，直至最后收敛，这就是EM算法的思想，下图给出了E-step和M-step。

下面计算引用【机器学习】EM——期望最大（非常详细）

随机初始化${\theta }_A=0.6$ 和 ${\theta }_B=0.5$，那么第一轮循环中，如果Unknown的硬币是A，得出5H5T的概率就是$0.6^5\ast 0.4^5$；如果是B，概率为$0.5^5\ast 0.5^5$。那么就有，Unknown为硬币A和B的概率分别是
$$\begin{gathered} P_A=\frac{0.6^5\ast 0.4^5}{0.6^5\ast 0.4^5+0.5^5\ast 0.5^5}=0.45 \\ {P}_B=\frac{0.5^5\ast 0.5^5}{0.6^5\ast 0.4^5+0.5^5\ast 0.5^5}=0.55 \end{gathered}$$

5个轮次我们都依照上面进行计算，可以得到如下情况

Coin A	Coin B
0.45	0.55
0.80	0.20
0.73	0.27
0.35	0.65
0.65	0.35

这一步对应了E-step

结合硬币A的概率和上面的投掷结果情况，我们可以利用期望求出硬币A和硬币B的贡献。

Coin	Statistics
Unknown	5H,5T
Unknown	9H,1T
Unknown	8H,2T
Unknown	4H,6T
Unknown	7H,3T

Coin A	Coin B
0.45	0.55
0.80	0.20
0.73	0.27
0.35	0.65
0.65	0.35

以第一轮硬币A为例，计算方式为：
$$H:0.45\ast 5=2.25T:0.45\ast 5=2.25$$
上面的计算代表着，硬币A对于投掷结果是5H5T的贡献情况。于是可以得到下面的贡献分布情况。

Coin A	Coin B
2.2H,2.2T	2.8H,2.8T
7.2H,0.8T	1.8H,0.2T
5.9H,1.5T	2.1H,0.5T
1.4H,2.1T	2.6H,3.9T
4.5H,1.9T	2.5H,1.1T

总计有

21.3H,8.6T	11.7H,8.4T

然后使用极大似然来估计${\theta }_A$和${\theta }_B$
$$\begin{gathered} {\theta }_A=\frac{21.3}{21.3+8.6}=0.71 \\ {\theta }_B=\frac{11.7}{11.7+8.4}=0.58 \end{gathered}$$
这一步就是M-step，如此反复迭代，可以得到最终的${\theta }_A$和${\theta }_B$。

原理

有一个样本集 { x 1 . . . , x m } \{x_1...,x_m\} {x1​...,xm​}，我们假设样本间互相独立，想要拟合模型$p(x;\theta )$到数据的参数。想要找到每个样本隐含的类别z，使得p(x,z)最大，那么可以得到下面的极大似然估计：
$$\begin{gathered} L(\theta )=\sum _{i=1}^m\log p(x_i;\theta ) \\ =\sum _{i=1}^{m}log\sum _{z}p(x_i,z;\theta ). \end{gathered}$$

第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。所以我们需要用EM算法来求z。对于参数估计，我们本质上还是想获得一个使似然函数最大化的那个参数θ，现在与极大似然不同的只是似然函数式中多了一个未知的变量z。

事实上，隐变量估计问题也可以通过梯度下降等优化算法，但事实由于求和项将随着隐变量的数目以指数级上升，会给梯度计算带来麻烦；而 EM 算法则可看作一种非梯度优化方法。

对于每个样本，我们使用$Q_i(z)$表示样本i隐含变量z的某种分布，且$Q_i(z)$满足：${\sum }_z^Z{Q}_i(z)=1,Q_i(z)\ge 0$，那么我们的目标最后就变成了找到合适的$\theta$和z使得$L(\theta )$最大。

我们将式子做如下变换：
$$\begin{gathered} \sum _i^n\log p(x_i;\theta )=\sum _i^n\log \sum _{z}p(x_i,z;\theta ) \\ =\sum _i^n\log \sum _z^Z{Q}_i(z)\frac{p(x_i,z;\theta )}{Q_i(z)} \\ \ge \sum _i^n\sum _z^Z{Q}_i(z)\log \frac{p(x_i,z;\theta }{Q_i(z)} \end{gathered}$$

第一步是求和每个样本的所有可能的类别 z 的联合概率密度函数，但是这一步直接求导非常困难，所以将其分母都乘以函数 $Q_i(z)$ ，转换到第二步。从第二步到第三步是利用 Jensen 不等式。

Jensen不等式：

• 如果f是凸函数，X是随机变量，则 $E[f(X)]\ge f(E[X])$ ，当f是严格凸函数时，则$E[f(X)]>f(E[X])$
• 如果f是凹函数，X是随机变量，则$E[f(X)]\le f(E[X])$ ，当f是严格凹函数时，则$E[f(X)]<f(E[X])$。
• 当$X=E[X]$时，即为常数时等式成立。
  

我们把第一步中的$\log$函数看成一个整体，由于$\log (x)$的二阶导数小于0，所以原函数为凹函数。我们把$Q_i(z)$看成一个概率$p_z$，把$\frac{p(x_i,z;\theta )}{Q_i(z)}$看成 z的函数$g(z)$。根据期望公式有：
$$E(z)=p_{z}g(z)=\sum _z^Z{Q}_i(z)[\frac{p(x_i,z;\theta )}{Q_i(z)}]$$

根据Jensen不等式：
$$\begin{gathered} f(\sum _z^Z{Q}_i(z)[\frac{p(x_i,z;\theta )}{Q_i(z)}]=f(E[z])\ge E[f(z)] \\ =\sum _z^Z{Q}_i(z)f([\frac{p(x_i,z;\theta )}{Q_i(z)}]) \end{gathered}$$

通过上面我们得到了：$L(\theta )\ge J(z,Q)$ 的形式（z 为隐变量），那么我们就可以通过不断最大化$J(z,Q)$的下界来使得$L(\theta )$不断提高。

图片来自 【机器学习】EM——期望最大（非常详细）

或许我们可以直接分别对$\theta$和z求偏导数，然后求得极值点，也就是最大化$L(\theta )$，但是式子中含有$\log {\sum }_{z}p$的形式，求导后就会变得非常复杂，因此很难求得我们需要的z和$\theta$。

也就是说，EM 算法通过引入隐含变量，使用 MLE（极大似然估计）进行迭代求解参数。通常引入隐含变量后会有两个参数，EM 算法首先会固定其中的第一个参数，然后使用 MLE 计算第二个变量值；接着通过固定第二个变量，再使用 MLE 估测第一个变量值，依次迭代，直至收敛到局部最优解。

最后可以总结一下EM算法的步骤：
第一步，初始化分布参数$\theta$；
第二步，重复E-step 和 M-step直到收敛：

• E-step:根据参数的初始值或上一次迭代的模型参数来计算出的隐性变量的后验概率（条件概率），其实就是隐性变量的期望值。作为隐藏变量的现有估计值：
  $$Q_i(z_j)=p(z_j|x_i;\theta )$$
• M-step:最大化似然函数从而获得新的参数值：
  $$\theta ={agrmax}_{\theta }\sum _i\sum _j{Q}_i(z_j)\log \frac{P(x_i,z_j;\theta )}{Q_i(z_j)}$$

代码

import numpy as np
from scipy.stats import multivariate_normal  #多元正态分布
from numpy import genfromtxt    #数据处理函数

def Estep(X, prior, mu, sigma):
    N, D = X.shape         #N表示数据的条数，D表示数据的维数
    K = len(prior)          #计算随机变量prior的长度
    gama_mat = np.zeros((N,K))       #初始化gama_mat为N行K列个0
    for i in range(0,N):        #遍历N条数据
        xi = X[i, :]        #对于每一行数据
        sum = 0
        for k in range(0, K):     #对K个prior求概率密度函数
            p = multivariate_normal.pdf(xi, mu[k,:], sigma[k,:,:])         #计算第k个prior的概率密度p
            sum += prior[k] * p      #求K个prior的总和
        for k in range(0,K):
            gama_mat[i, k] = prior[k] * multivariate_normal.pdf(xi, mu[k,:], sigma[k,:,:]) / sum   #求gama_mat概率矩阵
    return gama_mat   
       
def Mstep(X, gama_mat):        #将Estep得到的gama_mat，作为Mstep的输入参数
    (N, D) = X.shape
    K = np.size(gama_mat, 1)     #返回gama_mat的列数
    mu = np.zeros((K, D))        #给mu初始化为K行D列个0 
    sigma = np.zeros((K, D, D))   #给sigma初始化为K个D行D列的矩阵（方阵）
    prior = np.zeros(K)         #初始化prior为K个0
    for k in range(0, K):
        N_k = np.sum(gama_mat, 0)[k]    
        for i in range(0,N):
            mu[k] += gama_mat[i, k] * X[i]         
        mu[k] /= N_k             #得到均值矩阵
        for i in range(0, N):
            left = np.reshape((X[i] - mu[k]), (2,1))
            right = np.reshape((X[i] - mu[k]), (1,2))
            sigma[k] += gama_mat[i,k] * left * right     #得到sigma矩阵
        sigma[k] /= N_k
        prior[k] = N_k/N
    return mu, sigma, prior
    
def logLike(X, prior, Mu, Sigma):          #对数似然函数
    K = len(prior)                         
    N, M = np.shape(X)
    
    P = np.zeros([N, K])      # 初始化概率矩阵P
    for k in range(K):
        for i in range(N):
            P[i, k] = multivariate_normal.pdf(X[i], Mu[k, :], Sigma[k, :, :])  #P是一个NxK矩阵，其中(i,k)个元素表示第i个数据点在第j个簇中的可能性
    return np.sum(np.log(P.dot(prior)))     #返回似然函数值
    
def main():
    # Reading the data file
    X = genfromtxt('.\data\TrainingData_GMM.csv', delimiter=',')   #读取数据
    print('training data shape:', X.shape)    #打印数据的shape（5000,2）

    N, D = X.shape     #把5000、2分别赋予N，D表示有5000条数据，每条数据有两个特征
    K = 4       

    # initialization
    mu = np.zeros((K, D))      #给mu初始化为K行D列个0
    sigma = np.zeros((K, D, D))     #给sigma初始化为K个D行D列的矩阵（方阵）
     #mu[0] = [-0.5, -0.5]
     #mu[1] = [0.3, 0.3]
     #mu[2] = [-0.3, 0.3]
     #mu[3] = [1.3, -1.3]

    mu[0] = [1, 0]
    mu[1] = [0, 1]
    mu[2] = [0, 0]
    mu[3] = [1, 1]

    for k in range(0, K):            #K个二维单位矩阵作为sigma
        sigma[k] = [[1, 0], [0, 1]]

    prior = np.ones(K) / K        #先验概率为K个1/k
    iter = 0
    prevll = -999999
    ll_evol = []      #定义一个空数组ll_evol用于后面装log likelihood值

    print('initialization of the params:')
    print('mu:\n', mu)
    print('sigma:\n', sigma)
    print('prior:', prior)

    # Iterate with E-Step and M-step
    while (True):
        W = Estep(X, prior, mu, sigma)   #调用Estep
        mu, sigma, prior = Mstep(X, W)    #通过Mstep更新初始的参数mu、sigma、prior
        ll_train = logLike(X, prior, mu, sigma)   #调用对数似然函数
        print('iter:',iter, 'log likelihood:',ll_train)   #返回第几次迭代和对应的似然值
        ll_evol.append(ll_train)

        iter = iter + 1
        if (iter > 150 or abs(ll_train - prevll) < 0.01):   #当迭代次数大于150次或ll_train的值接近-999999时结束操作
            break

        abs(ll_train - prevll)   
        prevll = ll_train     #更新prevll的值
  
    import pickle    #pickle模块能将对象序列化
    with open('results.pkl', 'wb') as f:     
        pickle.dump([prior, mu, sigma, ll_evol], f)    #把prior, mu, sigma, ll_evol数据以二进制形式写入results.pkl文件中   
        
if __name__ == '__main__':
    main()
    

training data shape: (5000, 2)
initialization of the params:
mu:
 [[1. 0.]
 [0. 1.]
 [0. 0.]
 [1. 1.]]
sigma:
 [[[1. 0.]
  [0. 1.]]

 [[1. 0.]
  [0. 1.]]

 [[1. 0.]
  [0. 1.]]

 [[1. 0.]
  [0. 1.]]]
prior: [0.25 0.25 0.25 0.25]
iter: 0 log likelihood: -7758.98161924987
iter: 1 log likelihood: -6842.325796847348
iter: 2 log likelihood: -5877.092929285351
iter: 3 log likelihood: -5368.007026989181
iter: 4 log likelihood: -5158.951024247597
iter: 5 log likelihood: -4711.421034949499
iter: 6 log likelihood: -3814.307460091173
iter: 7 log likelihood: -2974.4639436766956
iter: 8 log likelihood: -2696.7255823784203
iter: 9 log likelihood: -2428.8597242586584
iter: 10 log likelihood: -2182.0334455177353
iter: 11 log likelihood: -2023.8042796345283
iter: 12 log likelihood: -1970.824061153165
iter: 13 log likelihood: -1956.739299533951
iter: 14 log likelihood: -1951.841638776625
iter: 15 log likelihood: -1949.6197599504249
iter: 16 log likelihood: -1948.411570496814
iter: 17 log likelihood: -1947.6739124172018
iter: 18 log likelihood: -1947.1944657348135
iter: 19 log likelihood: -1946.87392793565
iter: 20 log likelihood: -1946.6573830294146
iter: 21 log likelihood: -1946.5107251106433
iter: 22 log likelihood: -1946.4114561101228
iter: 23 log likelihood: -1946.3443661539402
iter: 24 log likelihood: -1946.299098150319
iter: 25 log likelihood: -1946.2685983080376
iter: 26 log likelihood: -1946.248073038022
iter: 27 log likelihood: -1946.2342732733175
iter: 28 log likelihood: -1946.225002083594

import numpy as np
from scipy.stats import multivariate_normal    #导入多元正太分布函数
from numpy import genfromtxt  
import matplotlib.pyplot as pyplot
import pickle 
with open('results.pkl', 'rb') as f:
    [prior, mu, sigma, ll_evol] = pickle.load(f)   #反序列化对象，将文件中的数据解析为python对象

pyplot.plot(ll_evol, 'o')  #画出29个训练数据的ll_evol值
pyplot.show() 

print('prior:',prior)       #进过训练后的prior
print('mu:', mu)            #进过训练后的mu
print('sigma:', sigma)      #进过训练后的sigma

X = genfromtxt('.\data\TrainingData_GMM.csv', delimiter=',')
print('data shape:', X.shape)

sel_num = 500
X_sel = []
sel_idxs = []
while len(sel_idxs) < sel_num:
    idx = np.random.randint(0, 5000, 1)    #如果len(sel_idxs)小于500，则从0到4999中返回一个随机整数给idx
    while idx in sel_idxs:                     #如果返回的idx在sel_idxs中已有，则继续从0到4999中返回一个随机整数给idx
        idx = np.random.randint(0, 5000, 1)
    sel_idxs.append(idx[0])            #把idx放入数组sel_idxs中
X_sel = X[sel_idxs]  

def get_label(x, prior, mu, sigma):
    K = len(prior)
    p = np.zeros(K)
    for k in range(0,K):
        p[k] = prior[k] * multivariate_normal.pdf(x, mu[k,:], sigma[k,:,:])     #求k的概率p
    label = np.argmax(p)      #输出最大p对应的索引值（即label）
    return label

lbs = []
for i in range(0, sel_num):
    lb = get_label(X_sel[i], prior, mu, sigma)      #调用get_label函数传参后，得到相应的分类标签
    lbs.append(lb)

# plot
pyplot.scatter(X_sel[:,0], X_sel[:,1], marker='o', c=lbs)      #做出lbs的散点图
pyplot.show()

prior: [0.23       0.27274177 0.26683872 0.23041951]
mu: [[ 1.20354296 -1.19685788]
 [ 0.14438214  0.14615911]
 [-0.44133345 -0.45061832]
 [-0.40658589  0.32258308]]
sigma: [[[ 0.02258855 -0.00761021]
  [-0.00761021  0.02361336]]

 [[ 0.00885564  0.00186778]
  [ 0.00186778  0.00880396]]

 [[ 0.07027983  0.0373684 ]
  [ 0.0373684   0.06508573]]

 [[ 0.03447208 -0.01299104]
  [-0.01299104  0.03454537]]]
data shape: (5000, 2)

参考

【机器学习】EM——期望最大（非常详细）

EM算法详解+通俗例子理解

(EM算法)The EM Algorithm

What is the expectation maximization
algorithm?

Python——EM(期望极大算法)实战(附详细代码与注解)(一)