机器学习课程笔记【十】- 混合高斯与期望最大化算法

本节为吴恩达教授机器学习课程笔记第十部分，混合高斯与期望最大化算法，主要包括：混合高斯模型与无监督的参数估计方法EM算法，EM算法即期望最大化算法的一般形式，一般形式EM算法重新对混合高斯模型参数拟合进行分析。最后附上EM算法拟合混合高斯分布参数的python代码。

1. 混合高斯与期望最大化算法用于密度估计

  给定训练集 { x ( 1 ) , . . . , x ( m ) } \{x^{(1)},...,x^{(m)}\} {x(1),...,x(m)}，因为我们讨论的是无监督学习算法，所以并没有对应的标签数据。
  我们希望用一个联合分布$p(x^{(i)},z^{(i)})=p(x^{(i)}|z^{(i)})p(z^{(i)})$，其中$z^{(i)}$是一个多项式分布，其中多项式的参数中有${\phi }_j\ge 0,{\sum }_{j=1}^k{\phi }_j=1并且{\phi }_j指出了p(z^{(i)}=j)$，并且$x^{(i)}|z^{(i)}=j\sim N({\mu }_j,{\Sigma }_j)$。令$k$表示$z^{(i)}$可以取的值的数量。也就是说，我们的模型假定每一个$x^{(i)}$都是由从 { 1 , 2 , . . . , k } \{1,2,...,k\} {1,2,...,k}中随机取值的$z^{(i)}$产生的，也即$x^{(i)}$是一个在$z^{(i)}$上的高斯分布中的$k$个值中的一个，这就是混合高斯模型。同时假定$z^{(i)}$是潜在的随机变量，也就是其取值是隐藏的或者未被观测到，这也将会增加估计问题的难度。

实际上混合高斯就是通过求解两个高斯模型并通过一定得权重将两个高斯模型融合成一个模型。

  可以看到，我们的模型的参数是$\phi 、\mu 和\Sigma$，为了对参数进行估计，我们根据数据写出极大似然函数：

  然而令其导数为0，我们无法找出参数的极大似然估计，随机变量$z^{(i)}$表明了每个$x^{(i)}$所属的$k$个高斯分布的值，如果我们直到$z^{(i)}$，那么最大似然问题就很容易解，此时对数似然函数可以写为：

  函数最大化，得到参数的值分别为：

  事实上，如果$z^{(i)}$已知，那么最大似然函数变的和高斯判别分析相似，唯一的不同在于这里$z^{(i)}$扮演了类别标签的角色。然而在密度估计问题中该值未知，这就要使用EM即期望最大化算法了。
  EM算法是一个有两个步骤的迭代算法，在我们的密度估计问题中应用，它的第一步即"E-step"尝试去猜测$z^{(i)}$的值，第二部即"M-step"根据这个猜测来更新模型中的参数。由于假设第一步的猜测是正确的，所以最大化的过程就会很简单，算法步骤如下所示：

  在第一步中，我们计算根据已有的参数来计算$z^{(i)}$在给定$x^{(i)}$时的后验概率，比如说使用贝叶斯法则，有：

  其中:

  是通过对一个均值${\mu }_j$和协方差${\Sigma }_j$的高斯分布的密度及进行估计得到的，${\omega }_j^{(i)}$时对$z^{(i)}$的一个弱估计。{若估计代表我们的猜测时概率从[0,1]区间取值；相反强估计即从{1,…,k}中取值}。
  EM算法和K均值算法相比，这里采用一个弱估计而不是一个实在的聚类依附关系。二者相同的地方在于他们都可能收敛到局部最优，所以可以用不同的初始参数重复运行几次算法。
  但是这里并没有给出一个严格的证明，来保证算法确实收敛。

2. 期望最大化算法全貌

  第一小节给出了应用EM算法对高斯混合分布的参数进行拟合，该节将给出EM算法的全貌，并且应用在更一般的带有隐含参数的估计问题中。

2.1 Jensen不等式

  令$f$为定义域为实数的函数，之前提到，如果如果$f^{\prime \prime }(x)\ge 0$那么$f$是凸函数，如果对于所有的$x$都有$f^{\prime \prime }(x)>0$我们说函数是严格凸函数。这里我们的函数$f$以向量为输入，推广之后即海森矩阵$H$是正的半正定矩阵($H\ge 0$)则函数为凸函数(严格凸对应于海森矩阵为严格的半正定的正值矩阵，H>0)，Jensen不等式定义如下：

  如果$f$是严格凸的，那么当且仅当$X=E[X]的概率为1时，E[(X)]=f(EX)$。
  对定理的理解，如下图所示：

2.2 EM算法

  假设我们现在面临一个估计问题，训练集 { x ( 1 ) , . . . , x ( m ) } \{x^{(1)},...,x^{(m)}\} {x(1),...,x(m)}包含$m$个独立的样本，我们的目标是对模型$p(x,z)$进行参数拟合，对数似然函数如下：

  因为直接求参数$\theta$的极大似然估计非常困难，所以这里引入一个潜在的随机变量$z^{(i)}$【通常这个潜在的随机变量的分布很容易观测到】，这样就简化了极大似然估计问题。
  在这样的设定中，EM算法能够给出一个更高效的极大似然估计。直接最大化$l(\theta )$不现实，我们策略是不断构造$l$的下界，然后再M步中对这个下界进行优化，不断提升这个下界，直到收敛。对于每一个$i$，令$Q_i$表示$z$的某个分布（${\Sigma }_z{Q}_i(z)=1,Q_i(z)\ge 0$），考虑：

  最后一步就是用Jensen不等式推到的，注意这里f(x)=logx是一个凹函数，公式的这一项：

  表示量：

  的期望，其中$z^{(i)}$由分布$Q_i$给出，通过Jensen不等式，有：

  这样，我们就给出$l(\theta )$的一个较小的下界，然后通过最大化这个下界得到一个较好的解，那么新的问题就出现了，如何选择$Q_i$呢？
  为了使得我们估计的下界对于特定的$\theta$值足够紧确，我们需要把上面式子中涉及Jensen不等式的推导步骤变为等式，为了保持等号成立，随机变量需要从一组常量中取值才能保证期望与本身相等，也就是说我们需要：

  而且这个$c$不能依赖$z^{(i)}$，即选择一个$Q_i$与后面的部分具有线性关系：

  实际上因为$Q_i$表示的是某种分布，所以我们已知：

  这样我们将$p(x^{(i)},z^{(i)};\theta )$归一化，得下面得式子：

  这样我们简单地把$Q_i$视为给定$x^{(i)}$和参数$\theta$时$z^{(i)}$的后验分布。
  这样给定了$Q_i$，根据上面的推导，我们就得到了我们想要最大化的对数似然函数的一个下界，这就是E-Step。之后的M步我们就根据E-Step的假设来拟合参数，步骤如下所示：

  至此我们还是没有证明EM算法确实收敛，现在假定${\theta }^{(t)}$和${\theta }^{(t+1)}$为EM算法两次成功的迭代中的参数，我们只需证明$l({\theta }^{(t)})\le l({\theta }^{(t+1)})$就可以说明E算法在不断提升这个下界，这个证明的关键使$Q_i$的选择，如果某一次迭代完成后参数为${\theta }^{(t)}$，那么我们可以令：

  这样有：

  这样${\theta }^{(t+1)}$可以通过最大化上式右边的部分来获得，因此有：

  第一行的不等式成立是因为：

  第二行的不等式成立是因为${\theta }^{(t+1)}$是通过如下过程选取的：

  至此我们证明了EM算法会使得似然函数持续收敛，那EM算法终止的条件就是某两次迭代之间变化的范围已经小于某个可以接受的值，那么认为算法已经收敛。
  最后值得一提的是，如果我们定义：

  则有：

  EM算法也可以看作对于$J$的坐标下降过程，其中E步关于$Q$最大化$J$，M步关于$\theta$最大化$J$。

2.3 重新审视混合高斯

  学习了一般EM算法的定义后，我们重新审视拟合混合高斯的参数$\theta 、\mu 、\Sigma$的例子，为了方便这里只作M步更新前两个参数的假设。
  这里E步非常简单，我们可以这样计算：

  这里$Q_i(z^{(i)}=j)$表示了$z^{(i)}$在服从分布$Q_i$的约束下，取值为$j$的概率。
  在M步我们最大化这个式子，来对参数进行拟合：

  首先关于${\mu }_l$来最大化这个式子，即这个式子关于${\mu }_l$求偏导，得：

  令导数为0，得到：

  这与第一小节中得结论相同。再来看如何拟合参数${\phi }_j$，通过观察发现，式子中与${\phi }_j$有关得部分很少，我们只需要最大化这个即可：

  然而有一个额外得约束条件即所有得${\phi }_j$和为1，因为这个参数表示一个概率值即：

  为了解决这个带有限制条件得优化问题，我们构建拉格朗日乘子式：

  其中$\beta$是拉格朗日乘子，对这个乘子式关于${\phi }_j$求导得：

  令导数为0，有：

  因为：

  且：

  则有：

  这个推导又用到了${\omega }_j^{(i)}=Q_i(z^{(i)}=j)$且概率和必为1，即${\Sigma }_j{\omega }_j^{(i)}=1$,这样：

最后附上代码，EM算法拟合混合高斯分布参数的python代码。

import numpy as np
import random
import math
import time

def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
    '''
    初始化数据集
    :return: 混合了两个高斯分布的数据
    '''
    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):
    '''
    根据高斯密度函数计算值
    :param dataSetArr: 可观测数据集
    :param mu: 均值
    :param sigmod: 方差
    :return: 整个可观测数据集的高斯分布密度（向量形式）
    '''
    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的响应度
    '''
    # 计算y0的响应度
    # 先计算模型0的响应度的分子
    gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
    # 模型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):
    mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)
    mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)

    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算法进行参数估计
    '''
    # 将可观测数据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()
    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)

---

欢迎扫描二维码关注微信公众号 深度学习与数学   [每天获取免费的大数据、AI等相关的学习资源、经典和最新的深度学习相关的论文研读，算法和其他互联网技能的学习，概率论、线性代数等高等数学知识的回顾]