机器学习之EM算法 学习笔记-Datawhale Task03

Author	Bryce230
e-mail	2540892461@qq.com
Software	win10，Pycharm2019.3.3，Python3.7.7

1 前言

概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。

EM算法是机器学习十大算法之一，主要分为E（期望）步和M（极大值）步。涉及到的数学知识有深度，推导比较繁琐。本文主要参考“如何通俗理解EM算法”一文以及自己的理解。如果文中有些图片显示不了，可以点这个：七月在线-如何通俗理解EM算法

2 EM算法的通俗理解

2.1 抛硬币

Nature Biotech在他的文章《Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature biotechnology, 26(8), 897.》中用了一个抛硬币的例子来讲EM算法的思想，如下图。
a 表示知道抛的是硬币A还是B，b 表示不知道抛的是哪枚硬币。

看不懂？我也看不懂，那就看看v_JULY_v大佬的解释。

有两枚硬币A和B，假定随机抛掷后正面朝上概率分别为PA，PB。为了估计这两个硬币朝上的概率，我们轮流抛硬币A和B，每一轮都连续抛5次，总共5轮：

硬币	结果	统计
A	正正反正反	3正-2反
B	反反正正反	2正-3反
A	正反反反反	1正-4反
B	正反反正正	3正-2反
A	反正正反反	2正-3反

硬币A被抛了15次，在第一轮、第三轮、第五轮分别出现了3次正、1次正、2次正，所以很容易估计出PA：
PA = （3+1+2）/ 15 = 0.4
类似的，PB也很容易计算出来：
PB= （2+3）/10 = 0.5

这好像有点简单，那么问题来了？如果我们不知道抛的硬币是A还是B呢（即硬币种类是隐变量），然后再轮流抛五轮，得到如下结果：

硬币	结果	统计
Unknown	正正反正反	3正-2反
Unknown	反反正正反	2正-3反
Unknown	正反反反反	1正-4反
Unknown	正反反正正	3正-2反
Unknown	反正正反反	2正-3反

现在我们的目标没变，还是估计PA和PB，需要怎么做呢？
显然，此时我们多了一个硬币种类的隐变量，设为z，可以把它认为是一个5维的向量（z1, z2, z3, z4, z5)，代表每次投掷时所使用的硬币，比如z1，就代表第一轮投掷时使用的硬币是A还是B。

但是，这个变量z不知道，就无法去估计PA和PB，所以，我们必须先估计出z，然后才能进一步估计PA和PB。
可要估计z，我们又得知道PA和PB，这样我们才能用极大似然概率法则去估计z，这不是鸡生蛋和蛋生鸡的问题吗，如何破？

答案就是我们初始化的PA和PB，按照最大似然概率就可以估计出z，然后基于z，按照最大似然概率可以反过来估计出P1和P2，当与我们初始化的PA和PB一样时，说明是P1和P2很有可能就是真实的值。如果新估计出来的PA和PB和我们初始化的值差别很大，怎么办呢？就是继续用新的P1和P2迭代，直至收敛。

先随便给PA和PB赋一个值，比如：

硬币A正面朝上的概率PA = 0.2
硬币B正面朝上的概率PB = 0.7

然后，我们看看第一轮抛掷最可能是哪个硬币。
如果是硬币A，得出3正2反的概率为 0.2 * 0.2 * 0.2 * 0.8 * 0.8 = 0.00512
如果是硬币B，得出3正2反的概率为 0.7 * 0.7 * 0.7 * 0.3 * 0.3 = 0.03087
然后依次求出其他4轮中的相应概率。做成表格如下（标粗表示其概率更大）：

轮数	若是硬币A	若是硬币B
1	0.00512，即0.2 0.2 0.2 0.8 0.8，3正-2反	0.03087，3正-2反
2	0.02048，即0.2 0.2 0.8 0.8 0.8，2正-3反	0.01323，2正-3反
3	0.08192，即0.2 0.8 0.8 0.8 0.8，1正-4反	0.00567，1正-4反
4	0.00512，即0.2 0.2 0.2 0.8 0.8，3正-2反	0.03087，3正-2反
5	0.02048，即0.2 0.2 0.8 0.8 0.8，2正-3反	0.01323，2正-3反

按照最大似然法则：
第1轮中最有可能的是硬币B
第2轮中最有可能的是硬币A
第3轮中最有可能的是硬币A
第4轮中最有可能的是硬币B
第5轮中最有可能的是硬币A

我们便可以按照最大似然概率法则来估计新的PA和PB：
PA = （2+1+2）/15 = 0.33
PB =（3+3）/10 = 0.6

设想我们是全知的神，知道每轮抛掷时的硬币就是如本文本节开头标示的那样，那么，PA和PB的最大似然估计就是0.4和0.5（下文中将这两个值称为PA和PB的真实值）。那么对比下我们初始化的PA和PB和新估计出的PA和PB：

初始化的PA	估计出的PA	真实的PA	初始化的PB	估计出的PB	真实的PB
0.2	0.33	0.4	0.7	0.6	0.5

我们估计的PA和PB相比于它们的初始值，更接近它们的真实值了！就这样，不断迭代，不断接近真实值，这就是EM算法的奇妙之处。

可以期待，我们继续按照上面的思路，用估计出的PA和PB再来估计z，再用z来估计新的PA和PB，反复迭代下去，就可以最终得到PA = 0.4，PB=0.5，此时无论怎样迭代，PA和PB的值都会保持0.4和0.5不变，于是乎，我们就找到了PA和PB的最大似然估计。

2.2 盛菜

你现在需要把锅里的一份菜平均分到两个盘子，该怎么分？
平均分配这个结果是“观测数据z”，为实现平均分配而给每个盘子分配多少菜是“待求参数θ”，分配菜的手感就是“概率分布”。

EM算法的思想就是：
1）给θ自主规定个初值（既然我不知道想实现“两个碟子平均分配锅里的菜”的话每个碟子需要有多少菜，那我就先估计个值）；
2） 根据给定观测数据和当前的参数θ，求未观测数据z的条件概率分布的期望（在上一步中，已经根据手感将菜倒进了两个碟子，然后这一步根据“两个碟子里都有菜”和“当前两个碟子都有多少菜”来判断自己倒菜的手感）；
3）上一步中z已经求出来了，于是根据极大似然估计求最优的θ’（手感已经有了，那就根据手感判断下盘子里应该有多少菜，然后把菜匀匀）；
4） 因为第二步和第三步的结果可能不是最优的，所以重复第二步和第三步，直到收敛（重复多次匀匀的过程，直到两个碟子中菜的量大致一样）。

上面的第二步被称作E步（求期望），第三步被称作M步（求极大化），从而不断的重复E、M，直至收敛，也即两份菜一样多。

3 EM算法

3.1 EM算法的流程

推导逼近和证明收敛这里不做展示，具体可参考Datawhale: Task3 EM.ipynb

EM算法的整个过程可见下图，我们固定θ，调整Q(z)使下界J(z,Q)上升至与L(θ)在此点θ处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（θt到θt+1），然后再固定θ，调整Q(z)……直到收敛到似然函数L(θ)的最大值处的θ*。

3.2 EM算法的另一种理解

如果我们定义

EM可以看做是J的坐标上升法（Coordinate ascent），E步固定$\theta$，优化Q，M步固定Q，优化$\theta$。

图中的直线式迭代优化的路径，可以看到每一步都会向最优值前进一步，而且前进路线是平行于坐标轴的，因为每一步只优化一个变量。

对应到EM上，就是：E步：固定θ，优化Q；M步：固定Q，优化θ；交替将极值推向最大。

3.3 EM算法的应用

EM算法广泛应用于GMM混合高斯模型、聚类、HMM等等。比如聚类

4 高斯混合分布

高斯混合模型，顾名思义，就是由多个高斯模型组合在一起的混合模型（可以理解为多个高斯分布函数的线性组合，理论上高斯混合模型是可以拟合任意类型的分布），例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的：
两个高斯模型可以拟合数据集，如图所示：

5 代码实现

5.1 混合高斯分布模型

EM算法更多是一种思想，用概率来解决问题的一种方法，具体的代码看自己选用模型，所以并没有通用的模型，本此代码主要是讲解混合高斯分布模型的求解过程。

```python
import numpy as np
import random
import math
import time

'''
数据集：伪造数据集（两个高斯分布混合）
数据集长度：1000
'''

def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
    '''
    初始化数据集
    这里通过服从高斯分布的随机函数来伪造数据集
    :param mu0: 高斯0的均值
    :param sigma0: 高斯0的方差
    :param mu1: 高斯1的均值
    :param sigma1: 高斯1的方差
    :param alpha0: 高斯0的系数
    :param alpha1: 高斯1的系数
    :return: 混合了两个高斯分布的数据
    '''
    # 定义数据集长度为1000
    length = 1000

    # 初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，以此来满足alpha的作用
    data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
    # 第二个高斯分布的数据
    data1 = np.random.normal(mu1, sigma1, int(length * alpha1))

    # 初始化总数据集
    # 两个高斯分布的数据混合后会放在该数据集中返回
    dataSet = []
    # 将第一个数据集的内容添加进去
    dataSet.extend(data0)
    # 添加第二个数据集的数据
    dataSet.extend(data1)
    # 对总的数据集进行打乱（其实不打乱也没事，只不过打乱一下直观上让人感觉已经混合了
    # 读者可以将下面这句话屏蔽以后看看效果是否有差别）
    random.shuffle(dataSet)

    #返回伪造好的数据集
    return dataSet

# 高斯分布公式，没有什么特殊的
def calcGauss(dataSetArr, mu, sigmod):
    '''
    根据高斯密度函数计算值
    依据：“9.3.1 高斯混合模型” 式9.25
    注：在公式中y是一个实数，但是在EM算法中(见算法9.2的E步)，需要对每个j
    都求一次yjk，在本实例中有1000个可观测数据，因此需要计算1000次。考虑到
    在E步时进行1000次高斯计算，程序上比较不简洁，因此这里的y是向量，在numpy
    的exp中如果exp内部值为向量，则对向量中每个值进行exp，输出仍是向量的形式。
    所以使用向量的形式1次计算即可将所有计算结果得出，程序上较为简洁

    :param dataSetArr: 可观测数据集
    :param mu: 均值
    :param sigmod: 方差
    :return: 整个可观测数据集的高斯分布密度（向量形式）
    '''
    # 计算过程就是依据式9.25写的，没有别的花样
    result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))
    # 返回结果
    return result


def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
    '''
    EM算法中的E步
    依据当前模型参数，计算分模型k对观数据y的响应度
    :param dataSetArr: 可观测数据y
    :param alpha0: 高斯模型0的系数
    :param mu0: 高斯模型0的均值
    :param sigmod0: 高斯模型0的方差
    :param alpha1: 高斯模型1的系数
    :param mu1: 高斯模型1的均值
    :param sigmod1: 高斯模型1的方差
    :return: 两个模型各自的响应度
    '''
    # 计算y0的响应度
    # 先计算模型0的响应度的分子
    gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
    #print("gamma0=",gamma0.shape) # 1000, 维向量
    # 模型1响应度的分子
    gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)

    # 两者相加为E步中的分布
    sum = gamma0 + gamma1
    # 各自相除，得到两个模型的响应度
    gamma0 = gamma0 / sum
    gamma1 = gamma1 / sum

    # 返回两个模型响应度
    return gamma0, gamma1

def M_step(muo, mu1, gamma0, gamma1, dataSetArr):
    # 依据算法9.2计算各个值
    # 这里没什么花样，对照书本公式看看这里就好了

    # np.dot 点积：[1,2] [2,3] = [2,6]
    mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)
    mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)

    # math.sqrt  平方根
    sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))
    sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))

    alpha0_new = np.sum(gamma0) / len(gamma0)
    alpha1_new = np.sum(gamma1) / len(gamma1)

    # 将更新的值返回
    return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new


## 训练主函数
def EM_Train(dataSetList, iter=500):
    '''
    根据EM算法进行参数估计
    算法依据“9.3.2 高斯混合模型参数估计的EM算法” 算法9.2
    :param dataSetList:数据集（可观测数据）
    :param iter: 迭代次数
    :return: 估计的参数
    '''
    # 将可观测数据y转换为数组形式，主要是为了方便后续运算
    dataSetArr = np.array(dataSetList)

    # 步骤1：对参数取初值，开始迭代
    alpha0 = 0.5
    mu0 = 0
    sigmod0 = 1
    alpha1 = 0.5
    mu1 = 1
    sigmod1 = 1

    # 开始迭代
    step = 0
    while (step < iter):
        # 每次进入一次迭代后迭代次数加1
        step += 1
        # 步骤2：E步：依据当前模型参数，计算分模型k对观测数据y的响应度
        gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
        # 步骤3：M步
        mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)

    # 迭代结束后将更新后的各参数返回
    return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1


if __name__ == '__main__':
    start = time.time()

    # 设置两个高斯模型进行混合，这里是初始化两个模型各自的参数
    # 见“9.3 EM算法在高斯混合模型学习中的应用”
    # alpha是“9.3.1 高斯混合模型” 定义9.2中的系数α
    # mu0是均值μ
    # sigmod是方差σ
    # 在设置上两个alpha的和必须为1，其他没有什么具体要求，符合高斯定义就可以

    alpha0 = 0.3  # 系数α
    mu0 = -2  # 均值μ
    sigmod0 = 0.5  # 方差σ

    alpha1 = 0.7  # 系数α
    mu1 = 0.5  # 均值μ
    sigmod1 = 1  # 方差σ

    # 初始化数据集
    dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)

    #打印设置的参数
    print('---------------------------')
    print('the Parameters set is:')
    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
        alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
    ))

    # 开始EM算法，进行参数估计
    alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)

    # 打印参数预测结果
    print('----------------------------')
    print('the Parameters predict is:')
    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
        alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
    ))

    # 打印时间
    print('----------------------------')
    print('time span:', time.time() - start)

运行结果为：

---------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
time span: 0.06682229042053223

5.2 绘图表示高斯混合模型

import math
import copy
import numpy as np
import matplotlib.pyplot as plt

#生成随机数据，4个高斯模型
def generate_data(sigma,N,mu1,mu2,mu3,mu4,alpha):
    global X                  #可观测数据集
    X = np.zeros((N, 2))       # 初始化X，2行N列。2维数据，N个样本
    X=np.mat(X)
    global mu                 #随机初始化mu1，mu2，mu3，mu4
    mu = np.random.random((4,2))
    mu=np.mat(mu)
    global excep              #期望第i个样本属于第j个模型的概率的期望
    excep=np.zeros((N,4))
    global alpha_             #初始化混合项系数
    alpha_=[0.25,0.25,0.25,0.25]
    for i in range(N):
        if np.random.random(1) < 0.1:  # 生成0-1之间随机数
            X[i,:]  = np.random.multivariate_normal(mu1, sigma, 1)     #用第一个高斯模型生成2维数据
        elif 0.1 <= np.random.random(1) < 0.3:
            X[i,:] = np.random.multivariate_normal(mu2, sigma, 1)      #用第二个高斯模型生成2维数据
        elif 0.3 <= np.random.random(1) < 0.6:
            X[i,:] = np.random.multivariate_normal(mu3, sigma, 1)      #用第三个高斯模型生成2维数据
        else:
            X[i,:] = np.random.multivariate_normal(mu4, sigma, 1)      #用第四个高斯模型生成2维数据

    print("可观测数据：\n",X)       #输出可观测样本
    print("初始化的mu1，mu2，mu3，mu4：",mu)      #输出初始化的mu


# E 期望
#  \hat{\gamma_{jk}}
def e_step(sigma,k,N):
    global X
    global mu
    global excep
    global alpha_
    for i in range(N):
        denom=0
        for j in range(0,k):
            #  sigma.I 表示矩阵的逆矩阵
            # np.transpose ：矩阵转置   np.linalg.det():矩阵求行列式
            denom += alpha_[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(0,k):
            numer = math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/np.sqrt(np.linalg.det(sigma))        #分子
            excep[i,j]=alpha_[j]*numer/denom      #求期望
    print("隐藏变量：\n",excep)


def m_step(k,N):
    global excep
    global X
    global alpha_
    for j in range(0,k):
        denom=0   #分母
        numer=0   #分子
        for i in range(N):
            numer += excep[i,j]*X[i,:]
            denom += excep[i,j]
        mu[j,:] = numer/denom    #求均值
        alpha_[j]=denom/N        #求混合项系数

        #     #可视化结果
def plotShow():
    # 画生成的原始数据
    plt.subplot(221)
    plt.scatter(X[:,0].tolist(), X[:,1].tolist(),c='b',s=25,alpha=0.4,marker='o')    #T散点颜色，s散点大小，alpha透明度，marker散点形状
    plt.title('random generated data')
    #画分类好的数据
    plt.subplot(222)
    plt.title('classified data through EM')
    order=np.zeros(N)
    color=['b','r','k','y']
    for i in range(N):
        for j in range(k):
            if excep[i,j]==max(excep[i,:]):
                order[i]=j     #选出X[i,:]属于第几个高斯模型
            probility[i] += alpha_[int(order[i])]*math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/(np.sqrt(np.linalg.det(sigma))*2*np.pi)    #计算混合高斯分布
        plt.scatter(X[i, 0], X[i, 1], c=color[int(order[i])], s=25, alpha=0.4, marker='o')      #绘制分类后的散点图
    #绘制三维图像
    ax = plt.subplot(223, projection='3d')
    plt.title('3d view')
    for i in range(N):
        ax.scatter(X[i, 0], X[i, 1], probility[i], c=color[int(order[i])])
    plt.show()


if __name__ == '__main__':
    iter_num=1000  #迭代次数
    N=500         #样本数目
    k=4            #高斯模型数
    probility = np.zeros(N)    #混合高斯分布
    u1=[5,35]
    u2=[30,40]
    u3=[20,20]
    u4=[45,15]
    sigma=np.mat([[30, 0], [0, 30]])               #协方差矩阵
    alpha=[0.1,0.2,0.3,0.4]         #混合项系数
    generate_data(sigma,N,u1,u2,u3,u4,alpha)     #生成数据
    #迭代计算
    for i in range(iter_num):
        err=0     #均值误差
        err_alpha=0    #混合项系数误差
        Old_mu = copy.deepcopy(mu)
        Old_alpha = copy.deepcopy(alpha_)

        e_step(sigma,k,N)     # E步
        m_step(k,N)           # M步

        print("迭代次数:",i+1)
        print("估计的均值:",mu)
        print("估计的混合项系数:",alpha_)
        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]-alpha_[z])
        if (err<=0.001) and (err_alpha<0.001):     #达到精度退出迭代
            print(err,err_alpha)
            break

    # 画图
    plotShow()

运行后的结果为：

6 参考文献

[1] 如何通俗理解EM算法。如果有些图片显示不了，可以点这个：七月在线-如何通俗理解EM算法
[2] EM算法 - 期望极大算法
[3] Datawhale: Task3 EM.ipynb
[4] 怎么通俗易懂地解释EM算法并且举个例子?