使用訓練的高斯混合模型標記新數據

Question

我不確定如何使用訓練的高斯混合模型（GMM）對一些新數據進行預測。例如，我已經從3個不同的類（集羣）中獲得了一些標記的數據。對於每一類數據點，我適合一個GMM（gm1，gm2和gm3）。假設我們知道每個類的高斯混合數（例如k1 = 2，k2 = 1和k3 = 3），或者可以使用Akaike信息準則（AIC）估計（優化）它的數量。然後，當我得到一些新的數據集時，我怎麼能知道它是否更可能屬於第1,2或3類？

一些MATLAB腳本顯示了我的意思是：

clc; clf; clear all; close all; 

%% Create some artificial training data 

% 1. Cluster 1 with two mixture of Gaussian (k1 = 2) 
rng default; % For reproducibility 
mu1     = [1 2]; 
sigma1    = [3 .2; .2 2]; 
mu2     = [-1 -2]; 
sigma2    = [2 0; 0 1]; 
X1     = [mvnrnd(mu1,sigma1,200); mvnrnd(mu2,sigma2,100)]; 

options1   = statset('Display', 'final'); 
k1     = 2; 
gm1     = fitgmdist(X1, k1, 'Options', options1); 


% 2. Cluster 2 with one mixture of Gaussian (k2 = 1) 
mu3     = [6 4]; 
sigma3    = [3 .1; .1 4]; 
X2     = mvnrnd(mu3,sigma3,300); 

options2   = statset('Display', 'final'); 
k2     = 1; 
gm2     = fitgmdist(X2, k2, 'Options', options2); 

% 3. Cluster 3 with three mixture of Gaussian (k3 = 3) 
mu4     = [-5 -6]; 
sigma4    = [1 .1; .1 1]; 
mu5     = [-5 -10]; 
sigma5    = [6 .1; .1 1]; 
mu6     = [-2 -15]; 
sigma6    = [8 .1; .1 4]; 
X3     = [mvnrnd(mu4,sigma4,200); mvnrnd(mu5,sigma5,300); mvnrnd(mu6,sigma6,100)]; 

options3   = statset('Display', 'final'); 
k3     = 3; 
gm3     = fitgmdist(X3, k3, 'Options', options3); 

% Display 
figure, 
scatter(X1(:,1),X1(:,2),10,'ko'); hold on; 
ezcontour(@(x,y)pdf(gm1, [x y]), [-12 12], [-12 12]); 
scatter(X2(:,1),X2(:,2),10,'ko'); 
ezcontour(@(x,y)pdf(gm2, [x y]), [-12 12], [-12 12]); 
scatter(X3(:,1),X3(:,2),10,'ko'); 
ezcontour(@(x,y)pdf(gm3, [x y]), [-12 12], [-12 12]); hold off; 

我們可以得到圖中：

然後我們有例如一些新的測試數據：

%% Create some artificial testing data 
mut1    = [6.1 3.8]; 
sigmat1    = [3.1 .1; .1 4.2]; 
mut2    = [5.8 4.5]; 
sigmat2    = [2.8 .1; .1 3.8]; 
Xt1     = [mvnrnd(mut1,sigmat1,500); mvnrnd(mut2,sigmat2,100)]; 

figure, 
scatter(Xt1(:,1),Xt1(:,2),10,'ko'); 
xlim([-12 12]); ylim([-12 12]); 

我故意將測試數據類似於Cluster 2數據。在我們使用GMM進行訓練之後，我們能否以某種方式預測新測試數據的標籤？對於每個類的預測，是否有可能得到一些概率（p1 = 18％，p2 = 80％和p3 = 2％）。因爲我們已經得到了p2 = 80％，所以我們可以有一個硬分類，即新的測試數據被標記爲簇2.

p.s .:我找到了這篇文章，但對我來說似乎理論上(A similar post)。如果你可以請在你的回覆中加入一些簡單的Matlab腳本。

非常感謝。 A.

---

編輯：

至於荷銀回答該問題的解決方案，我有更多的問題。

1. 荷銀創造了一個新的GMM使用整個數據集的一些初始化：

   % initial parameters of the new GMM (combination of the previous three) 
   % (note PComponents is normalized according to proportion of data in each subset) 
   S = struct('mu',[gm1.mu; gm2.mu; gm3.mu], ... 
       'Sigma',cat(3, gm1.Sigma, gm2.Sigma, gm3.Sigma), ... 
       'PComponents',[gm1.PComponents*n1, gm2.PComponents*n2, gm3.PComponents*n3]./n); 
   
   % train the final model over all instances 
   opts = statset('MaxIter',1000, 'Display','final'); 
   gmm = fitgmdist(X, k, 'Options',opts, 'Start',S); 
   

   什麼荷銀得到的是類似下面

   

   這可能不適合我因爲它將我標記的cluster1和cluster2與cluster1的一部分混合在一起。這是我想要避免的。

   這裏我介紹的是一個人工數值例子;然而，在我的實際應用中，它處理圖像分割問題（​​例如，cluster1是我的背景圖像，cluster2是我想分離的對象）。然後，我嘗試以某種方式「強制」單獨的GMM以適應不同的課程。如果兩個簇很遠（例如本例中的簇1和簇3），那麼使用Amro的方法來組合所有數據然後進行GMM擬合就沒有問題。但是，當我們對圖像數據進行訓練時，由於分辨率的限制（由於部分音量效應），將背景與對象分離絕不是完美的;因此，如圖所示，很可能我們有cluster1與cluster2重疊的情況。我想也許混合所有的數據，然後做擬合會導致一些問題，以進一步預測新的數據，我說得對嗎？

   然而，思維一點點之後，我想做現在要做的是：

   % Combine the mixture of Gaussian and form a new gmdistribution 
   muAll    = [gm1.mu; gm2.mu; gm3.mu]; 
   sigmaAll   = cat(3, gm1.Sigma, gm2.Sigma, gm3.Sigma); 
   
   gmAll    = gmdistribution(muAll, sigmaAll); 
   
   pt1     = posterior(gmAll, Xt1); 
   

   你們有什麼覺得？或者它相當於Amro的方法？如果是這樣，有沒有辦法強制我的訓練GMM分開？

2. 此外，我有關於使用posterior函數的基本原理的問題。基本上，我想估計給定GMM擬合的測試數據的可能性。那麼爲什麼我們現在計算後驗概率呢？或者這只是一個命名問題（換句話說，'後驗概率'='可能性'）？

3. 據我所知，GMM一直被用作無監督方法。有人甚至向我提到GMM是k均值聚類的概率版本。它有資格以這種'監督'的風格使用它嗎？任何推薦的論文或參考？

非常感謝您的回覆！ 答：

Answer 1

實際上，您已經培訓了三個GMM模型，而不是一個，每個模型都是混合物本身。通常，您將創建一個GMM具有多個組件，每個組件代表一個集羣...

所以我會做什麼，你的情況是創建培訓了整個數據集新的GMM模型（X1，X2和X3）元件數量等於三個GMM中所有元件的總和（即2+1+3 = 6高斯混合）。該模型將使用單獨訓練的參數進行初始化。

這裏的代碼來說明（我使用的是你在你的例子中創建相同的變量）：

% number of instances in each data subset 
n1 = size(X1,1); 
n2 = size(X2,1); 
n3 = size(X3,1); 

% the entire dataset 
X = [X1; X2; X3]; 
n = n1 + n2 + n3; 
k = k1 + k2 + k3; 

% initial parameters of the new GMM (combination of the previous three) 
% (note PComponents is normalized according to proportion of data in each subset) 
S = struct('mu',[gm1.mu; gm2.mu; gm3.mu], ... 
    'Sigma',cat(3, gm1.Sigma, gm2.Sigma, gm3.Sigma), ... 
    'PComponents',[gm1.PComponents*n1, gm2.PComponents*n2, gm3.PComponents*n3]./n); 

% train the final model over all instances 
opts = statset('MaxIter',1000, 'Display','final'); 
gmm = fitgmdist(X, k, 'Options',opts, 'Start',S); 

% display GMM density function over training data 
line(X(:,1), X(:,2), 'LineStyle','none', ... 
    'Marker','o', 'MarkerSize',1, 'Color','k') 
hold on 
ezcontour(@(x,y) pdf(gmm,[x y]), xlim(), ylim()) 
hold off 
title(sprintf('GMM over %d training instances',n)) 

---

現在，我們已經培養了GMM模型在整個訓練數據集（與k=6混合），我們可以用它來聚類新的數據實例：

cIdx = cluster(gmm, Xt1); 

這是作爲手動計算組件的後驗概率，並考慮該組件具有最大概率作爲簇指數相同：

pr = posterior(gmm, Xt1); 
[~,cIdx] = max(pr,[],2); 

作爲測試數據的預期幾乎95％的集羣作爲屬於同一成分：

>> tabulate(cIdx) 
    Value Count Percent 
     1  27  4.50% 
     2  0  0.00% 
     3  573  95.50% 

這裏是匹配高斯參數：

>> idx = 3; 
>> gmm.mu(idx,:) 
ans = 
    5.7779 4.1731 
>> gmm.Sigma(:,:,idx) 
ans = 
    2.9504 0.0801 
    0.0801 4.0907 

這確實對應於組件位於上圖右上角。

同樣，如果您檢查其他組件idx=1，它將是前一個左邊的那個，它解釋了600個測試實例中的27個是如何「錯誤分類」的......下面是GMM如何自信是在那些實例：

>> pr(cIdx==1,:) 
ans = 
    0.9813 0.0001 0.0186 0.0000 0.0000 0.0000 
    0.6926 0.0000 0.3074 0.0000 0.0000 0.0000 
    0.5069 0.0000 0.4931 0.0000 0.0000 0.0000 
    0.6904 0.0018 0.3078 0.0000 0.0000 0.0000 
    0.6954 0.0000 0.3046 0.0000 0.0000 0.0000 
    <... output truncated ...> 
    0.5077 0.0000 0.4923 0.0000 0.0000 0.0000 
    0.6859 0.0001 0.3141 0.0000 0.0000 0.0000 
    0.8481 0.0000 0.1519 0.0000 0.0000 0.0000 

這裏是覆蓋在上圖的頂部的測試實例：

hold on 
gscatter(Xt1(:,1), Xt1(:,2), cIdx) 
hold off 
title('clustered test instances') 

---

編輯：

我上面的例子是爲了展示如何使用GMM來對數據進行聚類（無監督學習）。根據我現在所理解的，你想要的是用現有的訓練模型對數據進行分類（監督學習）。我想我很困惑你使用的「集羣」一詞:)

無論如何，它應該是現在輕鬆;只需使用每個模型計算測試數據的類別 - 條件概率密度函數，並選擇具有最高可能性的模型作爲類別標籤（不需要將模型合併爲一個）。

所以繼續您的初始代碼，這將簡單地是：

p = [pdf(gm1,Xt), pdf(gm2,Xt), pdf(gm3,Xt)]; % P(x|model_i) 
[,cIdx] = max(p,[],2);       % argmax_i P(x|model_i) 

cIdx是在測試數據的每個實例的類的預測（1，2或3）。