高斯混合聚类(GMM)及代码实现
高斯混合聚类(GMM)及代码实现
by 2017-03-20 郭昱良 机器学习算法与Python学习
通过学习概率密度函数的Gaussian Mixture Model (GMM) 与 k-means 类似，不过 GMM 除了用在 clustering 上之外，还经常被用于 density estimation。对于二者的区别而言简单地说，k-means 的结果是每个数据点被 assign 到其中某一个 cluster ，而 GMM 则给出这些数据点被 assign 到每个 cluster 的概率。
作为一个流行的算法，GMM 肯定有它自己的一个相当体面的归纳偏执了。其实它的假设非常简单，顾名思义，Gaussian Mixture Model ，就是假设数据服从 Mixture Gaussian Distribution ，换句话说，数据可以看作是从数个 Gaussian Distribution 中生成出来的。实际上，一般在K-means 和 K-medoids 中所用的例子就是由多个Gaussian 分布随机生成的。实际上，从中心极限定理可以看出，Gaussian 分布这个假设其实是比较合理的，除此之外，Gaussian 分布在计算上也有一些很好的性质，所以，虽然我们可以用不同的分布来随意地构造 XX Mixture Model ，但是还是 GMM 最为流行。另外，Mixture Model 本身其实也是可以变得任意复杂的，通过增加 Model 的个数，我们可以任意地逼近任何连续的概率密分布。
每个 GMM 由 K 个 Gaussian 分布组成，每个 Gaussian 称为一个“Component”，这些 Component 线性加成在一起就组成了 GMM 的概率密度函数：
那么如何用 GMM 来做 clustering 呢？其实很简单，假定它们是由 GMM 生成出来的，那么我们只要根据数据推出 GMM 的概率分布来就可以了，然后 GMM 的 K 个 Component 实际上就对应了 K 个 cluster 。根据数据来推算概率密度通常被称作 density estimation 。
现在，假设我们有 N 个数据点，并假设它们服从某个分布（记作 p(x)），现在要确定里面的一些参数的值。 我们的想法是，找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大，而这个概率实际上就等于p(x)之和，我们把这个乘积称作似然函数 (Likelihood Function)。通常单个点的概率都很小，许多很小的数字相乘起来在计算机里很容易造成浮点数下溢，因此我们通常会对其取对数，把乘积变成加和 ，得到 log-likelihood function 。接下来我们只要将这个函数最大化（通常的做法是求导并令导数等于零，然后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，我们就认为这是最合适的参数，这样就完成了参数估计的过程。
下面让我们来看一看 GMM 的 log-likelihood function ：
求解最大似然估计过程如下：
- 估计数据由每个 Component 生成的概率（并不是每个 Component 被选中的概率）：对于每个数据 xi 来说，它由第 k 个 Component 生成的概率为
由于式子里的 方差 和 均值 也是需要我们估计的值，我们采用迭代法，在计算时将取上一次迭代所得的值（或者初始值）。
-
估计每个 Component 的参数：现在我们假设上一步中得到的 r(i,k) 就是正确的“数据 xi 由 Component k 生成的概率”。集中考虑所有的数据点，现在实际上可以看作 Component 生成了所有的 这些点。由于每个 Component 都是一个标准的 Gaussian 分布，可以很容易分布求出最大似然所对应的参数值：
-
重复迭代前面两步，直到似然函数的值收敛为止。
采用MATLAB自带的kmeansdata数据集进行验证仿真，具体代码如下所示。
'
%% 导入数据
load('kmeansdata')
%% 初始化混合模型参数
K = 3;
% 随机初始化均值和协方差
means = randn(K,2);
for k = 1:K
covs(:,:,k) = randeye(2);
end
priors = repmat(1/K,1,K); % 初始化，假设隐含变量服从先验均匀分布
%% 主算法
MaxIts = 100; % 最大迭代次数
N = size(X,1); % 样本数
q = zeros(N,K); % 后验概率
D = size(X,2); % 维数
cols = {'r','g','b'};
plotpoints = [1:1:10,12:2:30 40 50];
B(1) = -inf;
converged = 0;
it = 0;
tol = 1e-2;
while 1
it = it + 1;
% 把乘除化为对数加减运算，防止乘积结果过于接近于0
for k = 1:K
const = -(D/2)log(2pi) - 0.5log(det(covs(:,:,k)));
Xm = X - repmat(means(k,:),N,1);
temp(:,k) = const - 0.5 * diag(Xminv(covs(:,:,k))Xm'); end
% 计算似然下界
if it > 1
B(it) = sum(sum(q.log(repmat(priors,N,1))))+ sum(sum(q.temp)) - sum(sum(q.*log(q)));
if abs(B(it)-B(it-1))<tol
converged = 1;
end
end
if converged == 1 || it > MaxIts
break
end
% 计算每个样本属于第k类的后验概率
temp = temp + repmat(priors,N,1);
q = exp(temp - repmat(max(temp,[],2),1,K));
q(q < 1e-60) = 1e-60;
q(q > (1-(1e-60))) = 1-(1e-60);
q = q./repmat(sum(q,2),1,K); % 更新先验分布
priors = mean(q,1); % 更新均值
for k = 1:K
means(k,:) = sum(X.repmat(q(:,k),1,D),1)./sum(q(:,k)); end
% 更新方差
for k = 1:K
Xm = X - repmat(means(k,:),N,1);
covs(:,:,k) = (Xm.repmat(q(:,k),1,D))'*Xm;
covs(:,:,k) = covs(:,:,k)./sum(q(:,k));
end
end
%% plot the data
figure(1);hold on;
plot(X(:,1),X(:,2),'ko')
;for k = 1:K
plot_2D_gauss(means(k,:),covs(:,:,k), -2:0.1:5,-6:0.1:6);
end
ti = sprintf('After %g iterations',it);
title(ti)
%% 绘制似然下界迭代过程图
figure(2);
hold off
plot(2:length(B),B(2:end),'k');
xlabel('Iterations');
ylabel('Bound');
'
