高斯混合模型

使用單高斯模型來模組化有一些限制，例如，它一定隻有一個衆數，它一定對稱的。舉個例子，如果我們對下面的分布建立單高斯模型，會得到顯然相差很多的模型：

于是，我們引入混合高斯模型（Gaussian Mixture Model，GMM）。高斯混合模型就是多個單高斯模型的和。它的表達能力十分強，任何分布都可以用GMM來表示。例如，在下面這個圖中，彩色的線表示一個一個的單高斯模型，黑色的線是它們的和，一個高斯混合模型：

在小球檢測的栗子中，我們試圖對紅色小球建立單高斯模型（紅和綠這二進制），會發現紅色小球的觀測值不是很符合所建立的模型，如下圖：

此時，如果我們采取高斯混合模型（兩個二進制高斯分布），會發現效果好了很多，如下圖：

下面，我們來詳細地介紹一下高斯混合模型，高斯混合模型的數學形式如下：

其中，

是均值為

，協方差矩陣為

的單高斯模型，

是

的權重系數（

），

是單高斯模型的個數。

前面我們提到了GMM的優點（能夠表示任何分布），當然GMM也有缺點。其一，參數太多了，每一個單高斯模型都有均值、協方差矩陣和權重系數，另外還有單高斯模型的個數

也是其中一個參數，這使得求解GMM十分複雜。其二，有時候所建立的高斯混合模型對觀測值的模組化效果很好，但是其他值可能效果不好，我們稱這個缺點為過度拟合（Overfitting）。

2、求解高斯混合模型：EM算法

在求解單高斯分布的時候，我們用最大似然估計（MLE）的方法得到了理論上的最優解。當我們使用相同的方法試圖求解高斯混合模型的時候，會卡在中間步驟上（具體來說，是單高斯分布求和出現在了對數函數裡面）。索性我們可以用疊代的方法來求解GMM，具體來說，最大期望算法（Expectation Maximization algorith，EM）。

上一節提到，我們想要求解GMM如下：

這其中需要求解的變量很多：

、

。為了簡化問題，我們假定

是預先設定好的，并且每個單高斯分布的權重相等，即假定：

這樣，我們需要确定的就隻剩下每個單高斯分布的參數（

）了。

前面提到，這個模型是沒有解析解（理論最優）的。取而代之，我們采取疊代的方法逼近最優解（大家可以回想一下牛頓法求近似求解方程，或者遺傳算法）。

在牛頓法中，我們需要預先猜測一個初始解，同樣在EM算法中，我們也需要預先猜測一個初始解（

）。

在EM算法中，我們引入一個中間變量

，其中

是單高斯分布的個數，

是觀測值的數目。

直覺地可以這樣了解：在确定了所有單高斯分布之後，可以計算觀測值

發生的機率（分母部分）和第

個單高斯分布

下觀測值

發生的機率（分子部分），這樣

就可以了解為第

個高斯分布對觀測值

發生的機率的貢獻了。如下表所示：

表中最後一行

可以了解為第

個單高斯分布對整個GMM的貢獻。

例如，當GMM中

時（即由兩個單高斯分布組成）：

在這個圖中，對于觀測值

，

是第一個單高斯分布

下

發生的機率，

是第二個單高斯分布

是在整個GMM下

發生的機率；

表示

對

發生機率的貢獻，

發生機率的貢獻。

當我們算出了所有的

之後，我們再反過來用它更新每一個單高斯分布。更新的公式為：

大家可以想一想，對于第

個高斯分布

個觀測值

的貢獻，而

的平均值，用

來更新

就很好了解了，進而再更新協方差矩陣

。

這樣，兩組值互相更新，當相鄰兩個步驟值的差小于事先設定一個門檻值（threshold）時（收斂），我們停止循環，輸出

和

為所求。

值得一提的是，EM算法的結果和初值有很大關系，不同初值會得到不同的解。（想象一下GMM模型其實是多峰的，不同的初值可能最終收斂到不同的峰值，有的初值會收斂到全局最優，有的初值會收斂到局部最優）

實戰：小球檢測

具體任務在這個網址：小球檢測

簡單來說，給你一些圖檔作為訓練集，讓你對黃色小球進行模組化，模組化完成後，讀入最左邊的原始圖檔，要求能像最右邊的圖檔一樣辨別出小球的位置：

（注：我們使用matlab完成所有代碼）

第一步：建立顔色模型

前面幾節我們提到，要建立高斯分布模型，首先要有很多觀測值，是以，我們将在這一步裡面采集很多觀測值，再建立高斯模型。課程裡面提供了19張圖檔作為訓練集供采樣：

我們編寫代碼mark.m，作用為讀入一張一張圖檔，然後手動圈出圖檔中的黃色小球，将圈中部分的像素RGB值存入觀測值。

imagepath = './train';
Samples = [];
for k=1:19
    I = imread(sprintf('%s/%03d.png',imagepath,k)); % read files into RGB
    R = I(:,:,1);
    G = I(:,:,2);
    B = I(:,:,3);
      

</span>%<span style="color: #000000;"> Collect samples 
disp(</span><span style="color: #800000;">''</span><span style="color: #000000;">);
disp(</span><span style="color: #800000;">'</span><span style="color: #800000;">INTRUCTION: Click along the boundary of the ball. Double-click when you get back to the initial point.</span><span style="color: #800000;">'</span><span style="color: #000000;">)
disp(</span><span style="color: #800000;">'</span><span style="color: #800000;">INTRUCTION: You can maximize the window size of the figure for precise clicks.</span><span style="color: #800000;">'</span><span style="color: #000000;">)
figure(</span><span style="color: #800080;">1</span>), mask =<span style="color: #000000;"> roipoly(I); 
figure(</span><span style="color: #800080;">2</span>), imshow(mask); title(<span style="color: #800000;">'</span><span style="color: #800000;">Mask</span><span style="color: #800000;">'</span><span style="color: #000000;">);

sample_ind </span>= find(mask &gt; <span style="color: #800080;">0</span>); % <span style="color: #0000ff;">select</span><span style="color: #000000;"> marked pixels
R </span>=<span style="color: #000000;"> R(sample_ind);
G </span>=<span style="color: #000000;"> G(sample_ind);
B </span>=<span style="color: #000000;"> B(sample_ind);

Samples </span>= [Samples; [R G B]]; %<span style="color: #000000;"> insert selected pixels into samples

disp(</span><span style="color: #800000;">'</span><span style="color: #800000;">INTRUCTION: Press any key to continue. (Ctrl+c to exit)</span><span style="color: #800000;">'</span><span style="color: #000000;">)
pause
           

end

save('Samples.mat', 'Samples'); % save the samples to file

figure,

scatter3(Samples(:,1),Samples(:,2),Samples(:,3),'.');

title('Pixel Color Distribubtion');

xlabel('Red');

ylabel('Green');

zlabel('Blue');

結果我們采集了6699個觀測值，散點圖如下：

我們選用多元高斯分布作為我們的模型，應用第四節中的結論，編寫代碼coefficient.m求出平均值

和協方差矩陣

并存儲到本地以供後續使用。

mu = mean(Samples); % mu
      

sig=zeros(3,3); % sigma

for i=1:N

data=double(Samples(i,:));

sig=sig+(data-mu)'*(data-mu)/N;

% save the coefficients to files

save('mu.mat', 'mu');

save('sig.mat', 'sig');

接下來，我們需要編寫函數detecBall，以圖像為參數，以劃分好的黑白圖像和小球中心為輸出，如下，I是輸入的圖像，segI是劃分好的黑白圖像，loc是球的中心點。detecBall.m的具體代碼如下：

function [segI, loc] = detectBall(I)
load('mu.mat');
load('sig.mat');
      

Id=double(I); % array in size of (row, col, 3)

row=size(Id,1);

col=size(Id,2);

% x_i - mu

for i=1:3

Id(:,:,i) = Id(:,:,i) - mu(i);

% reshape the image to a matrix in size of (rowcol, 3)

Id=reshape(Id,rowcol,3);

% calc possibility using gaussian distribution

% be careful of using * and .* in matrix multiply

Id = exp(-0.5* sum(Idinv(sig).Id, 2)) ./ (2*pi)^1.5 ./ det(sig)^0.5;

% reshape back, now each pixels is with the value of the possibility

Id=reshape(Id,row,col);

% set threshold

thr=8e-06;

% binary image about if each pixel 'is ball'

Id=Id>thr;

% find the biggest ball area

segI = false(size(Id));

CC = bwconncomp(Id);

numPixels = cellfun(@numel,CC.PixelIdxList);

[biggest,idx] = max(numPixels);

segI(CC.PixelIdxList{idx}) = true;

%figure, imshow(segI); hold on;

S = regionprops(CC,'Centroid');

loc = S(idx).Centroid;

%plot(loc(1), loc(2),'r+');

需要注意的一點是，求像素屬于小球的機率的時候，盡量用矩陣運算，而不要用循環，否則會效率低下，請讀者仔細揣摩計算技巧，靈活運用矩陣的

運算。

接下來，我們就可以測試小球檢測的效果啦！

我們編寫test.m測試訓練集的19張圖檔：

imagepath = './train';
for k=1:19
    I = imread(sprintf('%s/%03d.png',imagepath,k));
    [segI, loc] = detectBall(I);
    figure, imshow(segI); hold on; 
    plot(loc(1), loc(2), '+b','MarkerSize',7); 
    disp('Press any key to continue. (Ctrl+c to exit)')
    pause
end