EM算法理论与Python代码实践

本次学习主要参考李航老师《统计学习方法》第九章EM算法及其推广和俞栋老师《语音识别实践》第一部分传统声学模型中隐马尔可夫模型及其变体

EM算法的引入
如果概率模型的变量都是观测变量，那么给定数据就可以直接用极大似然估计法或者贝叶斯估计法估计模型参数，但是当模型含有隐变量的时候就不能简单的使用这些估计方法。EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。 EM算法的每次迭代由两部组成：E步，求期望；M步，求极大。

EM算法
输入：观测变量数据Y，隐变量数据Z，联合分布$P=(Y,Z|\theta )$ ,条件分布$P=(Z|Y,\theta )$ ;
输出：模型参数$\theta$
（1）选择参数的初值${\theta }^{(0)}$ ，开始迭代；
（2）E步：记${\theta }^{(i)}$ 为第i次迭代参数$\theta$ 的估计值，在第i+1次迭代的E步，计算$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$
这里$P(Y,Z|{\theta }^{(i)})$ 是在给定观测数据Y和当前的参数估计[\theta ^{(i)}]下隐变量数据Z的条件概率分布；
（3）M步：求使$Q(\theta ,{\theta }^{(i)})$ 极大化的$\theta$ ，确定第i+1次迭代的参数估计值${\theta }^{(i+1)}$
${\theta }^{(i+1)}=argmaxQ(\theta ,{\theta }^{(i)})$
（4）重复（2）（3），直到收敛
其中： $Q(\theta ,{\theta }^{(i)})$ 是EM算法的核心，成为Q函数
给出停止迭代的条件，一般是对比较小的正数${\xi }_1,{\xi }_2$
若满足$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{\xi }_1$ 或$||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{\xi }_2$ 则停止迭代

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)}]$
注意，第一个变元表示要极大化的参数，第2个变元表示参数当前估计值，每次迭代实际都在求Q函数及其极大

EM.py

import numpy as np
import random
import math
import time


# 初始化数据集 通过高斯分布的随机函数来伪造数据集
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:混合了两个高斯分布的数据
    """
    #定义数据集长度为2000
    length = 20
    #初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，一次来满足alpha的作用
    data0 = np.random.normal(mu0,sigma0,int(length * alpha0)) # np.random.normal()的意思是一个正态分布
    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):
    """
     :param dataSetArr: 可观测数据集
    :param mu: 均值
    :param sigmod: 方差
    :return: 整个可观测数据集的高斯分布密度（向量形式）
    """
    result = (1/(math.sqrt(2*math.pi)*sigmod**2))*np.exp(-1*(dataSetArr - mu)**2 / (2*sigmod**2))
    return result

def E_step(dataSetArr,alpha0,mu0,sigmod0,alpha1,mu1,sigmod1):
    """
    根据当前模型参数，计算模型k对观测数据y的响应度
    :param dataSetArr: 可观测数据y
    :param alpha0: 高斯模型0的系数
    :param mu0: 高斯模型0的均值
    :param sigmod0: 高斯模型0的方差
    :param alpha1: 高斯模型1的系数
    :param mu1: 高斯模型1的均值
    :param sigmod1: 高斯模型1的方差
    :return: 两个模型各自的响应度
    """
    gamma0 = alpha0 * calcGauss(dataSetArr,mu0,sigmod0)
    gamma1 = alpha1 * calcGauss(dataSetArr,mu1,sigmod1)


    gamma0 = gamma0 /(gamma0 + gamma1)
    gamma1 = gamma1 / (gamma0 + gamma1)

    return gamma0,gamma1

def M_step(mu0,mu1,gamma0,gamma1,dataSetArr):
    mu0_new = np.dot(gamma0,dataSetArr)
    mu1_new = np.dot(gamma1,dataSetArr)

    sigmod0_new = math.sqrt(np.dot(gamma0,(dataSetArr - mu0)**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算法进行参数估计
    :param dataSetList: 可观测数据的数据集
    :param iter: 迭代次数
    :return: 估计的参数
    """
    dataSetArr = np.array(dataSetList)
    #步骤1 对参数选取初值，开始迭代
    alpha0 = 0.5
    mu0 = 0
    sigmod0 =1
    alpha1 = 0.5
    mu1 = 1
    sigmod1 = 1
    #开始迭代
    step = 0
    while(step<iter):
        step+=1
        gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
        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()

    #设置两个高斯模型进行混合，初始化两个模型各自的参数
    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(dataSetList)

    # 打印设置的参数
    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)