混合高斯模型（matlab）

推荐博客：http://blog.csdn.net/crzy_sparrow/article/details/7413019

背景模型有很多种，其中很多方法对光照的的突变和其它因素的适应能力不够，而高斯混合模型是最为成功的一种背景建模方法。

高斯背景模型是由Stauffer等人提出的经典的自适应混合高斯背景提取方法，是一种基于背景建模的方法，它是根据视频中的每个像素在时域上的分布情况来构建各个像素的颜色分布模型，依次来达到背景建模的目的。混合高斯背景模型是有限个高斯函数的加权和，它能描述像素的多峰状态，适用于对光照渐变、树木摇摆等复杂背景进行准确建模。此后经过很多研究人员的不断改进，该方法目前已经成为比较常用的背景提取方法。

混合高斯模型的原理说白了就一句话：任意形状的概率分布都可以用多个高斯分布函数去近似。换句话说，高斯混合模型（GMM）可以平滑任意形状的概率分布。其参数求解方法一般使用极大似然估计求解，但使用极大似然估计法往往不能获得完整数据（比如样本已知，但样本类别（属于哪个高斯分布）未知），于是出现了EM（最大期望值）求解方法。

虽然上面说的简单，但是混合高斯模型和EM求解的理论还是比较复杂的，我把我所找到的我认为能够快速掌握高斯混合模型的资料打包到了附件中，大家可以去下载，了解混合高斯模型以及EM的完整推导过程。

附件下载地址：

http://download.csdn.net/detail/crzy_sparrow/4187994

1）任意数据分布可用高斯混合模型（M个单高斯）表示（（1）式）

其中

2）N个样本集X的log似然函数如下：

（1）初始参数

（ 2）E步（求期望）

求取：

其实，这里最贴切的式子应该是log似然函数的期望表达式

事实上，3）中参数的求取也是用log似然函数的期望表达式求偏导等于0解得的简化后的式子，而这些式子至于（7）式有关，于是E步可以只求（7）式，以此简化计算，不需要每次都求偏导了。

（ 3）M步（最大化步骤）

（4）截止条件

不断迭代EM步骤，更新参数，直到似然函数前后差值小于一个阈值，或者参数前后之间的差（一般选择欧式距离）小于某个阈值，终止迭代，这里选择前者。

MATLAB对应程序如下：

（1）初始参数

    function [pMiu pPi pSigma] = init_params()%初始化参数pMiu = centroids;% K-by-D matrixpPi = zeros(1, K);%1-by-K matrixpSigma = zeros(D, D, K);%% hard assign x to each centroidsdistmat = repmat(sum(X.*X, 2), 1, K) + ... % X is a N-by-D data matrix.repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...% X->K列 U->N行 XU^T is N-by-K2*X*pMiu';%计算每个点到K个中心的距离[~, labels] = min(distmat, [], 2);%找到离X最近的pMiu，[C,I] labels代表这个最小值是从那列选出来的for k=1:KXk = X(labels == k, : );% Xk是所有被归到K类的X向量构成的矩阵pPi(k) = size(Xk, 1)/N;% 数一数几个归到K类的pSigma(:, :, k) = cov(Xk); %???计算协方差矩阵，D-by-D matrix,最小方差无偏估计endend

其中， pMiu 为k类的类中心，K*D,D为样本X的列数，样本X中一列极为一个样本，（比如空间中的点x，y，z）

pPi 为每类的比例数目，1*K,比如总共是个样本点，其中属于第一类的占0.2，属于低二类的占0.1，属于第k类的……

pSigma为每类的协方差矩阵，D*D*K,表示第k类的协方差矩阵为pSigma（:,:,k）

</pre></p><p><span style="font-size: 16px; line-height: 26px; font-family: Arial; background-color: rgb(255, 255, 255);"><span style="color:#ff0000;"><span style="color: rgb(255, 0, 0); font-family: Arial; font-size: 16px; line-height: 26px;"><span style="color: rgb(51, 51, 51); font-family: Arial; font-size: 16px; line-height: 26px;">       </span><span style="font-family: Arial; font-size: 16px; line-height: 26px; color: rgb(255, 0, 0);">（</span><span style="font-family: Arial; line-height: 26px; color: rgb(255, 0, 0); font-size: 18px;">2）计算N</span></span></span></span><span style="color: rgb(255, 0, 0); font-family: Arial; font-size: 18px; line-height: 26px;">         </span></p><p><span style="font-family: Arial; font-size: 18px; line-height: 26px;"></span><pre name="code" class="cpp">    function Px = calc_prob()Px = zeros(N, K);for k = 1:KXshift = X-repmat(pMiu(k, : ), N, 1);%x-ulemda=1e-5;conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%这里处理singular问题，为协方差矩阵加上一个很小lemda*Iinv_pSigma = inv(conv);%协方差的逆tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);%(X-U_k)sigma.*(X-U_k),tmp是个N*1的向量coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));%前面的参数Px(:, k) = coef * exp(-0.5*tmp);%把数据点 x 带入到 Gaussian model 里得到的值endend

这个N极为公式（1）中的N

（3）计算和更新 pPi，pMiu ， pSigma ；

Lprev = -inf;while truePx = calc_prob();%计算N(x|mu,sigma)% new value for pGammapGamma = Px .* repmat(pPi, N, 1);%估计 gamma 是个N*K的矩阵pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);%每个样本又第k类聚类，又称componen的生成概率% new value for parameters of each ComponentNk = sum(pGamma, 1);%N_KpMiu = diag(1./Nk) * pGamma' * X; % 数字 *( K-by-N * N-by-D)加个括号有助理解pPi = Nk/N;for kk = 1:KXshift = X-repmat(pMiu(kk, : ), N, 1);%x-upSigma(:, :, kk) = (Xshift' * ...(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%更新sigmaend% check for convergenceL = sum(log(Px*pPi'));if L-Lprev < thresholdbreak;endLprev = L;end

上式中的Px即为公式（1）

上式中的pGamma即为公式（7）

上式中的pPi，pMiu ， pSigma 即为按照M步进行更新。

上式中的停止准则按照公式（4）来计算

所有的程序如下：

gmm.m文件

function varargout = gmm( X,K_or_centroids )
nargout=0;
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
% ============================================================
%转载自http://blog.pluskid.org/?p=39
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
%  - X: N-by-D data matrix.%需要注意的是这里的X包括了全部
%  - K_OR_CENTROIDS: either K indicating the number of
%       components or a K-by-D matrix indicating the
%       choosing of the initial K centroids.
%
%  - PX: N-by-K matrix indicating the probability of each
%       component generating each point.
%  - MODEL: a structure containing the parameters for a GMM:
%       MODEL.Miu: a K-by-D matrix.
%       MODEL.Sigma: a D-by-D-by-K matrix.
%       MODEL.Pi: a 1-by-K vector.
% ============================================================
threshold = 1e-15;
[N,D] = size(X);%行和列
if isscalar(K_or_centroids)%标量返回TRUE否则返回falseK=K_or_centroids;rndp = randperm(N);%1到n这些数随机打乱得到的一个数字序列centroids = X(rndp(1:K),:);
elseK = size(K_or_centroids,1);%行数centroids = K_or_centroids;
end%初始化
[pMiu pPi pSigma] = init_params();Lprev = -inf;while truePx = calc_prob();%计算N(x|mu,sigma)% new value for pGammapGamma = Px .* repmat(pPi, N, 1);%估计 gamma 是个N*K的矩阵pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);%%求每个样本由第K个聚类，也叫“component“生成的概率% new value for parameters of each ComponentNk = sum(pGamma, 1);%N_KpMiu = diag(1./Nk) * pGamma' * X; % 数字 *( K-by-N * N-by-D)加个括号有助理解pPi = Nk/N;for kk = 1:KXshift = X-repmat(pMiu(kk, : ), N, 1);%x-upSigma(:, :, kk) = (Xshift' * ...(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%更新sigmaend% check for convergenceL = sum(log(Px*pPi'));if L-Lprev < thresholdbreak;endLprev = L;endif nargout == 1varargout = {Px};elsemodel = [];model.Miu = pMiu;model.Sigma = pSigma;model.Pi = pPi;varargout = {pGamma, model};%注意！！！！！这里和大神代码不同，他返回的是px，而我是 pGammaendfunction [pMiu pPi pSigma] = init_params()%初始化参数pMiu = centroids;% K-by-D matrixpPi = zeros(1, K);%1-by-K matrixpSigma = zeros(D, D, K);%% hard assign x to each centroidsdistmat = repmat(sum(X.*X, 2), 1, K) + ... % X is a N-by-D data matrix.repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...% X->K列 U->N行 XU^T is N-by-K2*X*pMiu';%计算每个点到K个中心的距离[~, labels] = min(distmat, [], 2);%找到离X最近的pMiu，[C,I] labels代表这个最小值是从那列选出来的for k=1:KXk = X(labels == k, : );% Xk是所有被归到K类的X向量构成的矩阵pPi(k) = size(Xk, 1)/N;% 数一数几个归到K类的pSigma(:, :, k) = cov(Xk); %???计算协方差矩阵，D-by-D matrix,最小方差无偏估计endendfunction Px = calc_prob()Px = zeros(N, K);for k = 1:KXshift = X-repmat(pMiu(k, : ), N, 1);%x-ulemda=1e-5;conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%这里处理singular问题，为协方差矩阵加上一个很小lemda*Iinv_pSigma = inv(conv);%协方差的逆tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);%(X-U_k)sigma.*(X-U_k),tmp是个N*1的向量coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));%前面的参数Px(:, k) = coef * exp(-0.5*tmp);%把数据点 x 带入到 Gaussian model 里得到的值endend
end
%repmat 通过拓展向量到矩阵
%inv 求逆
%min 求矩阵最小值，可以返回标签
%X(labels == k, : ) 对行做筛选
% size(Xk, 1) 求矩阵的长或宽
%scatter 对二维向量绘图

hungauss.m文件

X=[...
1 1 1.1;
1 2 1;
2 1 1;
2 2 1;
5 5 2;
5 6 2;
6 5 2;
6 6 2.1]%数据X自己改
K=2
%以上设置数据点
[Px,model]=gmm(X,K);%这里得到结果
[~,belong]=max(Px,[],2)%选概率最大的那个数，输出聚类结果
%以下绘图
z=zeros(size(X,1),3);
if isscalar(K)%获得类别个数K_number=K;
elseK_number = size(K, 1);
endfor k=1:K_numbercolor=[rand rand rand];%建立颜色矩阵，随机给个颜色csize=size(z(belong==k,:),1);%数一数有几行z(belong==k,:)=repmat(color,csize,1);%对属于某一类下的点染色
end
if size(X,2)==2 %二维或三维的可以画一画
figure('color','w');%把背景改成白的
scatter(X(:,1),X(:,2),30,z)
axis off;%关掉坐标系显示
elseif size(X,2)==3 %3维的情况figure('color','w');%把背景改成白的scatter3(X(:,1),X(:,2),X(:,3),30,z,'filled')end
end

运行结果如下：

混合高斯模型（matlab）相关推荐

基于混合高斯模型与帧差法结合的目标跟踪算法matlab仿真
目录一.理论基础二.核心程序三.仿真测试结果一.理论基础目标检测:混合高斯模型与帧差法结合的算法,与单独的混合高斯模型算法作对比,体现前者的优越性 3.要求和结果:对比改进前后的算法,可以非 ...
混合高斯模型介绍以及应用
混合高斯模型 1. 单一的高斯模型(Gaussian single model, GSM) 2. 混合高斯模型(GMM模型) 2.1 混合高斯模型直观上的理解和描述 2.2 极大似然估计(Maximu ...
【转载】混合高斯模型（Mixtures of Gaussians）和EM算法
混合高斯模型(Mixtures of Gaussians)和EM算法这篇讨论使用期望最大化算法(Expectation-Maximization)来进行密度估计(density estimation ...
混合高斯模型（Mixtures of Gaussians）和EM算法
与k-means一样,给定的训练样本是,我们将隐含类别标签用表示.与k-means的硬指定不同,我们首先认为是满足一定的概率分布的,这里我们认为满足多项式分布,,其中,有k个值{1,-,k}可以选取. ...
聚类(1)——混合高斯模型 Gaussian Mixture Model
聚类系列: 聚类(序)----监督学习与无监督学习聚类(1)----混合高斯模型 Gaussian Mixture Model 聚类(2)----层次聚类 Hierarchical Clusteri ...
混合高斯模型_EM算法求解高斯混合模型(GMM)
单维的高斯模型: 求解一维的高斯模型参数时,我们可以使用最大似然法,其对数似然函数的表达式(log-likelyhood)如下: 对均值和方差求偏导可以求的高斯分布中的在混合高斯模型中, 记其对数 ...
python gmm em算法 2维数据_AI大语音（六）——混合高斯模型（GMM）（深度解析）...
1 GMM基础高斯混合模型(GMM)指的是多个高斯分布函数的线性组合,理论上GMM可以拟合出任意类型的分布,通常用于解决同一集合下的数据包含多个不同的分布的情况. 灵魂的拷问:为什么GMM可以拟合出 ...
ASR 混合高斯模型GMM的理解
混合高斯模型(GMM)是使用非常广泛的统计模型,一种非常高调的说法是,混合高斯模型能拟合一切数据.虽然实际还是受到很多限制,比如混合高斯分布数量需要确定等等,不难看出其强大指出.此文包含以下内容: G ...
混合高斯模型（GMM）推导及实现
作者:桂. 时间:2017-03-20 06:20:54 链接:http://www.cnblogs.com/xingshansi/p/6584555.html 前言本文是曲线拟合与分布拟合系列的 ...

混合高斯模型（matlab）

混合高斯模型（matlab）相关推荐

最新文章

热门文章