6. EM算法-高斯混合模型GMM+Lasso详细代码实现

1. 前言

我们之前有介绍过4. EM算法-高斯混合模型GMM详细代码实现,在那片博文里面把GMM说涉及到的过程,可能会遇到的问题,基本讲了。今天我们升级下,主要一起解析下EM算法中GMM(搞事混合模型)带惩罚项的详细代码实现。

2. 原理

由于我们的极大似然公式加上了惩罚项,所以整个推算的过程在几个地方需要修改下。

在带penality的GMM中,我们假设协方差是一个对角矩阵,这样的话,我们计算高斯密度函数的时候,只需要把样本各个维度与对应的\(\mu_k\)和\(\sigma_k\)计算一维高斯分布,再相加即可。不需要通过多维高斯进行计算,也不需要协方差矩阵是半正定的要求。

我们给上面的(1)式加入一个惩罚项,
\[
\lambda\sum_{k=1}^K\sum_{j=1}^P\frac{|\mu_k-\bar{x}_j|}{s_j}
\]
其中的\(P\)是样本的维度。\(\bar{x}_j\)表示每个维度的平均值,\(s_j\)表示每个维度的标准差。这个penality是一个L1范式,对\(\mu_k\)进行约束。

加入penality后(1)变为
\[
L(\theta,\theta^{(j)})=\sum_{k=1}^Kn_k[log\pi_k-\frac{1}{2}(log(\boldsymbol{\Sigma_k})+\frac{{(x_i-\boldsymbol{\mu}_k})^2}{\boldsymbol{\Sigma}_k})] - \lambda\sum_{k=1}^K\sum_{j=1}^P\frac{|\mu_k-\bar{x}_j|}{s_j}
\]

这里需要注意的一点是,因为penality有一个绝对值,所以在对\(\mu_k\)求导的时候,需要分情况。于是(2)变成了
\[
\mu_k=\frac{1}{n_k}\sum_{i=1}^N\gamma_{ik}x_i
\]
\[
\mu_k=
\left \{\begin{array}{cc}
\frac{1}{n_k}(\sum_{i=1}^N\gamma_{ik}x_i - \frac{\lambda\sigma^2}{s_j}), & \mu_k >= \bar{x}_j\\frac{1}{n_k}(\sum_{i=1}^N\gamma_{ik}x_i + \frac{\lambda\sigma^2}{s_j}), & \mu_k < \bar{x}_j
\end{array}\right.
\]

3. 算法实现

  • 和不带惩罚项的GMM不同的是,我们GMM+LASSO的计算高斯密度函数有所变化。
#计算高斯密度概率函数,样本的高斯概率密度函数,其实就是每个一维mu,sigma的高斯的和
def log_prob(self, X, mu, sigma):
    N, D = X.shape
    logRes = np.zeros(N)
    for i in range(N):
        a = norm.logpdf(X[i,:], loc=mu, scale=sigma)
        logRes[i] = np.sum(a)
    return logRes
  • 在m-step中计算\(\mu_{k+1}\)的公式需要变化,先通过比较\(\mu_{kj}\)和\(means_{kj}\)的大小,来确定绝对值shift的符号。
def m_step(self, step):
    gammaNorm = np.array(np.sum(self.gamma, axis=0)).reshape(self.K, 1)
    self.alpha = gammaNorm / np.sum(gammaNorm)
    for k in range(self.K):
        Nk = gammaNorm[k]
        if Nk == 0:
            continue
        for j in range(self.D):
            if step >= self.beginPenaltyTime:
                # 算出penality的偏移量shift,通过当前维度的mu和样本均值比较,确定shift的符号,相当于把lasso的绝对值拆开了
                shift = np.square(self.sigma[k, j]) * self.penalty / (self.std[j] * Nk)
                if self.mu[k, j] >= self.means[j]:
                    shift = shift
                else:
                    shift = -shift
            else:
                shift = 0
            self.mu[k, j] = np.dot(self.gamma[:, k].T, self.X[:, j]) / Nk - shift
            self.sigma[k, j] = np.sqrt(np.sum(np.multiply(self.gamma[:, k], np.square(self.X[:, j] - self.mu[k, j]))) / Nk)
  • 最后需要修改loglikelihood的计算公式
def GMM_EM(self):
    self.init_paras()
    for i in range(self.times):
        #m step
        self.m_step(i)
        # e step
        logGammaNorm, self.gamma= self.e_step(self.X)
        #loglikelihood
        loglike = self.logLikelihood(logGammaNorm)
        #penalty
        pen = 0
        if i >= self.beginPenaltyTime:
            for j in range(self.D):
                pen += self.penalty * np.sum(abs(self.mu[:,j] - self.means[j])) / self.std[j]

        # print("step = %s, alpha = %s, loglike = %s"%(i, [round(p[0], 5) for p in self.alpha.tolist()], round(loglike - pen, 5)))
        # if abs(self.loglike - loglike) < self.tol:
        #     break
        # else:

        self.loglike = loglike - pen

4. GMM算法实现结果

用我实现的GMM+LASSO算法,对多个penality进行计算,选出loglikelihood最大的k和penality,与sklearn的结果比较。

fileName = amix1-est.dat, k = 2, penalty = 0 alpha = [0.52838, 0.47162], loglike = -693.34677
fileName = amix1-est.dat, k = 2, penalty = 0 alpha = [0.52838, 0.47162], loglike = -693.34677
fileName = amix1-est.dat, k = 2, penalty = 1 alpha = [0.52789, 0.47211], loglike = -695.26835
fileName = amix1-est.dat, k = 2, penalty = 1 alpha = [0.52789, 0.47211], loglike = -695.26835
fileName = amix1-est.dat, k = 2, penalty = 2 alpha = [0.52736, 0.47264], loglike = -697.17009
fileName = amix1-est.dat, k = 2, penalty = 2 alpha = [0.52736, 0.47264], loglike = -697.17009
myself GMM alpha = [0.52838, 0.47162], loglikelihood = -693.34677, bestP = 0
sklearn GMM alpha = [0.53372, 0.46628], loglikelihood = -176.73112
succ = 299/300
succ = 0.9966666666666667
[0 1 0 0 1 1 0 1 1 1 0 0 1 0 0 1 0 0 0 1]
[0 1 0 0 1 0 0 1 1 1 0 0 1 0 0 1 0 0 0 1]
fileName = amix1-tst.dat, loglike = -2389.1852339407087
fileName = amix1-val.dat, loglike = -358.1157431278091
fileName = amix2-est.dat, k = 2, penalty = 0 alpha = [0.56, 0.44], loglike = 53804.54265
fileName = amix2-est.dat, k = 2, penalty = 0 alpha = [0.82, 0.18], loglike = 24902.5522
fileName = amix2-est.dat, k = 2, penalty = 1 alpha = [0.82, 0.18], loglike = 23902.65183
fileName = amix2-est.dat, k = 2, penalty = 1 alpha = [0.56, 0.44], loglike = 52929.96459
fileName = amix2-est.dat, k = 2, penalty = 2 alpha = [0.82, 0.18], loglike = 22907.40397
fileName = amix2-est.dat, k = 2, penalty = 2 alpha = [0.82, 0.18], loglike = 22907.40397
myself GMM alpha = [0.56, 0.44], loglikelihood = 53804.54265, bestP = 0
sklearn GMM alpha = [0.56217, 0.43783], loglikelihood = 11738677.90164
succ = 200/200
succ = 1.0
[0 1 0 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 0 1]
[0 1 0 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 0 1]
fileName = amix2-tst.dat, loglike = 51502.878096147084
fileName = amix2-val.dat, loglike = 6071.217012747491
fileName = golub-est.dat, k = 2, penalty = 0 alpha = [0.575, 0.425], loglike = -24790.19895
fileName = golub-est.dat, k = 2, penalty = 0 alpha = [0.525, 0.475], loglike = -24440.82743
fileName = golub-est.dat, k = 2, penalty = 1 alpha = [0.55, 0.45], loglike = -25582.27485
fileName = golub-est.dat, k = 2, penalty = 1 alpha = [0.6, 0.4], loglike = -26137.97508
fileName = golub-est.dat, k = 2, penalty = 2 alpha = [0.55, 0.45], loglike = -26686.02411
fileName = golub-est.dat, k = 2, penalty = 2 alpha = [0.55, 0.45], loglike = -26941.68964
myself GMM alpha = [0.525, 0.475], loglikelihood = -24440.82743, bestP = 0
sklearn GMM alpha = [0.5119, 0.4881], loglikelihood = 13627728.10766
succ = 29/40
succ = 0.725
[0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 1]
[0 1 0 1 1 1 0 0 1 1 0 1 0 1 0 0 1 1 0 0]
fileName = golub-tst.dat, loglike = -12949.606698037718
fileName = golub-val.dat, loglike = -11131.35137056415

5. 总结

通过一番改造,实现了GMM+LASSO的代码,如果读者有什么好的改进方法,或者我有什么错误的地方,希望多多指教。

原文地址:https://www.cnblogs.com/huangyc/p/10279240.html

时间: 2024-10-25 06:24:02

6. EM算法-高斯混合模型GMM+Lasso详细代码实现的相关文章

贝叶斯来理解高斯混合模型GMM

最近学习基础算法<统计学习方法>,看到利用EM算法估计高斯混合模型(GMM)的时候,发现利用贝叶斯的来理解高斯混合模型的应用其实非常合适. 首先,假设我们对于贝叶斯比较熟悉,对高斯分布也熟悉.本文将GMM用于聚类来举例. 除了简单的高斯分布,理论上通过组合多个不同的高斯分布可以构成任意复杂的分布函数.如下图所示: 在最大似然,贝叶斯方法与朴素贝叶斯分类中,2.1中提到高斯概率密度用来计算连续变量情况下的朴素贝叶斯概率.该情况下的高斯分布是训练已知,然后对于输入变量求取其概率密度,结合类别的先验

基本算法思想Java实现的详细代码

基本算法思想Java实现的详细代码 算法是一个程序的灵魂,一个好的算法往往可以化繁为简,高效的求解问题.在程序设计中算法是独立于语言的,无论使用哪一种语言都可以使用这些算法,本文笔者将以Java语言为例介绍一些常用的算法思想. 分类 穷举算法思想 递推算法思想 递归算法思想 分治算法思想 概率算法思想  穷举算法思想 穷举算法的基本思想 从所有可能情况中搜索正确答案 1. 对于一种可能情况,计算其结果. 2. 判断结果是否满足,如不能满足者执行第一步来搜索下一个可能的情况:如满足则表示选找到一个

高斯混合模型(GMM)

1. 有时候单一高斯分布不能很好的描述分布 image.png 上图左面用单一高斯分布去描述,显然没有右图用两个高斯分布去描述的效果好. 2. 引入混合高斯分 这里插一句,为什么是“高斯混合模型”,而不是别的混合模型,因为从中心极限定理知,只要K足够大,模型足够复杂,样本量足够多,每一块小区域就可以用高斯分布描述.而且高斯函数具有良好的计算性能,所GMM被广泛地应用. 单一高斯分布公式 image.png 混合高斯分布 每个GMM由K个高斯分布组成,每个高斯分布称为一个组件(Component)

高斯混合模型GMM

1.简介 众所周知,GMM(Gaussian Mixture-based Model)是用来分离场景中前景和背景的,或者叫做背景扣除,那么什么叫做背景扣除(Background Subtraction)呢?我们知道所谓的监控系统中,通常都是利用静态相机来捕捉场景的,因此其中比较具有挑战的一步就是如何检测出场景中的突然闯入者,传统的应用中都会假设场景中没有这样的闯入者,而在实际监控场景中,这种情况却是实实在在存在发生的,所以我们不能忽略.对于闯入者的检测,如果我们知道一个场景的统计模型,那些不服从

EM算法原理以及高斯混合模型实践

EM算法有很多的应用: 最广泛的就是GMM混合高斯模型.聚类.HMM等等. The EM Algorithm 高斯混合模型(Mixtures of Gaussians)和EM算法 EM算法 求最大似然函数估计值的一般步骤: (1)写出似然函数: (2)对似然函数取对数,并整理: (3)求导数,令导数为0,得到似然方程: (4)解似然方程,得到的参数即为所求. 期望最大化算法(EM算法): 优点: 1. 简单稳定: 2. 通过E步骤和M步骤使得期望最大化,是自收敛的分类算法,既不需要事先设定类别也

机器学习笔记(十)EM算法及实践(以混合高斯模型(GMM)为例来次完整的EM)

今天要来讨论的是EM算法.第一眼看到EM我就想到了我大枫哥,EM Master,千里马,RUA!!!不知道看这个博客的人有没有懂这个梗的.好的,言归正传,今天要讲的EM算法,全称是Expectation maximization,期望最大化.怎么个意思呢,就是给你一堆观测样本,让你给出这个模型的参数估计.我靠,这套路我们前面讨论各种回归的时候不是已经用烂了吗?求期望,求对数期望,求导为0,得到参数估计值,这套路我懂啊,MLE!但问题在于,如果这个问题存在中间的隐变量呢?会不会把我们的套路给带崩呢

EM算法 - 2 - EM算法在高斯混合模型学习中的应用

声明: 1,本篇为个人对<2012.李航.统计学习方法.pdf>的学习总结,不得用作商用,欢迎转载,但请注明出处(即:本帖地址). 2,由于本人在学习初始时有很多数学知识都已忘记,所以为了弄懂其中的内容查阅了很多资料,所以里面应该会有引用其他帖子的小部分内容,如果原作者看到可以私信我,我会将您的帖子的地址付到下面. 3,如果有内容错误或不准确欢迎大家指正. 4,如果能帮到你,那真是太好了. 在开始讲解之前,我要先给看这篇文章的你道个歉,因为<2012.李航.统计学习方法.pdf>中

GMM学习笔记(EM算法求解)

提出混合模型主要是为了能更好地近似一些较复杂的样本分布,通过不断增加component个数,可以任意地逼近任何连续的概率分布,所以我们认为任何样本分布都可以用混合模型来建模.因为高斯函数具有一些很实用的性质,所以高斯混合模型被广泛地使用. GMM与kmeans类似,也是属于clustering,不同的是,kmeans是把每个样本点聚到其中一个cluster,而GMM是给出这些样本点到每个cluster的概率,每个component就是一个聚类中心. GMM(Gaussian Mixture Mo

K-Means(K均值)、GMM(高斯混合模型),通俗易懂,先收藏了!

1. 聚类算法都是无监督学习吗? 什么是聚类算法?聚类是一种机器学习技术,它涉及到数据点的分组.给定一组数据点,我们可以使用聚类算法将每个数据点划分为一个特定的组.理论上,同一组中的数据点应该具有相似的属性和/或特征,而不同组中的数据点应该具有高度不同的属性和/或特征.聚类是一种无监督学习的方法,是许多领域中常用的统计数据分析技术. 常用的算法包括K-MEANS.高斯混合模型(Gaussian Mixed Model,GMM).自组织映射神经网络(Self-Organizing Map,SOM)