EM算法及GMM实现

1、引言

E，expectation（期望）；M，maximization（极大化）；
EM算法，又称期望极大算法。

EM已知的是观察数据，未知的是隐含数据和模型参数，在E步，我们所做的事情是固定模型参数的值，优化隐含数据的分布，而在M步，我们所做的事情是固定隐含数据分布，优化模型参数的值。

为什么使用EM 算法？
EM算法使用启发式的迭代方法，先固定模型参数的值，猜想模型的隐含数据；然后极大化观测数据和隐含数据的似然函数，更新模型参数，不断迭代。
所以，EM算法适用于无监督样本及含有隐含层的样本。

总的来说，如果样本没有隐含层（隐含层是指，只能采样得到结果，不知道分布 ），可以使用的算法有，贝叶斯估计法，极大似然估计法来获得模型参数。但是，如果含有隐含层，那么，EM算法是一种不错的选择。

EM 算法希望解决的问题：


估计含隐变量的联合概率密度。

1.1 Jensen不等式

1.2 EM

1.3 多维高斯分布函数

2、深入了解EM算法

考虑 观测值 V, q(Z) 是所有的隐变量。 独立采样 q(Z) ， 我们得到以下观测数据的似然函数。

$$\ln p(V|\theta )=\sum _{Z}q(Z)\ln p(V|\theta )=\sum _{Z}q(Z)\ln \frac{p(V,Z)|\theta }{p(Z|V,\theta )}$$

$$\ln p(V|\theta )=\sum _{Z}q(Z)\ln \frac{p(V,Z)|\theta }{q(Z)}-\sum _{Z}q(Z)\ln \frac{p(Z|V,\theta )}{q(Z)}(1)$$

$$\ln p(V|\theta )=\mathcal{L}(q,\theta )+\mathrm{K}\mathrm{L}(q\Vert p)$$

其中，式中 KL 散度衡量两个分布的距离，变量离散和连续分别定义如下，
$$\begin{array}{rl}{D}_{KL}(P\Vert Q)& =\sum _{i}P(i)\log \frac{P(i)}{Q(i)}\\ D_{KL}(P\Vert Q)& =\int P(x)\log \frac{P(x)}{Q(x)}dx\end{array}$$
且具有恒大于等于0的性质。

$\mathcal{L}(q,\theta )$ 包含了 观测数据 V 和隐变量 Z 的联合概率分布。
而且 $\mathcal{L}(q,\theta )\le \ln p(V|\theta )$.

也就是 $\mathcal{L}(q,\theta )$ 是观测数据似然函数的下界。

而我们的EM算法就是最大化似然函数。

E-step：
假设我们初始参数$\theta ={\theta }^{old}$, 公式（1）告诉我们需要最大化下界 $\mathcal{L}(q,{\theta }^{old})$,并把${\theta }^{old}$固定，只考虑隐变量分布 q.
同时，似然函数 $\ln p(V|\theta )$ 和 q(Z) 无关， 最大化 $\mathcal{L}(q,{\theta }^{old})$ 相当于 最小化 KL散度，当且仅当 $q(Z)=p(Z|V,{\theta }^{old})$ , (这里是指隐变量的先验 == 隐变量的后验分布)，KL = 0 ， 最小。

E-step 解决的就是评估 $p(Z|V,{\theta }^{old})$.

M-step:
把 q(Z) 固定，并把 E-step中等式 $q(z)=p(Z/V,{\theta }^{old})$ 代入下界$\mathcal{L}(q,\theta )$.

有
$$L(q,\theta )=\sum _{Z}p\left(Z|V,{\theta }^{\text{ old }}\right)\ln p(V,Z|\theta )-\sum _{Z}p\left(Z|V,{\theta }^{\text{ old }}\right)\ln p\left(Z|V,{\theta }^{\text{ old }}\right)$$

$$L(q,\theta )=\mathcal{Q}\left(\theta ,{\theta }^{\text{ old }}\right)+const$$

这里稍微解释下第二项为什么是常数，其实，第二项是 q 分布的负熵，与 $\theta$独立，故为常数。

M-step 考虑 $\theta$ 下，最大化下界，给出新值 ${\theta }^{new}$ 。

重复两个步骤。

一下是两个步骤的简单示意图。

2.1 EM算法的收敛性

EM算法为什么能近似实现对观测数据的极大似然估计呢？

EM算法是通过迭代逐步近似极大化 $L(\theta )$. 我们希望新的估计值能使 $L(\theta )$ 增加。考虑前后 $L(\theta )$ 的差：

利用Jensen不等式，得到下界：

令：

则：

函数$B(\theta ,{\theta }^{(i)})$是$L(\theta )$的一个下界。

任何使得 B 增大的 $\theta$ 都可以使得 $L(\theta )$ 增大。


上式等价于EM算法的一次迭代，即Q函数的极大化。EM算法不断求解下界的极大化来逼近求解对数似然函数的极大化。
在LV额个过程中，对数似然函数 $L(\theta )$不断增大，但EM算法不能保证找到全局最优值。

3、应用
高斯混合模型的参数估计

1、明确隐变量；
2、E：确定Q函数；
3、M：求新一轮迭代的模型参数；

4、EM算法的推广

F函数：

GEM算法,广义期望极大, GEM算法的特点是每次迭代增加F函数值(并不一定是极大化F函数)，从而增加似然函数值。

算法1：极大化F函数

有的时候，并不是那么容易求F函数的极大化，于是，有了算法2。

算法2:使Q函数收敛

通过找θ（i+1）,满足（3）式，来求解极大化。

算法3：条件极大化d

主要进行d次条件极大化，核心思想是，通过保留一个变量θi，固定其他变量θj，来求出极大的θi.然后，通过这个方法，可以求出θ的所有极大值，实现d次条件极大化。

---

5. 实践

5.1 GMM 组分个数确定

常用划分验证集，选取验证集似然函数最高的组分个数；或者带有模型个数惩罚 的BIC， AIC准则。

BIC

5.2 GMM 参数

5.3 GMM
根据李航《统计学习方法》的例子，写EM算法。

import numpy as np
import math
pA, pB, pC = 0.5,0.5,0.5

class EM(object):
    '''
    根据抛掷硬币的结果，利用EM 算法估计三块硬币正面朝上的概率pA,pB,pC三个参数。
    '''
    def __init__(self, prob):
        '''
        初始化三个概率参数。
        '''
        self.pA,self.pB,self.pC = prob
        
    def eStep(self,i):
        '''
        E步，估计观测数据来自硬币B的概率。
        '''
        uB = self.pA*pow(self.pB, data[i])*pow((1 - self.pB),1-data[i])
        uC = (1-self.pA)*pow(self.pC, data[i])*pow((1-self.pC),(1-data[i]))
        return (uB)/(uB + uC)
    
    def mStep(self, data,epoch):
        '''
        M步，更新参数
        '''
        n = len(data)
        for i in range(epoch):
            uBs = np.array([self.eStep(k) for k in range(n)])
            self.pA = 1/n*sum(uBs)
            self.pB = sum(uBs*data)/sum(uBs)
            self.pC = sum((1-uBs)*data)/sum(1-uBs)
            print('epoch {},pA:{},pB:{},pc:{}'.format(i,self.pA,self.pB,self.pC))

        return self.pA,self.pB,self.pC
    
if __name__ == '__main__':
    '''
    为什么迭代1次就收敛了？
    '''
    data = [1,1,0,1,0,0,1,0,1,1]
    em= EM(prob = [0.4,0.6,0.7])
    f = em.mStep(data,epoch = 3)

5.4 使用EM算法来进行GMM聚类：

import numpy as np
import matplotlib.pyplot as plt
import random

def generateData(mu1 = [170, 70], sigma1 = [[4,0],[0,2]], mu2 = [155,52], 
                 sigma2 = [[2,0],[0,1.5]], SEED = 1,  ):
    '''
    产生数据，两种高斯成分的数据，每个样本100个，总共两百个样本
    '''
    np.random.seed(SEED)
#    s1 = np.random.normal(mu1, sigma1, 100)
    s1 = np.random.multivariate_normal(mu1, sigma1, 100)
    np.random.seed(SEED+1)
#    s2 = np.random.normal(mu2, sigma2, 100)
    s2 = np.random.multivariate_normal(mu2, sigma2, 100)
    return s1, s2

s1, s2 = generateData()
data = np.concatenate((s1,s2), axis = 0)

#from numpy.random import multivariate_normal
from scipy.stats import multivariate_normal

def phi(Y, mu_k, cov_k):
    norm = multivariate_normal(mean= mu_k, cov= cov_k)
    return norm.pdf(Y)

def Estep(Y, mu, cov, alpha):
    '''
    计算数据属于每种模型的概率，
    输入： 数据np.array Y; mu:每个模型对于数据每列的均值； cov,每个模型对于拟合数据Y的协方差矩阵，
    输出：gamma 矩阵，代表每个样本属于模型1，2，...的概率
    '''
    N= Y.shape[0]
    K = alpha.shape[0]
    assert N >1
    assert K >1
    gamma = np.mat(np.zeros((N, K)))
    
    prob = np.zeros((N, K))
    for k in range(K):
        prob[:, k] = phi(Y, mu[k], cov[k])
    prob = np.mat(prob)
    #属于每个模型的概率
    for k in range(K):
        gamma[:, k] = alpha[k]*prob[:, k]
    #在模型维度标准化
    for i in range(N):
        gamma[i, :] /= np.sum(gamma[i, :])
    return gamma

def Mstep(Y, gamma):
    '''
    根据给定的样本属于每个样本的概率，优化模型参数 均值和协方差。
    输入：数据np.array Y； 概率矩阵 np.array gamma
    输出：更新的 mu, cov
    '''
    N, D = Y.shape
    K = gamma.shape[1]
    mu = np.zeros((K, D))
    alpha = np.zeros(K)
    cov = []
    
    for k in range(K):
        Nk = np.sum(gamma[:, k])
        #更新均值
        for d in range(D):
            mu[k, d] = np.sum(np.multiply(gamma[:, k], Y[:, d].reshape(-1,1)))/ Nk
        #更新协方差
        cov_k = np.mat(np.zeros((D, D)))
        for i in range(N):
            cov_k += gamma[i, k] * (Y[i] - mu[k]).T * (Y[i] - mu[k])/Nk
        cov.append(cov_k)
        alpha[k] = Nk/N
    cov = np.array(cov)
    return mu, cov, alpha

def scale_data(Y):
	'''
	标准化数据；
	输入：数据Y， Y应该至少具有两维；
	输出：标准化的Y
	'''
	n = len(Y[0])
    for i in range(n):
        max_ = Y[:, i].max()
        min_ = Y[:, i].min()
        Y[:, i] = (Y[:, i] - min_)/ (max_ - min_)
    return Y

def init_params(shape, K, SEED = 1):
	'''
	输出：模型初始化参数 mu, cov, alpha
	'''
    N, D = shape
    np.random.seed(SEED)
    mu = np.random.rand(K, D)
    cov = np.array([np.eye(D)] * K)
    #属于每种模型的概率
    alpha =np.array([1.0 / K] * K)
    return mu, cov, alpha

def GMM(Y, K, times):
	'''
	模型GMM
	'''
    Y = scale_data(Y)
    mu, cov, alpha = init_params(Y.shape, K)
    for i in range(times):
        gamma = Estep(Y, mu, cov, alpha)
        mu, cov, alpha = Mstep(Y, gamma)
        if i%20 == 0:
            print('mu:{0}, cov:{1}, alpha:{2}'.format(mu, cov, alpha))
    return mu, cov, alpha

#模型个数，人为设定，可以根据实验的效果尝试。
K = 2
mu, cov, alpha = GMM(data, K, 1)
# 根据 GMM 模型，对样本数据进行聚类，一个模型对应一个类别
N = data.shape[0]
# 求当前模型参数下，各模型对样本的响应度矩阵
gamma = Estep(data, mu, cov, alpha)
# 对每个样本，求响应度最大的模型下标，作为其类别标识
category = gamma.argmax(axis=1).flatten().tolist()[0]
# 将每个样本放入对应类别的列表中
class1 = np.array([data[i] for i in range(N) if category[i] == 0])
class2 = np.array([data[i] for i in range(N) if category[i] == 1])

# 绘制聚类结果
plt.plot(class1[:, 0], class1[:,1], 'rs', label="class1")
plt.plot(class2[:, 0], class2[:, 1], 'bo', label="class2")
plt.legend(loc="best")
plt.title("GMM Clustering By EM Algorithm")
plt.show()

结果图：

---

最近开通了个公众号，主要分享推荐系统，风控等算法相关的内容，感兴趣的伙伴可以关注下。

---

reference：

1. 黄海广 GitHub;
2. [李航 统计学习方法];
3. 吴恩达课件；
4. 吴恩达课件中文版；
5. EM算法实践GMM python实现;
6. scipy.stats.multivariate_normal
7. EM towards 解释；
8. 斯坦福大学 EM算法解释；
9. Bayesian information criterion;
10. BIC推导过程 ;
11. GMM github;
12. (推荐)EM算法原理总结;