GMM（高斯混合模型）以及简单实现

本文主要是学习比较，主要参考这几篇博文，写的非常通俗易懂：

http://blog.pluskid.org/?p=39

http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/

http://blog.csdn.net/gugugujiawei/article/details/45583051

http://blog.csdn.net/hevc_cjl/article/details/9733945

http://wenku.baidu.com/link?url=YQZJEErAFb4yDLfkVqMYhOtcl88BVyo9oLrB0aquT3tv-v5lYRo_J3Qpzg-vQ4laxHSGIFGP9eDjbuIHUyeUK1tKpR8hdhHXJD3oSdTqtqe

拖了两天了，今天唯一的任务就是把GMM搞清楚，然后把博文写出来， 嗯嗯，就这样。

先看一下单高斯模型（SGM）：

公式都写成代码的形式了，工具还不太会用。

正态分布的概率密度是：

X~ N(u, sigma^2) –> pdf(x) = 1./(math.sqrt(2 * math.pi) * sigma) * exp( -(x-u)^2 / (2. * sigma^2) )

多维单高斯（正态）分布的概率密度是（在sre或者asr中一般去mfcc特征，x=(x1,x2,x3,…,xd). T的维度d一般为39或者42, 注意x为列向量， 实际计算中一般是行向量，便于计算）：

pdf(x) = 1./( (math.sqrt((2 * math.pi)^d) * np.linalg.det(C)) ) * exp(-0.5 * (x-u).T * np.linalg.inv(C) * (x-u))

其中C为x1,x2,..,xd的协方差矩阵,

u为x的期望(u的维数是xi的维数，比如这个地方是39，可以看作u1是每个x的第一个数做均值)，

np.linalg.det(C)求协方差矩阵的行列式

np.linalg.inv(C)是C逆，

由于数据类型都是numpy，直接用 (x-u).T表示转置。

对于SGM,可以用于模式识别中的二分类问题，比如有单高斯模型K, 判断一个训练样本是不是属于K类，求 N（x , c）,

其中 c ＝ （u , C）, u为mean vector , C为 covariance matrix.

N（x , c）= pdf = 1./( (math.sqrt((2 * math.pi)^d * np.linalg.det(C) ))) * exp(-0.5 * (x-u).T * np.linalg.inv(C) * (x-u))

如果pdf大于阈值t, 则表示属于C类， 否则不属于。其中阈值t的选择可以是经验值或者是实验的结果或者是其他的策略，这里不展开，下面代码中使用一个经验值，比如0.7

NOTE： SGM的几何意义，注意通过集合意义与GMM区分，理解内在的含义，二维点在平面上近似椭圆，三维点在空间中近似椭球，但是无论如何都只是一个超平面，但是GMM是多个不同的SGM,通过平移来平滑的描述任何的图形。

高斯混合模型 GMM：

例如有一批观察数据X=(X1, X2, X3, …, Xn), 样本个数为N，每一个样本的维度为D, 在空间的分布不是椭球壮，那么就不适合SGM的pdf来描述。 此时我们假设每个点均由一个SGM来描述，而这一批数据由M个SGM生成，SGM的均值和方差都未知，xi属于那个SGM未知且每个SGM在GMM的比例PIj未知 ， 将来自不同分布的点混在一起，该分布称为混合高斯分布。

P(Xi) = sum ( [PI[j] * SGM(X[i], Miu[j], Sigma[j]) for j in range(M)] ) （一）， 其中 sum(PI[j] for j in range(M) ) = 1,

SGM函数为N(Xi, uj, Sigmaj), Sigmaj是第j个SGM的协方差矩阵。其中

SGM(X[i], Miu[j], Sigma[j])：

d = len(X[i])

a = 1./( (math.sqrt((2 * math.pi)^d * np.linalg.det(Sigma[j]) )) )

b = exp( -0.5 * (X[i] - Miu[j]).T * np.linalg.inv(Sigma[j]) * (X[i] - Miu[j]) )

1./( (math.sqrt((2 * math.pi)^d) * np.linalg.det(C)) ) * exp(-0.5 * (x-u).T * np.linalg.inv(C) * (x-u)) (二)

其中式一表示的是 Xi 在所有SGM下的概率之和， 式二表示 Xi在第j个SGM下的概率密度函数

令lamda[j] ＝ （PI[j], Miu[j], Sigma[j]）， GMM 共有 M 个SGM, 现在我们需要通过样本集X来训练GMM所有的参数Lamda[0:M]。

我们要求样本X的概率最大来训练参数，用极大似然估计法来求，设极大似然函数为：

L(X, Lamda) :

for i in range(N):

for j in range(M):

temp += PI[j] * SGM(X[i] , Miu[j], Sigma[j])

loss += math.log10(temp)

return loss

我们有了目标函数，有了要训练的参数， 有了样本数据，下面我们看是训练GMM参数：

第一步

初始值的确定：

方案一： 协方差矩阵为单位矩阵， PI[j] = 1./M ; Miu[j]为随机数。 (本实验取方案一)

方案二： 用K-mean算法对样本进行巨雷，利用各类样本的均值作为Miu, 各类样本的协方差作为Sigma ， 各类样本占样本总数的比例为PI.

训练的过程用EM(Expectation Maximization )算法来训练参数。

第二步

E step：

这里由于SGM的值可能比较小，为了使值能相对均匀的分布，我们将Xi在第j个高斯分布下的概率定义为：

Gamma(i, j) = PI[j] * SGM(X[i], Miu[j], Sigma[j])

Gamma(i, j) = Gamma(i, j) / sum(Gamma(i, j) for j in range(M))

由此我们可以得到所有点在所有SGM下的pdf, 构成一个N * M 的matrix:

j = 1, i = 1…N , 得到N * 1 的vector

再令j = 1..M 得到N ＊ M的matrix, 矩阵的每一列为所有点在该SGM下的PDF。

sum = zeros(M)

for i in range(N):

for j in range(M):

sum[i] += PI[j] * SGM(X[i], Miu[j], Sigma[j])

pdf_matrix = zeroslike((N,M))

for i in range(N):

for j in range(M):

pdf_matrix(i, j) = PI[j] * SGM(X[i], Miu[j], Sigma[j]) / sum[i] #有些语法错误是因为这个markdown还不太会用

loss ＝ sum(pdf_matrix)

第三步：

M step, 更新参数，直到收敛

PI’[j] = sum( Gamma(i,j) for i in range(N) ) / N

Miu’[j] = sum( Gamma(i,j) * X[i] for i in range(N) ) / sum( Gamma(i,j) for i in range(N) ) #用PI’[j]

Sigma’[j] = sum( Gamma(i,j) for i in range(N) ) / sum( Gamma(i,j) * (X[i]-Miu[j])*(X[i]-Miu[j]).T for i in range(N) ) #用PI’[j], Miu’[j]

收敛条件： 不断的迭代E step and M step, 重复更新Lamda, 直到L(X, lamba) - L(X, lamba’) < delta , 其中delta = 10^(-5)

完整代码上面的链接里边也有， 不过貌似有点问题，有时间来填坑