クラスター分析
この例では、Statistics and Machine Learning Toolbox™ のクラスター解析を使用して観測値やオブジェクトの類似度と非類似度を調べる方法を示します。多くの場合、データは観測値のグループ (クラスター) に自然に分類されます。同じクラスター内のオブジェクトは特徴が似ており、異なるクラスター間ではオブジェクトの特徴が異なります。
k-means クラスタリングと階層クラスタリング
The Statistics and Machine Learning Toolbox には、k-means クラスタリングおよび階層クラスタリングを実行するための関数が含まれています。
k-means クラスタリングとは、データ内の観測を位置と互いからの距離をもつオブジェクトとして扱う分割法の一種です。k-means クラスタリングでは、オブジェクトを k 個の互いに排他的なクラスターに分割します。同じクラスター内のオブジェクトはできるだけ近くに、別のクラスター内のオブジェクトからはできるだけ遠くにあるように分割します。各クラスターは重心、つまり中心となる点で特徴付けられます。もちろん、クラスタリングで使用される距離が空間的距離を表すことはまれです。
階層クラスタリングは、クラスター木を作成することにより、さまざまな距離のスケールで同時にデータのグループ化を調べる方法です。この木は、k-means とは異なり、複数のクラスターで構成される 1 つの集合ではなく、ある 1 つのレベルに存在する複数のクラスターがまとまって次の 1 段高いレベルでクラスターとなる複数レベルの階層です。これにより、クラスタリングのどのレベルまたはスケールがアプリケーションにおいて最適であるかを決定することが可能になります。
この例で使用する関数の一部では、MATLAB® に組み込まれている乱数発生関数を呼び出します。この例で示された結果を再現するには、次のコマンドを実行して乱数発生器を既知状態に設定する必要があります。状態を設定しないと、些細な点で結果が異なる可能性があります。たとえば、クラスター番号が変わってしまうかもしれません。また、準クラスター解を得ることになる可能性もあります (この例では、最適ではない解の内容やその回避方法などについての説明もします)。
rng(6,'twister')
フィッシャーのアヤメのデータ
1920 年代に植物学者たちは、アヤメの標本 150 個 (3 種について 50 個ずつ) のがく片の長さと幅、花弁の長さと幅に関する測定値を収集しました。この測定値は、「フィッシャーのアヤメのデータ セット」として知られるようになりました。
このデータ セットの各観測値は既知種に由来するので、データをグループ化するための明らかな方法が既に存在しています。差し当たり、種情報は無視して、生の測定値だけを使用してデータをクラスター化することにします。完了したら、結果として得られたクラスターと実際の種を比較し、3 種類のアヤメに明確に識別できる特徴があるかどうかを調べます。
k-means クラスタリングを使用するフィッシャーのアヤメのデータのクラスタリング
関数kmeans
は k-means クラスタリングを実行します。このクラスタリングでは、すべてのクラスターに対して、各オブジェクトからそのオブジェクトが属するクラスターの重心までの距離の総和が最小になるようにオブジェクトをクラスターに割り当てる反復的なアルゴリズムを使用します。この関数をフィッシャーのアヤメのデータに対して実行すると、アヤメの標本ががく片と花弁の測定値に基づいていくつかのグループに自然に分けられます。k-means クラスタリングでは、作成するクラスターの数を指定しなければなりません。
最初に、データを読み込み、関数 kmeans
を呼び出します。その際、目的のクラスター数を 2 に設定し、ユークリッド距離の 2 乗を使用します。結果のクラスターがいかにうまく分離されるかについて知るには、シルエット プロットを作成します。シルエット プロットは、1 個のクラスター中の各点が近隣のクラスター中の点にどれくらい接近しているかについて基準を表示します。
load fisheriris [cidx2,cmeans2] = kmeans(meas,2,'dist','sqeuclidean'); [silh2,h] = silhouette(meas,cidx2,'sqeuclidean');
シルエット プロットから、両方のクラスターの大部分の点が 0.8 よりも大きい大きなシルエットの値をもつことがわかります。これは、それらの点が近隣するクラスターから十分に離れていることを示します。ただし、各クラスターにはシルエット値が低い点もいくつか含まれています。これは、それらの点が他のクラスターの点に近いことを示します。
このデータの 4 番目の測定値である花弁の幅は、3 番目の測定値である花弁の長さと高い相関関係を示しています。したがって、最初の 3 つの測定値の 3 次元プロットは、次元数を 4 にしなくてもデータをうまく表現することになります。kmeans
によって作成された各クラスターに対して異なる記号を使用してデータをプロットすると、シルエット値が小さい点は他のクラスターの点に近い点であることがわかります。
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}); hold on end plot3(cmeans2(:,1),cmeans2(:,2),cmeans2(:,3),'ko'); plot3(cmeans2(:,1),cmeans2(:,2),cmeans2(:,3),'kx'); hold off xlabel('Sepal Length'); ylabel('Sepal Width'); zlabel('Petal Length'); view(-137,10); grid on
各クラスターの重心は、円で囲んだ X でプロットされます。下のクラスター (三角形でプロット) の点のうち 3 つは、上のクラスター (正方形でプロット) の点に非常に近くなっています。上のクラスターは広い範囲に広がっているので、これらの 3 点は、ギャップにより自身のクラスター内の点の塊から離れていますが、上のクラスターの重心よりも下のクラスターの重心に近い位置にあります。k-means クラスタリングでは、距離だけを考慮し密度は考慮されないので、このような結果が生じることがあります。
クラスター数を増やして、関数 kmeans
がデータ内でグループ化構造をさらに見つけることができるかどうかを調べることができます。今回は、オプションの名前と値のペアの引数 'Display'
を使用して、クラスタリング アルゴリズムの各反復に関する情報を出力します。
[cidx3,cmeans3] = kmeans(meas,3,'Display','iter');
iter phase num sum 1 1 150 146.424 2 1 5 144.333 3 1 4 143.924 4 1 3 143.61 5 1 1 143.542 6 1 2 143.414 7 1 2 143.023 8 1 2 142.823 9 1 1 142.786 10 1 1 142.754 Best total sum of distances = 142.754
kmeans
アルゴリズム (アルゴリズムを参照) では、各反復で点から重心までの距離の総和が少なくなるクラスターに点を再度割り当ててから、新しいクラスター割り当てに対してクラスターの重心を計算します。距離の総和と再割り当て回数は、アルゴリズムが最小値に到達するまで、反復ごとに減少していくことに注意してください。関数 kmeans
で使用されるアルゴリズムは、2 段階で構成されます。前の例では、アルゴリズムの第 2 段階で再割り当てがまったく行われませんでした。これは、2、3 回ほどの反復後に第 1 段階で最小値に達したことを示しています。
既定の設定では、関数 kmeans
は無作為に選択された初期の重心位置の集合を使用して処理を開始します。kmeans
アルゴリズムは、局所的最小値である解に収束する可能性があります。つまり、kmeans
は、任意の 1 点を異なるクラスターに移動すると距離の総和が大きくなるようにデータを分割する可能性があります。ただし、他の多くのタイプの数値的最小化と同様に、kmeans
によって到達する解は、開始点に依存する場合があります。したがって、そのデータについて、距離の総和がより小さい他の解 (局所的最小値) が存在する可能性があります。オプションの名前と値のペアの引数 'Replicates'
を使用すると、異なる解をテストできます。複数の複製を指定すると、kmeans
は、無作為に選択した重心から始まるクラスタリング プロセスを各複製について繰り返します。その後、kmeans
は、すべての複製の中で距離の総和が最小である解を返します。
[cidx3,cmeans3,sumd3] = kmeans(meas,3,'replicates',5,'display','final');
Replicate 1, 9 iterations, total sum of distances = 78.8557. Replicate 2, 10 iterations, total sum of distances = 78.8557. Replicate 3, 8 iterations, total sum of distances = 78.8557. Replicate 4, 8 iterations, total sum of distances = 78.8557. Replicate 5, 1 iterations, total sum of distances = 78.8514. Best total sum of distances = 78.8514
出力は、この比較的簡単な問題に対してさえ、大域的でない最小が存在することを示します。この 5 回の繰り返しのそれぞれは、別の初期重心の集合から開始されました。出発点に応じて、関数 kmeans
は 2 つの異なる解のどちらかに到達しました。しかし、kmeans
によって返される最終的な解は、すべての複製のうちで距離の総和が最小値となる解になります。3 番目の出力引数には、その最適解について各クラスター内での距離の総和が格納されます。
sum(sumd3)
ans = 78.8514
この 3 クラスター解のシルエット プロットは、十分に離れたクラスターが 1 つ存在するが、他の 2 つのクラスターはあまり離れていないことを示しています。
[silh3,h] = silhouette(meas,cidx3,'sqeuclidean');
ここで再度、生データをプロットして、関数 kmeans
が点をクラスターに割り当てる様子を確認します。
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
関数 kmeans
が上のクラスターを 2 クラスター解から分離したこと、およびこの 2 つのクラスターが互いに非常に近い位置にあることがわかります。クラスタリング後にこのデータで何を行うのかに応じて、この 3 クラスター解が前の 2 クラスター解より便利なのかそれとも不便なのかが決まります。関数 silhouette
の最初の出力引数には、各点のシルエット値が格納されます。これを使用して 2 つの解を定量的に比較することができます。2 クラスター解の平均シルエット値は大きかったのですが、これは、離れたクラスターを作成するという観点から純粋に見て、より優れた解であることを示しています。
[mean(silh2) mean(silh3)]
ans = 0.8504 0.7357
これらのデータを別の距離でクラスター化することもできます。これらのデータでは、コサイン距離を使用するのが理にかなうかもしれません。測定値の絶対サイズは無視され、相対サイズだけが考慮されるからです。したがって、がく片と花弁のサイズは異なっていても形状が似ていた 2 本の花は、ユークリッド距離の 2 乗に関しては近くないかもしれませんが、コサイン距離に関しては近いと言えます。
[cidxCos,cmeansCos] = kmeans(meas,3,'dist','cos');
シルエット プロットから考えて、これらのクラスターは、ユークリッド距離の 2 乗を使用して見つかったクラスターよりもわずかに離れているように思われます。
[silhCos,h] = silhouette(meas,cidxCos,'cos');
[mean(silh2) mean(silh3) mean(silhCos)]
ans = 0.8504 0.7357 0.7491
クラスターの順序が前のシルエット プロットとは異なることに注意してください。これは、関数 kmeans
が初期クラスター割り当てを無作為に選択するからです。
生データをプロットすることにより、2 つの異なる距離を使用して作成されたクラスター形状の違いを確認することができます。この 2 つの解は似ていますが、2 つの上のクラスターはコサイン距離を使用するときの原点の向きに延長されます。
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]); h_gca = gca; h_gca.XTick = 1:4; h_gca.XTickLabel = names; end
このプロットから、3 つのクラスターそれぞれの標本のがく片と花弁は、平均して相対サイズが明らかに異なります。最初のクラスターには、花弁より明らかに小さいがく片があります。2 番目の 2 つのクラスターのがく片と花弁は、サイズに関して部分的に重なり合いますが、3 番目のクラスターのがく片と花弁の方が、2 番目より大きく重なり合っています。また、2 番目と 3 番目のクラスターには互いに非常に似通った標本がいくつかあるということもわかります。
データ内の各観測の種を知っているので、関数 kmeans
が発見したクラスターと実際の種を比較して、3 種に見分けがつくほど固有の物理的な特徴があるかどうかを確認することができます。実際のところ、次のプロットが示すとおり、コサイン距離を使用して作成されたクラスターは、種グループとわずかに 5 本の花しか違いません。この 5 つの点は星型でプロットされており、すべて上の 2 つのクラスターの境界付近にあります。
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
階層クラスタリングを使用するフィッシャーのアヤメのデータのクラスタリング
k-means クラスタリングでは、アヤメのデータの分割が 1 つしか生成されませんでしたが、データのグループ化のスケールをいろいろに変えて調べてみたいと思うかもしれません。階層クラスタリングを使用すると、複数のクラスターで構成される階層木を 1 つ作成することによりそうすることができます。
最初に、アヤメのデータでの観測間の距離を使用してクラスター木を作成します。ユークリッド距離を使用して開始します。
eucD = pdist(meas,'euclidean'); clustTreeEuc = linkage(eucD,'average');
共表形相関は、クラスター木が元の距離に一致していることを確認するための方法の 1 つです。値が大きい場合、それは観測間の 2 つで 1 ペアの関連性が実際のペアワイズ距離と相関しているという意味で、木が距離をよく近似していることを示しています。この木は、距離をかなりよく近似しているようです。
cophenet(clustTreeEuc,eucD)
ans = 0.8770
クラスターの階層を可視化するには、系統樹をプロットします。
[h,nodes] = dendrogram(clustTreeEuc,0);
h_gca = gca;
h_gca.TickDir = 'out';
h_gca.TickLength = [.002 0];
h_gca.XTickLabel = [];
この木のルート ノードは、残りのノードよりはるかに高い位置にあります。これは、k-means クラスタリングでの結果を裏付けています。つまり、2 つの大きい異なる観測グループが存在するということです。この 2 つのグループそれぞれで、距離のスケールを下げるにつれて、低水準のグループが現れてくるのがわかります。レベル、サイズ、および相違度がさまざまに異なるグループがたくさんあります。
k-means クラスタリングの結果に基づいて考えると、距離計量としてコサインも有望な選択肢かもしれません。結果として得られる階層木はまったく異なるものであり、これはアヤメのデータのグループ構造を見るまったく別の方法があることを示唆しています。
cosD = pdist(meas,'cosine'); clustTreeCos = linkage(cosD,'average'); cophenet(clustTreeCos,cosD)
ans = 0.9360
[h,nodes] = dendrogram(clustTreeCos,0);
h_gca = gca;
h_gca.TickDir = 'out';
h_gca.TickLength = [.002 0];
h_gca.XTickLabel = [];
この木の最高水準により、アヤメの標本が 2 つの非常に異なるグループに分けられます。この系統樹には、コサイン距離に関して、グループ間の違いと比較してグループ内の違いがユークリッド距離の場合よりはるかに小さいことが示されています。これこそまさに、このデータでは当然の結果と言えます。コサイン距離では原点から同じ "向き" にあるオブジェクトのゼロ対距離を計算するからです。
150 個の観測値では、プロットが乱れますが、木の最低水準を表示しない単純化された系統樹を作成することができます。
[h,nodes] = dendrogram(clustTreeCos,12);
この木にある最上位の 3 つのノードによって、3 つのサイズが等しいグループおよび他のいずれにも近くない 1 つの標本 (ラベル表示は葉ノード 5) に分けられます。
[sum(ismember(nodes,[11 12 9 10])) sum(ismember(nodes,[6 7 8])) ...
sum(ismember(nodes,[1 2 4 3])) sum(nodes==5)]
ans = 54 46 49 1
たいていの場合には、系統樹が十分な結果となります。もっとも、さらに 1 歩進んで関数 cluster
を使用して、木を切り分け、k-means の場合のように観測を特定のクラスターに明示的に分割することもできます。クラスターを作成するためのコサイン距離に基づく階層を使用しながら、3 つの最も高いノードの下で木を切る関連性の高さを指定し、クラスターを 4 つ作成し、クラスター化された生データをプロットします。
hidx = cluster(clustTreeCos,'criterion','distance','cutoff',.006); for i = 1:5 clust = find(hidx==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
このプロットは、コサイン距離を使用する階層クラスタリングの結果が、k-means クラスタリングの結果と定性的に類似していることを 3 つのクラスターで示しています。ただし、階層的クラスター木を作成することにより、k-means クラスタリングにおいて K の値をさまざまに変えて実験をかなり繰り返し実行しなければならないようにする要因を一度に可視化することができます。
階層クラスタリングでは、関連性をさまざまに変えて実験することもできます。たとえば、アヤメのデータを 1 つの関連性でクラスター化する場合、平均距離より離れているオブジェクトが連結されがちですが、データ構造のかなり異なる解釈を得ることができます。
clustTreeSng = linkage(eucD,'single'); [h,nodes] = dendrogram(clustTreeSng,0); h_gca = gca; h_gca.TickDir = 'out'; h_gca.TickLength = [.002 0]; h_gca.XTickLabel = [];