Catalog
  1. 1. 前情提要
    1. 1.1 为什么提出 —— 感性认知
      1. 以一个例子说明:
    2. 1.2 数学基础
  2. 2. 导出EM算法
  3. 3. EM算法的合理性
    1. 1. 证明单调递增
    2. 2.证明有上界
    3. 3.稳定到哪个值
  • 4. GMM中应用
  • 5. 代码
    1. EM 算法的实现
  • 《统计学习方法》(CH09 EM算法_GMM模型)

    EM算法(Expectation-Maximization algorithm)不是模型,只是一种用于计算模型中参数的方法。

    1. 前情提要

    1.1 为什么提出 —— 感性认知

    优秀的解说:如何感性地理解EM算法?

    以一个例子说明

    • 正常概率题:

      $A,B,C$三枚质地不均匀的硬币,正面朝上的概率分别为$P(A)=0.3,P(B)=0.8,P(C)=0.6$。先投掷A硬币,若A正面朝上,投掷$B$硬币,否则投$C$硬币。$B,C$的投掷结果为最终结果,求投掷出正面的概率。

      解答:

      $P($正面$)= P(A)P(B)+[1-P(A)]P(C)=0.30.8+0.70.6=0.66$

    • 变态概率题:

      现在不知道$P(A),P(B),P(C)$,知道100次投掷中$P(正)$出现了66次,求$P(A),P(B),P(C)$。

      思路:

      1. $P(正面)$还是等于$ P(A)P(B)+[1-P(A)]P(C)$;
      2. 先假设$P(A)=0.5,P(B)=0.5,P(C)=0.5$,带进去计算,$P(正)=0.5$,不对怎么办?
      3. 用P(正)===>估算B、A、C=====>代入对比P(正),再循环这个步骤

    1.2 数学基础

    1. 极大似然估计$(Maximum\ Likelihood\ Estimation)$

      已知观测数据,求常数参数

    2. $Q$函数

      定义:完全数据的对数似然函数关于在给定观测数据$Y$和当前参数下对未观测数据$Z$的条件概率分布$P(Z|Y,)$的期望为$Q$函数.

    3. $Jensen$不等式

      对于凸函数而言

    4. 拉格朗日乘子法$(Lagrange\ Multiplier)$

      拉格朗日乘法在$M-step$求导之后会使用

    5. 混合高斯模型$(Gaussian\ mixture\ model)$

      $EM$算法就是提出来解决$GMM$中因为含有隐变量,导致$log$极大似然函数的求解中$log$后面有加法,从而无法求解$MLE$

    2. 导出EM算法

    3. EM算法的合理性

    EM算法:

    $$ \widehat{\theta }=\arg \max _{\theta }\sum _{z}\log P\left( Y,Z|\theta \right) P\left( Z|Y,\theta ^{(i)}\right) dz$$

    需要证明以下几个点:

    1. 经过EM算法一定能保证极大似然函数估计值不断变大吗?=====>证明单调递增
    2. 这个估计值会收敛到哪一值吗?========================> 证明有上界
    3. 这个上界是似然函数的最大值或者一个极大值吗?===========>得出这个值

    1. 证明单调递增

    $$ \log P\left( Y|\theta \right) =\log P\left( Y,Z|\theta \right) -\log P\left( Z|Y,\theta \right) $$

    两边同时积分,因为左边没有跟Z相关的函数,所以还是$\log P( Y|\theta ) $

    $$ \log P\left( Y|\theta \right) =\sum _{Z}\log P\left( Y,Z|\theta \right) P\left( Z|Y,\theta ^{\left( i\right) }\right) -\sum _{Z}\log P\left( Z|Y,\theta \right) P\left( Z|Y,\theta ^{\left( i\right) }\right) $$

    令:$$Q\left( \theta ,\theta ^{(i)}\right) = \sum _{Z}\log P\left( Y,Z|\theta \right) P\left( Z|Y,\theta ^{\left( i\right) }\right) $$

    令:$$H\left( \theta ,\theta ^{(i)}\right)=\sum _{Z}\log P\left( Z|Y,\theta \right) P\left( Z|Y,\theta ^{\left( i\right) }\right)$$

    则:

    $$\begin{align} \log P\left( Y|\theta \right) =& Q\left( \theta ,\theta ^{(i)}\right) - H\left( \theta ,\theta ^{(i)}\right) \\log P\left( Y|\theta^{(i+1)} \right)- \log P\left( Y|\theta^{(i)} \right)=& Q\left( \theta^{(i+1)},\theta ^{(i)}\right) - Q\left( \theta ^{(i)},\theta ^{(i)}\right)-[H\left( \theta ^{(i+1)},\theta ^{(i+1)}\right)-H\left( \theta^{(i)} ,\theta ^{(i)}\right)] \end{align}$$

    1. $\theta ^{(i+1)}$是极大,所以$Q\left( \theta ^{(i+1)},\theta ^{(i)}\right)$ 是肯定大于$ Q\left( \theta^{(i)} ,\theta ^{(i)}\right)$

    2. 由$Jensen$不等式,$$H\left( \theta ^{(i+1)},\theta ^{(i)}\right)-H\left( \theta^{(i)} ,\theta ^{(i)}\right)=0 $$

    所以,$\log P\left( Y|\theta^{(i+1)} \right)- \log P\left( Y|\theta^{(i)} \right)\geq 0$$,$$ \log P\left( Y|\theta \right)$ 是一个单调递增的函数, $P( Y|\theta ^{(i+1)})\geq P( Y|\theta ^{(i)})$

    2.证明有上界

    因为 $P( Y|\theta )$有界,所以必定有界。

    3.稳定到哪个值

    4. GMM中应用

    高斯模型密度函数:$$f\left( x\right) =\dfrac{1}{\sqrt{2\pi }\sigma}e^{\left( -\dfrac{\left( x-\mu \right) ^{2}}{2\sigma^{2}}\right) }$$

    高斯混合模型优秀讲解:徐亦达老师讲解

    高斯混合模型概率分布模型:

    1. 明确隐函数 $Z$:写出要求的 $ \log P( Y|\theta ) $
    1. $E-Step$ :确定 $Q(\theta,\theta^{(i)})$
    1. $M-Step$:求最大值

    5. 代码

    EM 算法的实现

    1. 计算E(expectation),期望
    2. 计算M(maximization),最大值
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    import numpy as np
    import matplotlib.mlab as mlab
    import matplotlib.pyplot as plt
    import random
    import math
    import time
    def set_GM(mu, sigma, sample_num, alpha):
    '''
    :func 生成单个高斯分布
    :param mu: 高斯的均值
    :param sigma: 高斯的方差
    :param sample_num: 数据量
    '''
    ## 1. 用生成随机数的方式生成
    # 用numpy.random.normal() 高斯分布随机数
    data = np.random.normal(mu, sigma, int(sample_num * alpha))
    #
    # 画图
    plt.hist(data, bins=100, normed=True)
    plt.show()
    ## 这种方式生成的数很随机,
    # 虽然形状有点像正态分布,不能算是严格意义上的正态分布。
    # 代码中设定的smaple_num 加大,相当于增加样本数量,
    # 那么整个图像就会更加接近正态分布的形状。

    return data

    def EM(dataSet, EM_iter):
    '''
    算法依据“9.3.2 高斯混合模型参数估计的EM算法” 算法9.2
    :param dataSet: 数据集
    :param EM_iter: 迭代次数
    :return: 估计的参数
    '''
    dataSet = np.array(dataSet)
    # 1.取theta初值
    alpha0 = 0.5; mu0 = 0; sigmod0 = 1
    alpha1 = 0.5; mu1 = 1; sigmod1 = 1
    # 2.迭代
    step = 0
    while (step < EM_iter):
    # 3. E()计算期望
    gamma0, gamma1 = E(dataSet, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
    # 4. M()计算最大值
    mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = \
    M(mu0, mu1, gamma0, gamma1, dataSet)
    # 迭代次数加一
    step += 1
    # 返回theta_N
    return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1

    def E(dataSet, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
    '''
    计算分模型对数据的响应度:也就是每个数据在每个模型上的概率是多少

    每个模型都是一个高斯分布,就是求每个数据在每个高斯上的期望
    '''
    gamma_1 = alpha0 * E_Gauss(dataSet,mu0,sigmod0)
    gamma_2 = alpha1 * E_Gauss(dataSet,mu1,sigmod1)

    sum_pro = gamma_1 + gamma_2

    gamma_1 = gamma_1 / sum_pro
    gamma_2 = gamma_2 / sum_pro

    return gamma_1, gamma_2

    def M(mu0, mu1, gamma0, gamma1, dataSet):
    '''
    新一轮迭代的参数
    '''
    mu0_new = np.dot(gamma0, dataSet) / np.sum(gamma0)
    mu1_new = np.dot(gamma1, dataSet) / np.sum(gamma1)

    sigmod0_new = math.sqrt(np.dot(gamma0, (dataSet - mu0)**2) / np.sum(gamma0))
    sigmod1_new = math.sqrt(np.dot(gamma1, (dataSet - 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 E_Gauss(dataSet, mu, sigma):
    '''
    公式就是高斯分布密度,正态分布的密度公式
    '''
    Gauss_prob = (1 / (math.sqrt(2*math.pi) * sigma**2)) * np.exp(-1 * (dataSet - mu)**2) / (2 * sigma**2)

    return Gauss_prob

    if __name__ == '__main__':
    start = time.time()

    alpha0 = 0.3; mu0 = -2; sigmod0 = 0.5
    alpha1 = 0.7; mu1 = 0.5; sigmod1 = 1

    data_1 = set_GM(mu0,sigmod0,10000,alpha0)
    data_2 = set_GM(mu1,sigmod1,10000,alpha1)
    #初始化总数据集
    #两个高斯分布的数据混合后会放在该数据集中返回
    dataSet = []
    #将第一个数据集的内容添加进去
    dataSet.extend(data_1)
    #添加第二个数据集的数据
    dataSet.extend(data_2)
    #对总的数据集进行打乱(其实不打乱也没事,只不过打乱一下直观上让人感觉已经混合了
    # 读者可以将下面这句话屏蔽以后看看效果是否有差别)
    # random.shuffle(dataSet)

    plt.hist(dataSet, bins=100, normed=True)
    plt.show()

    alpha0_new, mu0_new, sigmod0_new, alpha1_new, mu1_new, sigmod1_new = EM(dataSet,1000)
    data_1 = set_GM(mu0_new,sigmod0_new,10000,alpha0_new)
    data_2 = set_GM(mu1_new,sigmod1_new,10000,alpha1_new)
    #初始化总数据集
    #两个高斯分布的数据混合后会放在该数据集中返回
    dataSet = []
    #将第一个数据集的内容添加进去
    dataSet.extend(data_1)
    #添加第二个数据集的数据
    dataSet.extend(data_2)
    #对总的数据集进行打乱(其实不打乱也没事,只不过打乱一下直观上让人感觉已经混合了
    # 读者可以将下面这句话屏蔽以后看看效果是否有差别)
    # random.shuffle(dataSet)

    plt.hist(dataSet, bins=100, normed=True)
    plt.show()
    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
    alpha0, mu0, sigmod0, alpha1, mu1, sigmod1
    ))
    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
    alpha0_new, mu0_new, sigmod0_new, alpha1_new, mu1_new, sigmod1_new
    ))

    png

    png

    png

    png

    png

    png

    alpha0:0.3, mu0:-2.0, sigmod0:0.5, alpha1:0.7, mu1:0.5, sigmod1:1.0
    alpha0:0.4, mu0:-1.6, sigmod0:0.8, alpha1:0.6, mu1:0.8, sigmod1:0.8
    1
    2
    if __name__ == '__main__':
    set_GM(mu,sigma)
    Donate
    • 微信
    • 支付寶

    Comment