Python实现EM

1.EM算法简介

EM算法也称期望最大化（Expectation-Maximum,简称EM）算法，如果概率模型的变量都是观测变量（数据中可见的变量），则可以直接用极大似然估计，或者用贝叶斯估计模型参数。但是，当模型含有隐变量（数据中看不到的变量）时，就不能简单地使用这些估计方法，而应该使用含有隐变量的概率模型参数的极大似然估计法，也即EM算法。
　EM算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含数据（EM算法的E步），接着基于观察数据和猜测的隐含数据一起来极大化对数似然，求解我们的模型参数（EM算法的M步)。由于我们之前的隐藏数据是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。不过没关系，我们基于当前得到的模型参数，继续猜测隐含数据（EM算法的E步），然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。
　从上面的描述可以看出，EM算法是迭代求解最大值的算法，同时算法在每一次迭代时分为两步，E步和M步。一轮轮迭代更新隐含数据和模型分布参数，直到收敛，即得到我们需要的模型参数。

2.EM算法的推导

对于m个样本观察数据 $x=(x^{(1)},x^{(2)},...x^{(m)})$x=(x(1),x(2),...x(m))中，找出样本的模型参数θ, 极大化模型分布的对数似然函数如下：
　
$\theta = arg \max \limits_{\theta}\sum\limits_{i=1}^m logP(x^{(i)}|\theta)$θ=argθmax​i=1∑m​logP(x(i)∣θ)

如果我们得到的观察数据有未观察到的隐含数据 $z=(z^{(1)},z^{(2)},...z^{(m)})$z=(z(1),z(2),...z(m))，此时我们的极大化模型分布的对数似然函数如下：

$\theta = arg \max \limits_{\theta}\sum\limits_{i=1}^m logP(x^{(i)}|\theta) = arg \max \limits_{\theta}\sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}， z^{(i)}|\theta)$θ=argθmax​i=1∑m​logP(x(i)∣θ)=argθmax​i=1∑m​logz(i)∑​P(x(i)，z(i)∣θ)

上面这个式子是没有 办法直接求出θ的。因此需要一些特殊的技巧，我们首先对这个式子进行缩放如下：

$\sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}， z^{(i)}|\theta) = \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)}， z^{(i)}|\theta)}{Q_i(z^{(i)})} \\ \geq \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}， z^{(i)}|\theta)}{Q_i(z^{(i)})}$i=1∑m​logz(i)∑​P(x(i)，z(i)∣θ)=i=1∑m​logz(i)∑​Qi​(z(i))Qi​(z(i))P(x(i)，z(i)∣θ)​≥i=1∑m​z(i)∑​Qi​(z(i))logQi​(z(i))P(x(i)，z(i)∣θ)​

上面第(1)式引入了一个未知的新的分布 $Q_i(z^{(i)})$Qi​(z(i))，第(2)式用到了Jensen不等式：

$log\sum\limits_j\lambda_jy_j \geq \sum\limits_j\lambda_jlogy_j\;\;, \lambda_j \geq 0, \sum\limits_j\lambda_j =1$logj∑​λj​yj​≥j∑​λj​logyj​,λj​≥0,j∑​λj​=1

或者说由于对数函数是凹函数，所以有: $f(E(x))≥E(f(x))如果f(x)是凹函数$f(E(x))≥E(f(x))如果f(x)是凹函数，此时如果要满足Jensen不等式的等号，则有： $\frac{P(x^{(i)}， z^{(i)}|\theta)}{Q_i(z^{(i)})} =c, c为常数$Qi​(z(i))P(x(i)，z(i)∣θ)​=c,c为常数

由于 $Q_i(z^{(i)})$Qi​(z(i))是一个分布，所以满足： $\sum\limits_{z}Q_i(z^{(i)}) =1$z∑​Qi​(z(i))=1

从上面两式，我们可以得到：

$Q_i(z^{(i)}) = \frac{P(x^{(i)}， z^{(i)}|\theta)}{\sum\limits_{z}P(x^{(i)}， z^{(i)}|\theta)} = \frac{P(x^{(i)}， z^{(i)}|\theta)}{P(x^{(i)}|\theta)} = P( z^{(i)}|x^{(i)}，\theta))$Qi​(z(i))=z∑​P(x(i)，z(i)∣θ)P(x(i)，z(i)∣θ)​=P(x(i)∣θ)P(x(i)，z(i)∣θ)​=P(z(i)∣x(i)，θ))

如果 $Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)}，\theta))$Qi​(z(i))=P(z(i)∣x(i)，θ)), 则第(2)式是我们的包含隐藏数据的对数似然的一个下界。如果我们能极大化这个下界，则也在尝试极大化我们的对数似然。即我们需要最大化下式：

$arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}， z^{(i)}|\theta)}{Q_i(z^{(i)})}$argθmax​i=1∑m​z(i)∑​Qi​(z(i))logQi​(z(i))P(x(i)，z(i)∣θ)​

去掉上式中为常数的部分，则我们需要极大化的对数似然下界为：

$arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}， z^{(i)}|\theta)}$argθmax​i=1∑m​z(i)∑​Qi​(z(i))logP(x(i)，z(i)∣θ)

上式也就是我们的EM算法的M步，那E步呢？注意到上式中 $Q_i(z^{(i)})$Qi​(z(i))是一个分布，因此 $\sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}， z^{(i)}|\theta)}$z(i)∑​Qi​(z(i))logP(x(i)，z(i)∣θ)可以理解为 $logP(x^{(i)}， z^{(i)}|\theta)$logP(x(i)，z(i)∣θ)基于条件概率分布 $Q_i(z^{(i)})$Qi​(z(i))的期望。
至此，我们理解了EM算法中E步和M步的具体数学含义。

3.EM算法步骤

输入：观测变量数据Y，隐变量数据Z,联合分布 $P(Y,Z |\theta)$P(Y,Z∣θ)，条件分布 $P(Z |Y,\theta)$P(Z∣Y,θ)；

输出：模型参数θ。

(1) 选择参数的初值 $θ^{(0)}$θ(0),开始迭代；

(2) E步：记 $θ^{(i)}$θ(i)为第i次迭代参数θ的估计值，在第i+1次迭代的E步，计算

$Q(θ,θ^{(i)})=E_z[logP(Y,Z|θ)|Y,θ^{(i)})]\\=\sum\limits_{Z}logP(Y,Z|θ)P(Z|Y,θ^{(i)})$Q(θ,θ(i))=Ez​[logP(Y,Z∣θ)∣Y,θ(i))]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))

这里， $P(Z|Y,θ^{(i)}$P(Z∣Y,θ(i)是在给定观测数据Y和当前的参数估计 $θ^{(i)}$θ(i)下隐变量数据z的条件概率分布；

(3) M步：求使屏幕 $Q^{(i+1)}$Q(i+1)极大化的θ，确定第i+1次迭代的参数的估计值 $θ^{(i+1)}$θ(i+1)

$θ^{(i+1)}=arg \max\limits_{θ}Q(θ,θ^{(i+1)})$θ(i+1)=argθmax​Q(θ,θ(i+1))

(4）重复第(2)步和第(3)步，直到收敛。

式 $Q(θ,θ^{(i)})=E_z[logP(Y,Z|θ)|Y,θ^{(i)})]\\=\sum\limits_{Z}logP(Y,Z|θ)P(Z|Y,θ^{(i)})$Q(θ,θ(i))=Ez​[logP(Y,Z∣θ)∣Y,θ(i))]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))的函数 $Q(θ,θ^{(i+1)})$Q(θ,θ(i+1))是EM算法的核心，称为Q函数（Q function)。

4.EM算法中的Q函数

定义(Q函数）完全数据的对数似然函数 $logP(Y,Z|θ)$logP(Y,Z∣θ)关于在给定观测数据Y和当前参数 $θ^{(i)}$θ(i)下对未观测数据Z的条件概率分布 $P(Z |Y,\theta^{(i)})$P(Z∣Y,θ(i))的期望称为Q函数，即

$Q(θ,θ^{(i)})=E_z[logP(Y,Z|θ)|Y,θ^{(i)})]$Q(θ,θ(i))=Ez​[logP(Y,Z∣θ)∣Y,θ(i))]

下面关于EM算法作几点说明：

步骤(1)参数的初值可以任意选择，但需注意EM算法对初值是敏感的。

步骤(2)E步求 $Q(θ,θ^{(i)})$Q(θ,θ(i))。Q函数式中Z是未观测数据，Y是观测数据。注意， $Q(θ,θ^{(i)})$Q(θ,θ(i))的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。

步骤(3)M步求 $Q(θ,θ^{(i)})$Q(θ,θ(i))的极大化，得到 $θ^{(i+1)}$θ(i+1)，完成一次迭代 $θ^{(i)}\rightarrowθ^{(i+1)}$θ(i)→θ(i+1)。后面将证明每次迭代使似然函数增大或达到局部极值。

步骤(4)给出停止迭代的条件，一般是对较小的正数 $ε_1,ε_2$ε1​,ε2​,若满足

$||θ^{(i+1)}-θ^{(i)}||&lt;ε_1 或||Q(θ^{(i+1)},θ^{(i)})-Q(θ^{(i)},θ^{(i)})||&lt;ε_2$∣∣θ(i+1)−θ(i)∣∣<ε1​或∣∣Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∣∣<ε2​

则停止迭代。

5.GMM算法使用Python实现

EM算法的一个重要应用是高斯混合模型(GMM)的参数估计。在许多情况下，EM算法是学习高斯混合模型的有效方法，敲公式太麻烦了，这里直接放图了，邹博-机器学习。邹博老师的公式比李航老师的要更好理解一些，公式原理是一致的。

E步计算gama值

tau1 = pi * multivariate_normal.pdf(data,mu1,sigma1)      # 概率密度，tau1，例如男性的概率，分子
tau2 = (1 - pi) * multivariate_normal.pdf(data,mu2,sigma2)# 概率密度，tau2，例如女性的概率
gamma = tau1 / (tau1 + tau2)  # 共有m个，每个样本属于tau1类的概率,Nk=sum(γ(i,k))

M步计算mu、sima、pi

mu1 = np.dot(gamma, data) / np.sum(gamma)                            # μ=γ(i,k)*x/Nk
mu2 = np.dot((1 - gamma), data) / np.sum((1 - gamma))     
sigma1 = np.dot(gamma * (data - mu1).T, data - mu1) / np.sum(gamma)  # ∑=γ(i,k)*(x-μ)^T*(x-μ)/Nk
sigma2 = np.dot((1 - gamma) * (data - mu2).T, data - mu2) / np.sum(1 - gamma)
pi = np.sum(gamma) / m                                    # π=Nk/N
print( i,'行均值,','mu1:',mu1,'mu2:',mu2)

使用python代码实现出来是不是非常简单呢，我们只是用了python中的scipy.stats计算了两个高斯分布的概率密度，其实自己编程实现高斯分布的概率密度函数也不复杂的。根据数据分布不同我们完全可以带入二项分布等公式计算概率，EM算法公式推导及公式记号显得极其复杂，在实际工程中使用起来极其简答，EM算法隐变量的求解思想与SMO算法，坐标轴下降法等类似。

EM算法运行结果如下：


以上所有源代码如下：

# -*- coding: utf-8 -*-
""" @Time : 2018/12/13 13:24 @Author : hanzi5 @Email : [email protected] @File : EM.py @Software: PyCharm """
import numpy as np
import scipy as sc
#import pandas as pd
from scipy.stats import multivariate_normal
import matplotlib
import matplotlib.pyplot as plt
from sklearn.metrics.pairwise import pairwise_distances_argmin

matplotlib.rcParams['font.family']='SimHei'  # 用来正常显示中文
plt.rcParams['axes.unicode_minus']=False     # 用来正常显示负号

if __name__ == '__main__':
    print('\n1、EM，开始')
    np.random.seed(100)       # 设置随机数种子
    mu1_fact = (0, 0)         # 设置第一组数据均值mu，两个维度均值都是0
    cov_fact = np.identity(2) # 设置协方差矩阵，单位阵
    data1 = np.random.multivariate_normal(mu1_fact, cov_fact, 400) # 随机产生400条符合mu1_fact、cov_fact高斯分布的数据
    mu2_fact = (5, 5)         # 设置第而组数据均值mu，两个维度均值都是5，与上一组数据分的更开
    data2 = np.random.multivariate_normal(mu2_fact, cov_fact, 100) # 随机产生100条符合mu2_fact、cov_fact高斯分布的数据
    data = np.vstack((data1, data2))           # 两组数据合并，总500条
    y = np.array([True] * 400 + [False] * 100) # 设置y数据，前400为true后100为false
    
    max_loop = 1000         # EM算法循环迭代最大次数
    m, n = data.shape       #数据行、列数
    # 方法一，随机指定mu1及mu2
    # mu1 = np.random.standard_normal(n)
    # mu2 = np.random.standard_normal(n)

    # 方法而，不随机产生mu1及mu2
    mu1 = data.min(axis=0)  # 第一组数据区数据中最小值作为初始值
    mu2 = data.max(axis=0)  # 第二组数据区数据中最大值作为初始值
    sigma1 = np.identity(n) # 使用单位阵作为初始值
    sigma2 = np.identity(n) # 使用单位阵作为初始值
    pi = 0.5                # 每组的概率　，EM算法对初值是
    
    # EM算法
    for i in range(max_loop):
        # E Step
        mu1_old=mu1.copy()  # 记录上一轮的mu1，用于判断跳出循环
        mu2_old=mu2.copy()  # 记录上一轮的mu2，用于判断跳出循环
        tau1 = pi * multivariate_normal.pdf(data,mu1,sigma1)      # 概率密度，tau1，例如男性的概率，分子
        tau2 = (1 - pi) * multivariate_normal.pdf(data,mu2,sigma2)# 概率密度，tau2，例如女性的概率
        gamma = tau1 / (tau1 + tau2)  # 共有m个，每个样本属于tau1类的概率,Nk=sum(γ(i,k))
        
        # M Step
        mu1 = np.dot(gamma, data) / np.sum(gamma)                            # μ=γ(i,k)*x/Nk
        mu2 = np.dot((1 - gamma), data) / np.sum((1 - gamma))     
        sigma1 = np.dot(gamma * (data - mu1).T, data - mu1) / np.sum(gamma)  # ∑=γ(i,k)*(x-μ)^T*(x-μ)/Nk
        sigma2 = np.dot((1 - gamma) * (data - mu2).T, data - mu2) / np.sum(1 - gamma)
        pi = np.sum(gamma) / m                                    # π=Nk/N
        print( i,'行均值,','mu1:',mu1,'mu2:',mu2)
        
        # 判断前后两次均值mu变化情况，变化非常小时，停止迭代
        epsilon = 0.00001
        if (((mu1-mu1_old)<epsilon).all()) and (((mu2-mu2_old)<epsilon).all()):
            break
        
    print( '类别概率:\t', pi)
    print( '均值:\t', mu1, mu2)
    print( '方差:\n', sigma1, '\n\n', sigma2, '\n')
    
    
    # 预测分类
    norm1 = multivariate_normal(mu1, sigma1)
    norm2 = multivariate_normal(mu2, sigma2)
    tau1 = norm1.pdf(data)
    tau2 = norm2.pdf(data)

    fig = plt.figure(figsize=(13, 7), facecolor='w')
    ax = fig.add_subplot(121)
    ax.scatter(data[:, 0], data[:, 1], c='b', s=30, marker='o')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_title('原始数据', fontsize=18)
    ax = fig.add_subplot(122)
    
    
    order = pairwise_distances_argmin([mu1_fact, mu2_fact], [mu1, mu2], metric='euclidean')
    if order[0] == 0:
        c1 = tau1 > tau2
    else:
        c1 = tau1 < tau2
    c2 = ~c1
    acc = np.mean(y == c1)
    print( '准确率：%.2f%%' % (100*acc))
    ax.scatter(data[c1, 0], data[c1, 1], c='r', s=30, marker='o')
    ax.scatter(data[c2, 0], data[c2, 1], c='g', s=30, marker='^')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_title('EM算法分类', fontsize=18)
    plt.tight_layout()
    plt.show()

参考资料：
1、《机器学习实战》Peter Harrington著
2、《机器学习》西瓜书，周志华著
3、 斯坦福大学公开课 ：机器学习课程 4、机器学习视频，邹博 5、《统计学习方法》李航