高斯混合模型（理論+opencv實現）

時間 2019-11-09

標籤高斯混合模型理論 opencv 實現简体版

原文原文鏈接

查資料的時候看了一個不文明的事情，轉載別人的東西而不標註出處，結果原創無人知曉，轉載不少人評論~~標註了轉載而不說出處這樣的人有點可恥！ html

寫在前面：ios

　　　　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.協方差ui

　　注：內容是http://pinkyjie.com/2010/08/31/covariance/（講的太好了，一會兒就理解了，對下面的高斯混合模型沒啥幫助，我覺得會用到的，哈哈）url

　　統計學的基本概念spa

　　　　學過幾率統計的孩子都知道，統計裏最基本的概念就是樣本的均值，方差，或者再加個標準差。首先咱們給你一個含有n個樣本的集合 $X = {X_{1}, \dots, X_{n}}$ code

$X = {X_{1}, \dots, X_{n}}$ orm

　　很顯然，均值描述的是樣本集合的中間點，它告訴咱們的信息是頗有限的，而標準差給咱們描述的則是樣本集合的各個樣本點到均值的距離之平均。以這兩個集合爲例，[0，8，12，20]和[8，9，11，12]，兩個集合的均值都是10，但顯然兩個集合差異是很大的，計算二者的標準差，前者是8.3，後者是1.8，顯而後者較爲集中，故其標準差小一些，標準差描述的就是這種「散佈度」。之因此除以n-1而不是除以n，是由於這樣能使咱們以較小的樣本集更好的逼近整體的標準差，即統計上所謂的「無偏估計」。而方差則僅僅是標準差的平方。

爲何須要協方差？

　　上面幾個統計量看似已經描述的差很少了，但咱們應該注意到，標準差和方差通常是用來描述一維數據的，但現實生活咱們經常遇到含有多維數據的數據集，最簡單的你們上學時免不了要統計多個學科的考試成績。面對這樣的數據集，咱們固然能夠按照每一維獨立的計算其方差，可是一般咱們還想了解更多，好比，一個男孩子的猥瑣程度跟他受女孩子歡迎程度是否存在一些聯繫啊，嘿嘿~協方差就是這樣一種用來度量兩個隨機變量關係的統計量，咱們能夠仿照方差的定義：

來度量各個維度偏離其均值的程度，標準差能夠這麼來定義：

　　協方差的結果有什麼意義呢？若是結果爲正值，則說明二者是正相關的(從協方差能夠引出「相關係數」的定義)，也就是說一我的越猥瑣就越受女孩子歡迎，嘿嘿，那必須的~結果爲負值就說明負相關的，越猥瑣女孩子越討厭，可能嗎？若是爲0，也是就是統計上說的「相互獨立」。

從協方差的定義上咱們也能夠看出一些顯而易見的性質，如：

1. c o v (X, X) = v a r (X)

2. c o v (X, Y) = c o v (Y, X)

協方差多了就是協方差矩陣

　　上一節提到的猥瑣和受歡迎的問題是典型二維問題，而協方差也只能處理二維問題，那維數多了天然就須要計算多個協方差，好比n維的數據集就須要計算 $\frac{n!}{(n - 2)! * 2}$

這個定義仍是很容易理解的，咱們能夠舉一個簡單的三維的例子，假設數據集有 ${x, y, z}$

可見，協方差矩陣是一個對稱的矩陣，並且對角線是各個維度上的方差。

Matlab協方差實戰

　　上面涉及的內容都比較容易，協方差矩陣彷佛也很簡單，但實戰起來就很容易讓人迷茫了。必需要明確一點，### 協方差矩陣計算的是不一樣維度之間的協方差，而不是不一樣樣本之間的。這個我將結合下面的例子說明，如下的演示將使用Matlab，爲了說明計算原理，不直接調用Matlab的cov函數。

首先，隨機產生一個10*3維的整數矩陣做爲樣本集，10爲樣本的個數，3爲樣本的維數。

1	MySample = fix(rand(10,3)*50)

　　根據公式，計算協方差須要計算均值，那是按行計算均值仍是按列呢，我一開始就總是困擾這個問題。前面咱們也特別強調了，協方差矩陣是計算不一樣維度間的協方差，要時刻牢記這一點。樣本矩陣的每行是一個樣本，每列爲一個維度，因此咱們要### 按列計算均值。爲了描述方便，咱們先將三個維度的數據分別賦值：

 
     dim1 = MySample(:,1); 
    
     dim2 = MySample(:,2); 
    
     dim3 = MySample(:,3);

計算dim1與dim2，dim1與dim3，dim2與dim3的協方差：

 
     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

搞清楚了這個後面就容易多了，協方差矩陣的對角線就是各個維度上的方差，下面咱們依次計算：

 
     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 }