クラスタリングを手計算することは理解を深める上で有用です。Rの出力を検算に用いれば良いと思うのですが、ウォード法の出力は手計算と一致してくれません。理由は、ユークリッド距離との整合性を保つため、Rではウォード距離の2倍の平方根を出力しているようです。以下、調べたことを記します。
参考にさせていただいたのは、
東京工業大学情報理工学院 機械学習 CSC.T254 岡崎 直観 先生の教材です。.

参考サイト.

1 距離.

Rでクラスター分析を行う関数は “hclust” です。

   hclust(dist(wine_data2), method = “@”).

クラスター間距離構造の手法を指定し、階層的クラスター分析を行います。@には

  • ”ward.D2″(ウォード法)
  • “single”(最短距離法)
  • “complete”(最長距離法)
  • “average” (群平均法)
  • “centroid” (重心法)

のいずれかを選択して入れます。.

2 ウォード法.

ウォード法は凝集型クラスタリングの手順でクラスタ内平方和の最小化を目指す。

(クラスター間変動を最大化するともいえる).

2.1 クラスター間距離.

ウォード法によるクラスター間距離は、
  \(d(C_1,C_2)=L(C_1\cup C_2)-L(C_1)-L(C_2)\)  ただし \(L(C)=\displaystyle\sum_{x\in C}D(x,g_C)^2\).
クラスター間距離は群間変動を意味し、 \(L(C)=\displaystyle\sum_{x\in C}D(x,g_C)^2\) によって計算される。
距離の近いクラスターを結合する方針に変わりはない。

特に、2点間のウォード距離は.
\(d_{ward}(C_1,C_2)=\left(c_1-\dfrac{c_1+c_2}{2}\right)^2+\left(c_2-\dfrac{c_1+c_2}{2}\right)^2=\dfrac{1}{2}(c_1-c_2)^2=\dfrac{1}{2}C_1C_2~^2\).
   \(\ C_1C_2=\sqrt{2\ d_{ward}(C_1,C_2)}\):ウォード距離の2倍の平方根がユークリッド距離.

ウォード距離が最短ならユークリッド距離も最短であるから、1番目に結合される元はどのメソッドでも同じである。
コーフェン行列は個体間の距離を要素に持つ行列であるが、ウォード距離のかわりに\(\sqrt{2\ d_{ward}(C_1,C_2)}\) を記入することにすれば、1番目の距離との整合性が保たれる。

ウォード法によるクラスター間距離は.
   \(d(C_1,C_2)=\dfrac{n_1n_2}{n_1+n_2}D(g_{C_1},g_{C_2})^2\).
で計算される。
これは、
   \(V(X)=p_1\sigma_1^2 +p_2\sigma_2^2+p_1p_2(\mu_1-\mu_2)^2\)  から導出される。.

実際、
   \(p_1p_2(\mu_1-\mu_2)^2=\dfrac{n_1n_2}{(n_1+n_2)^2}(g_{C_1}-g_{C_2})^2\).
の両辺に\((n_1+n_2)\)を乗ずれば偏差平方和になり、ウォード法によるクラスター間距離となる。

2.2 クラスター距離の更新.

堀田啓介先生「クラスター分析ノート」 によれば、
クラスター \(p,~q\) を結合したクラスターを \(t\)とする。
   データの個数を \(n_p,n_q\) とすると、\(n_t=n_p+n_q\)  \(\therefore~\dfrac{n_p}{n_t}+\dfrac{n_q}{n_t}=1\)

   郡内平均は \(\bar{x}_t=\dfrac{n_p}{n_t}\bar{x}_p+\dfrac{n_q}{n_t}\bar{x}_q\) (加重平均)

クラスター \(r\)\(t\) の距離 \(S_{tr}\)\(S_{pr}~,~S_{qr}~,~S_{pq}\) を用いて求めると.
   \(S_{tr} = \dfrac {1}{n_{p} + n_{q} + n_{r}} \{(n_{p} + n_{r}){S_{pr}} + (n_{q} + n_{r})S_{qr} - n_{r}{S_{pq}}\}\).

注意
次元が増えれば、各成分の平方和を合計することになり、1次元で成立した公式は次元が増えても成立する。
(導出の詳細は、当サイト「直前問題 クラスター分析」参照).

3 人文2018問題4.

3.1 データファイル.

par(family = "HiraKakuProN-W3") #日本語フォントの指定
sake<-matrix(c(51,6.7,9.8,9.8,2.4,45.5,6.3,8.2,1.8,0.9,40.7,5,6.6,3.8,1.4,35.9,13.3,7.8,2.5,1.4,34.3,7,10.7,2.5,1.7,33.7,2.5,20.3,2.5,0.8),5,6,byrow = FALSE);sake
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 51.0 45.5 40.7 35.9 34.3 33.7
## [2,]  6.7  6.3  5.0 13.3  7.0  2.5
## [3,]  9.8  8.2  6.6  7.8 10.7 20.3
## [4,]  9.8  1.8  3.8  2.5  2.5  2.5
## [5,]  2.4  0.9  1.4  1.4  1.7  0.8
colnames(sake)<- c("東京","高知","大阪","新潟","青森","宮崎")
rownames(sake)<- c("ビール","日本酒","焼酎","ワイン","ウィスキー")
sake
##            東京 高知 大阪 新潟 青森 宮崎
## ビール     51.0 45.5 40.7 35.9 34.3 33.7
## 日本酒      6.7  6.3  5.0 13.3  7.0  2.5
## 焼酎        9.8  8.2  6.6  7.8 10.7 20.3
## ワイン      9.8  1.8  3.8  2.5  2.5  2.5
## ウィスキー  2.4  0.9  1.4  1.4  1.7  0.8
class(sake)
## [1] "matrix" "array"
sake<-t(sake)#転置
sake
##      ビール 日本酒 焼酎 ワイン ウィスキー
## 東京   51.0    6.7  9.8    9.8        2.4
## 高知   45.5    6.3  8.2    1.8        0.9
## 大阪   40.7    5.0  6.6    3.8        1.4
## 新潟   35.9   13.3  7.8    2.5        1.4
## 青森   34.3    7.0 10.7    2.5        1.7
## 宮崎   33.7    2.5 20.3    2.5        0.8

3.2 距離行列.

sake.d <- dist(sake)
sake.d
##           東京      高知      大阪      新潟      青森
## 高知  9.960924                                        
## 大阪 12.498800  5.616048                              
## 新潟 18.162048 11.918893  9.749872                    
## 青森 18.263899 11.545995  7.971825  7.123903          
## 宮崎 21.977943 17.337532 15.652156 16.676031 10.657392

3.3 デンドログラム.

par(family = "HiraKakuProN-W3") #日本語フォントの指定
sake.hc3<-hclust(sake.d,method="ward.D2")
round(cophenetic(sake.hc3) ,1)
##      東京 高知 大阪 新潟 青森
## 高知 12.6                    
## 大阪 12.6  5.6               
## 新潟 21.6 21.6 21.6          
## 青森 21.6 21.6 21.6  7.1     
## 宮崎 21.6 21.6 21.6 15.6 15.6
plot(sake.hc3,hang=-1,main="ウォード法")

3.4 クラスター距離の計算.

3.4.1 (東京,高知,大阪)と(新潟,青森,宮崎).

group <- sake
L12 <- sum(diag(var(group)))*5
group1 <- rbind(sake[1,],sake[2,],sake[3,])
L1 <- sum(diag(var(group1)))*2
group2 <- rbind(sake[4,],sake[5,],sake[6,])
L2 <-sum(diag(var(group2)))*2
L <- L12-L1-L2
sqrt(2*L)
## [1] 21.57638

クラス間変動から算出したward距離を2倍し平方根をとって整合性を持たせた値がRで出力されている。
距離の更新公式から算出する計算は、当サイト「直前問題 クラスター分析」参照。

4 人文2019問題2.

man<-matrix(c(1,3,2,2,3,4,7,9,8,7,9,9),6,2,byrow = TRUE);man
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    2
## [3,]    3    4
## [4,]    7    9
## [5,]    8    7
## [6,]    9    9
colnames(man)<- c("築年数","所要時間")
rownames(man)<- c("A","B","C","D","E","F")
man
##   築年数 所要時間
## A      1        3
## B      2        2
## C      3        4
## D      7        9
## E      8        7
## F      9        9
summary(man)#築年数も所要時間も分布はほぼ同じなので標準化の必要はない
##      築年数        所要時間    
##  Min.   :1.00   Min.   :2.000  
##  1st Qu.:2.25   1st Qu.:3.250  
##  Median :5.00   Median :5.500  
##  Mean   :5.00   Mean   :5.667  
##  3rd Qu.:7.75   3rd Qu.:8.500  
##  Max.   :9.00   Max.   :9.000
man.d<- dist(man);man.d
##           A         B         C         D         E
## B  1.414214                                        
## C  2.236068  2.236068                              
## D  8.485281  8.602325  6.403124                    
## E  8.062258  7.810250  5.830952  2.236068          
## F 10.000000  9.899495  7.810250  2.000000  2.236068
(man.hc<-hclust(man.d,method = "single"))#最短法、デフォルトはcomplete最遠隣法
## 
## Call:
## hclust(d = man.d, method = "single")
## 
## Cluster method   : single 
## Distance         : euclidean 
## Number of objects: 6
summary(man.hc)
##             Length Class  Mode     
## merge       10     -none- numeric  
## height       5     -none- numeric  
## order        6     -none- numeric  
## labels       6     -none- character
## method       1     -none- character
## call         3     -none- call     
## dist.method  1     -none- character
man.hc$merge
##      [,1] [,2]
## [1,]   -1   -2
## [2,]   -4   -6
## [3,]   -3    1
## [4,]   -5    2
## [5,]    3    4
man.hc$height
## [1] 1.414214 2.000000 2.236068 2.236068 5.830952
man.hc$order
## [1] 3 1 2 5 4 6
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(man.hc,hang=-1,main="最短法") #下を揃える

cophenetic(man.hc)#コーフェン行列
##          A        B        C        D        E
## B 1.414214                                    
## C 2.236068 2.236068                           
## D 5.830952 5.830952 5.830952                  
## E 5.830952 5.830952 5.830952 2.236068         
## F 5.830952 5.830952 5.830952 2.000000 2.236068

4.0.1 k-means法.

km <- kmeans(man,2)#3つのクラスター
result <- km$cluster
result
## A B C D E F 
## 1 1 1 2 2 2
library(cluster)
clusplot(man, km$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)

5 準1級2017_問題13.

5.1 データファイル.

dat <- data.frame(Math=c(10,10,9,7,6),Eng=c(10,8,6,8,7))
library(magrittr)
rownames(dat) <- c(1:5)
dat
Math Eng
10 10
10 8
9 6
7 8
6 7

5.2 距離行列.

dat.d <- dist(dat)
round(dat.d,1)
##     1   2   3   4
## 2 2.0            
## 3 4.1 2.2        
## 4 3.6 3.0 2.8    
## 5 5.0 4.1 3.2 1.4

5.3 クラスタリング.

併合によりクラスターができたら、新クラスターとの距離を新たに計算して更新する。
更新の計算方法が「最小距離法」「最大距離法」「群平均法」「ウォード法」で異なっている。その違いに注意してみていこう。

5.3.1 最短距離法.

dat.hcs <- hclust(dat.d, method = "single")
dat.hcs
## 
## Call:
## hclust(d = dat.d, method = "single")
## 
## Cluster method   : single 
## Distance         : euclidean 
## Number of objects: 5
round(cophenetic(dat.hcs) ,1)
##     1   2   3   4
## 2 2.0            
## 3 2.2 2.2        
## 4 2.8 2.8 2.8    
## 5 2.8 2.8 2.8 1.4
plot(dat.hcs, hang = -1)

5.3.2 k-means 法.

https://qiita.com/Haruka-Ogawa/items/fcda36cc9060ba851225

Toshihiro Kamishima 様

#kmeans(x, centers, iter.max = 10, nstart = 1,algorithm =c("Hartigan-Wong", "Lloyd", "Forgy", "MacQueen"))
km <- kmeans(dat,2)#3つのクラスター
result <- km$cluster
result
## 1 2 3 4 5 
## 1 1 2 2 2
library(cluster)
clusplot(dat, km$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)