混合高斯模型（matlab）
阿新 • • 發佈：2019-02-10
推薦部落格：http://blog.csdn.net/crzy_sparrow/article/details/7413019
背景模型有很多種，其中很多方法對光照的的突變和其它因素的適應能力不夠，而高斯混合模型是最為成功的一種背景建模方法。
高斯背景模型是由Stauffer等人提出的經典的自適應混合高斯背景提取方法，是一種基於背景建模的方法，它是根據視訊中的每個畫素在時域上的分佈情況來構建各個畫素的顏色分佈模型，依次來達到背景建模的目的。混合高斯背景模型是有限個高斯函式的加權和，它能描述畫素的多峰狀態，適用於對光照漸變、樹木搖擺等複雜背景進行準確建模。此後經過很多研究人員的不斷改進，該方法目前已經成為比較常用的背景提取方法。
混合高斯模型的原理說白了就一句話：任意形狀的概率分佈都可以用多個高斯分佈函式去近似。換句話說，高斯混合模型（GMM）可以平滑任意形狀的概率分佈。其引數求解方法一般使用極大似然估計求解，但使用極大似然估計法往往不能獲得完整資料（比如樣本已知，但樣本類別（屬於哪個高斯分佈）未知），於是出現了EM（最大期望值）求解方法。
雖然上面說的簡單，但是混合高斯模型和EM求解的理論還是比較複雜的，我把我所找到的我認為能夠快速掌握高斯混合模型的資料打包到了附件中，大家可以去下載，瞭解混合高斯模型以及EM的完整推導過程。
附件下載地址：
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 matrix
pPi = zeros(1, K);%1-by-K matrix
pSigma = zeros(D, D, K);%
% hard assign x to each centroids
distmat = 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-K
2*X*pMiu';%計算每個點到K箇中心的距離
[~, labels] = min(distmat, [], 2);%找到離X最近的pMiu，[C,I] labels代表這個最小值是從那列選出來的
for k=1:K
Xk = X(labels == k, : );% Xk是所有被歸到K類的X向量構成的矩陣
pPi(k) = size(Xk, 1)/N;% 數一數幾個歸到K類的
pSigma(:, :, k) = cov(Xk); %???計算協方差矩陣，D-by-D matrix,最小方差無偏估計
end
endpPi 為每類的比例數目，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:K
Xshift = X-repmat(pMiu(k, : ), N, 1);%x-u
lemda=1e-5;
conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%這裡處理singular問題，為協方差矩陣加上一個很小lemda*I
inv_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 裡得到的值
end
end
這個N極為公式（1）中的N
（3）計算和更新 pPi，pMiu ，pSigma ；
Lprev = -inf;
while true
Px = calc_prob();%計算N(x|mu,sigma)
% new value for pGamma
pGamma = Px .* repmat(pPi, N, 1);%估計 gamma 是個N*K的矩陣
pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);%每個樣本又第k類聚類，又稱componen的生成概率
% new value for parameters of each Component
Nk = sum(pGamma, 1);%N_K
pMiu = diag(1./Nk) * pGamma' * X; % 數字 *( K-by-N * N-by-D)加個括號有助理解
pPi = Nk/N;
for kk = 1:K
Xshift = X-repmat(pMiu(kk, : ), N, 1);%x-u
pSigma(:, :, kk) = (Xshift' * ...
(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%更新sigma
end
% check for convergence
L = sum(log(Px*pPi'));
if L-Lprev < threshold
break;
end
Lprev = L;
end
上式中的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否則返回false
K=K_or_centroids;
rndp = randperm(N);%1到n這些數隨機打亂得到的一個數字序列
centroids = X(rndp(1:K),:);
else
K = size(K_or_centroids,1);%行數
centroids = K_or_centroids;
end
%初始化
[pMiu pPi pSigma] = init_params();
Lprev = -inf;
while true
Px = calc_prob();%計算N(x|mu,sigma)
% new value for pGamma
pGamma = Px .* repmat(pPi, N, 1);%估計 gamma 是個N*K的矩陣
pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);%%求每個樣本由第K個聚類，也叫“component“生成的概率
% new value for parameters of each Component
Nk = sum(pGamma, 1);%N_K
pMiu = diag(1./Nk) * pGamma' * X; % 數字 *( K-by-N * N-by-D)加個括號有助理解
pPi = Nk/N;
for kk = 1:K
Xshift = X-repmat(pMiu(kk, : ), N, 1);%x-u
pSigma(:, :, kk) = (Xshift' * ...
(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%更新sigma
end
% check for convergence
L = sum(log(Px*pPi'));
if L-Lprev < threshold
break;
end
Lprev = L;
end
if nargout == 1
varargout = {Px};
else
model = [];
model.Miu = pMiu;
model.Sigma = pSigma;
model.Pi = pPi;
varargout = {pGamma, model};%注意！！！！！這裡和大神程式碼不同，他返回的是px，而我是 pGamma
end
function [pMiu pPi pSigma] = init_params()%初始化引數
pMiu = centroids;% K-by-D matrix
pPi = zeros(1, K);%1-by-K matrix
pSigma = zeros(D, D, K);%
% hard assign x to each centroids
distmat = 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-K
2*X*pMiu';%計算每個點到K箇中心的距離
[~, labels] = min(distmat, [], 2);%找到離X最近的pMiu，[C,I] labels代表這個最小值是從那列選出來的
for k=1:K
Xk = X(labels == k, : );% Xk是所有被歸到K類的X向量構成的矩陣
pPi(k) = size(Xk, 1)/N;% 數一數幾個歸到K類的
pSigma(:, :, k) = cov(Xk); %???計算協方差矩陣，D-by-D matrix,最小方差無偏估計
end
end
function Px = calc_prob()
Px = zeros(N, K);
for k = 1:K
Xshift = X-repmat(pMiu(k, : ), N, 1);%x-u
lemda=1e-5;
conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%這裡處理singular問題，為協方差矩陣加上一個很小lemda*I
inv_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 裡得到的值
end
end
end
%repmat 通過拓展向量到矩陣
%inv 求逆
%min 求矩陣最小值，可以返回標籤
%X(labels == k, : ) 對行做篩選
% size(Xk, 1) 求矩陣的長或寬
%scatter 對二維向量繪圖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;
else
K_number = size(K, 1);
end
for k=1:K_number
color=[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;%關掉座標系顯示
else
if size(X,2)==3 %3維的情況
figure('color','w');%把背景改成白的
scatter3(X(:,1),X(:,2),X(:,3),30,z,'filled')
end
end
