混合高斯模型的EM求解(Mixtures of Gaussians)及Python实现源码

今天为大家带来混合高斯模型的EM推导求解过程。

全部代码如下!

def NDimensionGaussian(X_vector,U_Mean,CovarianceMatrix):
    #X=numpy.mat(X_vector)
    X=X_vector
    D=numpy.shape(X)[0]
    #U=numpy.mat(U_Mean)
    U=U_Mean
    #CM=numpy.mat(CovarianceMatrix)
    CM=CovarianceMatrix
    Y=X-U
    temp=Y.transpose() * CM.I * Y
    result=(1.0/((2*numpy.pi)**(D/2)))*(1.0/(numpy.linalg.det(CM)**0.5))*numpy.exp(-0.5*temp)
    return result

def CalMean(X):
    D,N=numpy.shape(X)
    MeanVector=numpy.mat(numpy.zeros((D,1)))
    for d in range(D):
        for n in range(N):
            MeanVector[d,0] += X[d,n]
        MeanVector[d,0] /= float(N)
    return MeanVector

def CalCovariance(X,MV):
    D,N=numpy.shape(X)
    CoV=numpy.mat(numpy.zeros((D,D)))
    for n in range(N):
        Temp=X[:,n]-MV
        CoV += Temp*Temp.transpose()
    CoV /= float(N)
    return CoV

def CalEnergy(Xn,Pik,Uk,Cov):
    D,N=numpy.shape(Xn)
    D_k,K=numpy.shape(Uk)
    if D!=D_k:
        print ('dimension not equal, break')
        return

    energy=0.0
    for n_iter in range(N):
        temp=0
        for k_iter in range(K):
            temp += Pik[0,k_iter] * NDimensionGaussian(Xn[:,n_iter],Uk[:,k_iter],Cov[k_iter])
        energy += numpy.log(temp)
    return float(energy)

def SequentialEMforMixGaussian(InputData,K):
    #初始化piK
    pi_Cof=numpy.mat(numpy.ones((1,K))*(1.0/float(K)))
    X=numpy.mat(InputData)
    X_mean=CalMean(X)
    print (X_mean)
    X_cov=CalCovariance(X,X_mean)
    print (X_cov)
    #初始化uK,其中第k列表示第k个高斯函数的均值向量
    #X为D维,N个样本点
    D,N=numpy.shape(X)
    print (D,N)
    UK=numpy.mat(numpy.zeros((D,K)))
    for d_iter in range(D):
        for k_iter in range(K):
            UK[d_iter,k_iter] = X_mean[d_iter,0] + (-1)**k_iter + (-1)**d_iter
    print (UK)
    #初始化k个协方差矩阵的列表
    List_cov=[]

    for k_iter in range(K):
        List_cov.append(numpy.mat(numpy.eye(X[:,0].size)))
    print (List_cov)

    List_cov_new=copy.deepcopy(List_cov)
    rZnk=numpy.mat(numpy.zeros((N,K)))
    denominator=numpy.mat(numpy.zeros((N,1)))
    rZnk_new=numpy.mat(numpy.zeros((N,K)))

    Nk=0.5*numpy.mat(numpy.ones((1,K)))
    print (Nk)
    Nk_new=numpy.mat(numpy.zeros((1,K)))
    UK_new=numpy.mat(numpy.zeros((D,K)))
    pi_Cof_new=numpy.mat(numpy.zeros((1,K)))

    for n_iter in range(1,N):
        #rZnk=pi_k*Gaussian(Xn|uk,Cov_k)/sum(pi_j*Gaussian(Xn|uj,Cov_j))
        for k_iter in range(K):
            rZnk_new[n_iter,k_iter] = pi_Cof[0,k_iter] * NDimensionGaussian(X[:,n_iter],UK[:,k_iter],List_cov[k_iter])
            denominator[n_iter,0] += rZnk_new[n_iter,k_iter]
        for k_iter in range(K):
            rZnk_new[n_iter,k_iter] /= denominator[n_iter,0]
            print ('rZnk_new', rZnk_new[n_iter,k_iter],'\n')
        for k_iter in range(K):
            Nk_new[0,k_iter] = Nk[0,k_iter] + rZnk_new[n_iter,k_iter] - rZnk[n_iter,k_iter]
            print ('Nk_new',Nk_new,'\n')
            ##############当前有(n_iter+1)样本###########################
            pi_Cof_new[0,k_iter] = Nk_new[0,k_iter]/float(n_iter+1)
            print ('pi_Cof_new',pi_Cof_new,'\n')
            UK_new[:,k_iter] = UK[:,k_iter] + ( (rZnk_new[n_iter,k_iter] - rZnk[n_iter,k_iter])/float(Nk_new[0,k_iter]) ) * (X[:,n_iter]-UK[:,k_iter])
            print ('UK_new',UK_new,'\n')
            Temp = X[:,n_iter] - UK_new[:,k_iter]
            List_cov_new[k_iter] = List_cov[k_iter] + ((rZnk_new[n_iter,k_iter] - rZnk[n_iter,k_iter])/float(Nk_new[0,k_iter]))*(Temp*Temp.transpose()-List_cov[k_iter])
            print ('List_cov_new',List_cov_new,'\n')

        rZnk=copy.deepcopy(rZnk_new)
        pi_Cof=copy.deepcopy(pi_Cof_new)
        UK_new=copy.deepcopy(UK)
        List_cov=copy.deepcopy(List_cov_new)
    print (pi_Cof,UK_new,List_cov)
    return pi_Cof,UK_new,List_cov

def BatchEMforMixGaussian(InputData,K,MaxIter):
    #初始化piK
    pi_Cof=numpy.mat(numpy.ones((1,K))*(1.0/float(K)))
    X=numpy.mat(InputData)
    X_mean=CalMean(X)
    print (X_mean)
    X_cov=CalCovariance(X,X_mean)
    print (X_cov)
    #初始化uK,其中第k列表示第k个高斯函数的均值向量
    #X为D维,N个样本点
    D,N=numpy.shape(X)
    print (D,N)
    UK=numpy.mat(numpy.zeros((D,K)))
    for d_iter in range(D):
        for k_iter in range(K):
            UK[d_iter,k_iter] = X_mean[d_iter,0] + (-1)**k_iter + (-1)**d_iter
    print (UK)
    #初始化k个协方差矩阵的列表
    List_cov=[]

    for k_iter in range(K):
        List_cov.append(numpy.mat(numpy.eye(X[:,0].size)))
    print (List_cov)

    energy_new=0
    energy_old=CalEnergy(X,pi_Cof,UK,List_cov)
    print (energy_old)
    currentIter=0
    while True:
        currentIter += 1

        List_cov_new=[]
        rZnk=numpy.mat(numpy.zeros((N,K)))
        denominator=numpy.mat(numpy.zeros((N,1)))
        Nk=numpy.mat(numpy.zeros((1,K)))
        UK_new=numpy.mat(numpy.zeros((D,K)))
        pi_new=numpy.mat(numpy.zeros((1,K)))

        #rZnk=pi_k*Gaussian(Xn|uk,Cov_k)/sum(pi_j*Gaussian(Xn|uj,Cov_j))
        for n_iter in range(N):
            for k_iter in range(K):
                rZnk[n_iter,k_iter] = pi_Cof[0,k_iter] * NDimensionGaussian(X[:,n_iter],UK[:,k_iter],List_cov[k_iter])
                denominator[n_iter,0] += rZnk[n_iter,k_iter]
            for k_iter in range(K):
                rZnk[n_iter,k_iter] /= denominator[n_iter,0]
                #print 'rZnk', rZnk[n_iter,k_iter]

        #pi_new=sum(rZnk)
        for k_iter in range(K):
            for n_iter in range(N):
                Nk[0,k_iter] += rZnk[n_iter,k_iter]
            pi_new[0,k_iter] = Nk[0,k_iter]/(float(N))
            #print 'pi_k_new',pi_new[0,k_iter]

        #uk_new= (1/sum(rZnk))*sum(rZnk*Xn)
        for k_iter in range(K):
            for n_iter in range(N):
                UK_new[:,k_iter] += (1.0/float(Nk[0,k_iter]))*rZnk[n_iter,k_iter]*X[:,n_iter]
            #print 'UK_new',UK_new[:,k_iter]

        for k_iter in range(K):
            X_cov_new=numpy.mat(numpy.zeros((D,D)))
            for n_iter in range(N):
                Temp = X[:,n_iter] - UK_new[:,k_iter]
                X_cov_new += (1.0/float(Nk[0,k_iter]))*rZnk[n_iter,k_iter] * Temp * Temp.transpose()
            #print 'X_cov_new',X_cov_new
            List_cov_new.append(X_cov_new)

        energy_new=CalEnergy(X,pi_new,UK_new,List_cov)
        print ('energy_new',energy_new)
        #print pi_new
        #print UK_new
        #print List_cov_new
        if energy_old>=energy_new or currentIter>MaxIter:
            UK=copy.deepcopy(UK_new)
            pi_Cof=copy.deepcopy(pi_new)
            List_cov=copy.deepcopy(List_cov_new)
            break
        else:
            UK=copy.deepcopy(UK_new)
            pi_Cof=copy.deepcopy(pi_new)
            List_cov=copy.deepcopy(List_cov_new)
            energy_old=energy_new

return pi_Cof,UK,List_cov
时间: 2024-10-12 07:35:58

混合高斯模型的EM求解(Mixtures of Gaussians)及Python实现源码的相关文章

混合高斯模型(Mixtures of Gaussians)和EM算法

混合高斯模型(Mixtures of Gaussians)和EM算法 主要内容: 1. 概率论预备知识 2. 单高斯模型 3. 混合高斯模型 4. EM算法 5. K-means聚类算法 一.概率论预备知识 1. 数学期望/均值.方差/标准差 设离散型随机变量X的分布律为 则称为X的数学期望或均值 设连续型随机变量X的概率密度函数(pdf)为 则其数学期望定义为: 随机变量X的方差: 随机变量X的标准差: 2. 正态分布.协方差 正态分布: 概率密度函数: 设(X,Y)为二维随机变量,若存在,则

【转载】混合高斯模型(Mixtures of Gaussians)和EM算法

混合高斯模型(Mixtures of Gaussians)和EM算法 这篇讨论使用期望最大化算法(Expectation-Maximization)来进行密度估计(density estimation). 与k-means一样,给定的训练样本是,我们将隐含类别标签用表示.与k-means的硬指定不同,我们首先认为是满足一定的概率分布的,这里我们认为满足多项式分布,,其中,有k个值{1,…,k}可以选取.而且我们认为在给定后,满足多值高斯分布,即.由此可以得到联合分布. 整个模型简单描述为对于每个

EM算法与混合高斯模型

很早就想看看EM算法,这个算法在HMM(隐马尔科夫模型)得到很好的应用.这个算法公式太多就手写了这部分主体部分. 好的参考博客:最大似然估计到EM,讲了具体例子通熟易懂. JerryLead博客很不错 混合高斯模型算法

混合高斯模型(Mixtures of Gaussians)

http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006924.html 这篇讨论使用期望最大化算法(Expectation-Maximization)来进行密度估计(density estimation). 与k-means一样,给定的训练样本是,我们将隐含类别标签用表示.与k-means的硬指定不同,我们首先认为是满足一定的概率分布的,这里我们认为满足多项式分布,,其中,有k个值{1,…,k}可以选取.而且我们认为在给定后,满足多值高斯

斯坦福ML公开课笔记13A——混合高斯模型、混合贝叶斯模型

本文对应公开课的第13个视频,这个视频仍然和EM算法非常相关,第12个视频讲解了EM算法的基础,本视频则是在讲EM算法的应用.本视频的主要内容包括混合高斯模型(Mixture of Gaussian, MoG)的EM推导.混合贝叶斯模型(Mixture of Naive Bayes,MoNB)的EM推导.因子分析模型(Factor Analysis Model)及其EM求解.由于本章内容较多,故而分为AB两篇,本篇介绍至混合模型的问题. 很久没有写这个系列的笔记了,各种事情加各种懒导致的.虽然慢

混合高斯模型算法(转)

下面介绍一下几种典型的机器算法 首先第一种是高斯混合模型算法: 高斯模型有单高斯模型(SGM)和混合高斯模型(GMM)两种. (1)单高斯模型: 为简单起见,阈值t的选取一般靠经验值来设定.通常意义下,我们一般取t=0.7-0.75之间. 二维情况如下所示: (2)混合高斯模型: 对于(b)图所示的情况,很明显,单高斯模型是无法解决的.为了解决这个问题,人们提出了高斯混合模型(GMM),顾名思义,就是数据可以看作是从数个高 斯分布中生成出来的.虽然我们可以用不同的分布来随意地构造 XX Mixt

聚类——混合高斯模型 Gaussian Mixture Model

转自: http://blog.csdn.net/jwh_bupt/article/details/7663885聚类系列: 聚类(序)----监督学习与无监督学习 聚类(1)----混合高斯模型 Gaussian Mixture Model 聚类(2)----层次聚类 Hierarchical Clustering 聚类(3)----谱聚类 Spectral Clustering -------------------------------- 聚类的方法有很多种,k-means要数最简单的一

混合高斯模型聚类

混合高斯模型简介 混合高斯模型基于多变量正态分布.类gmdistribution通过使用EM算法来拟合数据,它基于各观测量计算各成分密度的后验概率. 高斯混合模型常用于聚类,通过选择成分最大化后验概率来完成聚类.与k-means聚类相似,高斯混合模型也使用迭代算法计算,最终收敛到局部最优.高斯混合模型在各类尺寸不同.聚类间有相关关系的的时候可能比k-means聚类更合适.使用高斯混合模型的聚类属于软聚类方法(一个观测量按概率属于各个类,而不是完全属于某个类),各点的后验概率提示了各数据点属于各个

OpenCV运动目标检测——帧间差,混合高斯模型方法

一.简单的帧间差方法 帧差法是在连续的图像序列中两个或三个相邻帧间采用基于像素的时间差分并且闽值化来提取图像中的运动区域. 代码: int _tmain(int argc, _TCHAR* argv[]) { VideoCapture capture("bike.avi"); if(!capture.isOpened()) return -1; double rate = capture.get(CV_CAP_PROP_FPS); int delay = 1000/rate; Mat