Stanford機器學習課程筆記4-Kmeans與高斯混合模型
這一部分屬於無監督學習的內容，無監督學習內容主要包括：Kmeans聚類演算法、高斯混合模型及EM演算法、Factor Analysis、PCA、ICA等。本文是Kmeans聚類演算法、高斯混合模型的筆記，EM演算法是適用於存在latent/hidden變數的通用演算法，高斯混合模型僅僅是EM演算法的一種特殊情況，關於EM演算法的推到參見Andrew Ng講義。由於公式太多，最近時間又忙實習的事就簡單寫一些，回頭看時還得參考Ng的筆記和自己的列印Notes上的筆記，這裡的程式對理解可能能提供另外的一些幫助。
Kmeans聚類
其實聚類演算法除了Kmeans，還有其它的聚類方式，記得曾經接觸到的就有層次聚類、基於模糊等價關係的
- 針對浮點型資料的運算，為提高運算效率，將浮點型Round成整形（我操作的話一般先乘一個係數，以使資料分佈在整個整數空間，降低Round的資料損失）
- 層次化的Kmeans聚類(HIKM)
這兩種改進都能在VLFeat裡見到C程式碼實現。作為完善基礎知識，我還是願意先來複習下最基本的Kmeans聚類，後面再根據自己的經驗補充點HIKM的東西。這是Stanford機器學習課程中第一個非監督學習演算法，為嚴謹一些，先給出一般Kmeans聚類問題的描述：
- 給定資料集:
- 將資料group成K個clusters
請注意：相對於監督學習演算法，Kmeans雖然沒有label資訊，但聚類數K卻是已知的，也就是說：使用Kmeans我們得計劃好要聚成多少個cluster。當然也有一些Kmeans基礎上的改進方法可以通過閾值直接判斷聚類數目，這裡不討論。Kmeans演算法的流程是：
先隨機初始化K個聚類中心，計算樣本資料到聚類中心的距離，用聚類結果計算新的中心，如此迭代。提出兩個問題：
- 隨機初始化聚類中心怎麼個隨機法：這K個點是資料點還是非資料點。其實除了一般的隨機，還有一種叫K-means++(參考2)的初始化方式
凡是迭代演算法都有一個問題值得討論：收斂準則。白話描述Kmeans收斂準則就是：聚類中心u穩定，不再發生太大的變化。嚴謹一些，定義Distortion Function：
J表示每個訓練樣本到被分配到的聚類中心的距離之和。當J最小時則Kmeans演算法收斂。然而，由於J是非凸函式，利用座標下降演算法無法保證能收斂到全域性最優解(座標下降演算法在Andrew Ng在SVM講義部分中有討論，是梯度下降演算法的多維情況：先固定A變數，B變數用梯度下降迭代；再固定B變數，A變數用梯度下降迭代)。然而，在實際應用中，對譬如（1）存在多個區域性最有之間波動的情況很少見（2）即使收斂到區域性最有Kmeans演算法，對於實際問題這個區域性最優解也一般能夠接受。所以Kmeans依然是使用最廣泛的無監督聚類演算法。
Matlab的Statistics and Machine Learning Toolbox自帶kmeans演算法，更多內容可以參考一下Matlab中的doc kmeans或者http://cn.mathworks.com/help/stats/kmeans.html 。
混合高斯模型（Mixture of Gaussian）
尖子班和普通班的成績混在一起了，要把這兩個班的成績分開。只知道2個班的成績都滿足高斯分佈，但不知道每個班的高斯分佈的平均分和方差波動。如何通過無監督聚類的方法將這兩個班的資料分開？上面的例子就是高斯混合模型的一個例子。
回想GDA（高斯判別分析）的生成模型，當時建立的模型是：p(y)服從伯努利實驗，p(x|y)服從高斯分佈，即隱藏在資料背後的高斯模型差生了訓練資料，通過極大似然法使聯合概率p(x,y)最大，求解得到伯努利引數phi，高斯分佈引數u,sigma。混合高斯模型由於是無監督的資料，沒有y，所以建模時假定一個latent/hidden(隱藏)變數z（等價GDA中的y），利用GDA的完全一樣的極大似然法求解phi(z)、u(z)、sigma(z)，但此時的pht、u、sigma都是和latent變數z有關的不確定的數。這就需要一種迭代演算法：先估計z，更確切的說，這裡是要來估計z的概率；再用z的估計值計算phi、u和sigma；更新z的估計，再迭代。。。。
高斯混合模型是EM演算法的一個特例，迭代演算法分為E-Step（估計z的概率）和M-Step（似然函式最大化），
在E-Step中，估計的是z的後驗概率，可以先通過初始化的phi、u、sigma計算似然概率和先驗概率，再用Bayes Rule得到z的後驗估計。EM演算法與Kmeans演算法一樣可能收斂到區域性最優，有點不同的是EM演算法的聚類中心數是可以自動決定的而Kmeans是預先給定的。下面是從 http://www.mathworks.com/matlabcentral/fileexchange/26184-em-algorithm-for-gaussian-mixture-model 找到的一份高斯混合模型的EM程式碼，也可以下載完整的EM Example在Matlab上執行
function [label, model, llh] = emgm(X, init)
% Perform EM algorithm for fitting the Gaussian mixture model.
% X: d x n data matrix
% init: k (1 x 1) or label (1 x n, 1<=label(i)<=k) or center (d x k)
% Written by Michael Chen ([email protected]).
%% initialization
fprintf('EM for Gaussian mixture: running ... \n');
R = initialization(X,init);
[~,label(1,:)] = max(R,[],2);
R = R(:,unique(label));
tol = 1e-10;
maxiter = 500;
llh = -inf(1,maxiter);
converged = false;
t = 1;
while ~converged && t < maxiter
t = t+1;
model = maximization(X,R);
[R, llh(t)] = expectation(X,model);
[~,label(:)] = max(R,[],2);
u = unique(label); % non-empty components
if size(R,2) ~= size(u,2)
R = R(:,u); % remove empty components
else
converged = llh(t)-llh(t-1) < tol*abs(llh(t));
end
end
llh = llh(2:t);
if converged
fprintf('Converged in %d steps.\n',t-1);
else
fprintf('Not converged in %d steps.\n',maxiter);
end
function R = initialization(X, init)
[d,n] = size(X);
if isstruct(init) % initialize with a model
R = expectation(X,init);
elseif length(init) == 1 % random initialization
k = init;
idx = randsample(n,k);
m = X(:,idx);
[~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
[u,~,label] = unique(label);
while k ~= length(u)
idx = randsample(n,k);
m = X(:,idx);
[~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
[u,~,label] = unique(label);
end
R = full(sparse(1:n,label,1,n,k,n));
elseif size(init,1) == 1 && size(init,2) == n % initialize with labels
label = init;
k = max(label);
R = full(sparse(1:n,label,1,n,k,n));
elseif size(init,1) == d %initialize with only centers
k = size(init,2);
m = init;
[~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
R = full(sparse(1:n,label,1,n,k,n));
else
error('ERROR: init is not valid.');
end
function [R, llh] = expectation(X, model)
mu = model.mu;
Sigma = model.Sigma;
w = model.weight;
n = size(X,2);
k = size(mu,2);
logRho = zeros(n,k);
for i = 1:k
logRho(:,i) = loggausspdf(X,mu(:,i),Sigma(:,:,i));
end
logRho = bsxfun(@plus,logRho,log(w));
T = logsumexp(logRho,2);
llh = sum(T)/n; % loglikelihood
logR = bsxfun(@minus,logRho,T);
R = exp(logR);
function model = maximization(X, R)
[d,n] = size(X);
k = size(R,2);
nk = sum(R,1);
w = nk/n;
mu = bsxfun(@times, X*R, 1./nk);
Sigma = zeros(d,d,k);
sqrtR = sqrt(R);
for i = 1:k
Xo = bsxfun(@minus,X,mu(:,i));
Xo = bsxfun(@times,Xo,sqrtR(:,i)');
Sigma(:,:,i) = Xo*Xo'/nk(i);
Sigma(:,:,i) = Sigma(:,:,i)+eye(d)*(1e-6); % add a prior for numerical stability
end
model.mu = mu;
model.Sigma = Sigma;
model.weight = w;
function y = loggausspdf(X, mu, Sigma)
d = size(X,1);
X = bsxfun(@minus,X,mu);
[U,p]= chol(Sigma);
if p ~= 0
error('ERROR: Sigma is not PD.');
end
Q = U'\X;
q = dot(Q,Q,1); % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U))); % normalization constant
y = -(c+q)/2;參考
- Andrew Ng Lecture Notes.
- D. Arthur and S. Vassilvitskii. k-means++: The advantages of careful seeding. In Proc. ACM-SIAM Symp. on Discrete Algorithms, 2007.