第一章引子

假设放在你面前有5篮子鸡蛋&＃xff0c;每个篮子有且仅有一种蛋&＃xff0c;这些蛋表面上一模一样&＃xff0c;就是每一种蛋涵盖有且只有一种维生素&＃xff0c;分别是A、B、C、D、E。这个时候&＃xff0c;你需要估计这五个篮子的鸡蛋的平均重量μ。 首先有个总的假设&＃xff1a; 假设每一种维生素的鸡蛋的重量都服从高斯分布。 这个时候&＃xff0c;因为每个篮子的鸡蛋包含有且只有一种&＃xff0c;并且彼此之间相同的维生素&＃xff0c;即每个篮子的鸡蛋都服从相同的分布&＃xff0c;这个时候可以用极大似然估计去估计每一种维生素鸡蛋的平均重量。

现在问题来了: 我把那5种鸡蛋混在一起&＃xff0c;这个时候要你去估计这5中鸡蛋的平均重量和方差&＃xff1f; 仍旧假设&＃xff1a;每一种维生素的鸡蛋的重量都服从高斯分布 这个时候有两个参数需要估计均值μ和方差σ。 你从5中鸡蛋里面去拿一种鸡蛋&＃xff0c;每种鸡蛋被拿到的概率是呈一定分布的(不一定就是均匀分布&＃xff0c;然后概率是1/5)。假设第j中鸡蛋被拿到的概率是φ_j。

因此你有三个未知量&＃xff0c;即隐含量: φ_j、均值μ和方差σ。

现在是如果你知道类别z^(j)&＃xff0c;你就能根据极大似然估计计算其他两个隐含量。令每个鸡蛋的重量用x⁽ⁱ⁾表示&＃xff0c;x⁽ⁱ⁾服从高斯分布。 。(x⁽ⁱ⁾|z⁽ⁱ⁾&＃61;j)~N(μ,σ²)

高斯混合模型(gaussian mixture model-GMM)就是上述问题的抽象模型&＃xff0c;即对于每一个样本x⁽ⁱ⁾&＃xff0c;先从k个类别中按某种分布抽取一个z⁽ⁱ⁾&＃xff0c;然后从这个类别下的高斯分布中生成一个样本x⁽ⁱ⁾

第二章 GMM数学推导

假设独立样本集是{x⁽¹⁾......x^(m)}&＃xff0c;这些样本是从k个高斯分布的数据里面抽取出来的。从k不同分布中抽取某一类别是呈某种分布,假设抽取到类别z⁽ⁱ⁾的概率是p(z⁽ⁱ⁾&＃61;j)&＃61; φ_j。因此这里面会有三个隐含变量: φ_j、均值μ和方差σ&＃xff0c;令这三个变量构成一个集合θ.

则利用极大似然估计θ就是我们非常熟悉的极大似然估计函数:

即:

就连乘变成连加&＃xff0c;取对数有:

接下来&＃xff0c;我们利用EM算法来估计这些参数。EM算法请参看(http://blog.csdn.net/u010866505/article/details/77877345).

1.       先假设我们知道样本的类别z⁽ⁱ⁾&＃xff0c;然后计算期望E&＃xff0c;即后验概率&＃xff1a;

2.       然后就是M-step&＃xff1a;

求偏导可以计算均值 μj:

在公式(1)中&＃xff0c;如果均值和方差固定的话,那么(1)式可以简化成:

、

为了计算优化(3)式&＃xff0c;而 又满足一定的条件&＃xff0c;即 &＃xff0c;如果你知道大名鼎鼎的支持向量机(SVM)的优化目标函数是如何得到的话&＃xff1f;这里也就明白了。拉格朗日乘子。构造拉格朗日函数如下:

求导:

令(5)式&＃61;0&＃xff0c;计算有:

上面(6)式两边做如下处理:

                                                                                

因此(6)式得到如下变换:

接下来就是计算方差:对(1)公式计算方差 的偏导:

接下来要计算方差的偏导&＃xff0c;因此(9)式中括号内的第一项和最后一项可以不用考虑。于是有:

令(10)式等于0&＃xff0c;可以计算得到方差的公式如下&＃xff1a;

(2)&＃43;(8)&＃43;(11)就是我们估计参数的最后的解。

 第三章 GMM代码-python实现
如下是GMM算法的简单实现&＃xff0c;

#! /usr/bin/env python
#! -*- coding&＃61;utf-8 -*-#模拟两个正态分布的参数from numpy import *
import numpy as np
import random
import copySIGMA &＃61; 6
EPS &＃61; 0.0001
#均值不同的样本
def generate_data(): Miu1 &＃61; 2Miu2 &＃61; 4sigma1 &＃61; 1sigma2 &＃61; 2alpha1 &＃61; 0.4alpha2 &＃61; 0.6N &＃61; 5000N1 &＃61; int(alpha1 * N)X &＃61; mat(zeros((N,1)))for i in range(N1):temp &＃61; random.uniform(0,0.5)X[i] &＃61; temp * sigma1 &＃43; Miu1for i in range(N-N1):temp &＃61; random.uniform(0,0.5)X[i&＃43;N1] &＃61; temp * sigma2 &＃43; Miu2return X#EM算法
def my_GMM(X):k &＃61; 2N &＃61; len(X)Miu &＃61; np.random.rand(k,1)Posterior &＃61; mat(zeros((N,k))) sigma &＃61; np.random.rand(k,1)sigma[0]&＃61;1#sigma[1]&＃61;2alpha &＃61; np.random.rand(k,1)alpha[0] &＃61; 0.1alpha[1] &＃61; 0.9dominator &＃61; 0numerator &＃61; 0#先求后验概率print sigmafor it in range(1000):for i in range(N):dominator &＃61; 0for j in range(k):dominator &＃61; dominator &＃43; np.exp(-1.0/(2.0*sigma[j]) * (X[i] - Miu[j])**2)#print -1.0/(2.0*sigma[j]),(X[i] - Miu[j])**2,-1.0/(2.0*sigma[j]) * (X[i] - Miu[j])**2,np.exp(-1.0/(2.0*sigma[j]) * (X[i] - Miu[j])**2)#returnfor j in range(k):numerator &＃61; np.exp(-1.0/(2.0*sigma[j]) * (X[i] - Miu[j])**2)Posterior[i,j] &＃61; numerator/dominator oldMiu &＃61; copy.deepcopy(Miu)oldalpha &＃61; copy.deepcopy(alpha)oldsigma &＃61; copy.deepcopy(sigma)#最大化 for j in range(k):numerator &＃61; 0dominator &＃61; 0for i in range(N):numerator &＃61; numerator &＃43; Posterior[i,j] * X[i]dominator &＃61; dominator &＃43; Posterior[i,j]Miu[j] &＃61; numerator/dominatoralpha[j] &＃61; dominator/Ntmp &＃61; 0for i in range(N):tmp &＃61; tmp &＃43; Posterior[i,j] * (X[i] - Miu[j])**2#print tmp,Posterior[i,j],(X[i] - Miu[j])**2 sigma[j] &＃61; tmp/dominatorprint tmp, dominator, sigma[j]if ((abs(Miu - oldMiu)).sum()

参考资料

[1]:http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html