【machine learning】GMM算法(Python版)

本文参考CSDN大神的博文,并在讲述中引入自己的理解,纯粹理清思路,并将代码改为了Python版本。(在更改的过程中,一方面理清自己对GMM的理解,一方面学习了numpy的应用,不过也许是Python粉指数超标才觉得有必要改(⊙o⊙))

一、GMM模型

事实上,GMM 和 k-means 很像,不过 GMM 是学习出一些概率密度函数来(所以 GMM 除了用在 clustering 上之外,还经常被用于 density estimation ),简单地说,k-means 的结果是每个数据点被 assign 到其中某一个 cluster 了,而 GMM 则给出这些数据点被 assign 到每个 cluster 的概率,又称作 soft assignment 。

得出一个概率有很多好处,因为它的信息量比简单的一个结果要多,比如,我可以把这个概率转换为一个 score ,表示算法对自己得出的这个结果的把握,参考pluskid大神博文

  • 也许我可以对同一个任务,用多个方法得到结果,最后选取“把握”最大的那个结果;
  • 另一个很常见的方法是在诸如疾病诊断之类的场所,机器对于那些很容易分辨的情况(患病或者不患病的概率很高)可以自动区分,而对于那种很难分辨的情况,比如,49% 的概率患病,51% 的概率正常,如果仅仅简单地使用 50% 的阈值将患者诊断为“正常”的话,风险是非常大的,因此,在机器对自己的结果把握很小的情况下,会“拒绝发表评论”,而把这个任务留给有经验的医生去解决。

每个GMM由K个Gaussian分布组成,每个Gaussian称为一个“Component”,这些Component 线性加成在一起就组成了GMM 的概率密度函数:

根据上面的式子,如果我们要从 GMM 的分布中随机地取一个点的话,实际上可以分为两步:首先随机地在这 K个Gaussian Component 之中选一个,每个 Component 被选中的概率实际上就是它的系数 pi(k) ,选中了 Component 之后,再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布,转化为了已知的问题。

那么如何用 GMM 来做 clustering 呢?其实很简单,现在我们有了数据,假定它们是由 GMM 生成出来的,那么我们只要根据数据推出 GMM 的概率分布来就可以了,然后 GMM 的 K 个 Component 实际上就对应了 K 个 cluster 了。根据数据来推算概率密度通常被称作 density estimation ,特别地,当我们在已知(或假定)了概率密度函数的形式,而要估计其中的参数的过程被称作“参数估计”。

二、参数与似然函数

现在假设我们有 N 个数据点,并假设它们服从某个分布(记作 p(x)),现在要确定里面的一些参数的值,例如,在 GMM 中,我们就需要确定 影响因子pi(k)、各类均值pMiu(k) 和 各类协方差pSigma(k) 这些参数。 我们的想法是,找到这样一组参数,它所确定的概率分布生成这些给定的数据点的概率最大,而这个概率实际上就等于 ,我们把这个乘积称作似然函数 (Likelihood Function)。通常单个点的概率都很小,许多很小的数字相乘起来在计算机里很容易造成浮点数下溢,因此我们通常会对其取对数,把乘积变成加和∑Ni=1logp(xi),得到log-likelihood function 。接下来我们只要将这个函数最大化(通常的做法是求导并令导数等于零,然后解方程),亦即找到这样一组参数值,它让似然函数取得最大值,我们就认为这是最合适的参数,这样就完成了参数估计的过程。

下面让我们来看一看 GMM 的 log-likelihood function :

由于在对数函数里面又有加和,我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题,我们采取之前从 GMM 中随机选点的办法:分成两步,实际上也就类似于K-means 的两步。

三、算法流程

1、估计数据由每个 Component 生成的概率(并不是每个 Component 被选中的概率):对于每个数据 xi 来说,它由第 k 个 Component 生成的概率为

其中N(xi|μk,Σk)就是后验概率

2、通过极大似然估计可以通过求到令参数=0得到参数pMiu,pSigma的值。

其中 Nk=∑Ni=1γ(i,k) ,并且 πk也顺理成章地可以估计为 Nk/N(这里的顺理成章要考证起来有很多要说)

3、重复迭代前面两步,直到似然函数的值收敛为止。

声明:这里完全可以对照EM算法的两步走,对应关系如下:

  • 均值和方差对应θ
  • Component 生成的概率为隐藏变量
  • 最大似然函数L(θ)=

四、GMM聚类代码

GMM.py

#! /usr/bin/env python
#coding=utf-8

from numpy import *
import pylab
import random,math

def loadDataSet(fileName):      #general function to parse tab -delimited floats
    dataMat = []                #assume last column is target value
    fr = open(fileName)
    for line in fr.readlines():
        curLine = line.strip().split(‘\t‘)
        fltLine = map(float,curLine) #map all elements to float()
        dataMat.append(fltLine)
    return dataMat

def gmm(file, K_or_centroids):
# ============================================================
# Expectation-Maximization iteration implementation of
# Gaussian Mixture Model.
#
# PX = GMM(X, K_OR_CENTROIDS)
# [PX MODEL] = GMM(X, K_OR_CENTROIDS)
#
#  - X: N-by-D data matrix.
#  - K_OR_CENTROIDS: either K indicating the number of
#       components or a K-by-D matrix indicating the
#       choosing of the initial K centroids.
#
#  - PX: N-by-K matrix indicating the probability of each
#       component generating each point.
#  - MODEL: a structure containing the parameters for a GMM:
#       MODEL.Miu: a K-by-D matrix.
#       MODEL.Sigma: a D-by-D-by-K matrix.
#       MODEL.Pi: a 1-by-K vector.
# ============================================================
    ## Generate Initial Centroids
    threshold = 1e-15
    dataMat = mat(loadDataSet(file))
    [N, D] = shape(dataMat)
    K_or_centroids = 2
    # K_or_centroids可以是一个整数,也可以是k个质心的二维列向量
    if shape(K_or_centroids)==(): #if K_or_centroid is a 1*1 number
        K = K_or_centroids
        Rn_index = range(N)
        random.shuffle(Rn_index) #random index N samples
        centroids = dataMat[Rn_index[0:K], :]; #generate K random centroid
    else: # K_or_centroid is a initial K centroid
        K = size(K_or_centroids)[0];
        centroids = K_or_centroids;  

    ## initial values
    [pMiu,pPi,pSigma] = init_params(dataMat,centroids,K,N,D)
    Lprev = -inf #上一次聚类的误差  

    # EM Algorithm
    while True:
        # Estimation Step
        Px = calc_prob(pMiu,pSigma,dataMat,K,N,D)

        # new value for pGamma(N*k), pGamma(i,k) = Xi由第k个Gaussian生成的概率
        # 或者说xi中有pGamma(i,k)是由第k个Gaussian生成的
        pGamma = mat(array(Px) * array(tile(pPi, (N, 1))))  #分子 = pi(k) * N(xi | pMiu(k), pSigma(k))
        pGamma = pGamma / tile(sum(pGamma, 1), (1, K)) #分母 = pi(j) * N(xi | pMiu(j), pSigma(j))对所有j求和  

        ## Maximization Step - through Maximize likelihood Estimation
        #print ‘dtypeddddddddd:‘,pGamma.dtype
        Nk = sum(pGamma, 0) #Nk(1*k) = 第k个高斯生成每个样本的概率的和,所有Nk的总和为N。  

        # update pMiu
        pMiu = mat(diag((1/Nk).tolist()[0])) * (pGamma.T) * dataMat #update pMiu through MLE(通过令导数 = 0得到)
        pPi = Nk/N

        # update k个 pSigma
        print ‘kk=‘,K
        for kk in range(K):
            Xshift = dataMat-tile(pMiu[kk], (N, 1))  

            Xshift.T * mat(diag(pGamma[:, kk].T.tolist()[0])) *  Xshift / 2

            pSigma[:, :, kk] = (Xshift.T *                 mat(diag(pGamma[:, kk].T.tolist()[0])) * Xshift) / Nk[kk]

        # check for convergence
        L = sum(log(Px*(pPi.T)))
        if L-Lprev < threshold:
            break
        Lprev = L

    return Px

def init_params(X,centroids,K,N,D):
    pMiu = centroids #k*D, 即k类的中心点
    pPi = zeros([1, K]) #k类GMM所占权重(influence factor)
    pSigma = zeros([D, D, K]) #k类GMM的协方差矩阵,每个是D*D的  

    # 距离矩阵,计算N*K的矩阵(x-pMiu)^2 = x^2+pMiu^2-2*x*Miu
    #x^2, N*1的矩阵replicateK列\#pMiu^2,1*K的矩阵replicateN行
    distmat = tile(sum(power(X,2), 1),(1, K)) +         tile(transpose(sum(power(pMiu,2), 1)),(N, 1)) -          2*X*transpose(pMiu)
    labels = distmat.argmin(1) #Return the minimum from each row  

    # 获取k类的pPi和协方差矩阵
    for k in range(K):
        boolList = (labels==k).tolist()
        indexList = [boolList.index(i) for i in boolList if i==[True]]
        Xk = X[indexList, :]
        #print cov(Xk)
        # 也可以用shape(XK)[0]
        pPi[0][k] = float(size(Xk, 0))/N
        pSigma[:, :, k] = cov(transpose(Xk))  

    return pMiu,pPi,pSigma

# 计算每个数据由第k类生成的概率矩阵Px
def calc_prob(pMiu,pSigma,X,K,N,D):
    # Gaussian posterior probability
    # N(x|pMiu,pSigma) = 1/((2pi)^(D/2))*(1/(abs(sigma))^0.5)*exp(-1/2*(x-pMiu)‘pSigma^(-1)*(x-pMiu))
    Px = mat(zeros([N, K]))
    for k in range(K):
        Xshift = X-tile(pMiu[k, :],(N, 1)) #X-pMiu
        #inv_pSigma = mat(pSigma[:, :, k]).I
        inv_pSigma = linalg.pinv(mat(pSigma[:, :, k]))

        tmp = sum(array((Xshift*inv_pSigma)) * array(Xshift), 1) # 这里应变为一列数
        tmp = mat(tmp).T
        #print linalg.det(inv_pSigma),‘54545‘

        Sigema = linalg.det(mat(inv_pSigma))

        if Sigema < 0:
            Sigema=0

        coef = power((2*(math.pi)),(-D/2)) * sqrt(Sigema)
        Px[:, k] = coef * exp(-0.5*tmp)
    return Px

main.py

#! /usr/bin/env python
#coding=utf-8

import GMM
‘‘‘
def showFigure(dataMat,k,clusterAssment):

    tag=[‘go‘,‘or‘,‘yo‘,‘ko‘]
    for i in range(k):

        datalist = dataMat[nonzero(clusterAssment[:,0].A==i)[0]]
        pylab.plot(datalist[:,0],datalist[:,1],tag[i])
    pylab.show()
‘‘‘
if __name__ == ‘__main__‘:
    GMM.gmm(‘testSet.txt‘,2)

说明:

  • 程序由matlab代码改过了,利用了numpy库的函数,需要下载软件可参考链接
  • 其中的分析数据源我用了KMeans算法中的数据,不过遇到了奇异矩阵的问题(因为数据分布本来就有悖于高斯分布),即使用了伪逆也没解决问题,看之后学习能否回头再解决。也可查看博文评论围观别人的解决方法
  • matlab在矩阵的处理上的确优于Python太多了,一方面不用导入库,也没有存在array,mat等类的转换和众多函数的比较和考虑。不过Python综合应用多,上至Web下至PC,从测试到开发都有很多人用,相比matlab还是做测试多一点。
时间: 2024-10-30 15:00:02

【machine learning】GMM算法(Python版)的相关文章

机器学习系统设计(Building Machine Learning Systems with Python)- Willi Richert Luis Pedro Coelho

机器学习系统设计(Building Machine Learning Systems with Python)- Willi Richert Luis Pedro Coelho 总述 本书是 2014 的,看完以后才发现有第二版的更新,2016.建议阅读最新版,有能力的建议阅读英文版,中文翻译有些地方比较别扭(但英文版的书确实是有些贵). 我读书的目的:泛读主要是想窥视他人思考的方式. 作者写书的目标:面向初学者,但有时间看看也不错.作者说"我希望它能激发你的好奇心,并足以让你保持渴望,不断探索

Building.Machine.Learning.Systems.with.Python(2013.7).Willi.Richert.文字版

下载地址 近日网盘管得紧,pdf不能随便分享,下载后请去掉最后的字母 f 即可变成pdf文档了.

学习 Machine Learning Mastery With Python (1)

1 介绍 1.1 机器学习的错误的想法 一定要对python 编程和python语法非常了解 深入学习scikit learn使用的机器学习算法的理论和参数 避免或者不能接触实际项目中的其他部分. 对某些人可能是适用的, 但是对于很多人, 这样会让人觉得望而生畏,经过很长的准备过程学习理论和算法 , 才能开始实际的项目. 很遗憾的是, 这是市面上大部分书籍所采取的办法. 1.2 python 机器学习 这本书主要介绍机器学习的一个子领域叫做预测模型.这个是机器学习在工业界最有用的应用. 本书的组

Machine Learning:机器学习算法

原文链接:https://riboseyim.github.io/2018/02/10/Machine-Learning-Algorithms/ 摘要 机器学习算法分类:监督学习.半监督学习.无监督学习.强化学习 基本的机器学习算法:线性回归.支持向量机(SVM).最近邻居(KNN).逻辑回归.决策树.k平均.随机森林.朴素贝叶斯.降维.梯度增强 目录 监督学习(Supervised learning) 监督学习 原文地址:https://www.cnblogs.com/Antiver/p/99

Building Machine Learning Systems with Python 2

1>监督学习(分类):先让机器学习一下每种花朵的样本数据,然后让他根据这些信息,对未标志出花朵种类的图像进行分类. 2>特征:我们把数据中所有测量的结果都叫特征. 2>交叉验证:极端的叫去一法(leave-one-out)从训练集中拿出一个样本,并在缺少这个样本的数据上训练一个模型,然后看模型是否能够对这个样本正确分类 3>分类模型的组成: 模型结构:采用一个阀值在一个特征上进行划分. 搜素过程:尽可能多的尝试所有特征和阀值的组合. 损失函数:用他来确定哪些可能性不会太差. 4&g

Getting started with machine learning in Python

Getting started with machine learning in Python Machine learning is a field that uses algorithms to learn from data and make predictions. Practically, this means that we can feed data into an algorithm, and use it to make predictions about what might

Python (1) - 7 Steps to Mastering Machine Learning With Python

Step 1: Basic Python Skills install Anacondaincluding numpy, scikit-learn, and matplotlib Step 2: Foundational Machine Learning Skills Unofficial Andrew Ng course notes Tom Mitchell Machine Learning Lectures Step 3: Scientific Python Packages Overvie

Python -- machine learning, neural network -- PyBrain 机器学习 神经网络

I am using pybrain on my Linuxmint 13 x86_64 PC. As what it is described: PyBrain is a modular Machine Learning Library for Python. Its goal is to offer flexible, easy-to-use yet still powerful algorithms for Machine Learning Tasks and a variety of p

Awesome Machine Learning

Awesome Machine Learning  A curated list of awesome machine learning frameworks, libraries and software (by language). Inspired by awesome-php. If you want to contribute to this list (please do), send me a pull request or contact me @josephmisiti Als

机器学习---用python实现感知机算法和口袋算法(Machine Learning PLA Pocket Algorithm Application)

之前在<机器学习---感知机(Machine Learning Perceptron)>一文中介绍了感知机算法的理论知识,现在让我们来实践一下. 有两个数据文件:data1和data2,分别用于PLA和Pocket Algorithm.可在以下地址下载: 先回顾一下感知机算法: 1,初始化w 2,找出一个分类错误点 3,修正错误,假设迭代次数为t次(t=1,2,...),那么修正公式为: 4,直至没有分类错误点,返回最终的w 接下来让我们安照算法步骤,一步一步进行. 首先导入需要用到的库,其中