高斯混合模型

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

 

於是，我們引入混合高斯模型（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);
    
    % Collect samples 
    disp('');
    disp('INTRUCTION: Click along the boundary of the ball. Double-click when you get back to the initial point.')
    disp('INTRUCTION: You can maximize the window size of the figure for precise clicks.')
    figure(1), mask = roipoly(I); 
    figure(2), imshow(mask); title('Mask');

    sample_ind = find(mask > 0); % select marked pixels
    R = R(sample_ind);
    G = G(sample_ind);
    B = B(sample_ind);
    
    Samples = [Samples; [R G B]]; % insert selected pixels into samples
    
    disp('INTRUCTION: Press any key to continue. (Ctrl+c to exit)')
    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;
end

% 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);
end

% reshape the image to a matrix in size of (row*col, 3)
Id=reshape(Id,row*col,3);
 
% calc possibility using gaussian distribution
% be careful of using * and .* in matrix multiply
Id = exp(-0.5* sum(Id*inv(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

 

效果如下：