使用训练的高斯混合模型标记新数据

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）。