EM（Expection-Maximization）算法

最新推荐文章于 2021-07-27 07:28:19 发布

如风过境YD

最新推荐文章于 2021-07-27 07:28:19 发布

阅读量469

点赞数

分类专栏： Machine Learning

本文链接：https://blog.csdn.net/qq_35649945/article/details/98672662

版权

Machine Learning 专栏收录该内容

27 篇文章 5 订阅

订阅专栏

贝叶斯分类中我们考虑了样本的属性独立性假设，一对一和一对多的依赖关系，那么在实际情况中，当属性有缺失时应该怎么办呢？此时可以引入EM算法解决问题。
首先来了解一下EM算法，又称期望最大化算法，算法分为两步，E步和M步，算法主要解决概率模型中含有观测变量和隐变量的样本。观测变量即给定数据已知样本属性值，可以用极大似然法或贝叶斯估计来估计模型参数，隐变量即样本的某种属性值有未知(我们知道数据可以分成两类，但是随机抽取一个数据点，我们不知道这个数据点属于第一类还是第二类，它的归属我们观察不到，因此引入一个隐含变量 $\boldsymbol{z}$ 来描述这个现象。)此时不能简单的采用这些方法。

常用场景

观测数据的点 X 既可以是离散的也可以是连续的，但每个数据点相关联都是一组观测向量。
丢失值 Z 是离散的，是从一组固定值中抽取的结果，每个观测值都与一个潜在变量关联。
参数是连续的，但是有两种：一种是参数与所有的数据点相关；另一种是参数与一个特定的隐变量相关。

已知某属性的概率模型，Y表示观测变量的数据，其概率分布 $P(Y|\theta)$ , $\theta$ 是需要估计的的模型参数， Z表示隐变量，其对数似然函数为 $logP(Y|\theta)$ ，j假设Y和Z的联合概率分布为 $P(Y,Z|\theta),\quad$ 则对于该属性的对数似然函数是 $logP(Y,Z|\theta)$ 。

EM算法通过迭代求 $L(\theta)=logP(Y|\theta)=log\sum_zP(Y,Z|\theta)$ 的极大似然估计，每次迭代包含E步：求期望；M步：求极大化；
首先选取参数的初值，记作 $\theta$ ,然后通过以下步骤迭代计算参数的估计值，直至收敛，第 ${i+1}$ 次的迭代过程如下：
E步
计算在模型参数 $\theta^i$ 下观测数据Y的期望
$\begin{aligned} Q\left(\theta, \theta^{(i)}\right) &=E_{z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right] \\ &=\sum_{Z} \log P(Y, Z | \theta) P\left(Z | Y, \theta^{(i)}\right) \end{aligned}$
其中，P(Z|Y, $\theta^i$ )是在给定数据Y和当前参数估计 $\theta^i$ 下隐变量数据Z的条件概率分布。
M步
求使得 $Q(\theta,\theta^i)$ 极大化的 $\theta$ ,确定第i+1次的参数的估计值 $\theta^{i+1}$ ,
$\theta^{i+1}=argmax_{\theta} Q(\theta,\theta^i)$
重复上面的步骤直至收敛。
注：

参数的初值可随意选择，但EM算法对初值是敏感的，
$Q(\theta,\theta^i)$ 中 $\theta$ 是要求的极大化的参数， $\theta^i$ 是当前的参数估计值。
Q函数是完全数据（Y+Z）的对数似然函数 $logP(Y,Z|\theta）$ 关于在给定观测数据Y和当前参数 $\theta^i$ 下对未观测数据Z的条件概率分布 $P(Z|Y,\theta^i)$ 的期望。
每次迭代是在求 $Q$ 函数及其极大。停止迭代的条件一般是对于较小的正数 $\epsilon_1,\epsilon_2$ ,满足：
$||\theta^{i+1}-\theta^i||<\epsilon_1\quad或\quad||Q(\theta^{i+1},\theta^i)-Q(\theta^i,\theta^i)||<\epsilon_2$
则停止迭代。

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

高斯混合模型即存在一个样本集中，不同的类符合不同参数的高斯分布，假设有N个样本，有K个类，有高斯分布 $N_1(\mu_1,\sigma_1^2),N_2(\mu_2,\sigma_2^2),...,N_k(\mu_k,\sigma_k^2)$ 。取第k个模型的概率为 $\alpha_k$ , $y_i$ 表示观测数据。
高斯混合模型一般是如下形式的概率分布模型：
$P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k)$
其中K表示分模型的个数， $\alpha_k>0,\sum_1^K\alpha_k=1,\theta_k=(\mu_k,\sigma_k^2),$
$\phi\left(y | \theta_{k}\right)=\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)$

E步
估计样本数据由每个模型（类）生成的概率，对于每个观测数据 $y_i$ (已知数据),它由第k个类生成的概率：
$\gamma(i,k)=\frac{\alpha_k*\phi(x_i|\mu_k,\sigma_k)}{\sum_{j=1}^K\phi(x_i|\mu_j,\sigma_k)}=\frac{\alpha_k*f_{\mu_k,\sigma_k}(x_i)}{\sum_{j=1}^Kf_{\mu_k,\sigma_k}(x_i)},表示x_i属于第k个类（高斯分模型）的概率。$
上式中 $\mu，\sigma$ 同样也是估计值，因此采用迭代法,在计算 $\gamma(i,k)$ 时假定 $\mu,\sigma$ 已知（即设置初始值）。
M步
由此可以算出任意一个 $x_i$ 的值属于第k个模型的概率。然后使用 $x_i*\gamma(i,k),k={1,2,..,K}$ ,
然后将这些数据重新计算 $\mu_i,\sigma_i$ , （ $\hat{\gamma}_{j k}$ 是迭代过程中在当前模型参数下的第j个观测数据来自第k个分模型的概率，称为观测数据 $y_j$ 的响应度。与上面的 $\gamma(i,k)$ 意义相同）。在参数最大化的过程中使用对数似然估计，由于极大似然估计可推得极大似然估计的结果与矩估计的结果近似相同，因此可以采用下面的方法进行参数的估计。
$\begin{array}{c}{\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}* y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K}\\ \\ {\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K}\end{array}\\ \hat{\alpha}_k=\frac{\sum_{j=1}^N\hat{\gamma}_{jk}}{N}$ ，
其中 $\sum_{j=1}^N\gamma_{jk}$ 表示所有样本属于第k个分模型的概率相加之和。
然后在迭代计算 $\gamma(i,k)$ 继续执行E步，M步，直至收敛（参数不再变化）。

举例：
假设有N=10000个样本，代表有10000个学生，已知每个样本有属性身高h和性别sex，假设样本身高分布符合高斯分布，显然男女符合不同的高斯分布。假设男生符合分布 $N_1(\mu_1,\sigma_1)$ ，女生符合分布 $N_2(\mu_2,\sigma_2)$ 。已知先验概率P(x=男)=P(x=女)= $\pi$ =0.5;
根据EM算法：

首先为 $\mu_1,\mu_2,\sigma_1,\sigma_2$ 赋初值。
假设有样本h=1.80,根据初始模型参数值计算该样本为男或女的概率：
$P(x=男|h=1.80)=P(x=男)*f_{\mu_1,\sigma_1}(h=1.80)=m$ ;
$P(x=女|h=1.80)=P(x=女)*f_{\mu_1,\sigma_1}(h=1.80)=n$ ;
同理，计算其他已知h的样本;
然后得出分模型对于观测数据（已知数据） $y_j$ 的响应度
$\hat{\gamma}_{j k}=\frac{\alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K$
即： $\hat{\gamma}_{(j ，男)}=\frac{\pi*m}{\pi*m+\pi*n}， \hat{\gamma}_{(j ,女)}=\frac{\pi*n}{\pi*m+\pi*n}$
然后更新数据，例如 $\hat{\gamma}_{(j,女)}$ , $\quad1.80*\hat{\gamma}_{(j，男)}$ ,将新得到的数据按照k=男，女的分模型分开，然后求参。
然后求参数最大化，即：
$\begin{array}{c}{\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2}\\ \\ {\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2,}\end{array}\\ \hat{\alpha}_k=\frac{\sum_{j=1}^N\hat{\gamma}_{jk}}{N}，相当于求P(x=男)，P(x=女) 上面的1,2分别代表男女。$
然后重复上面的步骤，直到收敛。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
DEBUG=True
def debug(*args,**kwargs):
    global DEBUG
    if DEBUG:
        print(*args,**kwargs)
    
def phi(Y,mu_k,cov_k):
    norm=multivariate_normal(mean=mu_k,cov=cov_k)#高斯分布函数
    return norm.pdf(Y)
#E步计算期望（响应度，当前参数下样本属于第k个模型的概率）
def getExpectation(Y,mu,cov,alpha):
    N=Y.shape[0]
    K=alpha.shape[0]
    #判断是否只有一个模型或样本
    assert N>1,"must be more than one sample"
    assert K>1,"must be more than one model"
    #gamma响应度矩阵
    gamma=np.mat(np.zeros((N,K)))
    #prob各模型中所有样本出现的概率，gamma=alpha*prob
    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
#M步，计算均值，协方差矩阵，alpha各模型占总体的概率（初始设为均匀分布）
def maximize(Y,gamma):
    N,D=Y.shape
    #K表示模型数
    K=gamma.shape[1]
    mu=np.zeros((K,D))
    cov=[]
    alpha=np.zeros(K)
    for k in range(K):
        #Nk第k个模型对所有样本的响应度之和
        Nk=np.sum(gamma[:,k])
        #更新mu,mu[k,d]表示均值矩阵mu的第k各模型第d维的均值
        for d in range(D):
            mu[k,d]=np.sum(np.multiply(gamma[:,k],Y[:,d]))/Nk
        #更新cov  
        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
        alpha[k]=Nk/N
    cov=np.array(cov)
    return mu,cov,alpha
def scale_data(Y):
    for i in range(Y.shape[1]):
        max_=Y[:,i].max()
        min_=Y[:,i].min()
        Y[:,i]=(Y[:,i]-min_)/(max_-min_)
    debug("Data scaled")
    return Y
def init_params(shape,K):
    N,D=shape
    mu=np.random.rand(K,D)
    cov=np.array([np.eye(D)]*K)
    alpha=np.array([1.0/K]*K)
    debug("Parameters initalized")
    debug('mu:',mu,'cov:',cov,'alpha:',alpha,sep='\n')
    return mu,cov,alpha

def GMM_EM(Y,K,times):
    Y=scale_data(Y)
    #print(Y)
    Y1=np.array(Y)
    mu,cov,alpha=init_params(Y.shape,K)
    for i in range(times):
        gamma=getExpectation(Y,mu,cov,alpha)
        mu,cov,alpha=maximize(Y,gamma)
        if i in (0,20,50):
            category=gamma.argmax(axis=1).flatten().tolist()[0]
            class1=np.array([Y1[i] for i in range(N) if category[i]==0])
            class2=np.array([Y1[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,itercount={}'.format(i))
            plt.show()
    debug('{sep}Result {sep}'.format(sep='-'*20))
    debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n")
    return mu,cov,alpha

import matplotlib.pyplot as plt
DEBUG=True
Y=np.loadtxt('gmm.data')
matY=np.matrix(Y,copy=True)
K=2
mu,cov,alpha=GMM_EM(matY,K,100)

N=Y.shape[0]
gamma=getExpectation(matY,mu,cov,alpha)
category=gamma.argmax(axis=1).flatten().tolist()[0]
class1=np.array([Y[i] for i in range(N) if category[i]==0])
class2=np.array([Y[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()

在这里插入图片描述

3.600000 79.000000
1.800000 54.000000
3.333000 74.000000
2.283000 62.000000
4.533000 85.000000
2.883000 55.000000
4.700000 88.000000
3.600000 85.000000
1.950000 51.000000
4.350000 85.000000
1.833000 54.000000
3.917000 84.000000
4.200000 78.000000
1.750000 47.000000
4.700000 83.000000
2.167000 52.000000
1.750000 62.000000
4.800000 84.000000
1.600000 52.000000
4.250000 79.000000
1.800000 51.000000
1.750000 47.000000
3.450000 78.000000
3.067000 69.000000
4.533000 74.000000
3.600000 83.000000
1.967000 55.000000
4.083000 76.000000
3.850000 78.000000
4.433000 79.000000
4.300000 73.000000
4.467000 77.000000
3.367000 66.000000
4.033000 80.000000
3.833000 74.000000
2.017000 52.000000
1.867000 48.000000
4.833000 80.000000
1.833000 59.000000
4.783000 90.000000
4.350000 80.000000
1.883000 58.000000
4.567000 84.000000
1.750000 58.000000
4.533000 73.000000
3.317000 83.000000
3.833000 64.000000
2.100000 53.000000
4.633000 82.000000
2.000000 59.000000
4.800000 75.000000
4.716000 90.000000
1.833000 54.000000
4.833000 80.000000
1.733000 54.000000
4.883000 83.000000
3.717000 71.000000
1.667000 64.000000
4.567000 77.000000
4.317000 81.000000
2.233000 59.000000
4.500000 84.000000
1.750000 48.000000
4.800000 82.000000
1.817000 60.000000
4.400000 92.000000
4.167000 78.000000
4.700000 78.000000
2.067000 65.000000
4.700000 73.000000
4.033000 82.000000
1.967000 56.000000
4.500000 79.000000
4.000000 71.000000
1.983000 62.000000
5.067000 76.000000
2.017000 60.000000
4.567000 78.000000
3.883000 76.000000
3.600000 83.000000
4.133000 75.000000
4.333000 82.000000
4.100000 70.000000
2.633000 65.000000
4.067000 73.000000
4.933000 88.000000
3.950000 76.000000
4.517000 80.000000
2.167000 48.000000
4.000000 86.000000
2.200000 60.000000
4.333000 90.000000
1.867000 50.000000
4.817000 78.000000
1.833000 63.000000
4.300000 72.000000
4.667000 84.000000
3.750000 75.000000
1.867000 51.000000
4.900000 82.000000
2.483000 62.000000
4.367000 88.000000
2.100000 49.000000
4.500000 83.000000
4.050000 81.000000
1.867000 47.000000
4.700000 84.000000
1.783000 52.000000
4.850000 86.000000
3.683000 81.000000
4.733000 75.000000
2.300000 59.000000
4.900000 89.000000
4.417000 79.000000
1.700000 59.000000
4.633000 81.000000
2.317000 50.000000
4.600000 85.000000
1.817000 59.000000
4.417000 87.000000
2.617000 53.000000
4.067000 69.000000
4.250000 77.000000
1.967000 56.000000
4.600000 88.000000
3.767000 81.000000
1.917000 45.000000
4.500000 82.000000
2.267000 55.000000
4.650000 90.000000
1.867000 45.000000
4.167000 83.000000
2.800000 56.000000
4.333000 89.000000
1.833000 46.000000
4.383000 82.000000
1.883000 51.000000
4.933000 86.000000
2.033000 53.000000
3.733000 79.000000
4.233000 81.000000
2.233000 60.000000
4.533000 82.000000
4.817000 77.000000
4.333000 76.000000
1.983000 59.000000
4.633000 80.000000
2.017000 49.000000
5.100000 96.000000
1.800000 53.000000
5.033000 77.000000
4.000000 77.000000
2.400000 65.000000
4.600000 81.000000
3.567000 71.000000
4.000000 70.000000
4.500000 81.000000
4.083000 93.000000
1.800000 53.000000
3.967000 89.000000
2.200000 45.000000
4.150000 86.000000
2.000000 58.000000
3.833000 78.000000
3.500000 66.000000
4.583000 76.000000
2.367000 63.000000
5.000000 88.000000
1.933000 52.000000
4.617000 93.000000
1.917000 49.000000
2.083000 57.000000
4.583000 77.000000
3.333000 68.000000
4.167000 81.000000
4.333000 81.000000
4.500000 73.000000
2.417000 50.000000
4.000000 85.000000
4.167000 74.000000
1.883000 55.000000
4.583000 77.000000
4.250000 83.000000
3.767000 83.000000
2.033000 51.000000
4.433000 78.000000
4.083000 84.000000
1.833000 46.000000
4.417000 83.000000
2.183000 55.000000
4.800000 81.000000
1.833000 57.000000
4.800000 76.000000
4.100000 84.000000
3.966000 77.000000
4.233000 81.000000
3.500000 87.000000
4.366000 77.000000
2.250000 51.000000
4.667000 78.000000
2.100000 60.000000
4.350000 82.000000
4.133000 91.000000
1.867000 53.000000
4.600000 78.000000
1.783000 46.000000
4.367000 77.000000
3.850000 84.000000
1.933000 49.000000
4.500000 83.000000
2.383000 71.000000
4.700000 80.000000
1.867000 49.000000
3.833000 75.000000
3.417000 64.000000
4.233000 76.000000
2.400000 53.000000
4.800000 94.000000
2.000000 55.000000
4.150000 76.000000
1.867000 50.000000
4.267000 82.000000
1.750000 54.000000
4.483000 75.000000
4.000000 78.000000
4.117000 79.000000
4.083000 78.000000
4.267000 78.000000
3.917000 70.000000
4.550000 79.000000
4.083000 70.000000
2.417000 54.000000
4.183000 86.000000
2.217000 50.000000
4.450000 90.000000
1.883000 54.000000
1.850000 54.000000
4.283000 77.000000
3.950000 79.000000
2.333000 64.000000
4.150000 75.000000
2.350000 47.000000
4.933000 86.000000
2.900000 63.000000
4.583000 85.000000
3.833000 82.000000
2.083000 57.000000
4.367000 82.000000
2.133000 67.000000
4.350000 74.000000
2.200000 54.000000
4.450000 83.000000
3.567000 73.000000
4.500000 73.000000
4.150000 88.000000
3.817000 80.000000
3.917000 71.000000
4.450000 83.000000
2.000000 56.000000
4.283000 79.000000
4.767000 78.000000
4.533000 84.000000
1.850000 58.000000
4.250000 83.000000
1.983000 43.000000
2.250000 60.000000
4.750000 75.000000
4.117000 81.000000
2.150000 46.000000
4.417000 90.000000
1.817000 46.000000
4.467000 74.000000