机器学习（一） 极大似然估计与EM算法

最近在复习机器学习的知识，正好每天整理一点。

这些都是根据网上的资料整理而成，仅供自己学习，如有侵权，请联系我删除，谢谢！

一.极大似然估计

1.1概念

极大似然估计实际上就是根据已知的样本结果信息，反推最具有可能导致这结果出现（最大概率）的模型参数值。

换句话说，即模型已定，参数未知。我们要去推导参数。

由一个实际问题理解：

假如我们需要调查学校的男生和女生的身高分布 ，我们抽取100个男生和100个女生，将他们按照性别划分为两组。然后，统计抽样得到100个男生的身高数据和100个女生的身高数据。如果我们知道他们的身高服从正态分布，但是这个分布的均值  和方差  是不知道，这两个参数就是我们需要估计的。

我们可以把它转化成数学模型：

问题数学化：样本集  。概率密度是： 抽到第i个男生身高的概率。由于100个样本之间独立同分布，所以同时抽到这100个男生的概率是它们各自概率的乘积，就是从分布是p(X|θ)的总体样本中抽取到这100个样本的概率，也就是样本集X中各个样本的联合概率，用下式表示：

这就反应了在参数取Θ时，出现这些样本数据的概率。

所以我们需要找到一个参数θ，使得抽到X这组样本的概率最大，也就是说需要其对应的似然函数L(θ)最大。满足条件的θ叫做θ的最大似然估计值，记为：

 

1.2最大似然函数估计值的求解步骤

• 首先，写出似然函数：

• 其次，对似然函数取对数：

• 然后，对上式求导，另导数为0，得到似然方程。
• 最后，解似然方程，得到的参数值即为所求。

多数情况下，我们是根据已知条件来推算结果，而极大似然估计是已知结果，寻求使该结果出现的可能性最大的条件，以此作为估计值。

 

凸函数定义:设f是定义域为实数的函数，如果对所有的实数x，f(x)的二阶导数都大于0，那么f是凸函数。

Jensen不等式定义如下：

如果f是凸函数，X是随机变量，那么：  。当且仅当X是常量时，该式取等号。其中，E(X)表示X的数学期望。

注：Jensen不等式应用于凹函数时，不等号方向反向。当且仅当x是常量时，该不等式取等号。

 

二.EM算法

2.1问题描述

我们目前有100个男生和100个女生的身高，但是我们不知道这200个数据中哪个是男生的身高，哪个是女生的身高，即抽取得到的每个样本都不知道是从哪个分布中抽取的。这个时候，对于每个样本，就有两个未知量需要估计：

（1）这个身高数据是来自于男生数据集合还是来自于女生？

（2）男生、女生身高数据集的正态分布的参数分别是多少？

这个问题用EM算法求解，则可分为以下步骤:

（1）初始化参数：先初始化男生身高的正态分布的参数：如均值=1.65，方差=0.15

（2）计算每一个人更可能属于男生分布或者女生分布；

（3）通过分为男生的n个人来重新估计男生身高分布的参数（最大似然估计），女生分布也按照相同的方式估计出来，更新分布。

（4）这时候两个分布的概率也变了，然后重复步骤（1）至（3），直到参数不发生变化为止。

 

2.2证明：

暂且略过，可看文末的参考链接

2.3 EM算法流程:

输入：观察到的数据，联合分布  ，条件分布  ，最大迭代次数J。

算法步骤：

（1）随机初始化模型参数θ的初值  。

（2）j=1,2,...,J 开始EM算法迭代：

• E步：计算联合分布的条件概率期望：

• M步：极大化  ,得到  :

• 如果 已经收敛，则算法结束。否则继续进行E步和M步进行迭代。

输出：模型参数θ。

代码示例：

# -*- coding: utf-8 -*-

import numpy as np
import math
import copy
import matplotlib.pyplot as plt

isdebug = True


# 指定k个高斯分布参数，这里指定k=2。注意2个高斯分布具有相同均方差Sigma，均值分别为Mu1,Mu2。
def init_data(Sigma, Mu1, Mu2, k, N):
    global X
    global Mu
    global Expectations
    X = np.zeros((1, N))
    Mu = np.random.random(k)
    Expectations = np.zeros((N, k))
    for i in range(0, N):
        if np.random.random(1) > 0.5:
            X[0, i] = np.random.normal(Mu1, Sigma)
        else:
            X[0, i] = np.random.normal(Mu2, Sigma)
    if isdebug:
        print("***********")
        print("初始观测数据X：")
        print(X)


# EM算法：步骤1，计算E[zij]
def e_step(Sigma, k, N):
    global Expectations
    global Mu
    global X
    for i in range(0, N):
        Denom = 0
        Numer = [0.0] * k
        for j in range(0, k):
            Numer[j] = math.exp((-1 / (2 * (float(Sigma ** 2)))) * (float(X[0, i] - Mu[j])) ** 2)
            Denom += Numer[j]
        for j in range(0, k):
            Expectations[i, j] = Numer[j] / Denom
    if isdebug:
        print("***********")
        print("隐藏变量E（Z）：")
        print(Expectations)


# EM算法：步骤2，求最大化E[zij]的参数Mu
def m_step(k, N):
    global Expectations
    global X
    for j in range(0, k):
        Numer = 0
        Denom = 0
        for i in range(0, N):
            Numer += Expectations[i, j] * X[0, i]
            Denom += Expectations[i, j]
        Mu[j] = Numer / Denom


# 算法迭代iter_num次，或达到精度Epsilon停止迭代
def run(Sigma, Mu1, Mu2, k, N, iter_num, Epsilon):
    init_data(Sigma, Mu1, Mu2, k, N)
    print("初始<u1,u2>:", Mu)
    for i in range(iter_num):
        Old_Mu = copy.deepcopy(Mu)
        e_step(Sigma, k, N)
        m_step(k, N)
        print(i, Mu)
        if sum(abs(Mu - Old_Mu)) < Epsilon:
            break


if __name__ == '__main__':
    sigma = 6  # 高斯分布具有相同的方差
    mu1 = 40  # 第一个高斯分布的均值 用于产生样本
    mu2 = 20  # 第二个高斯分布的均值 用于产生样本
    k = 2  # 高斯分布的个数
    N = 1000  # 样本个数
    iter_num = 1000  # 最大迭代次数
    epsilon = 0.0001  # 当两次误差小于这个时退出
    run(sigma, mu1, mu2, k, N, iter_num, epsilon)

    plt.hist(X[0, :], 50)
    plt.show()

 

 

 

 

参考链接:

EM算法详解 - 知乎 (zhihu.com)

(142条消息) 机器学习——EM算法及代码实现_程旭员的博客-CSDN博客_em算法代码