從高斯分布、機器人誤差、EM演算法到小球檢測

01-24

Coursera上的課程（Robotics: Estimation and Learning），感覺講得特別棒，寫下自己的理解。

高斯分布被廣泛應用於對機器人誤差的建模。在這篇筆記中，我們將會：

介紹如何使用一元高斯分布、多元高斯分布和高斯混合模型對機器人誤差進行建模。
介紹求解這些高斯分布的演算法。
以小球檢測這一實際應用來實踐我們的模型和演算法。

1. 一元高斯分布

在這一節我們將介紹如何使用一元高斯分布對機器人誤差進行建模。

在進入正題之前，我們需要首先了解為什麼學習高斯分布？什麼使得高斯分布有效並且重要？為什麼使用高斯分布對雜訊和不確定性進行建模（而不是其他概率分布模型，比如均勻分布）？

高斯分布使得只需要兩個參數就能確定一個連續的概率分布。這兩個參數分別是均值和方差，這很簡潔。這也使得高斯分布易於計算和推斷。
高斯分布有一些較好的數學性質。例如，兩個高斯分布的積還是高斯分布。這使得你在對高斯分布做運算的時候不需要考慮其他的概率分布。
理論上來說，根據中心極限定理，大量相互獨立的隨機變數，其均值的分布以高斯分布為極限。而雜訊和不確定性就是包含大量相互獨立的隨機變數的數據，用高斯分布對它們進行建模是不二選擇。

說了這麼多，我們來舉個栗子吧，怎樣使用高斯分布來對一個目標顏色進行建模？

下圖是足球機器人通過頭頂攝像機拍攝到的照片。

顯然，圖片中有兩種顏色的球，黃球和紅球。我們設想這個足球機器人想檢測圖片中的黃球。對於人類來說，一眼就能分辨出黃球在哪，但對一個機器人來說這並不簡單。機器人讀入的圖片是一個一個像素的數值，它需要建立從像素的數值到「黃色」或者「紅色」的映射。不妨讓我們來看看「黃球」裡面的像素數值是怎樣的。我們記錄黃球的每個像素的色相值（Hue，HSV模型中表示色彩的基本屬性，就是平常所說的顏色名稱，如紅色、黃色等），然後畫出每一個色相值出現次數的柱狀圖如右下。

我們可以看到，」黃色「這一種顏色可以有很多個色相值（換成RGB模型也是一樣的），這就是我們所說的雜訊或者不確定性。從分布上來看，」黃色「的色相值以大約53為中心然後向兩邊擴散。為了記錄「黃色」，一種方法是記錄每一個「黃色」像素的色相值，然而這將十分消耗存儲空間，確切來說，有多少像素就消耗多大的空間。

另一種簡潔有效的方法是使用高斯分布。

讓我們先來溫故一下高斯分布的數學表達形式再返回來解釋這個例子吧。

高斯分布的數學形式如下：

$p(x) = frac{1}{sqrt{2pi} sigma}e^{-frac{(x-mu)^2}{2sigma^2}}$

其中， $x$ 是變數， $p(x)$ 表示 $x$ 的出現的概率， $mu$ 是均值， $sigma$ 是標準差（ $sigma^2$ 是方差）。 $mu$ 和 $sigma$ 是之前我們提到的表示一個高斯分布的兩個參數。

在我們上面的栗子中，假設我們已經完成了對黃色的高斯分布建模， $x$ 是一個採樣像素的色相值，那麼 $p(x)$ 就表示該像素屬於「黃色」的概率。

我們繼續介紹高斯分布，為了更加直觀地了解高斯分布長什麼樣子，我們把高斯分布畫出來。下圖是 $mu=0,sigma=1$ 的情況下的高斯分布圖（標準高斯分布）。

當均值變化而方差不變的時候，可以看作將圖像進行了平移。

當均值不變，而標準差變化的時候，可以看作將圖像進行了「擠」或「壓」（方差越大越扁）。

總結來說， $mu$ 決定了高斯分布的對稱軸， $sigma$ 決定了高斯分布的擴散程度。

讓我們返回到小球檢測的栗子。我們只使用一個高斯分布（即兩個值）就可以近似地存儲所有的「黃色」像素，如下圖。

那麼如何通過數據來求 $mu$ 和 $sigma$ 這兩個參數呢？我們將在下一節中介紹最大似然估計方法求解一元高斯分布。

2. 求解一元高斯分布：最大似然估計

在這一節，我們將介紹使用最大似然估計（Maximum Likelihood Estimation，MLE）的方法來求解觀測數據的一元高斯分布。

首先，我們引入似然函數（或者說似然，Likelihood）的概念。

似然指的是在給定模型參數的情況下觀測值出現的概率。它被定義成了一個條件概率 $p({x_i}|mu,sigma)$ 。比如，在上一節的小球例子中， ${x_i}$ 表示所有「黃色」像素點的色相值，如果給定了高斯分布的 $mu$ 和 $sigma$ ，我們就能算出 ${x_i}$ 出現的概率。

那麼我們的問題就可以表述為，給定觀測值 ${x_i}$ ，求 $mu$ 和 $sigma$ ，使得似然函數最大： $hat{mu}, hat{sigma}= ext{arg} underset{mu,sigma}{max} {p({x_i}|mu,sigma)}$

其中， $hat{mu}$ 和 $hat{sigma}$ 表示估計值。

這裡，觀測值 ${x_i}$ 是所有的觀測值的集合，我們假設這些觀測值兩兩相互獨立（在我們的小球檢測栗子中也確實是這樣的）。那麼：

$p({x_i}|mu,sigma) = prod_{i=1}^{N}p(x_i|mu,sigma)$

這裡， $N$ 表示觀測值集合的大小為 $N$ 。那麼我們的目標就變成了：

$hat{mu}, hat{sigma}= ext{arg} underset{mu,sigma}{max} {prod_{i=1}^Np(x_i|mu,sigma)}$

根據對數函數的單調性( $x_1<x_2equiv ln x_1<ln x_2$ )，取對數再求最大值等價於直接求最大值，即：

$ext{arg} underset{mu,sigma}{max} {prod_{i=1}^Np(x_i|mu,sigma)} = ext{arg} underset{mu,sigma}{max}ln{ {prod_{i=1}^Np(x_i|mu,sigma)} }$

那麼我們的目標就變成了：

$egin{align}hat{mu},hat{sigma} &= ext{arg} underset{mu,sigma}{max}ln{ {prod_{i=1}^Np(x_i|mu,sigma)} }\ &= ext{arg} underset{mu,sigma}{max}sum_{i=1}^{N}ln p(x_i|mu,sigma)end{align}$

而將 $ln p(x_i | mu,sigma)$ 帶入高斯分布我們有：

$egin{align}ln p(x_i|mu,sigma) &= ln ( frac{1}{sqrt{2pi}sigma}e^{-frac{(x_i-mu)^2}{2sigma^2}} ) \&= -frac{(x_i-mu)^2}{2sigma^2}-lnsigma-lnsqrt{2pi}end{align}$

那麼我們的目標就變成了：

$egin{align}hat{mu},hat{sigma} &= ext{arg} underset{mu,sigma}{max}sum_{i=1}^{N} (-frac{(x_i-mu)^2}{2sigma^2}-lnsigma-lnsqrt{2pi}) \&= ext{arg} underset{mu,sigma}{max}sum_{i=1}^{N} (-frac{(x_i-mu)^2}{2sigma^2}-lnsigma) \&= ext{arg} underset{mu,sigma}{min}sum_{i=1}^{N} (frac{(x_i-mu)^2}{2sigma^2}+lnsigma)end{align}$

讓我們給後面的目標函數做個記號，令

$J(mu,sigma) = sum_{i=1}^{N} (frac{(x_i-mu)^2}{2sigma^2}+lnsigma)$

我們的目標就是使得 $J(mu,sigma)$ 最小，而函數 $J(mu,sigma)$ 取最小值時其關於 $mu$ 的偏導數等於0：

$egin{align}frac{partial J}{partial mu}&= sum_{i=1}^Nfrac{-2(x_i-mu)}{2sigma^2}\&=frac{1}{sigma^2}(Nmu-sum_{i=1}^Nx_i)\&=0end{align}$

求得：

$hatmu = frac{1}{N}sum_{i=1}^{N}x_i$

將 $mu=hatmu$ 帶入函數 $J(mu,sigma)$ 得到函數 $J$ ，其取最小值時關於 $sigma$ 的導數等於0：

$egin{align}frac{d J$

求得：

$hatsigma^2 = frac{1}{N}sum_{i=1}^N (x_i-hatmu)^2$

總結來說，最後的答案為：

$hatmu = frac{1}{N}sum_{i=1}^{N}x_i$ $hatsigma^2 = frac{1}{N}sum_{i=1}^N (x_i-hatmu)^2$

回到我們小球檢測的栗子，我們可以求出的高斯分布如下圖：

在這一節中，我們介紹了怎樣使用MLE來求解一元高斯分布。然而實際情況中，可能需要提取的特徵值不止色相值這一個（比如RGB三個值）。這樣一元高斯分布不再適用。為了解決這個問題，我們將在下一節中介紹多元高斯分布。

3. 多元高斯分布

讓我們回到小球檢測的栗子，在一元高斯分布下，我們只使用了色相值這一個性質。然而，顏色其實是用多個維度來定義的。比如，在HSV模型下，除了色相值還有飽和度（Saturation）和亮度（Value）。而我們通常使用的三原色光模式（RGB模型）將顏色表示成紅色（R）、綠色（G）和藍色（B）的疊加。如果我們用RGB值來表示一個顏色，怎樣表示我們栗子中的小球呢？我們將圖片中所有像素點的RGB值用散點圖的形式畫出來可以得到下面的圖：

那我們怎樣對這種圖形進行建模呢？如這一節的題目所說，我們將一元高斯分布擴展到多元高斯分布並對RGB值進行建模。

讓我們首先來介紹多元高斯分布的數學形式吧：

$p(oldsymbol{x}) = frac{1}{(2pi)^{D/2}|Sigma|^{1/2}}exp{-frac{1}{2}(oldsymbol{x}-oldsymbolmu)^TSigma^{-1}(oldsymbol{x}-oldsymbolmu)}$

多元高斯分布和一元高斯分布是十分相似的，我們用加粗的 $oldsymbol{x}$ 來表示變數（一個向量）， $D$ 表示維度（元的數目），加粗的 $oldsymbolmu$ 表示平均向量，大寫的 $Sigma$ 表示協方差矩陣（Covariance Matrix，是一個方陣）， $|Sigma|$ 表示 $Sigma$ 的行列式值， $(oldsymbol{x}-oldsymbolmu)^T$ 表示矩陣 $(oldsymbol{x}-oldsymbolmu)$ 的轉置。

值得一提的是協方差矩陣，它由兩部分組成，方差（Variance）和相關性（Correlation），對角線上的值表示方差，非對角線上的值表示維度之間的相關性。拿一個二維協方差矩陣作栗子：

其中，對角線上的 $sigma^2_{x_1}$ 和 $sigma^2_{x_2}$ 分別表示變數 $x_1$ 和 $x_2$ 的獨立方差，非對角線上的 $sigma_{x_1}sigma_{x_2}$ 表示兩個變數之間的相關性（注意 $sigma_{x_1}sigma_{x_2}$ 和 $sigma_{x_2}sigma_{x_1}$ 是相等的）。

回到小球檢測的栗子，我們考慮用RGB來對「紅色」小球進行多元高斯分布的建模，那麼各個參數就如下圖所示了：

我們來看一下標準二元高斯分布圖：

各個參數的變化對多元高斯分布的影響請參考課程內容本身。值得一提的是，協方差矩陣中的相關性參數是一元高斯分布中沒有的。那麼怎樣求解多元高斯分布呢？確切地說，怎樣求解多元高斯分布中的平均向量 $oldsymbolmu$ 和協方差矩陣 $Sigma$ 呢？

4. 求解多元高斯分布：最大似然估計

和求解一元高斯分布類似，我們將問題描述為：給定觀測值 ${oldsymbol{x_i}}$ ，求 $oldsymbolmu$ 和 $Sigma$ ，使得似然函數最大：

$hat{oldsymbolmu}, hat{Sigma}= ext{arg} underset{oldsymbolmu,Sigma}{max} {p({oldsymbol{x_i}}|oldsymbolmu,Sigma)}$

同樣，假設觀測值兩兩相互獨立，根據獨立概率公式，我們有：

$hat{oldsymbolmu}, hat{Sigma}= ext{arg} underset{oldsymbolmu,Sigma}{max} {prod_{i=1}^N p(oldsymbol{x_i}|oldsymbolmu,Sigma)}$

同樣（1）取對數，（2）將多元高斯分布的形式帶入，我們有：

$egin{align}hat{oldsymbolmu}, hat{Sigma}&= ext{arg} underset{oldsymbolmu,Sigma}{max} ln prod_{i=1}^N p(oldsymbol{x_i}|oldsymbolmu,Sigma)\&= ext{arg} underset{oldsymbolmu,Sigma}{max} sum_{i=1}^N ln p(oldsymbol{x_i} | oldsymbolmu, Sigma) \&= ext{arg} underset{oldsymbolmu,Sigma}{max} sum_{i=1}^N ( -frac{1}{2}(oldsymbol{x_i}-oldsymbolmu)^TSigma^{-1}(oldsymbol{x_i}-oldsymbolmu) - frac{1}{2}ln |Sigma| - frac{D}{2}ln (2pi) ) \&= ext{arg} underset{oldsymbolmu,Sigma}{min} sum_{i=1}^N ( frac{1}{2}(oldsymbol{x_i}-oldsymbolmu)^TSigma^{-1}(oldsymbol{x_i}-oldsymbolmu) + frac{1}{2}ln |Sigma| )end{align}$

我們給目標函數做個記號，令

$J(oldsymbolmu, Sigma) = sum_{i=1}^N ( frac{1}{2}(oldsymbol{x_i}-oldsymbolmu)^TSigma^{-1}(oldsymbol{x_i}-oldsymbolmu) + frac{1}{2}ln |Sigma| )$

我們仍然分別對 $oldsymbolmu$ 和 $Sigma$ 求偏導來計算 $hat{oldsymbolmu}$ 和 $hat{Sigma}$ 。（這裡需要矩陣求導的知識，可以參考Matrix Calculus Manual）

$egin{align}frac{partial J}{partial oldsymbolmu}&= frac{partial}{partialoldsymbolmu}sum_{i=1}^N(frac{1}{2}(oldsymbol{x_i}-oldsymbolmu)^TSigma^{-1}(oldsymbol{x_i}-oldsymbolmu)+frac{1}{2}|Sigma|)\&=frac{partial}{partialoldsymbolmu}frac{1}{2}sum_{i=1}^N(oldsymbol{x_i}^TSigma^{-1}oldsymbol{x_i}-oldsymbol{x_i}^TSigma^{-1}oldsymbolmu-oldsymbolmu^TSigma^{-1}oldsymbol{x_i}+oldsymbolmu^TSigma^{-1}oldsymbolmu)\&=frac{partial}{partialmu}sum_{i=1}^N(frac{1}{2}oldsymbolmu^TSigma^{-1}oldsymbolmu-oldsymbolmuSigma^{-1}oldsymbol{x_i})\&=Sigma^{-1}sum_{i=1}^N(oldsymbolmu-oldsymbol{x_i})\&=oldsymbol{0}end{align}$ $egin{align}frac{partial J}{partialSigma}&=frac{partial}{partialSigma}sum_{i=1}^N(frac{1}{2}(oldsymbol{x_i}-hat{oldsymbolmu})^TSigma^{-1}(oldsymbol{x_i}-hat{oldsymbolmu})+frac{1}{2}ln|Sigma|)\&=frac{1}{2}sum_{i=1}^N(-Sigma^{-1}(oldsymbol{x_i}-hat{oldsymbolmu})(oldsymbol{x_i}-hat{oldsymbolmu})^TSigma^{-1}+Sigma^{-1})\&=frac{1}{2}Sigma^{-1}(-(sum_{i=1}^N(oldsymbol{x_i}-hat{oldsymbolmu})(oldsymbol{x_i}-hat{oldsymbolmu})^T)Sigma^{-1}+Noldsymbol{I})\&=oldsymbol{0}end{align}$

求得，

$hat{oldsymbolmu}=frac{1}{N}sum_{i=1}^Noldsymbol{x_i}$

$hatSigma=frac{1}{N}sum_{i=1}^N(oldsymbol{x_i}-hat{oldsymbolmu})(oldsymbol{x_i}-hat{oldsymbolmu})^T$

5. 高斯混合模型

在第一節和第四節，我們介紹了一元高斯分布和多元高斯分布，他們都屬於單高斯模型（Single Gaussian Model）。使用單高斯模型來建模有一些限制，例如，它一定只有一個眾數，它一定對稱的。舉個例子，如果我們對下面的分布建立單高斯模型，會得到顯然相差很多的模型：

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

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

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

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

$p(oldsymbol{x})=sum_{k=1}^Kw_kg_k(oldsymbol{x}|oldsymbolmu_k,Sigma_k)$

其中， $g_k$ 是均值為 $oldsymbolmu_k$ ，協方差矩陣為 $Sigma_k$ 的單高斯模型， $w_k$ 是 $g_k$ 的權重係數（ $w_k>0, sum_{k=1}^Kw_k=1$ ）， $K$ 是單高斯模型的個數。

前面我們提到了GMM的優點（能夠表示任何分布），當然GMM也有缺點。其一，參數太多了，每一個單高斯模型都有均值、協方差矩陣和權重係數，另外還有單高斯模型的個數 $K$ 也是其中一個參數，這使得求解GMM十分複雜。其二，有時候所建立的高斯混合模型對觀測值的建模效果很好，但是其他值可能效果不好，我們稱這個缺點為過度擬合（Overfitting）。

這一節我們介紹了高斯混合模型，與前面類似，我們將在下節介紹求解混合高斯模型的方法。

6. 求解高斯混合模型：EM演算法

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

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

$p(oldsymbol{x})=sum_{k=1}^Kw_kg_k(oldsymbol{x}|oldsymbolmu_k,Sigma_k)$

這其中需要求解的變數很多： $K$ 、 $w_k$ 、 $oldsymbolmu_k$ 、 $Sigma_k$ 。為了簡化問題，我們假定 $K$ 是預先設定好的，並且每個單高斯分布的權重相等，即假定：

$w_k = frac{1}{K},k=1,2,cdots,K$

這樣，我們需要確定的就只剩下每個單高斯分布的參數（ $oldsymbolmu_k$ 、 $Sigma_k$ ）了。

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

在牛頓法中，我們需要預先猜測一個初始解，同樣在EM演算法中，我們也需要預先猜測一個初始解（ $oldsymbolmu_1,Sigma_1,oldsymbolmu_2,Sigma_2,cdots,oldsymbolmu_K,Sigma_K$ ）。

在EM演算法中，我們引入一個中間變數 $z_k^i = frac{g_k(oldsymbol{x_i}|oldsymbolmu_k,Sigma_k)}{sum_{i=1}^Kg_k(oldsymbol{x_i}|oldsymbolmu_k,Sigma_k)}(k=1,2,cdots,K; i=1,2,cdots,N)$ ，其中 $K$ 是單高斯分布的個數， $N$ 是觀測值的數目。

直觀地可以這樣理解：在確定了所有單高斯分布之後，可以計算觀測值 $oldsymbol{x_i}$ 發生的概率（分母部分）和第 $k$ 個單高斯分布 $g_k(oldsymbolmu_k,Sigma_k)$ 下觀測值 $x_i$ 發生的概率（分子部分），這樣 $z^i_k$ 就可以理解為第 $k$ 個高斯分布對觀測值 $x_i$ 發生的概率的貢獻了。如下表所示：

$egin{tabular}{|c|c|c|c|c|c|c|c|}hline& $frac{1}{K}$ & $frac{1}{K}$ & $cdots$ & $frac{1}{K}$ & $cdots$ & $frac{1}{K}$ \& $oldsymbolmu_1$ & $oldsymbolmu_2$ & $cdots$ & $oldsymbolmu_k$ & $cdots$ & $oldsymbolmu_K$ \& $Sigma_1$ & $Sigma_2$ & $cdots$ & $Sigma_k$ & $cdots$ & $Sigma_K$ \hline$x_1$ & $z^1_1$ & $z^1_2$ & $cdots$ & $z^1_k$ & $cdots$ & $z^1_K$ \hline$x_2$ & $z^2_1$ & $z^2_2$ & $cdots$ & $z^2_k$ & $cdots$ & $z^2_K$ \hline$vdots$ & $vdots$ & $vdots$ & & $vdots$ & & $vdots$ \hline$x_i$ & $z^i_1$ & $z^i_2$ & $cdots$ & $z^i_k$ & $cdots$ & $z^i_K$ \hline$vdots$ & $vdots$ & $vdots$ & & $vdots$ & & $vdots$ \hline$x_N$ & $z^N_1$ & $z^N_2$ & $cdots$ & $z^N_k$ & $cdots$ & $z^N_K$ \hline& $z_1$ & $z_2$ & $cdots$ & $z_k$ & $cdots$ & $z_K$ \hlineend{tabular}$

表中最後一行 $z_k=sum_{i=1}^Nz^i_k$ 可以理解為第 $k$ 個單高斯分布對整個GMM的貢獻。

例如，當GMM中 $K=2$ 時（即由兩個單高斯分布組成）：

在這個圖中，對於觀測值 $oldsymbol{x_i}$ ， $p_1/2$ 是第一個單高斯分布 $g_1$ 下 $oldsymbol{x_i}$ 發生的概率， $p_2/2$ 是第二個單高斯分布 $g_2$ 下 $oldsymbol{x_i}$ 發生的概率， $(p_1+p_2)/2$ 是在整個GMM下 $oldsymbol{x_i}$ 發生的概率； $z_1^i=frac{p_1}{p_1+p_2}$ 表示 $g_1$ 對 $oldsymbol{x_i}$ 發生概率的貢獻， $z_2^i=frac{p_2}{p_1+p_2}$ 表示 $g_2$ 對 $oldsymbol{x_i}$ 發生概率的貢獻。

當我們算出了所有的 $z_k^i$ 之後，我們再反過來用它更新每一個單高斯分布。更新的公式為： $oldsymbolmu_k = frac{1}{z_k}sum_{i=1}^Nz_k^ioldsymbol{x_i}$

$Sigma_k = frac{1}{z_k}sum_{i=1}^Nz_k^i(oldsymbol{x_i}-oldsymbolmu_k)(oldsymbol{x_i}-oldsymbolmu_k)^T$

大家可以想一想，對於第 $k$ 個高斯分布 $g_k$ ， $z_k^i/z_k$ 可以理解為第 $i$ 個觀測值 $oldsymbol{x_i}$ 的貢獻，而 $oldsymbolmu_k$ 表示 $g_k$ 的平均值，用 $frac{1}{z_k}sum_{i=1}^Nz_k^ioldsymbol{x_i}$ 來更新 $oldsymbolmu_k$ 就很好理解了，進而再更新協方差矩陣 $Sigma_k$ 。

這樣，兩組值互相更新，當相鄰兩個步驟值的差小於事先設定一個閾值（threshold）時（收斂），我們停止循環，輸出 $oldsymbolmu_k$ 和 $Sigma_k$ 為所求。

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

7. 實戰：小球檢測

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

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

（註：我們使用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)") pauseendsave("Samples.mat", "Samples"); % save the samples to filefigure, scatter3(Samples(:,1),Samples(:,2),Samples(:,3),".");title("Pixel Color Distribubtion");xlabel("Red");ylabel("Green");zlabel("Blue");

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

我們選用多元高斯分布作為我們的模型，應用第四節中的結論，編寫代碼coefficient.m求出平均值 $oldsymbolmu$ 和協方差矩陣 $Sigma$ 並存儲到本地以供後續使用。

mu = mean(Samples); % musig=zeros(3,3); % sigmafor i=1:N data=double(Samples(i,:)); sig=sig+(data-mu)"*(data-mu)/N;end% save the coefficients to filessave("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 - mufor 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 multiplyId = 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 possibilityId=reshape(Id,row,col);% set thresholdthr=8e-06;% binary image about if each pixel "is ball"Id=Id>thr;% find the biggest ball areasegI = 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)") pauseend

效果如下：

成功！當然，我們也可以用其他種類的高斯模型對黃色小球進行建模，這裡就不做示範啦！

感謝閱讀！