使用EM算法优化的GMM

先上代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# -*- 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 ,从而得到一个(或者无穷多个)“超级函数”,能够拟合这个领域内所有的问题。

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

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



本文链接: http://home.meng.uno/articles/177fbbcc/ 欢迎转载!

© 2018.02.08 - 2020.06.02 Mengmeng Kuang  保留所有权利!

UV : | PV :

:D 获取中...

Creative Commons License