查資料的時候看了一個不文明的事情，轉載別人的東西而不標注出處，結果原創無人知曉，轉載很多人評論~~標注了轉載而不說出處這樣的人有點可恥！
寫在前面：
Gaussian Mixture Model (GMM)。事實上，GMM 和 k-means 很像，不過 GMM 是學習出一些概率密度函數來（所以 GMM 除了用在 clustering 上之外，還經常被用於 density estimation ），簡單地說，k-means 的結果是每個數據點被 assign 到其中某一個 cluster 了，而 GMM 則給出這些數據點被 assign 到每個 cluster 的概率，又稱作 soft assignment 。
得出一個概率有很多好處，因為它的信息量比簡單的一個結果要多，比如，我可以把這個概率轉換為一個 score ，表示算法對自己得出的這個結果的把握。也許我可以對同一個任務，用多個方法得到結果，最后選取“把握”最大的那個結果；另一個很常見的方法是在諸如疾病診斷之類的場所，機器對於那些很容易分辨的情況（患病或者不患病的概率很高）可以自動區分，而對於那種很難分辨的情況，比如，49% 的概率患病，51% 的概率正常，如果僅僅簡單地使用 50% 的閾值將患者診斷為“正常”的話，風險是非常大的，因此，在機器對自己的結果把握很小的情況下，會“拒絕發表評論”，而把這個任務留給有經驗的醫生去解決。
准備階段
1.協方差
注：內容是http://pinkyjie.com/2010/08/31/covariance/（講的太好了，一下子就理解了，對下面的高斯混合模型沒啥幫助，我以為會用到的，哈哈）
統計學的基本概念
學過概率統計的孩子都知道，統計里最基本的概念就是樣本的均值，方差，或者再加個標准差。首先我們給你一個含有n個樣本的集合X={X1,…,Xn}X={X1,…,Xn}，依次給出這些概念的公式描述，
這些高中學過數學的孩子都應該知道吧，一帶而過。
很顯然，均值描述的是樣本集合的中間點，它告訴我們的信息是很有限的，而標准差給我們描述的則是樣本集合的各個樣本點到均值的距離之平均。以這兩個集合為例，[0，8，12，20]和[8，9，11，12]，兩個集合的均值都是10，但顯然兩個集合差別是很大的，計算兩者的標准差，前者是8.3，后者是1.8，顯然后者較為集中，故其標准差小一些，標准差描述的就是這種“散布度”。之所以除以n-1而不是除以n，是因為這樣能使我們以較小的樣本集更好的逼近總體的標准差，即統計上所謂的“無偏估計”。而方差則僅僅是標准差的平方。
為什么需要協方差？
上面幾個統計量看似已經描述的差不多了，但我們應該注意到，標准差和方差一般是用來描述一維數據的，但現實生活我們常常遇到含有多維數據的數據集，最簡單的大家上學時免不了要統計多個學科的考試成績。面對這樣的數據集，我們當然可以按照每一維獨立的計算其方差，但是通常我們還想了解更多，比如，一個男孩子的猥瑣程度跟他受女孩子歡迎程度是否存在一些聯系啊，嘿嘿~協方差就是這樣一種用來度量兩個隨機變量關系的統計量，我們可以仿照方差的定義：
來度量各個維度偏離其均值的程度，標准差可以這么來定義：
協方差的結果有什么意義呢？如果結果為正值，則說明兩者是正相關的(從協方差可以引出“相關系數”的定義)，也就是說一個人越猥瑣就越受女孩子歡迎，嘿嘿，那必須的~結果為負值就說明負相關的，越猥瑣女孩子越討厭，可能嗎？如果為0，也是就是統計上說的“相互獨立”。
從協方差的定義上我們也可以看出一些顯而易見的性質，如：
協方差多了就是協方差矩陣
上一節提到的猥瑣和受歡迎的問題是典型二維問題，而協方差也只能處理二維問題，那維數多了自然就需要計算多個協方差，比如n維的數據集就需要計算n!(n−2)!∗2n!(n−2)!∗2個協方差，那自然而然的我們會想到使用矩陣來組織這些數據。給出協方差矩陣的定義：
這個定義還是很容易理解的，我們可以舉一個簡單的三維的例子，假設數據集有{x,y,z}{x,y,z}三個維度，則協方差矩陣為
可見，協方差矩陣是一個對稱的矩陣，而且對角線是各個維度上的方差。
Matlab協方差實戰
上面涉及的內容都比較容易，協方差矩陣似乎也很簡單，但實戰起來就很容易讓人迷茫了。必須要明確一點，### 協方差矩陣計算的是不同維度之間的協方差，而不是不同樣本之間的。這個我將結合下面的例子說明，以下的演示將使用Matlab，為了說明計算原理，不直接調用Matlab的cov函數。
首先，隨機產生一個10*3維的整數矩陣作為樣本集，10為樣本的個數，3為樣本的維數。
|
1
|
MySample = fix(rand(10,3)*50)
|
根據公式，計算協方差需要計算均值，那是按行計算均值還是按列呢，我一開始就老是困擾這個問題。前面我們也特別強調了，協方差矩陣是計算不同維度間的協方差，要時刻牢記這一點。樣本矩陣的每行是一個樣本，每列為一個維度，所以我們要### 按列計算均值。為了描述方便，我們先將三個維度的數據分別賦值：
|
1
2
3
|
dim1 = MySample(:,1);
dim2 = MySample(:,2);
dim3 = MySample(:,3);
|
計算dim1與dim2，dim1與dim3，dim2與dim3的協方差：
|
1
2
3
|
sum( (dim1-mean(dim1)) .* (dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) % 得到 74.5333
sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10.0889
sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -106.4000
|
搞清楚了這個后面就容易多了，協方差矩陣的對角線就是各個維度上的方差，下面我們依次計算：
|
1
2
3
|
std(dim1)^2 % 得到 108.3222
std(dim2)^2 % 得到 260.6222
std(dim3)^2 % 得到 94.1778
|
這樣，我們就得到了計算協方差矩陣所需要的所有數據，調用Matlab自帶的cov函數進行驗證：
|
1
|
cov(MySample)
|
把我們計算的數據對號入座，是不是一摸一樣？
Update：
今天突然發現，原來協方差矩陣還可以這樣計算，先讓樣本矩陣中心化，即每一維度減去該維度的均值，使每一維度上的均值為0，然后直接用新的到的樣本矩陣乘上它的轉置，然后除以(N-1)即可。其實這種方法也是由前面的公式通道而來，只不過理解起來不是很直觀，但在抽象的公式推導時還是很常用的！同樣給出Matlab代碼實現：
|
1
2
|
X = MySample - repmat(mean(MySample),10,1); % 中心化樣本矩陣，使各維度均值為0
C = (X'*X)./(size(X,1)-1);
|
總結
理解協方差矩陣的關鍵就在於牢記它計算的是不同維度之間的協方差，而不是不同樣本之間，拿到一個樣本矩陣，我們最先要明確的就是一行是一個樣本還是一個維度，心中明確這個整個計算過程就會順流而下，這么一來就不會迷茫了~
2.高斯函數(模型)
記得上次看高斯函數還是在高斯濾波和雙邊濾波的時候，現在再來提到都已經忘得差不多了。。。。
知道下面幾個公式就可以了：
高斯一維函數：
高斯概率分布函數：
3.EM算法
本來打算在這里敘述的，感覺內容有點多，而且以后還會用到這個算法，還是另外開一篇博客敘述：http://www.cnblogs.com/wjy-lulu/p/7010258.html
4.混合高斯模型
注釋：這里是http://blog.sina.com.cn/s/blog_54d460e40101ec00.html博客的內容，看了很多篇文章，感覺還是這篇講的好，不是很深入但是把原理講的很清楚。
如果理解了我另一篇博客的EM算法，看這個混合高斯模型很簡單了。。。。
高斯混合模型--GMM（Gaussian Mixture Model）
統計學習的模型有兩種，一種是概率模型，一種是非概率模型。
所謂概率模型，是指訓練模型的形式是P(Y|X)。輸入是X，輸出是Y，訓練后模型得到的輸出不是一個具體的值，而是一系列的概率值（對應於分類問題來說，就是輸入X對應於各個不同Y（類）的概率），然后我們選取概率最大的那個類作為判決對象（軟分類--soft assignment）。所謂非概率模型，是指訓練模型是一個決策函數Y=f(X)，輸入數據X是多少就可以投影得到唯一的Y，即判決結果（硬分類--hard assignment）。
所謂混合高斯模型（GMM）就是指對樣本的概率密度分布進行估計，而估計采用的模型（訓練模型）是幾個高斯模型的加權和（具體是幾個要在模型訓練前建立好）。每個高斯模型就代表了一個類（一個Cluster）。對樣本中的數據分別在幾個高斯模型上投影，就會分別得到在各個類上的概率。然后我們可以選取概率最大的類所為判決結果。
從中心極限定理的角度上看，把混合模型假設為高斯的是比較合理的，當然，也可以根據實際數據定義成任何分布的Mixture Model,不過定義為高斯的在計算上有一些方便之處，另外，理論上可以通過增加Model的個數，用GMM近似任何概率分布。
混合高斯模型的定義為：
其中K 為模型的個數；πk為第k個高斯的權重；p（x / k） 則為第k個高斯概率密度，其均值為μk，方差為σk。對此概率密度的估計就是要求出πk、μk 和σk 各個變量。當求出p（x ）的表達式后，求和式的各項的結果就分別代表樣本x 屬於各個類的概率。
在做參數估計的時候，常采用的是最大似然方法。最大似然法就是使樣本點在估計的概率密度函數上的概率值最大。由於概率值一般都很小，N 很大的時候, 連乘的結果非常小，容易造成浮點數下溢。所以我們通常取log，將目標改寫成：
也就是最大化對數似然函數，完整形式為：
一般用來做參數估計的時候，我們都是通過對待求變量進行求導來求極值，在上式中，log函數中又有求和，你想用求導的方法算的話方程組將會非常復雜，沒有閉合解。可以采用的求解方法是EM算法——將求解分為兩步：第一步,假設知道各個高斯模型的參數（可以初始化一個，或者基於上一步迭代結果），去估計每個高斯模型的權值；第二步,基於估計的權值，回過頭再去確定高斯模型的參數。重復這兩個步驟，直到波動很小，近似達到極值（注意這里是極值不是最值，EM算法會陷入局部最優）。具體表達如下：
注釋：這是需要求得參數有兩類：（每個模型權重W）（每個模型的參數），這是按照EM思想，首先假定（參數）把每個樣本點分類計算（權重），然后通過分類的樣本點計算新的（參數），接着對比假設參數和計算的參數是否符合精度，符合結束，不符合就繼續上面的操作去迭代。
A、（E step）
對於第i個樣本xi 來說，它由第k 個model 生成的概率為：
在這一步，假設高斯模型的參數和是已知的（由上一步迭代而來或由初始值決定）。
B、（M step）
C、重復上述兩步驟直到算法收斂。
聚類實例：
1 #include <opencv2/opencv.hpp> 2 #include <iostream> 3 4 using namespace cv; 5 using namespace std; 6 7 int main(int argc, char** argv) { 8 Mat img(500, 500, CV_8UC3); 9 RNG rng(12345); 10 const int Max_nCluster = 5; 11 Scalar colorTab[] = { 12 Scalar(0, 0, 255), 13 Scalar(0, 255, 0), 14 Scalar(255, 0, 0), 15 Scalar(0, 255, 255), 16 Scalar(255, 0, 255) 17 }; 18 //int numCluster = rng.uniform(2, Max_nCluster + 1);//隨機類數 19 int numCluster = 3; 20 int sampleCount = rng.uniform(5, 1000);//樣本點數量 21 Mat matPoints(sampleCount, 2, CV_64FC1);//樣本點矩陣：sampleCount X 2 22 Mat labels; 23 Mat centers; 24 // 生成隨機數 25 for (int k = 0; k < numCluster; k++) { 26 Point center;//隨機產生中心點 27 center.x = rng.uniform(0, img.cols); 28 center.y = rng.uniform(0, img.rows); 29 Mat pointChunk = matPoints.rowRange(k*sampleCount / numCluster, 30 (k + 1)*sampleCount / numCluster); 31 //-----符合高斯分布的隨機高斯 32 rng.fill(pointChunk, RNG::NORMAL, Scalar(center.x, center.y, 0, 0), Scalar(img.cols*0.05, img.rows*0.05, 0, 0)); 33 } 34 randShuffle(matPoints, 1, &rng);//打亂高斯生成的數據點順序 35 // EM Cluster Train 36 Ptr<ml::EM> em_model = ml::EM::create(); 37 em_model->setClustersNumber(numCluster); 38 em_model->setCovarianceMatrixType(ml::EM::COV_MAT_SPHERICAL); 39 em_model->setTermCriteria(TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 100, 0.1)); 40 em_model->trainEM(matPoints, noArray(), labels, noArray()); 41 42 // 用不同顏色顯示分類 43 img = Scalar::all(255); 44 for (int i = 0; i < sampleCount; i++) { 45 int index = labels.at<int>(i); 46 Point p = Point(matPoints.at<double>(i,0), matPoints.at<double>(i,1)); 47 circle(img, p, 2, colorTab[index], -1, 8); 48 } 49 // classify every image pixels 50 Mat sample(1, 2, CV_32FC1); 51 for (int row = 0; row < img.rows; row++) { 52 for (int col = 0; col < img.cols; col++) { 53 sample.at<float>(0) = (float)col; 54 sample.at<float>(1) = (float)row; 55 int response = cvRound(em_model->predict2(sample, noArray())[1]); 56 Scalar c = colorTab[response]; 57 circle(img, Point(col, row), 1, c*0.75, -1); 58 } 59 } 60 imshow("Demo", img); 61 waitKey(0); 62 return 0; 63 }
圖像分割：
1 #include <opencv2/opencv.hpp> 2 #include <iostream> 3 4 using namespace cv; 5 using namespace cv::ml; 6 using namespace std; 7 8 int main(int argc, char** argv) { 9 Mat src = imread("toux.jpg"); 10 if (src.empty()) { 11 printf("could not load iamge...\n"); 12 return -1; 13 } 14 char* inputWinTitle = "input image"; 15 namedWindow(inputWinTitle, CV_WINDOW_AUTOSIZE); 16 imshow(inputWinTitle, src); 17 18 // 初始化 19 int numCluster = 3; 20 const Scalar colors[] = { 21 Scalar(255, 0, 0), 22 Scalar(0, 255, 0), 23 Scalar(0, 0, 255), 24 Scalar(255, 255, 0) 25 }; 26 27 int width = src.cols; 28 int height = src.rows; 29 int dims = src.channels(); 30 int nsamples = width*height; 31 Mat points(nsamples, dims, CV_64FC1); 32 Mat labels; 33 Mat result = Mat::zeros(src.size(), CV_8UC3); 34 35 // 圖像RGB像素數據轉換為樣本數據 36 int index = 0; 37 for (int row = 0; row < height; row++) { 38 for (int col = 0; col < width; col++) { 39 index = row*width + col; 40 Vec3b rgb = src.at<Vec3b>(row, col); 41 points.at<double>(index, 0) = static_cast<int>(rgb[0]); 42 points.at<double>(index, 1) = static_cast<int>(rgb[1]); 43 points.at<double>(index, 2) = static_cast<int>(rgb[2]); 44 } 45 } 46 47 // EM Cluster Train 48 Ptr<EM> em_model = EM::create(); 49 em_model->setClustersNumber(numCluster); 50 em_model->setCovarianceMatrixType(EM::COV_MAT_SPHERICAL); 51 em_model->setTermCriteria(TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 100, 0.1)); 52 em_model->trainEM(points, noArray(), labels, noArray()); 53 54 // 對每個像素標記顏色與顯示 55 Mat sample(1, dims, CV_64FC1);// 56 double time = getTickCount(); 57 int r = 0, g = 0, b = 0; 58 for (int row = 0; row < height; row++) { 59 for (int col = 0; col < width; col++) { 60 index = row*width + col; 61 62 b = src.at<Vec3b>(row, col)[0]; 63 g = src.at<Vec3b>(row, col)[1]; 64 r = src.at<Vec3b>(row, col)[2]; 65 sample.at<double>(0, 0) = static_cast<double>(b); 66 sample.at<double>(0, 1) = static_cast<double>(g); 67 sample.at<double>(0, 2) = static_cast<double>(r); 68 int response = cvRound(em_model->predict2(sample, noArray())[1]); 69 Scalar c = colors[response]; 70 result.at<Vec3b>(row, col)[0] = c[0]; 71 result.at<Vec3b>(row, col)[1] = c[1]; 72 result.at<Vec3b>(row, col)[2] = c[2]; 73 } 74 } 75 printf("execution time(ms) : %.2f\n", (getTickCount() - time) / getTickFrequency() * 1000); 76 imshow("EM-Segmentation", result); 77 78 waitKey(0); 79 return 0; 80 }
參考：
http://www.mamicode.com/info-detail-466443.html（用大白文敘述GMM，娛樂的看看還好）
http://blog.pluskid.org/?p=39（講的很好）
http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006924.html（高斯混合模型的EM算法推導公式，感興趣可以看看，反正我沒看）
代碼參考賈志剛老師的資料，他的代碼有的不行，自己改了一下就可以了。