文章目录

  • 0、前言
  • 1、EM算法引入
  • 2、具体的EM算法
  • 3、EM算法推导
  • 3.1 Jensen不等式
  • 3.2 EM推导
  • 3.3 EM算法的收敛性
  • 4、EM算法在高斯混合模型中的应用
  • 4.1 高斯混合模型
  • 4.2 混合高斯分布模型python实现
  • 5 EM算法在HMM模型中的应用
  • 5.1 HMM模型
  • 5.2 EM算法在HMM中应用
  • 参考


0、前言

EM算法,也叫最大期望算法,或者是期望最大化算法,是机器学习十大算法之一,它很简单,但是也同样很有深度。

  • 简单是因为它就分两步求解问题:
  • E步:求期望(expectation)
  • M步:求极大(maximization)
  • 深度在于它的数学推理涉及到比较繁杂的概率公式等。

所以本文会介绍很多概率方面的知识,不懂的同学可以先去了解一些知识,当然本文也会尽可能的讲解清楚这些知识,讲的不好的地方麻烦大家评论指出,后续不断改进完善。

1、EM算法引入

当模型含有隐变量时,就不能简单的使用这些方法,EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法,我们讨论极大似然估计,极大后验概率估计与其类似。

参考统计学习方法书中的一个例子来引入EM算法:

python EMD分解代码 python em算法_概率模型python EMD分解代码 python em算法_算法_02python EMD分解代码 python em算法_算法_03,进行如下实验:

  • 先掷硬币A,根据硬币A的结果选出硬币B和硬币C,正面选硬币B,反面选硬币C;
  • 通过选择出的硬币(B或者C),记录掷硬币的结果出现正面为1,反面为0。

  也就是说一共需要掷两次硬币,第一次为硬币A,第二次硬币需要根据A的结果来选择,只记录一次结果,就是第二次掷的硬币的正反面。把这整个过程看作一次实验,如此独立地重复n次实验。我们当前规定n=10,则10次的结果如下所示:

python EMD分解代码 python em算法_python EMD分解代码_04

  假设现在我们只知道到记录掷硬币的结果,也就是上面的10次结果,不知道每次实验掷硬币的过程,也就是不知道每次实验第一次掷硬币A的结果;现在需要估计每次实验结果出现正面的概率,怎么求呢?

我们来构建这样一个三硬币模型:

硬币A的结果

第二次的硬币

第二次的结果

概率

正面

硬币B

正面

python EMD分解代码 python em算法_算法_05

正面

硬币B

反面

python EMD分解代码 python em算法_机器学习_06

反面

硬币C

正面

python EMD分解代码 python em算法_python EMD分解代码_07

反面

硬币C

反面

python EMD分解代码 python em算法_机器学习_08

python EMD分解代码 python em算法_python EMD分解代码_09,则第二次硬币为B的概率为python EMD分解代码 python em算法_概率模型,则B的结果概率为python EMD分解代码 python em算法_算法_11,当出现正面,就令python EMD分解代码 python em算法_python EMD分解代码_12,概率就为python EMD分解代码 python em算法_python EMD分解代码_13。对C同理。计算模型的概率python EMD分解代码 python em算法_python_14,其中python EMD分解代码 python em算法_python EMD分解代码_15是模型的参数,python EMD分解代码 python em算法_python EMD分解代码_09受到硬币A的结果的影响,记硬币A的结果为python EMD分解代码 python em算法_概率模型_17,则python EMD分解代码 python em算法_python_14是联合概率python EMD分解代码 python em算法_概率模型_19的边际概率,即python EMD分解代码 python em算法_python EMD分解代码_20,所以:

python EMD分解代码 python em算法_python EMD分解代码_21

  • python EMD分解代码 python em算法_算法_22,表示这此看到的是正面,这个正面有可能是B的正面,也可能是C的正面,则python EMD分解代码 python em算法_概率模型_23
  • python EMD分解代码 python em算法_python_24,则python EMD分解代码 python em算法_python_25

python EMD分解代码 python em算法_python EMD分解代码_09是观测变量,表示一次观测结果是python EMD分解代码 python em算法_算法_27python EMD分解代码 python em算法_python EMD分解代码_28python EMD分解代码 python em算法_概率模型_17是隐藏变量,表示掷硬币A的结果,这个是观测不到结果的,python EMD分解代码 python em算法_python EMD分解代码_15表示模型参数。若考虑python EMD分解代码 python em算法_概率模型_31次实验,将观测数据表示为python EMD分解代码 python em算法_算法_32,未观测的数据表示为python EMD分解代码 python em算法_机器学习_33,则观测函数的似然函数是:

python EMD分解代码 python em算法_机器学习_34

考虑求模型参数python EMD分解代码 python em算法_python EMD分解代码_15的极大似然估计,即:
python EMD分解代码 python em算法_机器学习_36
  这个问题没有解析解,只有通过迭代方法来求解,EM算法就是可以用于求解这个问题的一种迭代算法,下面给出EM算法的迭代过程:

  • 首先选取初始值,记做python EMD分解代码 python em算法_概率模型_37,第python EMD分解代码 python em算法_机器学习_38次的迭代参数的估计值为python EMD分解代码 python em算法_python_39
  • E步:计算在模型参数python EMD分解代码 python em算法_机器学习_40下观测变量python EMD分解代码 python em算法_算法_41来源于硬币B的概率:
    python EMD分解代码 python em算法_机器学习_42
    备注一下:这个公式的分母是python EMD分解代码 python em算法_算法_43,分子表示结果来源于B硬币的概率python EMD分解代码 python em算法_python EMD分解代码_44
  • M步:计算模型参数的新估计值:

python EMD分解代码 python em算法_概率模型_45

因为第二次使用B硬币是基于A硬币出现正面的结果,所以A硬币出现正面的概率就是python EMD分解代码 python em算法_算法_46的平均值。

python EMD分解代码 python em算法_算法_47

分子乘以python EMD分解代码 python em算法_概率模型_48,所以其实是计算B硬币出现正面的概率。

python EMD分解代码 python em算法_python_49

其中,python EMD分解代码 python em算法_算法_50表示出现C硬币的概率,即A出现反面的概率。

python EMD分解代码 python em算法_算法_51python EMD分解代码 python em算法_python_52一个闭环流程,接下来可以通过迭代法来做完成。针对上述例子,我们假设初始值为python EMD分解代码 python em算法_概率模型_53,因为对python EMD分解代码 python em算法_算法_54python EMD分解代码 python em算法_概率模型_55均有python EMD分解代码 python em算法_概率模型_56,利用迭代公式计算得到python EMD分解代码 python em算法_python_57,继续迭代得到最终的参数:

python EMD分解代码 python em算法_算法_58

如果一开始初始值选择为:python EMD分解代码 python em算法_机器学习_59,那么得到的模型参数的极大似然估计是:

python EMD分解代码 python em算法_python_60

这说明EM算法与初值的选择有关,选择不同的初值可能得到不同的参数估计值

  这个例子中你只观察到了硬币抛完的结果,并不了解A硬币抛完之后,是选择了B硬币抛还是C硬币抛,这时候概率模型就存在着隐含变量

2、具体的EM算法

  • 输入:观测变量数据Y,隐变量数据Z,联合分布python EMD分解代码 python em算法_python EMD分解代码_44,条件分布python EMD分解代码 python em算法_机器学习_62
  • 输出:模型参数python EMD分解代码 python em算法_算法_63
  • (1) 选择参数的初值python EMD分解代码 python em算法_算法_64,开始迭代
  • (2) E步:记python EMD分解代码 python em算法_算法_65为第python EMD分解代码 python em算法_机器学习_38次迭代参数python EMD分解代码 python em算法_算法_63的估计值,在第python EMD分解代码 python em算法_python_68次迭代的E步,计算python EMD分解代码 python em算法_机器学习_69函数:

python EMD分解代码 python em算法_python_70

这里,python EMD分解代码 python em算法_python EMD分解代码_71是在给定观测数据python EMD分解代码 python em算法_机器学习_72和当前的参数估计python EMD分解代码 python em算法_python EMD分解代码_73下隐变量数据python EMD分解代码 python em算法_概率模型_74的条件概率分布;

  • (3) M步:求使python EMD分解代码 python em算法_python EMD分解代码_75极大化的python EMD分解代码 python em算法_机器学习_76,确定第i+1次迭代的参数的估计值python EMD分解代码 python em算法_python EMD分解代码_77

python EMD分解代码 python em算法_机器学习_78

其中,python EMD分解代码 python em算法_python_79是EM算法的核心,称为Q函数(Q function),这个是需要自己构造的,可以看出这是一个期望,而M步是极大化这个期望值,这也是EM(Expectation Maximization )算法名字的来源。

  • (4) 重复第(2)步和第(3)步,直到收敛,收敛条件:

python EMD分解代码 python em算法_算法_80
收敛迭代就结束了。我们来拆解一下这个M步骤,

3、EM算法推导

3.1 Jensen不等式

  Jensen不等式,这个公式在算法推导和证明收敛都能用到,Jensen不等式主要是如下的结论:

  • python EMD分解代码 python em算法_python EMD分解代码_81是凸函数

python EMD分解代码 python em算法_算法_82

  • python EMD分解代码 python em算法_python EMD分解代码_81是凹函数

python EMD分解代码 python em算法_概率模型_84

3.2 EM推导

python EMD分解代码 python em算法_算法_51最大,其中Y是观测数据,python EMD分解代码 python em算法_算法_86是模型参数,意思就是在参数python EMD分解代码 python em算法_算法_86的情况下,观测数据Y出现的概率最大。然而当我们面对一个含有隐变量的概率模型时,我们也想使在参数python EMD分解代码 python em算法_算法_86的情况下,观测数据Y出现的概率最大,即最大化python EMD分解代码 python em算法_算法_51。但是我们不能直接使用极大似然估计,因为参数python EMD分解代码 python em算法_算法_86包含有隐变量的信息,而Y与隐变量相关,估计Y需要知道Z,估计Z需要知道Y,形成一个死循环。可以将隐变量分离出来,对于每一个样本,我们需要确定Z来使联合概率python EMD分解代码 python em算法_python_91最大,由条件概率公式:python EMD分解代码 python em算法_python EMD分解代码_92

python EMD分解代码 python em算法_python_93

python EMD分解代码 python em算法_算法_51出现的概率最大,为了方便计算,对似然函数取对数。所以我们的目标是极大化观测数据Y关于参数python EMD分解代码 python em算法_算法_86的对数似然函数python EMD分解代码 python em算法_机器学习_96:

python EMD分解代码 python em算法_python EMD分解代码_97

python EMD分解代码 python em算法_机器学习_98,刚好可以利用Jensen不等式。

python EMD分解代码 python em算法_机器学习_96每一次更新python EMD分解代码 python em算法_机器学习_76时,都能使python EMD分解代码 python em算法_python_101增大,也就是假设当前第python EMD分解代码 python em算法_python_102次迭代的参数估计值为python EMD分解代码 python em算法_python EMD分解代码_73,则在这python EMD分解代码 python em算法_概率模型_104次迭代中,我们需要重新估计参数值,使python EMD分解代码 python em算法_机器学习_96增大,即python EMD分解代码 python em算法_python_106。则python EMD分解代码 python em算法_机器学习_107

python EMD分解代码 python em算法_python_108

这一个步骤就是采用了凸函数的Jensen不等式做转换。因为python EMD分解代码 python em算法_概率模型_74是隐藏变量,所以有:
python EMD分解代码 python em算法_机器学习_110

于是继续变:

python EMD分解代码 python em算法_机器学习_111
也就是:python EMD分解代码 python em算法_算法_112,有下界(不等式右边就是下界),任何使下界增大的python EMD分解代码 python em算法_算法_86都可以使python EMD分解代码 python em算法_机器学习_96增大,因此我们需要最大化下界,来得到近似的极大值。

python EMD分解代码 python em算法_概率模型_104次迭代的参数估计值python EMD分解代码 python em算法_python EMD分解代码_116,就是使下界最大化的python EMD分解代码 python em算法_算法_86

python EMD分解代码 python em算法_算法_118

其中python EMD分解代码 python em算法_算法_119分母提出来是关于python EMD分解代码 python em算法_概率模型_74python EMD分解代码 python em算法_机器学习_121,是一个常数可以去掉。其中python EMD分解代码 python em算法_python_122是EM算法中非常重要的Q函数,很多时候写成期望的形式:

python EMD分解代码 python em算法_机器学习_123

python EMD分解代码 python em算法_概率模型_124,也就是python EMD分解代码 python em算法_机器学习_96的下界,可以看出python EMD分解代码 python em算法_概率模型_126,当python EMD分解代码 python em算法_算法_127时,等号成立。函数python EMD分解代码 python em算法_机器学习_128的增加会导致python EMD分解代码 python em算法_机器学习_96在每次迭代都是增加的,我们想尽快迭代到最优值,就要求在每一次迭代的时候都取到函数python EMD分解代码 python em算法_机器学习_128的最大值,也就是python EMD分解代码 python em算法_python EMD分解代码_116,也就是Q函数的最大值。然后基于python EMD分解代码 python em算法_python EMD分解代码_116,重新计算python EMD分解代码 python em算法_python EMD分解代码_133,会得到类似的图,继续寻找最大值python EMD分解代码 python em算法_概率模型_134,以此类推下去。在这个过程中,对数似然函数python EMD分解代码 python em算法_机器学习_96不断增大,直到达到全局最优值。

python EMD分解代码 python em算法_python_136

3.3 EM算法的收敛性

  EM算法的最大优点是简单性和普适性,但是EM算法得到的估计序列是否收敛呢?是收敛到全局最大值还是局部最大值?

定理一:设python EMD分解代码 python em算法_算法_51是观测数据的似然函数,python EMD分解代码 python em算法_python EMD分解代码_138是EM算法得到的参数估计序列,python EMD分解代码 python em算法_机器学习_139是对应的似然函数序列,则python EMD分解代码 python em算法_概率模型_140单调递增的,即:

python EMD分解代码 python em算法_算法_141

  证明过程:
由条件概率python EMD分解代码 python em算法_概率模型_142可得:

python EMD分解代码 python em算法_机器学习_143

取对数有:

python EMD分解代码 python em算法_算法_144

由上文可知:

python EMD分解代码 python em算法_机器学习_145

令:

python EMD分解代码 python em算法_算法_146

python EMD分解代码 python em算法_算法_147

则:

python EMD分解代码 python em算法_python_148

要证明python EMD分解代码 python em算法_python_149,只需证明上式的右边是非负的。其中python EMD分解代码 python em算法_python EMD分解代码_150是必大于0的,因为python EMD分解代码 python em算法_python EMD分解代码_116python EMD分解代码 python em算法_机器学习_152的最大值。则:

python EMD分解代码 python em算法_算法_153

由此可以证明python EMD分解代码 python em算法_python_149

定理二:设python EMD分解代码 python em算法_机器学习_155为观测数据的对数似然函数,python EMD分解代码 python em算法_python EMD分解代码_138是EM算法得到的参数估计序列,python EMD分解代码 python em算法_概率模型_157是对应的对数似然函数序列。

  1. 如果python EMD分解代码 python em算法_机器学习_158有上界,则$python EMD分解代码 python em算法_算法_159收敛到某一值python EMD分解代码 python em算法_概率模型_160
  2. 在函数python EMD分解代码 python em算法_算法_161python EMD分解代码 python em算法_python_101满足一定条件下,由EM算法得到的参数估计序列python EMD分解代码 python em算法_概率模型_163的收敛值python EMD分解代码 python em算法_算法_164python EMD分解代码 python em算法_python_101的稳定点。

这里就不证明了

4、EM算法在高斯混合模型中的应用

4.1 高斯混合模型

  EM算法的一个重要应用场景就是高斯混合模型的参数估计。高斯混合模型就是由多个高斯模型组合在一起的混合模型(可以理解为多个高斯分布函数的线性组合,理论上高斯混合模型是可以拟合任意类型的分布),例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的:

python EMD分解代码 python em算法_python EMD分解代码_166

两个高斯模型可以拟合数据集,如图所示:

python EMD分解代码 python em算法_机器学习_167

如果有多个高斯模型,公式表示为:

python EMD分解代码 python em算法_算法_168

python EMD分解代码 python em算法_算法_169表示为第k个高斯分布密度模型,定义如上,其中python EMD分解代码 python em算法_python EMD分解代码_170表示被选中的概率。在本次模型python EMD分解代码 python em算法_python_14中,观测数据是已知的,而观测数据具体来自哪个模型是未知的,有点像之前提过的三硬币模型,我们来对比一下,A硬币就像是概率python EMD分解代码 python em算法_python EMD分解代码_170,用来表明具体的模型,而B、C硬币就是具体的模型,只不过这里有很多个模型,不仅仅是B、C这两个模型。我们用python EMD分解代码 python em算法_python_173来表示,则:

python EMD分解代码 python em算法_算法_174

所以一个观测数据python EMD分解代码 python em算法_机器学习_175的隐藏数据python EMD分解代码 python em算法_算法_176,那么完全似然函数就是:

python EMD分解代码 python em算法_python EMD分解代码_177

取对数之后等于:

python EMD分解代码 python em算法_python_178

  • E 步 :
    python EMD分解代码 python em算法_python_179
    其中我们定义python EMD分解代码 python em算法_python EMD分解代码_180
    python EMD分解代码 python em算法_算法_181
    于是化简得到:
    python EMD分解代码 python em算法_概率模型_182

E 步 在代码设计上只有python EMD分解代码 python em算法_算法_183有用,用于M步的计算。

  • M步,
    python EMD分解代码 python em算法_机器学习_184
    python EMD分解代码 python em算法_python EMD分解代码_75求导,得到每个未知量的偏导,使其偏导等于0,求解得到:
    python EMD分解代码 python em算法_概率模型_186
    给一个初始值,来回迭代就可以求得值内容。这一块主要用到了python EMD分解代码 python em算法_python EMD分解代码_187的导数,并且用到了E步的python EMD分解代码 python em算法_python EMD分解代码_180

4.2 混合高斯分布模型python实现

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

这其中的M步 完全按照了 公式来计算。

import numpy as np
import random
import math
import time
'''
数据集:伪造数据集(两个高斯分布混合)
数据集长度:1000
------------------------------
运行结果:
----------------------------
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
----------------------------
'''

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)

E步主要计算内容:
其中我们定义python EMD分解代码 python em算法_算法_183
python EMD分解代码 python em算法_概率模型_190

M步 主要计算内容
这一步骤主要是对Q函数求导后的数据进行计算,利用了 E 步 的python EMD分解代码 python em算法_算法_183

import math
import copy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
#生成随机数据,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.matrix(X)
    global mu                 #随机初始化mu1,mu2,mu3,mu4
    mu = np.random.random((4,2))
    mu=np.matrix(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.matrix([[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()

python EMD分解代码 python em算法_概率模型_192

5 EM算法在HMM模型中的应用

5.1 HMM模型

  HMM模型也是一个含有隐变量的模型,具体参考HMM模型详解

5.2 EM算法在HMM中应用

  HMM模型在学习参数的时候,可以使用极大似然估计,也可以使用EM算法(在HMM称为bamu-welch算法)

class HMM(object):
    def __init__(self, n, m, a=None, b=None, pi=None):
        # 可能的隐藏状态数
        self.N = n
        # 可能的观测数
        self.M = m
        # 状态转移概率矩阵
        self.A = a
        # 观测概率矩阵
        self.B = b
        # 初始状态概率矩阵
        self.Pi = pi
        # 观测序列
        self.X = None
        # 状态序列
        self.Y = None
        # 序列长度
        self.T = 0

        # 定义前向算法
        self.alpha = None
        # 定义后向算法
        self.beta = None
    def baum_welch(self, x, criterion=0.001):
        self.T = len(x)
        self.X = x

        while True:
            # 为了得到alpha和beta的矩阵
            _ = self.forward(self.X)
            _ = self.backward(self.X)
            xi = np.zeros((self.T - 1, self.N, self.N), dtype=float)
            for t in range(self.T - 1):
            # 笨办法
            # for i in range(self.N):
            # gamma[t][i] = self.calculate_gamma(t, i)
            #         for j in range(self.N):
            #             xi[t][i][j] = self.calculate_psi(t, i, j)
            # for i in range(self.N):
            #     gamma[self.T-1][i] = self.calculate_gamma(self.T-1, i)

            # numpy矩阵的办法
                denominator = np.sum(np.dot(self.alpha[t, :], self.A) *
                                     self.B[:, self.X[t + 1]] * self.beta[t + 1, :])
                for i in range(self.N):
                    molecular = self.alpha[t, i] * self.A[i, :] * self.B[:, self.X[t+1]]*self.beta[t+1, :]
                    xi[t, i, :] = molecular / denominator
            gamma = np.sum(xi, axis=2)
            prod = (self.alpha[self.T-1, :]*self.beta[self.T-1, :])
            gamma = np.vstack((gamma, prod / np.sum(prod)))

            new_pi = gamma[0, :]
            new_a = np.sum(xi, axis=0) / np.sum(gamma[:-1, :], axis=0).reshape(-1, 1)
            new_b = np.zeros(self.B.shape, dtype=float)

            for k in range(self.B.shape[1]):
                mask = self.X == k
                new_b[:, k] = np.sum(gamma[mask, :], axis=0) / np.sum(gamma, axis=0)

            if np.max(abs(self.Pi - new_pi)) < criterion and \
                    np.max(abs(self.A - new_a)) < criterion and \
                    np.max(abs(self.B - new_b)) < criterion:
                break

        self.A, self.B, self.Pi = new_a, new_b, new_pi

参考

主要参考统计学习方法

统计学习方法-代码解读

EM算法 - 期望极大算法