使用EM算法优化的GMM

先上代码：

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as mat
import matplotlib.mlab as mlab
#EM算法
def EM(dataSet,K):
    (N, M) = np.shape(dataSet)
    W = np.zeros([N, K])
    P= N/K
    for k in range(K):
        W[np.floor(k*P):np.floor((k+1)*P), k] = 1
    A,M,S = Mstep(dataSet,W)
    return W, A, M, S
#M的跨度
def Mstep(data,W):
    (N, M) = np.shape(data)
    K = np.size(W,1) 
    Nk = np.sum(W,0)
    A = Nk/np.sum(Nk)
    Mm = data.T.dot(W).dot(np.diag(np.reciprocal(Nk)))
    S = np.zeros([M,M,K])
    for k in range(K):
        datMean = data.T - Mm[0:,k][None].T.dot(np.ones([1,N]))
        S[:,:,k] = (datMean.dot(np.diag(W[0:,k])).dot(datMean.T))/Nk[k]
    return A,Mm,S
#E的跨度
def Estep(data,A,M,S):
    N = np.size(data,0)
    K = np.size(A)
    W = np.zeros([N,K])
    for k in range(K):
        for i in range(N):
            W[i,k] = A[k]*multivariate(data[i,:][None].T, \
                     M[:,k][None].T,S[:,:,k])
    W = W*np.reciprocal(np.sum(W,1)[None].T)
    return W
def multivariate(x, m, s):
    if len(x) == len(m) and (len(x), len(x)) == s.shape:
        det = np.linalg.det(s)
        const = 1.0/(np.math.pow((2*np.pi), float(len(x))/2) * np.math.pow(det, 1.0/2))
        x_m = np.matrix(x - m)
        inv_ = np.linalg.inv(s)
        result = np.math.pow(np.math.e, -0.5 * (x_m.T * inv_ * x_m))
        return const * result
    else:
        return -1
#GMM主程序
def GMM():
    # 加载文件
    input_file = open('points.dat')
    lines = input_file.readlines()
    Data = np.array([line.strip().split() for line in lines]).astype(np.float)
    (x, y) = np.shape(Data)
    mat.draw()
    mat.pause(0.01)
    mat.subplot(111)
    mat.plot(x, y, 'b*')
    learn = Data[np.math.ceil(x*0.8):x, 0:]
    train = Data[:np.math.floor(x*0.8), 0:]
    trainnum = 16
    (W, Alpha, Mu, Sigma) = EM(train,trainnum)
    m = np.arange(-4.0, 4.0, 0.1)
    n = np.arange(-4.0, 4.0, 0.1)
    ax, ay = np.meshgrid(m, n)
    i = 0
    prev = -9999
    mat.clf()
    while(True):    
        if(False):
            SigmaSum = np.sum(Sigma,2)
            for k in range(trainnum):
                Sigma[:,:,k] = SigmaSum
        W = Estep(train,Alpha,Mu,Sigma)
        Alpha,Mu,Sigma = Mstep(train,W)
       # trains = logLike(train,Alpha,Mu,Sigma)
        N,M = np.shape(train)
        P = np.zeros([N,len(Alpha)])
        for k in range(len(Alpha)):
            for j in range(N):
                P[j,k] = multivariate(train[j,:][None].T,Mu[0:,k][None].T,Sigma[:,:,k])
        trains = np.sum(np.log(P.dot(Alpha)))
        i = i + 1
        #画图，训练和测试样本
        mat.subplot(211)
        mat.scatter(train[0:,0],train[0:,1])
        mat.hold(True)
        for k in range(0, trainnum):
            az = mlab.bivariate_normal(ax, ay, Sigma[0, 0, k], Sigma[1, \                1, k], Mu[0,k], Mu[1,k], Sigma[1, 0, k])
            try:
                mat.contour(ax, ay, az)
            except:
                continue
        mat.hold(False)
        # Render these
        mat.draw()
        mat.pause(0.01)
        mat.subplot(212)
        mat.scatter(learn[0:,0],learn[0:,1])
        mat.hold(True)
        for k in range(0, trainnum):
            az = mlab.bivariate_normal(ax, ay, Sigma[0, 0, k], Sigma[1, \                1, k], Mu[0,k], Mu[1,k], Sigma[1, 0, k])
            try:
                mat.contour(ax, ay, az)
            except:
                continue
        mat.hold(False)
        if(i>150 or abs(trains - prev)< 0.01):
            break
        prev = trains
if __name__ == '__main__':
    GMM()

实验要求

实验目标

实现一个混合高斯模型，并且用EM算法估计模型中的参数。

实验过程

用混合高斯模型产生k个高斯分布的数据（其中参数自己设定），然后用你实现的EM算法估计参数，看看每次迭代后似然值变化情况，考察EM算法是否可以获得正确的结果（与你设定的结果比较）。
可以UCI上找一个简单问题数据，用你实现的GMM进行聚类。

算法原理

GMM算法

高斯模型就是用高斯概率密度函数（正态分布曲线）精确地量化事物，将一个事物分解为若干的基于高斯概率密度函数（正态分布曲线）形成的模型。

对图像背景建立高斯模型的原理及过程：图像灰度直方图反映的是图像中某个灰度值出现的频次，也可以认为是图像灰度概率密度的估计。如果图像所包含的目标区域和背景区域相比比较大，且背景区域和目标区域在灰度上有一定的差异，那么该图像的灰度直方图呈现双峰-谷形状，其中一个峰对应于目标，另一个峰对应于背景的中心灰度。对于复杂的图像，尤其是医学图像，一般是多峰的。通过将直方图的多峰特性看作是多个高斯分布的叠加，可以解决图像的分割问题。 在智能监控系统中，对于运动目标的检测是中心内容，而在运动目标检测提取中，背景目标对于目标的识别和跟踪至关重要。而建模正是背景目标提取的一个重要环节。

我们首先要提起背景和前景的概念，前景是指在假设背景为静止的情况下，任何有意义的运动物体即为前景。建模的基本思想是从当前帧中提取前景，其目的是使背景更接近当前视频帧的背景。即利用当前帧和视频序列中的当前背景帧进行加权平均来更新背景,但是由于光照突变以及其他外界环境的影响，一般的建模后的背景并非十分干净清晰，而高斯混合模型是是建模最为成功的方法之一。

混合高斯模型使用K（基本为3到5个）个高斯模型来表征图像中各个像素点的特征,在新一帧图像获得后更新混合高斯模型, 用当前图像中的每个像素点与混合高斯模型匹配,如果成功则判定该点为背景点, 否则为前景点。

通观整个高斯模型，主要是有方差和均值两个参数决定，对均值和方差的学习，采取不同的学习机制,将直接影响到模型的稳定性、精确性和收敛性 。由于我们是对运动目标的背景提取建模，因此需要对高斯模型中方差和均值两个参数实时更新。为提高模型的学习能力,改进方法对均值和方差的更新采用不同的学习率;为提高在繁忙的场景下,大而慢的运动目标的检测效果,引入权值均值的概念,建立背景图像并实时更新,然后结合权值、权值均值和背景图像对像素点进行前景和背景的分类。

具体实现过程：

1. 为图像的每个像素点指定一个初始的均值、标准差以及权重。
2. 收集N（一般取200以上，否则很难得到像样的结果）帧图像利用在线EM算法得到每个像 素点的均值、标准差以及权重。
3. 从N+1帧开始检测，检测的方法：
   对每个像素点：
   1. 将所有的高斯核按照 ω / σ 降序排序
   2. 选择满足下式的前M个高斯核：M = arg min(ω / σ > T)
   3. 如果当前像素点的像素值在中有一个满足：就可以认为其为背景点。
4. 更新背景图像，用EM算法。

EM算法

EM 算法是 Dempster，Laind，Rubin 于 1977 年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行 MLE 估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据，截尾数据，带有噪声等所谓的不完全数据(incomplete data)。

最大期望算法经过两个步骤交替进行计算：

• 第一步是计算期望（E），利用对隐藏变量的现有估计值，计算其最大似然估计值；
• 第二步是最大化（M），最大化在 E 步上求得的最大似然值来计算参数的值。

M 步上找到的参数估计值被用于下一个 E 步计算中，这个过程不断交替进行。

通过交替使用这两个步骤，EM算法逐步改进模型的参数，使参数和训练样本的似然概率逐渐增大，最后终止于一个极大点。直观地理解EM算法，它也可被看作为一个逐次逼近算法：事先并不知道模型的参数，可以随机的选择一套参数或者事先粗略地给定某个初始参数λ0 ，确定出对应于这组参数的最可能的状态，计算每个训练样本的可能结果的概率，在当前的状态下再由样本对参数修正，重新估计参数λ，并在新的参数下重新确定模型的状态，这样，通过多次的迭代，循环直至某个收敛条件满足为止，就可以使得模型的参数逐渐逼近真实参数。

具体实现步骤：

• 给出种类数k,初始化每一类高斯分布的均值μk，方差∑k以及每一类的概率πk；
• 执行EM；
• 计算似然，如果没有达到预期效果，则返回第2步；
• 计算每个数据点对于k个高斯分布的似然，选择似然最大的一类作为数据的最终分类。

其他的混合模型，例如朴素贝叶斯混合模型也是可以使用EM算法推出使用的，这一算法虽然在GMM中作为参数使用，但是其仍然可以单独发挥作用。我觉得EM算法就是相互迭代（毕竟其由E和M两部分组成嘛），求出一个稳定值，而这种相互迭代的方法用的范围挺广的，例如混合模型，K-means等都需要使用。

与K-means的区别

在上文我也已经提到了EM算法可以用在K-means等其他需要迭代的方法上的事实，其实，我觉得GMM 和 K-means 很像，只不过后者要简单，而且相对来说实现并不是很高效。不过 GMM 是学习出一些概率密度函数来（所以 GMM 除了用在 clustering 上之外，还经常被用于 density estimation ），简单地说，K-means 的结果是每个数据点被 assign 到其中某一个 cluster 了，而 GMM 则给出这些数据点被 assign 到每个 cluster 的概率，又称作 soft assignment。

实验环境

Spyder 作为Python开发的集成开发环境；
编程语言：Python 3.5；
操作系统：macOS Sierra。

小结

GMM算法作为EM算法族的一个例子，它指定了各个参与杂合的分布都是高斯分布，即分布参数表现为均值Mu和方差Sigma。通过EM算法作为计算使用的框架，迭代地算出各个高斯分布的参数。

GMM与K-means的思考

提到GMM不得不提K-means，总结了网上的资料以及老师上课的课件，我将两者的区别与联系陈述如下：

• 两者的联系:

都是迭代执行的算法，且迭代的策略也相同：算法开始执行时先对需要计算的参数赋初值，然后交替执行两个步骤，一个步骤是对数据的估计（k-means是估计每个点所属簇；GMM是计算隐含变量的期望）；第二步是用上一步算出的估计值重新计算参数值，更新目标参数（k-means是计算簇心位置；GMM是计算各个高斯分布的中心位置和协方差矩阵）

• 两者的区别:

首先，两者需要计算的参数不同：K-means是簇心位置；GMM是各个高斯分布的参数；其次，两者计算目标参数的方法不同：K-means是计算当前簇中所有元素的位置的均值；GMM是基于概率的算法，是通过计算似然函数的最大值实现分布参数的求解的。

关于GMM引发的过拟合的思考

首先我想提到这样的一个“人辨认其他生物（例如鱼）”的例子。当我们被告知水里游的那个生物是鱼之后，我们会使用“在同样的地方生活的是同一种东西”这类似的假设，归纳出“在水里游的都是鱼”这样一个结论。当然这个过程是完全“本能”的，如果不仔细去想，我们也不会了解自己是怎样“认识鱼”的。另一个值得注意的地方是这样的假设并不总是完全正确的，甚至可以说总是会有这样那样的缺陷的，因为我们有可能会把虾、龟、甚至是潜水员当做鱼。也许你觉得可以通过修改前提假设来解决这个问题，例如，基于“生活在同样的地方并且穿着同样衣服的是同一种东西”这个假设，你得出结论：在水里有并且身上长有鳞片的是鱼。可是这样还是有问题，因为有些没有长鳞片的鱼现在又被你排除在外了。

机器在识别方面面临着和人一样的问题，在机器学习中，一个学习算法也会有一个前提假设，这里被称作“归纳偏执”。例如线性回归，目的是要找一个函数尽可能好地拟合给定的数据点，它的归纳偏执就是“满足要求的函数必须是线性函数”。一个没有归纳偏执的学习算法从某种意义上来说毫无用处，就像一个完全没有归纳能力的人一样，在第一次看到鱼的时候有人告诉他那是鱼，下次看到另一条鱼了，他并不知道那也是鱼，因为两条鱼总有一些地方不一样的，或者就算是同一条鱼，在河里不同的地方看到，或者只是看到的时间不一样，也会被他认为是不同的，因为他无法归纳，无法提取主要矛盾、忽略次要因素，只好要求所有的条件都完全一样──然而哲学家已经告诉过我们了：世界上不会有任何样东西是完全一样的，所以这个人即使是有无比强悍的记忆力，也绝学不到任何一点知识。

于是有了上面的铺垫，我们就可以引出论题——“过拟合 ”，就像前面的回归的问题，如果去掉“线性函数”这个归纳偏执，因为对于 N 个点，我们总是可以构造一个 N-1 次多项式函数，让它完美地穿过所有的这 N 个点，或者如果我用任何大于 N-1 次的多项式函数的话，我们甚至可以构造出无穷多个满足条件的函数出来。如果假定特定领域里的问题所给定的数据个数总是有个上限的话，我可以取一个足够大的 N ，从而得到一个（或者无穷多个）“超级函数”，能够拟合这个领域内所有的问题。

没有归纳偏执或者归纳偏执太宽泛会导致过拟合 ，然而另一个极端──限制过大的归纳偏执也是有问题的：如果数据本身并不是线性的，强行用线性函数去做回归通常并不能得到好结果（例如我在“实验一：多项式拟合曲线”中就做过相应的测试）。难点正在于在这之间寻找一个临界点。不过我们在这里相对于机器来说有一个很大的优势：人通常不会孤立地用某一个独立的系统和模型去处理问题，一个人每天都会从各个来源获取大量的信息，并且通过各种手段进行整合处理，归纳所得的所有知识最终得以统一地存储起来，并能有机地组合起来去解决特定的问题。

以上就是我关于“过拟合”的一点不全面的思考！

The link of this page is https://blog.nooa.tech/articles/177fbbcc/ . Welcome to reproduce it!