05 EM算法 - 高斯混合模型 - GMM

04 EM算法 - EM算法收敛证明

GMM (Gaussian Mixture Model, 高斯混合模型)是指该算法由多个高斯模型线性叠加混合而成。每个高斯模型称之为component。

GMM算法 描述的是数据的本身存在的一种分布，即样本特征属性的分布，和预测值Y无关。显然GMM算法是无监督的算法，常用于聚类应用中，component的个数就可以认为是类别的数量。

回到昨天说的例子：随机选择1000名用户，测量用户的身高；若样本中存在男性和女性，身高分别服从高斯分布N(μ1,σ1)和N(μ2,σ2)的分布，试估计参数:μ1,σ1,μ2,σ2；

1、如果明确的知道样本的情况(即男性和女性数据是分开的)，那么我们使用极大似然估计来估计这个参数值。

2、如果样本是混合而成的，不能明确的区分开，那么就没法直接使用极大似然估计来进行参数的估计。

我们可以认为当前的1000条数据组成的集X，是由两个高斯分布叠加而成的(男性的分布和女性的分布)。

如果能找到一种办法把每一个高斯分布对应的参数π、 μ、σ求出来，那么对应的模型就求解出来了。

如果模型求解出来后，如何对数据进行聚类？

这个公式求出来的分别是男性和女性身高分布的概率密度，如果把π、 μ、σ都求出来，以后我们可以构建出一个 能够根据样本特征 计算出样本属于男性或女性的可能性。

实际做样本分类的时候，我们把样本X的特征x1~xn分别代入两个公式中，求出来的两个结果分别是：样本X的性别是男、是女的可能性。如果是男的可能性大于是女的可能性，我们就把样本X归入男性的分类。

假定 GMM 由k个Gaussian分布线性叠加而成，那么概率密度函数如下：

分析第1个等式：

p(x)： 概率密度函数，k个Gaussian分布线性叠加而成的概率密度函数。

∑p(k)p(x|k)： k个某种模型叠加的概率密度函数。

p(k)： 每个模型占的权重，即上面提到的π。

p(x|k)： 给定类别k后，对应的x的概率密度函数。

分析第2个等式： 目标 - 将公式写成高斯分布的样子。

π k ： 即p(k)

p(x;μ k ,∑ k )： 多元高斯(正态)分布。有了观测数据x后，在 给定了条件 下的高斯分布。这个 条件 是 1、第k个分类的均值μ k ; 2、第k个分类的方差∑ k ;

深入分析p(x;μ k ,∑ k )的参数：

如果样本有n个特征，所有的特征x1~xn一起服从一个多元的高斯分布(正态分布)，所有特征的均值应该是一个向量 (μ 1 ~μ n )；

μ k ： 第k个分类的情况下(第k个高斯分布的情况下对应的每一列的均值)；μ k = (μ k1 ~μ kn )

∑ k ： 协方差矩阵(对称阵)。现在有n个特征，协方差矩阵是一个n×n的矩阵。现在我们要算的是：

cov(x1,x1)，cov(x1,x2)，，cov(x1,xn)

cov(x2,x1)，cov(x2,x2)，，cov(x2,xn)

cov(xn,x1)，cov(x1,x2)，，cov(xn,xn)

其中， 对角线 cov(x1,x1)、cov(x2,x2)， ，cov(xn,xn)中，x1和x1的协方差 = x1的方差；即cov(x1,x1) = var(x1)；所以 对角线上两个特征的协方差 = 对应的特征的方差。

协方差 （Covariance）在 概率论 和 统计学 中用于衡量两个变量的总体 误差 。而 方差 是协方差的一种特殊情况，即当两个变量是相同的情况。

协方差表示的是两个变量的总体的 误差 ，这与只表示一个变量误差的 方差 不同。 如果两个 变量 的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。 如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。

理解了公式后，再来看看公式在图像上是如何体现的：

如果样本X只有一个特征x1，在二维的坐标系上的表示出来。特征x1是由n个单变量样本的高斯分布叠加而成的。向量x1 k = ∑ k (x1 (1) ,x1 (2) ,~,x1 (n) )，如k=(男、女)，累加男性分类下的特征高斯分布和女性分类下的高斯分布；

图中 红色曲线 表示原有数据的分布情况，我认为这个原有数据是由多个比较的高斯分布叠加而成的， 蓝色曲线 表示单个单个高斯分布的分布情况。向量x1 = (x1 (1) ,x1 (2) ,~,x1 (n) )；

PS： 蓝1+蓝2=红 体现的就是公式 p(x) = ∑πp(x;μ,∑k)；

在得知数据的特征 x=(x1~xn) 后，如果我们想把数据合理得聚类到一个分类中，我们该如何去计算呢？

既然我已经得到了k个高斯分布对应的概率密度函数(现在设k=3，共3个分类)，将当前特征的x=(x1~xn)代入我们的概率密度函数： p(x) = ∑πp(x;μ,∑k)；

我们分别计算p(蓝1)、p(蓝2)、p(蓝3)，蓝色三条线各对应k分类中的一个，哪个数大，我认为当前的样本该分到哪一类。

GMM算法的两个前提：

1、数据服从高斯分布；

2、我们人为定义了分类个数k。

问：我们人为假定了高斯分布的分类个数k，就类似于我们聚簇时分的聚簇中心个数一样。参数π、μ、σ该如何求出来

答：和K-Means算法一样，我们可以用 EM算法 来求解这个问题。 GMM也满足EM算法的聚类思想，首先人为得定义了聚类的个数k，从数据特征X中发掘潜在关系的一种模型。而且我还默认数据是服从多个高斯分布的。

GMM算法中的隐含条件是：第k个模型占的权重 - 、 第k个高斯分布的情况下对应的每一列的均值 - 、协方差矩阵 cov(xi,xj) - ；因为本质上我们是知道数据原有的分类状况的，只是无法观测到隐含在数据中的这些特性，使用EM的思想可以迭代得求解出这些隐含变量。

对联合概率密度函数求对数似然函数：

对联合概率密度函数求对数后，原本 连乘 的最大似然估计变成了 连加 的函数状态。

EM算法求解 - E步：

套用公式后，我们可以假定隐含变量z的分布：Q(z (i) = j)；

我们认为分布wj (i) = 第i个观测值对应的隐含分类第z (i) 类； = 以(看不见的参数π、μ、∑)为参数的情况下，输入第i观测值的特征x后得到的分类z (i) 类；

EM算法求解 - M步：

M步第1行就是上一章通过化简找到 下界 的那个函数：

1、对均值求偏导:

2、对方差求偏导:

3、对概率使用拉格朗日乘子法求解:

06 EM算法 - 案例一 - EM分类初识及GMM算法实现

　　混合高斯模型：

　　运动检测的一般方法

　　目前，运动物体检测的问题主要分为两类，摄像机固定和摄像机运动。对于摄像机运动的运动物体检测问题，比较著名的解决方案是光流法，通过求解偏微分方程求的图像序列的光流场，从而预测摄像机的运动状态。对于摄像机固定的情形，当然也可以用光流法，但是由于光流法的复杂性，往往难以实时的计算，所以我采用高斯背景模型。因为，在摄像机固定的情况下，背景的变化是缓慢的，而且大都是光照，风等等的影响，通过对背景建模，对一幅给定图像分离前景和背景，一般来说，前景就是运动物体，从而达到运动物体检测的目的。

　　单分布高斯背景模型

　　单分布高斯背景模型认为，对一个背景图像，特定像素亮度的分布满足高斯分布，即对背景图像B，(x,y)点的亮度满足：　IB(x,y) ~ N(u,d)　这样我们的背景模型的每个象素属性包括两个参数：平均值u 和 方差d。　对于一幅给定的图像G，如果 Exp(-(IG(x,y)-u(x,y))^2/(2d^2)) > T，认为(x,y)是背景点，反之是前景点。　同时，随着时间的变化，背景图像也会发生缓慢的变化，这时我们要不断更新每个象素点的参数　u(t+1,x,y) = au(t,x,y) + (1-a)I(x,y)　这里，a称为更新参数，表示背景变化的速度，一般情况下，我们不更新d（实验中发现更不更新d，效果变化不大）。　高斯混合模型是用于背景提取的方法，OpenCV的cvaux中cvbgfg_gaussmixcpp文件根据文献An improved adaptive background mixture model for real-time tracking with shadow中提供的方法编写了高斯混合模型函数。其中定义了CvGaussBGModel类用于存放高斯混合模型的各个参数。我用OpenCV使用高斯混合模型函数分以下几步：　1。需要用到icvUpdateGaussianBGModel，icvReleaseGaussianBGModel两个函数，但是源程序中将这两个函数定义为内部函数，需要做一些修改，首先将cvbgfg_gaussmixcpp中前面两个函数的声明static void CV_CDECL icvReleaseGaussianBGModel( CvGaussBGModel bg_model );　static int CV_CDECL icvUpdateGaussianBGModel( IplImage curr_frame, CvGaussBGModel bg_model );两行代码注释掉。然后在cvbgfg_gaussmixcpp中间部分两个函数的定义部分，函数头static int和static void改成CV_IMPL int 和CV_IMPL void。最后在cvauxh文件中CVAPI(CvBGStatModel) cvCreateGaussianBGModel( IplImage first_frame,　CvGaussBGStatModelParams parameters CV_DEFAULT(NULL));这句后面加上以下两句CVAPI(void) icvReleaseGaussianBGModel( CvGaussBGModel bg_model );　CVAPI(int) icvUpdateGaussianBGModel( IplImage curr_frame, CvGaussBGModel bg_model );　程序修改完毕，点rebuild all，全部重新编译。　2。在程序初始化部分定义高斯混合模型参数CvGaussBGModel bg_model=NULL;在读取第一帧图像（背景图像）时，进行高斯背景建模bg_model = (CvGaussBGModel)cvCreateGaussianBGModel(image, 0);image可以是灰度图象也可以是彩色图像。接下来再读取当前帧时，更新高斯模型　regioncount=icvUpdateGaussianBGModel(currframe, bg_model );regioncount的含义我不确定，我理解是代表背景中不同颜色区域的个数，这个参数我没有用到，它只是icvUpdateGaussianBGModel函数的返回值。　3。现在bg_model已经保存了经过高斯混合模型分类后的结果，bg_model->background保存了背景图像，bg_model->foreground保存了前景图像。　include <stdioh>　#include <cvh>　#include <cxcoreh>　#include <highguih>　#include <cvauxh>//必须引此头文件　int main( int argc, char argv )　{　IplImage pFrame = NULL;　IplImage pFrImg = NULL;　IplImage pBkImg = NULL;　CvCapture pCapture = NULL;　int nFrmNum = 0;　cvNamedWindow("video", 1);　cvNamedWindow("background",1);　cvNamedWindow("foreground",1);　cvMoveWindow("video", 30, 0);　cvMoveWindow("background", 360, 0);　cvMoveWindow("foreground", 690, 0);　if( argc > 2 )　{　fprintf(stderr, "Usage: bkgrd [video_file_name]\n");　return -1;　}　//打开视频文件　if(argc == 2)　if( !(pCapture = cvCaptureFromFile(argv[1])))　{　fprintf(stderr, "Can not open video file %s\n", argv[1]);　return -2;　}　//打开摄像头　if (argc == 1)　if( !(pCapture = cvCaptureFromCAM(-1)))　{　fprintf(stderr, "Can not open camera\n");　return -2;　}　//初始化高斯混合模型参数　CvGaussBGModel bg_model=NULL;　while(pFrame = cvQueryFrame( pCapture ))　{　nFrmNum++;　if(nFrmNum == 1)　{　pBkImg = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U,3);　pFrImg = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U,1);　//高斯背景建模，pFrame可以是多通道图像也可以是单通道图像　//cvCreateGaussianBGModel函数返回值为CvBGStatModel，　//需要强制转换成CvGaussBGModel　bg_model = (CvGaussBGModel)cvCreateGaussianBGModel(pFrame, 0);　}　else　{　//更新高斯模型　cvUpdateBGStatModel(pFrame, (CvBGStatModel )bg_model );　//pFrImg为前景图像，只能为单通道　//pBkImg为背景图像，可以为单通道或与pFrame通道数相同　cvCopy(bg_model->foreground,pFrImg,0);　cvCopy(bg_model->background,pBkImg,0);　//把图像正过来　pBkImg->origin=1;　pFrImg->origin=1;　cvShowImage("video", pFrame);　cvShowImage("background", pBkImg);　cvShowImage("foreground", pFrImg);　if( cvWaitKey(2) >= 0 )　break;　}　}　//释放高斯模型参数占用内存　cvReleaseBGStatModel((CvBGStatModel)&bg_model);　cvDestroyWindow("video");　cvDestroyWindow("background");　cvDestroyWindow("foreground");　cvReleaseImage(&pFrImg);　cvReleaseImage(&pBkImg);　cvReleaseCapture(&pCapture);　return 0;　}

线性核、多项式核、高斯核。

特征维数高选择线性核

样本数量可观、特征少选择高斯核（非线性核）

样本数量非常多选择线性核（避免造成庞大的计算量）

准确率。

因为将泰勒展开式代入高斯核,将会得到一个无穷维度的映射。

1）LR是参数模型，SVM是非参数模型。2）从目标函数来看，区别在于逻辑回归采用的是logistical loss，SVM采用的是hinge loss这两个损失函数的目的都是增加对分类影响较大的数据点的权重，减少与分类关系较小的数据点的权重。3）SVM的处理方法是只考虑support vectors，也就是和分类最相关的少数点，去学习分类器。而逻辑回归通过非线性映射，大大减小了离分类平面较远的点的权重，相对提升了与分类最相关的数据点的权重。4）逻辑回归相对来说模型更简单，好理解，特别是大规模线性分类时比较方便。而SVM的理解和优化相对来说复杂一些，SVM转化为对偶问题后,分类只需要计算与少数几个支持向量的距离,这个在进行复杂核函数计算时优势很明显,能够大大简化模型和计算。5）logic 能做的 svm能做，但可能在准确率上有问题，svm能做的logic有的做不了。

假如特征维数较多，可以选择SVM模型，如果样本数量较大可以选择LR模型，但是LR模型需要进行数据预处理；假如缺失值较多可以选择决策树。选定完模型后，相应的目标函数就确定了。还可以在考虑正负样例比比，通过上下集采样平衡正负样例比。

ID3决策树优先选择信息增益大的属性来对样本进行划分,但是这样的分裂节点方法有一个很大的缺点,当一个属性可取值数目较多时,可能在这个属性对应值下的样本只有一个或者很少个,此时它的信息增益将很高,ID3会认为这个属性很适合划分,但实际情况下叫多属性的取值会使模型的泛化能力较差,所以C45不采用信息增益作为划分依据,而是采用信息增益率作为划分依据。但是仍不能完全解决以上问题,而是有所改善,这个时候引入了CART树,它使用gini系数作为节点的分裂依据。

基本思想就是连续特征离散化。现将特征排序，m个样本的连续特征A有m个，从小到大排列为a1,a2,,am,则C45取相邻两样本值的平均数，一共取得m-1个划分点，其中第i个划分点 表示为

要注意的是，与离散属性不同的是，如果当前节点为连续属性，则该属性后面还可以参与子节点的产生选择过程。