em算法 c语言,EM算法原理与应用（附代码）

1.     定义

最大期望算法(Expectation-maximizationalgorithm，又译为期望最大化算法)，是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐性变量。2.     算法步骤

最大期望算法经过两个步骤交替进行计算：

第一步是计算期望(E)，利用对隐藏变量的现有估计值，计算其最大似然估计值；

第二步是最大化(M)，最大化在E步上求得的最大似然值来计算参数的值。M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行。3.     什么是似然函数？

在数理统计学中，似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。“似然性”与“或然性”或“概率”意思相近，都是指某种事件发生的可能性。

而极大似然就相当于最大可能的意思。

比如你一位同学和一位猎人一起外出打猎，一只野兔从前方窜过。只听一声枪响，野兔应声到下，如果要你推测，这一发命中的子弹是谁打的？你就会想，只发一枪便打中，由于猎人命中的概率一般大于你那位同学命中的概率，从而推断出这一枪应该是猎人射中的。

这个例子所作的推断就体现了最大似然法的基本思想。

多数情况下我们是根据已知条件来推算结果，而最大似然估计是已经知道了结果，然后寻求使该结果出现的可能性最大的条件，以此作为估计值。

概率是根据条件推测结果，而似然则是根据结果反推条件。在这种意义上，似然函数可以理解为条件概率的逆反。

假定已知某个参数B时，推测事件A会发生的概率写作：

通过贝叶斯公式，可以得出

求极大似然函数估计值的一般步骤：

[1]写出似然函数；

[2]对似然函数取对数，并整理；

[3]求导数，令导数为0，得到似然方程；

[4]解似然方程，得到的参数即为所求。

4.  EM算法的经典案例：抛硬币

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

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

PA =(3+1+2)/ 15 = 0.4

PB=(2+3)/10 = 0.5

问题来了，如果我们不知道抛的硬币是A还是B呢(即硬币种类是隐变量)，然后再轮流抛五轮，得到如下结果：

问题变得有意思了。现在我们的目标没变，还是估计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，还是按照最大似然概率法则去估计新的PA和PB，不断迭代，直至收敛。

我们不妨这样，先随便给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轮中的相应概率。做成表格如下(标粗表示其概率更大)：

按照最大似然法则：

第1轮中最有可能的是硬币B，第2轮中最有可能的是硬币A，第3轮中最有可能的是硬币A，第4轮中最有可能的是硬币B，第5轮中最有可能的是硬币A。

我们就把概率更大，即更可能是A的，即第2轮、第3轮、第5轮出现正的次数2、1、2相加，除以A被抛的总次数15(A抛了三轮，每轮5次)，作为z的估计值，B的计算方法类似。然后我们便可以按照最大似然概率法则来估计新的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。

就这样，不断迭代不断接近真实值，这就是EM算法的奇妙之处。

可以期待，我们继续按照上面的思路，用估计出的PA和PB再来估计z，再用z来估计新的PA和PB，反复迭代下去，找到保持不变的值。我们就找到了PA和PB的最大似然估计。

假设我们有一个样本集{x(1),…,x(m)}，包含m个独立的样本，现在我们想找到每个样本隐含的类别z，能使得p(x,z)最大。p(x,z)的极大似然估计如下：

第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。

对于参数估计，我们本质上还是想获得一个使似然函数最大化的那个参数θ，现在与极大似然不同的只是似然函数式中多了一个未知的变量z，见下式(1)。也就是说我们的目标是找到适合的θ和z，以让L(θ)最大。那我们也许会想，你就是多了一个未知的变量而已啊，我也可以分别对未知的θ和z分别求偏导，再令其等于0，求解出来不也一样吗？

本质上我们是需要最大化(1)式，也就是似然函数，我们把分子分母同乘以一个相等的函数(即隐变量Z的概率分布Qi(z(i))，其概率之和等于1，即，即得到上图中的(2)式，我们通过Jensen不等式可得到(3)式，此时(3)式变成了“对数的和”，如此求导就容易了。

当我们在求(2)式的极大值时，(2)式不太好计算，我们便想(2)式大于某个方便计算的下界(3)式，而我们尽可能让下界(3)式最大化，直到(3)式的计算结果等于(2)式，便也就间接求得了(2)式的极大值；

关于Jensen不等式，定义如下：如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])，通俗的说法是函数的期望大于等于期望的函数。特别地，如果f是严格凸函数，当且仅当P(X = EX) = 1，即X是常量时，上式取等号，即E[f(X)] = f(EX)。

5. EM算法应用

关于EM算法的应用其实很多，最广泛的就是GMM混合高斯模型、聚类、HMM等等。如下是一个在Python 3环境使用EM算法求解GMM的编程实现。

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import multivariate_normal

#构建测试数据

N = 200; pi1 = np.array([0.6, 0.3, 0.1])

mu1 = np.array([[0,4], [-2,0], [3,-3]])

sigma1 = np.array([[[3,0],[0,0.5]], [[1,0],[0,2]], [[.5,0],[0,.5]]])

gen = [np.random.multivariate_normal(mu, sigma, int(pi*N)) for mu,sigma, pi in zip(mu1, sigma1, pi1)]

X = np.concatenate(gen)

#初始化: mu, sigma, pi =均值,协方差矩阵,混合系数

theta = {}; param = {}

theta['pi'] = [1/3, 1/3, 1/3]           #均匀初始化

theta['mu'] = np.random.random((3, 2))  #随机初始化

theta['sigma'] = np.array([np.eye(2)]*3) #初始化为单位正定矩阵

param['k'] = len(pi1); param['N'] = X.shape[0]; param['dim'] =X.shape[1]

#定义函数

def GMM_component(X, theta, c):

'''

由联合正态分布计算GMM的单个成员

'''

returntheta['pi'][c]*multivariate_normal(theta['mu'][c], theta['sigma'][c,...]).pdf(X)

def E_step(theta, param):

'''

E步：更新隐变量概率分布q(Z)。

'''

q = np.zeros((param['k'],param['N']))

for i in range(param['k']):

q[i, :] = GMM_component(X,theta, i)

q /= q.sum(axis=0)

return q

def M_step(X, q, theta, param):

'''

M步：使用q(Z)更新GMM参数。

'''

pi_temp = q.sum(axis=1);pi_temp /= param['N'] #计算pi

mu_temp = q.dot(X); mu_temp /=q.sum(axis=1)[:, None] #计算mu

sigma_temp =np.zeros((param['k'], param['dim'], param['dim']))

for i in range(param['k']):

ys = X - mu_temp[i, :]

sigma_temp[i] =np.sum(q[i, :, None, None]*np.matmul(ys[..., None], ys[:, None, :]), axis=0)

sigma_temp /= np.sum(q,axis=1)[:, None, None] #计算sigma

theta['pi'] = pi_temp; theta['mu']= mu_temp; theta['sigma'] = sigma_temp

return theta

def likelihood(X, theta):

'''

计算GMM的对数似然。

'''

ll = 0

for i in range(param['k']):

ll += GMM_component(X,theta, i)

ll = np.log(ll).sum()

return ll

def EM_GMM(X, theta, param, eps=1e-5, max_iter=1000):

'''

高斯混合模型的EM算法求解

theta: GMM模型参数; param:其它系数; eps:计算精度; max_iter:最大迭代次数

返回对数似然和参数theta，theta是包含pi、mu、sigma的Python字典

'''

for i in range(max_iter):

ll_old = 0

# E-step

q = E_step(theta, param)

# M-step

theta = M_step(X, q,theta, param)

ll_new = likelihood(X,theta)

if np.abs(ll_new - ll_old)< eps:

break;

else:

ll_old = ll_new

return ll_new, theta

# EM算法求解GMM，最大迭代次数为1e5

ll, theta2 = EM_GMM(X, theta, param, max_iter=10000)

#由theta计算联合正态分布的概率密度

L = 100; Xlim = [-6, 6]; Ylim = [-6, 6]

XGrid, YGrid = np.meshgrid(np.linspace(Xlim[0], Xlim[1], L),np.linspace(Ylim[0], Ylim[1], L))

Xout = np.vstack([XGrid.ravel(), YGrid.ravel()]).T

MVN = np.zeros(L*L)

for i in range(param['k']):

MVN += GMM_component(Xout,theta, i)

MVN = MVN.reshape((L, L))

#绘制结果

plt.plot(X[:, 0], X[:, 1], 'x', c='gray', zorder=1)

plt.contour(XGrid, YGrid, MVN, 5, colors=('k',), linewidths=(2,))

参考文献：July在其CSDN博客上的EM笔记《如何通俗理解EM算法》。