EM算法及其应用

EM算法简介

首先上一段EM算法的wiki定义:

expectation–maximization (EM) algorithm is an iterative method to find maximum likelihood(MLE) or maximum a posteriori (MAP) estimates of parameters in statistical models, where the model depends on unobserved latent variables.

就是EM算法是: 一种迭代式的算法,用于含有隐变量的概率参数模型的最大似然估计或极大后验概率估计.

网上已经有很多很优秀的博客讲EM算法的了,再次就不赘述了,只复述一些关键性的步骤,相关链接见本文参考部分.

(1) 起因: 给定一系列样本,求解含有隐变量的极大似然估计(MLE)

其中z表示隐变量. 由于隐变量的存在,无法直接使用MLE去求解theta,EM的策略是先建立极大似然函数的下界(E-Step),然后去优化下界逼近原始的极大解(M-Step),不停迭代直到收敛到局部最优解.

(2) 求解: 算法推导

Qi表示隐变量z的分布,需要满足条件:

,比如要将班上学生聚类,假设隐藏变量z是身高,那么Qi就是连续的高斯分布,如果按照隐藏变量是男女,那么就是伯努利分布.

主要是公式2到公式3比较难懂,使用的是Jensen不等式,具体可以看这篇博客有详细的数学解释,此处不赘述.

(3) 结论: 算法总结

公式3表示是对极大似然函数求下界,此时我们假定theta已近给定,通过调整Qi的值使得下界不断的上升去逼近真实值. 当不等式变成等式的时候表示已经调整到和真实值一样的水平了,由Jensen不等式知道此时的随机变量是常数C,即:

进一步推导得到:

得到第一个重要的结论: theta已知的情况下,使得下界提升的Qi就是后验概率,解决了Qi如何选择的问题,其实这就是E-Step.

在找到使得下界提升的Qi之后,固定住Qi,M-Step就是使用MLE极大化此时的下界.

总结下就是:

套路就是: 首先猜下隐类别变量z,之后更新其它参数(theta)

图解就是:

当收敛到theta*时或者||theta(t+1)-theta(t)|| < thresh就可以迭代停止了. 从算法的过程来看,EM算法对初始值敏感同时不能保证收敛到全局最优解. 至于后续的证明EM算法的收敛性,大家看我参考处的相关博客链接或者李航博士的<<统计学习方法>>一书第9章有详细的证明.

EM算法的应用

GMM

GMM(Gaussian Mixture Model)就是指对样本的概率密度(density estimation)分布进行估计,而估计采用的模型是多个高斯模型的加权和,其中的每个高斯模型就代表了一个类(Cluster). 实际分布中可以把模型定义为任何分布的mixture model,为何是高斯混合模型呢? 原因如下两点:

  • 计算比较方便
  • 理论任意多的高斯分布可以近似任意概率分布

问题简化为:

随机变量X是由K个高斯分布混合而成,各个高斯分布的权重(概率)是Φi, 各个高斯分布的均值和方差为μi, ∑i. 观测到随机变量X的一系列样本,估计参数Φ, μ, ∑.

和EM算法之前的引入一样,隐含类别标签用Zi表示,表示样本属于类别Zi,可以假定Zi服从多项式分布,即:

换句话来说就是第j个cluster的权重是Φj.

假设有K个类别(cluster). 假定在给定Zi后,Xi服从高斯分布,即:

联合概率分布是:

故此时的极大似然函数是:

参考EM算法的套路,首先猜测隐类别变量z,然后更新其它参数(Φ, μ, ∑).

Wji表示第i个数据点属于第j个cluster的概率. 具体的Wji的计算可以使用贝叶斯公式:

sklearn中的GMM

API地址在这:GuassianMixture_API 官方的guide是这么介绍的:

The GaussianMixture object implements the expectation-maximization (EM) algorithm for fitting mixture-of-Gaussian models.

可以看出是用EM算法求解的GMM. 官方有个示例, 示例地址是使用EM算法来进行density estimation的.

代码直接粘贴来,如下:

import matplotlib as mpl
import matplotlib.pyplot as plt

import numpy as np

from sklearn import datasets
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import StratifiedKFold

print(__doc__)

colors = ['navy', 'turquoise', 'darkorange']


def make_ellipses(gmm, ax):
    for n, color in enumerate(colors):
        if gmm.covariance_type == 'full':
            covariances = gmm.covariances_[n][:2, :2]
        elif gmm.covariance_type == 'tied':
            covariances = gmm.covariances_[:2, :2]
        elif gmm.covariance_type == 'diag':
            covariances = np.diag(gmm.covariances_[n][:2])
        elif gmm.covariance_type == 'spherical':
            covariances = np.eye(gmm.means_.shape[1]) * gmm.covariances_[n]
        v, w = np.linalg.eigh(covariances)
        u = w[0] / np.linalg.norm(w[0])
        angle = np.arctan2(u[1], u[0])
        angle = 180 * angle / np.pi  # convert to degrees
        v = 2. * np.sqrt(2.) * np.sqrt(v)
        ell = mpl.patches.Ellipse(gmm.means_[n, :2], v[0], v[1],
                                  180 + angle, color=color)
        ell.set_clip_box(ax.bbox)
        ell.set_alpha(0.5)
        ax.add_artist(ell)

iris = datasets.load_iris()

# Break up the dataset into non-overlapping training (75%) and testing
# (25%) sets.
skf = StratifiedKFold(n_splits=4)
# Only take the first fold.
train_index, test_index = next(iter(skf.split(iris.data, iris.target)))


X_train = iris.data[train_index]
y_train = iris.target[train_index]
X_test = iris.data[test_index]
y_test = iris.target[test_index]

n_classes = len(np.unique(y_train))

# Try GMMs using different types of covariances.
estimators = dict((cov_type, GaussianMixture(n_components=n_classes,
                   covariance_type=cov_type, max_iter=20, random_state=0))
                  for cov_type in ['spherical', 'diag', 'tied', 'full'])

n_estimators = len(estimators)

plt.figure(figsize=(3 * n_estimators // 2, 6))
plt.subplots_adjust(bottom=.01, top=0.95, hspace=.15, wspace=.05,
                    left=.01, right=.99)


for index, (name, estimator) in enumerate(estimators.items()):
    # Since we have class labels for the training data, we can
    # initialize the GMM parameters in a supervised manner.
    estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)
                                    for i in range(n_classes)])

    # Train the other parameters using the EM algorithm.
    estimator.fit(X_train)

    h = plt.subplot(2, n_estimators // 2, index + 1)
    make_ellipses(estimator, h)

    for n, color in enumerate(colors):
        data = iris.data[iris.target == n]
        plt.scatter(data[:, 0], data[:, 1], s=0.8, color=color,
                    label=iris.target_names[n])
    # Plot the test data with crosses
    for n, color in enumerate(colors):
        data = X_test[y_test == n]
        plt.scatter(data[:, 0], data[:, 1], marker='x', color=color)

    y_train_pred = estimator.predict(X_train)
    train_accuracy = np.mean(y_train_pred.ravel() == y_train.ravel()) * 100
    plt.text(0.05, 0.9, 'Train accuracy: %.1f' % train_accuracy,
             transform=h.transAxes)

    y_test_pred = estimator.predict(X_test)
    test_accuracy = np.mean(y_test_pred.ravel() == y_test.ravel()) * 100
    plt.text(0.05, 0.8, 'Test accuracy: %.1f' % test_accuracy,
             transform=h.transAxes)

    plt.xticks(())
    plt.yticks(())
    plt.title(name)

plt.legend(scatterpoints=1, loc='lower right', prop=dict(size=12))


plt.show()

代码大意是使用不同的covariance类型({‘full’, ‘tied’, ‘diag’, ‘spherical’}),来观察GMM对iris数据集的聚类效果. iris数据集由150个样本组成,每个样本的特征是4维,3个类别(setosa,versicolor,virginica).

结果如下:

EM还有用在DGM(Bayesian network)中的,这些就比较高深了,暂时还没做了解,以后再补.

参考 1. EM算法在wiki上的解释 2. Jerry Lead的博客 3. zouxy09的博客 4. 一个EM算法的总结 5. GMM模型 6. sina博客介绍的GMM 7. scikit-learn中的GMM

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏数据科学学习手札

(数据科学学习手札09)系统聚类算法Python与R的比较

上一篇笔者以自己编写代码的方式实现了重心法下的系统聚类(又称层次聚类)算法,通过与Scipy和R中各自自带的系统聚类方法进行比较,显然这些权威的快捷方法更为高效...

2948
来自专栏云时之间

深度学习与神经网络:正则化小栗子(附代码)

在上一篇文章中我们简单说了说AIC,BIC和L1,L2正则化的东西,而今天这篇文章,我们将着重说说正则化.

3986
来自专栏林德熙的博客

卷积神经网络全面解析

本文转自: http://www.moonshile.com/post/juan-ji-shen-jing-wang-luo-quan-mian-jie-xi ...

1032
来自专栏人工智能LeadAI

主成分分析降维(MNIST数据集)

今天看了用主成分分析简化数据,就顺便用MNIST数据集做了下实验,想直观地看一下效果,并通过完成这个小demo深入理解下原理。 我发现“是什么、能做什么、怎么用...

3498
来自专栏数据派THU

一文读懂支持向量积核函数(附公式)

来源:jerrylead 本文通过多个例子为你介绍支持向量积核函数,助你更好地理解。 核函数(Kernels) 考虑我们最初在“线性回归”中提出的问题,特征是房...

47014
来自专栏人工智能

掌握机器学习数学基础之概率统计(三)

标题: 机器学习为什么要使用概率 概率学派和贝叶斯学派 何为随机变量和何又为概率分布? 条件概率,联合概率和全概率公式: 边缘概率 独立性和条件独立性 期望、方...

25510
来自专栏黄成甲

数据分析之因子分析

系统聚类分析可以对变量进行分类,但是难以判断变量分类结果的合理性。另外,如果要衡量每个变量对类别的贡献,也难以通过聚类分析来实现。这个时候就要采用因子分析来实现...

1154
来自专栏机器学习原理

我的机器学习概率论篇排列 组合古典概率联合概率条件概率全概率公式贝叶斯公式独立事件随机变量离散型随机变量连续型随机变量期望和方差三个基本定理参数估计

前言: 概率论的理解有些抽象,掌握概率论的方法,用实际样本去无限接近真实,熟练掌握并且使用一些最基本的概念是前提,比如,均值,方差 排列 组合 计算各种...

3516
来自专栏ATYUN订阅号

Python中的统计假设检验速查表

本文是一个机器学习项目中最流行的统计假设检验的速查表,包含使用Python接口的示例。

1406
来自专栏各种机器学习基础算法

K-NN算法与K-Means算法的原理与区别(附带源码示例)

KNN算法 K-Means算法 目标  确定某个元素所属的分类 将已存在的一系列元素分类 算法类别 监督的分类算法 无监督的聚类算法 ...

2877

扫码关注云+社区