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)})$中，找出样本的模型参数θ, 极大化模型分布的对数似然函数如下：
　
$\theta =arg\underset{\theta }{\max }\sum _{i=1}^{m}logP(x^{(i)}|\theta )$

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

$\theta =arg\underset{\theta }{\max }\sum _{i=1}^{m}logP(x^{(i)}|\theta )=arg\underset{\theta }{\max }\sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)}，z^{(i)}|\theta )$

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

$$\begin{gathered} \sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)}，z^{(i)}|\theta )=\sum _{i=1}^{m}log\sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{P(x^{(i)}，z^{(i)}|\theta )}{Q_i(z^{(i)})} \\ \ge \sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{P(x^{(i)}，z^{(i)}|\theta )}{Q_i(z^{(i)})} \end{gathered}$$

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

$log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log{y}_j,{\lambda }_j\ge 0,\sum _j{\lambda }_j=1$

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

由于$Q_i(z^{(i)})$是一个分布，所以满足：$\sum _z{Q}_i(z^{(i)})=1$

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

$Q_i(z^{(i)})=\frac{P(x^{(i)}，z^{(i)}|\theta )}{\sum _{z}P(x^{(i)}，z^{(i)}|\theta )}=\frac{P(x^{(i)}，z^{(i)}|\theta )}{P(x^{(i)}|\theta )}=P(z^{(i)}|x^{(i)}，\theta ))$

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

$arg\underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{P(x^{(i)}，z^{(i)}|\theta )}{Q_i(z^{(i)})}$

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

$arg\underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})logP(x^{(i)}，z^{(i)}|\theta )$

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

3.EM算法步骤

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

输出：模型参数θ。

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

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

$$\begin{gathered} Q(\theta ,{\theta }^{(i)})=E_z[logP(Y,Z|\theta )|Y,{\theta }^{(i)})] \\ =\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)}) \end{gathered}$$

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

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

${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i+1)})$

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

式$$\begin{gathered} Q(\theta ,{\theta }^{(i)})=E_z[logP(Y,Z|\theta )|Y,{\theta }^{(i)})] \\ =\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)}) \end{gathered}$$的函数$Q(\theta ,{\theta }^{(i+1)})$是EM算法的核心，称为Q函数（Q function)。

4.EM算法中的Q函数

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

$Q(\theta ,{\theta }^{(i)})=E_z[logP(Y,Z|\theta )|Y,{\theta }^{(i)})]$

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

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

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

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

步骤(4)给出停止迭代的条件，一般是对较小的正数${\epsilon }_1,{\epsilon }_2$,若满足

$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{\epsilon }_1或||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{\epsilon }_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   : hanzi5@yeah.net
 @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、《统计学习方法》李航