高斯混合聚類(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');
'