机器学习高斯混合模型(后篇):GMM求解完整代码实现

《实例》阐述算法,通俗易懂,助您对算法的理解达到一个新高度。包含但不限于:经典算法,机器学习,深度学习,LeetCode 题解,Kaggle 实战。期待您的到来!

01

回顾

前面推送中,我们介绍了高斯混合模型(GMM)的聚类原理,以及聚类求解的公式推导,如果您想了解这部分,请参考之前的推送:

机器学习高斯混合模型:聚类原理分析(前篇)

机器学习高斯混合模型(中篇):聚类求解

总结来说,GMM是非常好的聚类利器,它不光能给出样本所属的类别,还能给出属于每个类别的概率,进而转化成得分值,有时所属每个簇的得分值具有重要的意义(意义说明详见之前两篇的推送)。GMM求解的思路本质上是借助最大期望算法的思路来进行求解,关于最大期望算法的原理例子解析,请参考:

机器学习期望最大算法:实例解析

接下来,就到了GMM的EM算法求解的代码实现环节了,如果我们能把一种聚类算法的思路从原理,到公式,再到代码实现,都能走一遍,那么无疑可以表明您对本算法和这一类的算法都有一个全新的理解。手写不掉包代码实现算法的结果,如果能与sklean中的实现基本一致,那么说明才说明您对这个算法正真了解了,在这个编码的过程,将是您对python,Numpy等常用科学计算工具的实践过程,总之意义挺大,锻炼价值也很大。

废话少说,让我们开始GMM模型的EM算法的代码实现之旅吧!

02

GMM的EM求解之数据生成

我们先从一维的数据样本的聚类开始说起,先易后难。首先阐述下GMM的EM求解思路。

1. 数据准备

我们借助sklearn的API,生成3堆一维高斯分布的数据,一维在此处是指数据的特征只有一个。首先本次实验导入的所有库包括:

import numpy as np

import numpy.linalg as la

import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs

生成数据的过程如下:

#生成的簇,和对应的分类

#这是sklearn的聚类结果

#下面自己编码GMM实现聚类,看看与sklearn的结果是够一致

x,label = make_blobs(n_samples=150,n_features=1, centers=3,

cluster_std=[0.01,0.02,0.03],

center_box=(0.1,0.4, 2.0),

random_state=2)

plt.scatter(np.arange(0,150),x[:, 0], marker='o', c=label)

plt.title('GMM classification')

plt.xlabel('point Id')

plt.ylabel('x1 attribute')

plt.show()

sklearn生成的满足高斯分布的3簇:

那么,我们一直这些样本点,如何进行正确的聚类呢?也就是能聚类出和上图差不多的效果来。

03

EM求解代码解析

1 初始化参数

需要初始化的参数包括:

  • 每个簇的均值,数组的形状参考注释(K by D的意思是K行D列);
  • 协方差(这个需要特别注意,一维高斯是方差,二维以上是协方差,形状也需要特别注意:D by D by K);
  • 每个簇的影响系数

#初始化参数

def initParams(K,D):

#每个簇的中心值:K by D

aves = np.random.rand(K,D)

#每个簇的偏差 D by D by K

sigmas = np.zeros((D,D,K))

###D by D 必须是对称矩阵

sig = np.eye(D)

for k in np.arange(0,K):

sigmas[:,:,k] = sig

#每个簇的影响系数:1 by K

pPis = np.random.rand(1,K)

return aves,sigmas,pPis

2 样本点对GMM的贡献系数求解

求解的公式如下,关于这个公式的具体含义,请参考本文开头列出的推送文章。

#样本点对簇的贡献系数

#pPi : 1 by K

#px: N by K

# return value: N by K

def fgamma(px,pPi):

#gamma公式的分子部分

#fenzi: N by K

fenzi = pPi * px

#gamma公式的分母部分

#fenmu: N by 1

fenmu = np.sum(fenzi,axis=1).reshape(N,1)

return fenzi/fenmu

一点说明:

在用Numpy求解时,数组的运算可以省掉C++,Java等的很多for循环,可以看求解上面这个公式只需要短短3行代码,可以说说很简洁,但是对于以前使用Java,C++的小伙伴,上手Numpy需要做一个思维转化,同时也要注意标注每个数组的shape,这对于我们后续检查bug非常重要。

3 簇中的样本点的贡献和

从第2步中得出的每个样本点的贡献,然后累加即可:

# 每个簇中的样本点的贡献系数之和

# gam: N by K

# return value: 1 by K

def fNk(gam):

nk = np.sum(gam,axis=0)

return nk.reshape(1,K)

上面相当于EM算法的E步,下面总结M步,是利用最大似然估计各个簇的分布参数。

4 每个簇的均值和协方差求解

每个簇的样本和协方差的求解公式如下:

#每个簇的均值

# Nk: 1 by K

# gam: N by K

# x : N by D

#return value: K by D

def faverage(aves,Nk,gam,x):

#print(np.shape(Nk))

for k in np.arange(0,K):

# sum : D

sumd = np.sum((gam[:,k].reshape(N,1)) * x,axis=0)

aves[k,:] = sumd.reshape(1,D)/Nk[:,k]

#每个簇的方差

# Nk: 1 by K

# gam: N by K

# x : N by D

# aves: K by D

#return value: D by D by K

def fsigma(sigmas,Nk,gam,x,aves):

for k in np.arange(0,K):

#shift: N by DA

shift = x - aves[k,:]

#shift_gam: N by D

shift_gam = gam[:,k].reshape(N,1)*shift

#shift2 : D by D

shift2 = shift_gam.T.dot(shift)

sigmas[:,:,k] = shift2/Nk[:,k]

return sigmas

5 多维高斯分布的概率密度公式求解

多维高斯分布的概率密度公式见下,式子中 d 表示维数(也就是特征个数),求和符号指:协方差(二维及以上是个方阵)

# D-dimension prob density

# x : N by D

# aves : K by D

# sigmas: D by D by K

# return value: N by K

def fpx(x,aves,sigmas):

Px = np.zeros((N,K))

# coef1 : 1 by 1

coef1 = np.power((2*np.pi),(D/2.0))

for k in np.arange(0,K):

# coef2 : 1 by 1

coef2 = np.power((la.det(sigmas[:,:,k])),0.5)

coef3 = 1/(coef1 * coef2)

# shift: N by D

shift = x - aves[k,:]

# sigmaInv: D by D

sigmaInv = la.inv(sigmas[:,:,k])

epow = -0.5*(shift.dot(sigmaInv)*shift)

# epowsum : N

epowsum = np.sum(epow,axis=1)

Px[:,k] = coef3 * np.exp(epowsum)

return Px

6 迭代停止策略

各个样本点的最大似然估计值趋于稳定(小于某个阈值:比如:1e-15),最大似然估计的公式如下:

#迭代求解的停止策略

#px: N by K

#pPi: 1 by K

#Loss function 1 by 1

def fL(px, pPi):

# sub: N by 1

sub = np.sum(pPi*px,axis=1)

logsub = np.log(sub)

curL = np.sum(logsub)

return curL

# stop iterator strategy

def stop_iter(threshold,preL,curL):

return np.abs(curL-preL) < threshold

04

GMM聚类接口编写

有了以上EM算法的各个函数后,下面可以编写GMM聚类的对外接口了。

# GMM

# return value: N by K

def GMM(x,K):

#loss value initilize

preL = -np.inf;

# aves 每个簇的中心值:K by D

# sigmas 每个簇的偏差 D by D by K

# pPi 每个簇的影响系数:1 by K

aves,sigmas,pPi = initParams(K,D)

while True:

# px: 每个数据所属簇的概率 N by K

px = fpx(x,aves,sigmas)

#print(px)

# 贡献系数 N by K

gam = fgamma(px,pPi)

#每个簇中的样本点的贡献系数之和 1 by K

Nk = fNk(gam)

pPi = Nk/N

# 每个簇的均值 K by D

faverage(aves,Nk,gam,x)

#每个簇的方差 D by D by K

fsigma(sigmas,Nk,gam,x,aves)

# loss function

curL = fL(px, pPi)

#迭代求解的停止策略

if stop_iter(1e-15,preL, curL):

break

preL = curL

return px,aves,sigmas

#返回聚类的结果:N

def classifior(px):

rslt = []

for row in px:

rslt.append(np.where(row==np.max(row)))

return np.array(rslt).reshape(-1)

05

模拟一维高斯分布的聚类

# K: 簇的个数

# D: 数据的维数(特征数或属性数)

# N:样本点个数

K = 3

D = 1

N = 150

#一维特征的GMM聚类模拟

px,aves,sigmas =GMM(x,3)

mylabel = classifior(px)

#可以看到不掉包的实现与sklearn的模拟结果是基本一致的

plt.scatter(np.arange(0,150),x[:, 0], marker='o', c=mylabel)

plt.title('GMM classification')

plt.xlabel('x1 attribute')

plt.ylabel('x2 attribute')

plt.show()

可以看到不掉包的实现与sklearn的掉包实现结果是基本一致的。

一维高斯分布的协方差是方差,是一个数。虽然以上算法能实现多维的高斯分布的聚类,但是鉴于篇幅,明天推送关于多维的高斯分布的聚类的结果展示,协方差,概率密度图等都有着非常重要的应用,并且它们也是非常有意思的。

谢谢您的阅读!

辅助材料:

机器学习储备(13):概率密度和高斯分布例子解析

原文发布于微信公众号 - 算法channel(alg-channel)

原文发表时间:2017-12-03

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏新智元

计算机视觉识别简史:从 AlexNet、ResNet 到 Mask RCNN

【新智元导读】 Medium 用户 Đặng Hà Thế Hiển 制作了一张信息图示,用专业、简洁并且最有吸引力的方式——信息图示,讲述计算机视觉(CV)物...

48480
来自专栏Pytorch实践

Pytorch实现LSTM时间序列预测

摘要:本文主要基于Pytorch深度学习框架,实现LSTM神经网络模型,用于时间序列的预测。 开发环境说明: Python 35 Pytorch 0.2 CP...

1.6K70
来自专栏智能算法

主成分分析到底怎么分析?

PCA(Principal Component Analysis)是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提...

351100
来自专栏机器学习算法与Python学习

推荐 | 图解机器学习

关键字全网搜索最新排名 【机器学习算法】:排名第一 【机器学习】:排名第一 【Python】:排名第三 【算法】:排名第四 源 | ynaughty 每当提...

43050
来自专栏智能算法

降维方法(一):PCA原理

PCA(Principal Component Analysis)是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提...

44490
来自专栏数据派THU

一文图解机器学习的基本算法!

来源:软件定义世界 本文长度为2877字,建议阅读6分钟 本文为你分析如何选择机器学习的各种方法。 每当提到机器学习,大家总是被其中的各种各样的算法和方法搞晕,...

25450
来自专栏数据派THU

【一图看懂】计算机视觉识别简史:从 AlexNet、ResNet 到 Mask RCNN

原文:medium 来源:新智元 作者:Đặng Hà Thế Hiển 编译:新智元编辑部 本文长度为5000字,建议阅读8分钟 本文通过一张信息图示,讲述计...

36070
来自专栏刘笑江的专栏

HTML

21130
来自专栏机器学习算法与理论

《白话深度学习与Tensorflow》学习笔记(7)RBM限制玻尔兹曼机

不受限的:在实际工程中用处不大,所以目前只需要研究受限玻尔兹曼机。 一层是可视化层,一般来说是输入层,另一层是隐含层,也就是我们一般指的特征提取层。 RBM是可...

416110
来自专栏林欣哲

深度学习的训练

今天来聊聊深度学习的训练方法和注意事项 数据集的拆分: 首先要准备好已经处理好的数据集(注意数据集要满足独立同分布),分为训练集、验证集、测试集。可按80%,1...

48380

扫码关注云+社区

领取腾讯云代金券