可能是最容易理解的EM算法入门文章

可能是最容易理解的EM算法入门文章

学习过机器学习的同学一定听说过或使用过EM算法,不知道你们第一次见到这个算法是什么感觉,反正我第一次见表情就这个样子

这个推导啥子的也太难了把。不过经过我不停不停不停不停的看这个算法,到今天我突然觉得自己好像明白了,然后我决定把我的理解写成一篇文章,毕竟只有给别人讲明白了才能算自己真正的明白。那么就进入我们这篇文章的主题:EM算法。

我们先讲一下极大似然估计法,然后再引申出EM算法

1.极大似然估计法

假设我们有如下的一维高斯分布

X\sim N=(\mu, \sigma^{2})

X 的概率密度函数为:

f(x;\mu,\sigma^{2})=\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}

其似然函数为

L(\mu,\sigma^{^{2}})=\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}

求对数为

lnL=-\frac{n}{2}ln\pi-\frac{n}{2}ln\sigma^{2}-\frac{1}{2\sigma^{2}\sum_{i=1}^{n}{(x_{i}-\mu)^{2}}}

对其求导,可以得到如下似然方程组

\frac{\delta}{\delta\mu}lnL=\frac{1}{\sigma^{2}}[\sum_{i=1}^{n}x_{i}-n\mu]

\frac{\delta}{\delta\sigma^{2}}=-\frac{n}{2\sigma^{2}}+\frac{1}{2(\sigma^{2})^{2}}[\sum_{i=1}^{n}x_{i}-n\mu]

我们可以使用

  • 梯度下降法
  • 极大似然估计法

这两种方法来根据样本估计高斯分布的参数,具体代码如下:

def cal_guassian_theta():
    # u=1 siga=4
    y = 1 + 2*np.random.randn(1000,1)
    n,m = y.shape

    # 梯度下降
    u1 = np.random.randn(1)[0]
    siga1 = np.random.randn(1)[0]
    lr = 0.001
    for i in range(1000):
        lu = (np.sum(y,axis=0) - n*u1)/(siga1)
        lsiga = (np.sum((y-u1)**2, axis=0)/siga1 - n)/(2*siga1)
        u1 = u1+lr*lu
        siga1 = siga1+lr*lsiga
    print("u1:   %lf, siga1:   %lf"%(u1, siga1))

    # 解析解
    u2 = np.sum(y, axis = 0)/n
    siga2 = np.sum((y-u2)**2, axis = 0)/n
    print("u2:   %lf, siga2:   %lf"%(u2, siga2))

我们使用均值为1,标准差为2的高斯分布随机生成了1000个样本,然后分别使用梯度下降和极大似然估计法两种方式来估计参数,得到的参数如下:

两种方法得到的结果还是挺不错的。

2. EM算法

极大似然算法确实可以很方便的根据样本估算模型的参数,如果样本来自一个以上的模型,我们又不知道某个样本点到底是来自某个模型的,那么此时极大似然算法就无能为力了。

我们依旧用高斯分布来举例子,混合高斯分布的模型如下:

P(y|\theta)=\sum_{k=1}^{K}\alpha_{k}\phi(y|\theta_{k})

其中 \alpha_{k} 是系数, \alpha_{k}\geq0 ,同时 \sum_{k=1}^{K}a_{k}=1\phi(y|\theta_{k}) 是高斯分布密度函数, \theta_{k}=(u_{k}, \sigma_{k}^{2})

这个时候我们需要估计的参数有

\theta=(\alpha_{1}, ..., \alpha_{k};\theta_{1}, .....,\theta_{k})

此时阻挡我们使用极大似然法的原因就是:我们不知道到底哪些样本点由哪个模型生成。

现在假设我们有1000个样本点,由两个独立的高斯分布生成。我们知道其中第一类有300个,第二类有700个,那么我们就可以对两个高斯分布分别使用极大似然法估计他们的参数了。


但事实上我们知道的只有一堆样本点以及其可能的类别数 K ,至于某个样本到底属于那个模型我们是不知道。此时就要到EM算法登场的时候了,EM算法的主要思想如下:

  • E步:先随便设置一下各种参数\theta=(\alpha_{1}, ..., \alpha_{k};\theta_{1}, .....,\theta_{k}),然后再算一下在当前情况下每个样本点属于哪一个模型的概率值;
  • M步:此时我们知道了一个样本点属于某个模型的概率,然后再次计算各个模型的参数(具体计算方法在下面);然后返回上一步,直至算法收敛。

现在我们知道了EM算法的思想,那么EM算法是怎么在第二步估算出各个模型的参数呢。

我们先介绍一些概念:用 Y 来表示观测随机变量的数据, Z 表示隐随机变量的数据(比如上述混合高斯分布样本点属于某个模型的概率), YZ 连在一起称为完全数据(这个我们是没法知道的),观测数据 Y 也被称为不完全数据(这个我们知道), Z 被称为隐变量(我们不知道),假定给定观测数据 Y ,其概率分布是 P(Y|\theta) ,其中 \theta 是需要估计的模型参数,那么不完全数据 Y 的似然函数是 P(Y|\theta) ,对数似然函数为 L(\theta)=logP(Y|\theta)YZ 的联合概率分布是 P(Y,Z|\theta) ,其对数似然函数是 logP(Y,Z|\theta)

EM算法是通过迭代来求 L(\theta)=logP(Y|\theta)的极大似然估计,也就是在估计出来的参数条件下,模型产生给定样本点的概率最大~。因此我们要最大化下式。

L(\theta)=logP(Y|\theta)=log\sum_{Z}P(Y,Z|\theta)=log(\sum_{Z}P(Y|Z,\theta)P(Z|\theta))

上式中最右边的 P(Z|\theta) 指的的一个模型被选择的概率, P(Y|Z,\theta) 是指我们选定了一个模型,此模型产生这个样本点的概率

上图是一维高斯混合分布,黄色的那个高斯分布均值为0,方差为1,被选择的概率为0.3;红色的那个高斯分布均值为3,方差为4,被选择的概率为0.7。

(下面部分参考《李航统计学习》P.159,推导更详细了一点)

EM是通过迭代的方法来逐步逼近近似极大化 L(\theta) ,假设某一次我们得到了模型的参数估计值为 \theta^{(i)} (是一个我们知道的值),我们要求估计新的 \theta可以时 L(\theta) 增大,即 L(\theta )>L(\theta^{(i)}) 。我们计算两者的差

L(\theta)-L(\theta^{(i)})=log(\sum_{Z}P(Y|Z,\theta)P(Z|\theta))-logP(Y|\theta^{(i)})

利用Jenson不等式

L(\theta)-L(\theta^{i})=log(\sum_{Z}P(Y|Z, \theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z, \theta^{(i)})})-logP(Y|\theta^{(i)})\\ \geq \sum_{Z}P(Z|Y, \theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y, \theta^{(i)})}-logP(Y|\theta^{(i)})\\ =\sum_{Z}P(Z|Y, \theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y, \theta^{(i)})}-logP(Y|\theta^{(i)})\sum_{Z}P(Z|Y, \theta^{(i)})\\ =\sum_{Z}P(Z|Y, \theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y, \theta^{(i)})P(Y|\theta^{(i)})}


第一步到第二步除了使用了Jenson不等式,还使用了 P(X|Y)=\frac{P(Y|X)P(X)}{P(Y)} ,其中的 P(X),P(Y) 都被约去。

B(\theta,\theta^{(i)})=L(\theta^{(i)})+\sum_{Z}P(Z|Y, \theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y, \theta^{(i)})P(Y|\theta^{(i)})}

则有 L(\theta) \geq B(\theta, \theta^{(i)}) ,为了使 L(\theta) 尽可能的大,我们应该选择 \theta^{(i+1)} 使  B(\theta, \theta^{(i)}) 达到极大值。即 \theta^{(i+1)}=arg\max_{\theta}B(\theta,\theta^{(i)})

通过省去对 \theta 极大化是常数的项。

\theta^{(i+1)}=arg\max_{\theta}(L(\theta^{(i)})+\sum_{Z}P(Z|Y, \theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y, \theta^{(i)})P(Y|\theta^{(i)})})\\ =arg\max_{\theta}(\sum_{Z}P(Z|Y,\theta^{(i)})logP(Y|Z,\theta)P(Z|\theta))\\ =arg\max_{\theta}(\sum_{Z}P(Z|Y,\theta^{(i)})logP(Y,Z|\theta))

其中 P(Z|Y,\theta^{(i)}) 指的是当我们知道模型的参数和样本的分布情况时,此时隐变量的状态。如果模型为高斯分布,那么 P(Z|Y,\theta^{(i)}) 指的就是样本点属于某个模型的概率; logP(Y,Z|\theta) 就是我们要找到 \theta 使得在当前 Y,Z 的情况下,获得函数的一个极大值。

EM算法的流程如下:

(1)随机选择参数 \theta^{(0)} ,开始迭代

(2)E步:计算 Q(\theta,\theta^{(i)})=\sum_{Z}P(Z|Y,\theta^{(i)})logP(Y,Z|\theta)

(3)M步:最大化 Q(\theta,\theta^{(i)})

(4)重复(2),(3)步直到收敛

对高斯混合模型使用EM算法估计参数,其中第一个高斯分布均值为3,方差为4,系数为0.7;第二个高斯分布均值为0,方差为1,系数为0.3。

def cal_mix_guassian_theta():
    # u = 3, siga = 4, alpha = 0.7
    f1 = 3+2*np.random.randn(700,1)
    # u = 0, siga = 1, alpha = 0.3 
    f2 = 1+np.random.rand(300,1)
    f = np.concatenate((f1,f2),axis=0)
    n,m = f.shape

    alpha = np.random.rand(1,2)
    alpha = alpha/np.sum(alpha, axis=1)
    u = np.random.rand(1,2)
    siga = np.random.rand(1,2)

    # EM算法求解
    for i in range(100):
        #1.E步
        gamma1 = gaussian(f, u[0][0], siga[0][0])
        gamma2 = gaussian(f, u[0][1], siga[0][1])
        gamma = np.concatenate((gamma1, gamma2), axis = 1)
        gamma = alpha*gamma/np.sum(alpha*gamma, axis = 1, keepdims=True)
        #2.M步
        u = np.sum(gamma*f, axis = 0, keepdims=True)/np.sum(gamma, axis = 0,keepdims=True)
        siga = np.sum(gamma*((f-u)**2), axis = 0, keepdims=True)/np.sum(gamma, axis = 0,keepdims=True)
        siga = siga**(1/2)
        alpha = np.sum(gamma, axis = 0, keepdims=True)/n

    print("alpha1:  %lf,alpha2:    %lf"%(alpha[0][0], alpha[0][1]))
    print("u1:  %lf,u2:    %lf"%(u[0][0], u[0][1]))
    print("siga1:  %lf,siga2:    %lf"%(siga[0][0], siga[0][1]))

算法的估计值如下:

发布于 2019-04-09

文章被以下专栏收录

    建立这个公众号的目的是给在从事人工智能相关领域(包括视觉领域,机器人领域,语音处理合成领域,医学人工智能领域,优化领域以及搜索领域等)工作和学习的人士以及对人工智能感兴趣的兴趣团体提供交流和分析信息的平台。 群主会定期组织相关的人工智能相关的学习会和分享会。以及撰写相关的博文供大家学习和讨论。同时也会分享一些开源技术以供大家参考。群主的愿景是希望把世界各个地区的人工智能从业人工凝聚起来,组成一个完整的人工智能团体。

    机器学习爱好者,知识星球:92416895,微信公众号:机器学习初学者,www.ai-start.com,机器学习爱好者qq群:727137612