前言：這兩天着手做畢設，在今年的研究生數學建模的基礎上研究“大數據下多流形聚類分析”問題，導師要求我這周把每一個算法的實現對比一下效果展示給他看，表示今天google的搜索結果中沒有找到諸如SSC的函數教程，又養成了不copy代碼的習慣，那就自己從頭開始學唄，剛好mathworks上面提供一篇詳細的聚類分析的教程，特此翻譯一下，希望自己和讀者都能更好的咬文嚼字，以作為未來幾天高維度數據matlab聚類實現的熱身運動。
下面的例子將演示如何使用Matlab內的統計和機器學習工具箱中的聚類分析算法檢驗觀測值或對象之間的異同。俗話說的好：“物以類聚，人以群分”。不管是人還是物總是呈現出抱團或者集群的狀態；其中，相同的種類集群對象的特征是相似的，不同簇對象的特性不一致。
本文內容如下：
- K均值和層次聚類
- 費雪鳶尾花數據
- 對費雪鳶尾花數據進行K-均值聚類
K均值和層次聚類
統計和機器學習工具箱中有兩個聚類算法的函數實現，分別是：K均值和層次聚類算法。
K-均值聚類算法是會把數據中的觀測值作為擁有位置和觀測值之間的距離的對象來划分，也就是說該算法會根據對象的位置和對象間的距離來決定聚類結果。
層次聚類算法是一種對每一個觀察值進行分組觀察的算法，該算法會在各種距離的基礎上創建一個集群樹。這棵樹不同於單一類集合的K-均值聚類，它神奇的地方在於其“多層次結構”，也就是說某一個級別的不同類可能會在下一個更高級別中被聚為一類。意味着你可以根據你的應用選擇最適合的聚類規模——“聚類之樹”的大小和高度。
在本次示例使用的函數中會調用Matlab內置函數rng來生成隨機數。因此為了讓本次示例中的結果一致，你應該通過執行以下命令來設置隨機種子控制隨機數的產生。如果你不設置成相同的隨機狀態將會導致一些不必要的麻煩，例如，你可能看發現集群的編號順序不同；還有就是可能會導致一個非最優的聚類結果（本次示例包含非最優聚類結果探討以及避免該結果的方法）。
rng(14,'twister');
譯者注：由於我的Matlab屬於早期版本R2010a，因此沒有內置隨機數控制器rng函數；故各位使用早期版本的同志可以通過如下方式設置隨機種子，詳情見官網解釋：
rand('twister',14);
費雪鳶尾花數據
在20世紀20年代，植物學家埃德加·安德森集從鳶尾花樣本中測量了三個不同品種的鳶尾花萼片長寬和花瓣長寬的數據，其中每個種類分別有50個小樣本。該測量數據被稱為費雪鳶尾花數據或安德森鳶尾花數據。
在這些數據中的各個觀測值都來自於已知物種，因此能夠很清晰的將這些數據進行分組。現在，我們將忽略數據中的品種信息，只使用原始的測量數據來對觀測值進行聚類。當我們完成聚類的時候，就可以將所得的聚類結果與實際物種分類結果進行對比，以查看鳶尾花的三個種類是否表現出明顯的個性特征。
對費雪鳶尾花數據進行K-均值聚類
函數kmeans實現K-均值聚類，通過使用迭代算法讓所有觀測值分配到類中，使得每一個對象到對應類的質心的距離之和最小。該方法應用於本次的費雪鳶尾花數據時，將能夠很自然的根據數據中的萼片和花瓣的測量數據對鳶尾花樣本進行種類分組。使用K-均值聚類的時候，你必須要指定聚類數。
首先，加載數據然后調用函數kmeans，該函數指定了希望聚類的數量為2，並使用（平方）歐式距離計算點到質心的距離。你可以繪制輪廓圖來更好的知道聚類分離的效果，輪廓圖顯示的是觀察值到同一類內的點之間的距離與到其他類內的點之間的距離的對比。
load fisheriris % 加載數據到矩陣變量meas到向量species中 [cidx2,cmeans2] = kmeans(meas,2,'dist','sqeuclidean'); % 參數‘dist’的值為‘sqeuclidean’ [silh2,h] = silhouette(meas,cidx2,'sqeuclidean'); % 繪制輪廓圖並返回輪廓值到si1h2
從上面所示的輪廓圖可以你可以看到，兩類中大部分的點都有很高的輪廓值（大於0.8），說明這些點能夠從臨近的群中很好的分離出來，但是每一個類中同樣有少部分點的輪廓值是比較低的，表明他們很接近其他類的點，很容易與其他類的點混淆，其特征不夠明顯。
事實表明，數據中的第四個測量變量——花瓣寬度與第三個測量變量花瓣長度之間的關系相當緊密。因此，我們可以在不考慮第四個測量變量的情況下根據前三個變量數據繪制出一個展示效果良好的三維圖。如果根據k-均值聚類的結果給每一個類的數據用不同的符號表示繪制散點圖，你將會很方便的識別出那些輪廓值低的點，因為這些點與其他類的點靠在一起。
ptsymb = {'bs','r^','md','go','c+'};
for i = 1:2
clust = find(cidx2==i);
plot3(meas(clust,1),meas(clust,2),meas(clust,3),ptsymb{i}); % 使用藍色方塊和紅色三角形分別標記2類的觀察值
hold on
end
plot3(cmeans2(:,1),cmeans2(:,2),cmeans2(:,3),'ko');
plot3(cmeans2(:,1),cmeans2(:,2),cmeans2(:,3),'kx'); % 使用符號來標記兩個類的質心
hold off
xlabel('Sepal Length'); % 標記x軸
ylabel('Sepal Width');
zlabel('Petal Length');
view(-137,10); % 旋轉
grid on % 添加網格
上圖中，很明顯的有上下兩個類；且使用圓叉符號標記每一個類的質心。其中屬於下層的有3個使用三角形標記的觀察值非常靠近上層中使用正方形標記的群集。但是，實際上，由於上部分的類展開的比較分散，因此那3個點相對於上面的類而言跟接近下面的類的質心，即使他們與被划分的那一個類別之間被一條溝分離。由於k-均值聚類只考慮距離而不是密度，也就使得這樣的聚類結果產生。
你可以通過增加簇的數目來查看kmeans函數是否可以從數據中找到進一步的分組結構。這一次，我們使用可選參數“display”打印出聚類算法在每一次迭代過程中的詳細信息。
[cidx3,cmeans3] = kmeans(meas,3,'display','iter');
在每一次迭代中，k均值算法都會重新類內的點以降低點到質心之間的距離之和，然后為新的類的划分重新計算聚類中心。注意在每一次迭代后重新分配的觀察值的數目都會減少直到距離之和達到最小值的時候質心穩定，不再需要重新聚類。在本次示例中k均值算法由兩個階段組成，第二階段不再進行任何重新分配，表明第一階段經過3此迭代后達到最小值78.8514。
默認情況下，k均值算法會使用隨機選擇的初始質心位置處開始聚類過程。正如其他類型的數值最小化算法一樣，k均值算法要想達到距離最小值有時候取決於初始點的選定；並且有可能因此只會收斂到局部最小值，這就會造成分配點到某一個類的時候計算的該點到質心的聚類之和相對來說不是全局最小值，從而因此錯過更好的聚類結果。好消息是你可以指定可選參數“replicates”來解決這一問題。當你指定該參數為1以上的時候，k均值算法將會為每一次重復聚類選定不同的隨機產生的質心開始聚類過程。
[cidx3,cmeans3,sumd3] = kmeans(meas,3,'replicates',5,'display','final'); % 指定重復5次選取初始質心
這里告訴我們，隨機種子的設定非常重要，一旦隨機子不同，選取的隨機數就不同，造成不同的初始點選定。其中，在選定隨機種子為13的時候，在重復第5次聚類過程經過7次迭代，出現了局部最小值為142.754。因此，對於這個簡單的示例，非全局最小值仍然存在。其中，第三個輸出參數“sumd3”包含了最佳聚類方案的每一個類內距離之和數據。
sum(sumd3) % 對類內距離之和求和
為這個聚為3類的解決方案繪制輪廓圖，該圖表明有一個類能夠很好的從其他類中分離，但是另外兩個類區分的並不是很明顯。
[silh3,h] = silhouette(meas,cidx3,'sqeuclidean');
你還可以通過再次繪制原始數據的散點圖觀察點的聚類情況。
for i = 1:3
clust = find(cidx3==i);
plot3(meas(clust,1),meas(clust,2),meas(clust,3),ptsymb{i});
hold on
end
plot3(cmeans3(:,1),cmeans3(:,2),cmeans3(:,3),'ko');
plot3(cmeans3(:,1),cmeans3(:,2),cmeans3(:,3),'kx');
hold off
xlabel('Sepal Length');
ylabel('Sepal Width');
zlabel('Petal Length');
view(-137,10);
grid on
由上圖可以發現原來屬於上面的那個類被聚為兩個類別，並且這兩個類里面的觀察值彼此之間挨得很緊密。當然，這三個聚類結果相對於前面的兩個聚類結果而言是否對你有用取決於你的打算。函數“silhouette”的第一個輸出參數包含每一個點的輪廓值，根據輪廓值的均值結果對比可知輪廓均值較大的那個解決方案即聚為兩類相對來說更合理。
[mean(silh2) mean(silh3)]
你也可以使用不同的距離指標來對這些數據重新聚類。余弦距離來評估這些數據的相似性使用可能更有意義，因為它會只考慮數據的相對大小而忽略測量數據的絕對大小。因此，兩個花朵的大小不一，但擁有相同形狀的花瓣和萼片，使用歐式距離就無法決定它們相似但是使用余弦距離卻可以下出這兩朵花相似的結論。
[cidxCos,cmeansCos] = kmeans(meas,3,'dist','cos');
從剪影圖的現實結果來看，余弦距離相對歐式距離的聚類結果似乎只顯示出輕微的較好的分離效果。
[silhCos,h] = silhouette(meas,cidxCos,'cos'); [mean(silh2) mean(silh3) mean(silhCos)]
注意，本次聚類的類別編號和上一次聚類的編號不同，這是因為k均值算法隨機的選擇初始值。通過繪制原始數據散點圖，可以看出使用不同距離創建的群集形狀的差異。這兩種聚類的結果相似，但是使用余弦距離處理的處於坐標上方的兩個簇集合形狀彼此之間互相滲透，互相延伸。
for i = 1:3
clust = find(cidxCos==i);
plot3(meas(clust,1),meas(clust,2),meas(clust,3),ptsymb{i});
hold on
end
hold off
xlabel('Sepal Length');
ylabel('Sepal Width');
zlabel('Petal Length');
view(-137,10);
grid on
上圖並沒有畫出質心，因為余弦距離質心對應於原始數據空間中原點上的的射線（這句話目前不能理解）。然而，你可以在平行坐標軸上繪制標准化后的數據散點圖進而可視化聚類中心之間的差異。
lnsymb = {'b-','r-','m-'};
names = {'SL','SW','PL','PW'};
meas0 = meas ./ repmat(sqrt(sum(meas.^2,2)),1,4); % 將每一個數據除以數據進行平方求和開根號后的結果從而達到標准化
ymin = min(min(meas0));
ymax = max(max(meas0));
for i = 1:3
subplot(1,3,i);
plot(meas0(cidxCos==i,:)',lnsymb{i});
hold on;
plot(cmeansCos(i,:)','k-','LineWidth',2);
hold off;
title(sprintf('Cluster %d',i));
xlim([.9, 4.1]);
ylim([ymin, ymax]);
% 以下部分在R2010a版本中運行不出
h_gca = gca;
h_gca.XTick = 1:4;
h_gca.XTickLabel = names;
end
譯者注：可能是Matlab的版本原因，導致定義結構化數據h_gca的時候報錯，並且橫坐標的刻度的label也沒有成功標記，因此該部分代碼修改如下：
set(gca,'XTick',1:4); set(gca,'XTickLabel',names);
上圖清晰的表明，三個簇中平均每一個花朵樣品的花瓣和萼片大小都不同。第一簇中的花瓣比他們的萼片嚴格小（嚴格遞減）。后面兩類的花瓣和萼片的大小都有重疊部分，也就是說后面兩類的部分花瓣和萼片大小差不多，但是第三類比第二類重疊的部分更多。你也可以看出第二類和大三類的樣品的屬性非常相似。
因為實際上的數據樣本中每一個觀察值所屬的種類都是已知的，因此可以通過將k均值聚類結果和實際分類結果對比發現這三個物種之間是否存在可以用來區分的物理特性。事實上，正如下圖所示，使用余弦距離創建的類於實際分組不同的觀察值只有5個，這5個用星星標記的點都處於坐標系上方兩個類之間的邊界附近。
subplot(1,1,1);
for i = 1:3
clust = find(cidxCos==i);
plot3(meas(clust,1),meas(clust,2),meas(clust,3),ptsymb{i});
hold on
end
xlabel('Sepal Length');
ylabel('Sepal Width');
zlabel('Petal Length');
view(-137,10);
grid on
sidx = grp2idx(species);
miss = find(cidxCos ~= sidx);
plot3(meas(miss,1),meas(miss,2),meas(miss,3),'k*');
legend({'setosa','versicolor','virginica'});
hold off
譯者注：原網站的部分圖片出現錯誤，本次翻譯都已經改正；所有代碼均在Matlab2010a版本上運行無誤，如有翻譯不當地方，歡迎指正。關於第二個層次聚類方法，后面將會繼續更——使用費雪鳶尾花數據實現層次聚類。