集群分析

Allan

2018-09-11

Coding

一般集群分析

要件：

Distance：計算點與點以及群與群之間的相似程度，常見的演算法有：
- 歐式距離
- 曼哈頓距離
Centroid：決定每群中心點的方法

大致可分為：

階層式分群（Hierarchical Clustering）
- 聚合式 Agglomerative（AGNES）
- 分割式 Divisive（DIANA）
- 其他：CURE、BIRCH、CHAMELEON
切割式分群（Partitional Clustering）
- K-means
- K-medoid／PAM（Patition around medoid）

Distance

歐式距離

在 2 維空間(平面)中，點 x=(x1, x2)與點 y=(y1, y2)之間的歐式距離為

一般通式寫成：

曼哈頓距離（絕對誤差）

在 2 維空間(平面)中，點 x=(x1, x2)與點 y=(y1, y2)之間的曼哈頓距離為

一般通式寫成：

切割式分群（Partitional Clustering）

K-means

原理

給定一個資料集，按照樣本之間的距離大小，將資料集切割成 k 群，讓群內的點盡量緊密地連在一起（越小越好），而群與群之間的距離越大越好

流程

給定群數 k，隨機產生 k 個群聚中心點（資料空間內任意值）
計算每個資料點與每個中心點的距離
資料點與其最相近的中心點會被劃分成一群，形成 k 群
利用目前得到的分群重新計算中心點（平均每個點的值）
重複 2~4 的步驟直到收斂（達到最大迭代次數 or 群中心點每次更換的變動距離最小）

圖片來源：http://www.cnblogs.com/skyEva/p/5630820.html

範例（鳶尾花資料集）

鳶尾花資料集(iris dataset)
- Sepal.length（花萼長度）
- Sepal.width（花萼寬度）
- Petal.length（花瓣長度）
- Petal.width（花瓣寬度）
- Species（物種）：山鳶尾(setosa)、變色鳶尾(versicolor)、維吉尼亞鳶尾(virginica)
套件讀取

1 2	library(factoextra) library(cluster)

資料處理

# 讀取資料
data("iris")
# 將特徵標準化以利分群
iris.x <- apply(iris[, -5], 2, function(x){(x-mean(x))/sd(x)})

k-means

# k-means
km.iris <- kmeans(iris.x, centers = 3)
# 分群視覺化
fviz_cluster(km.iris, data = iris.x, geom = "point",
             stand = T, ellipse.type = "norm")

K-medoid / PAM（Patition around medoid）

原理

與 K-means 的做法類似，差別在於==k-means==所選定的中心點為群內所有資料點之平均值，而==k-medoid==則是以群內某個樣本點作為中心點

流程

給定群數 k，隨機選取 k 個群聚中心點（必須是樣本點）
計算每個資料點到每個中心點的距離
資料點與其最相近的中心點會被劃分成一群，形成 k 群
利用目前得到的分群重新計算中心點
1. 計算群內所有樣本點至某一個樣本點的曼哈頓距離和（絕對誤差）
2. 選出使群內絕對誤差距離最小之樣本點作為新的中心點
重複 2~4 的步驟直到收斂（達到最大迭代次數 or 群中心點每次更換的變動距離最小）

範例（鳶尾花資料集）

# K-medoid
pam.iris <- pam(iris.x, k = 3)

# 分群視覺化
fviz_cluster(pam.iris, stand = T, geom = "point",
             ellipse.type = "norm")

階層式分群（Hierarchical Clustering）

原理

透過一種階層架構的方式，將資料層層聚合或分裂，以產生樹狀結構
常見的方式有：

聚合（AGNES）：由樹的底層開始，由下而上將資料或群集合併
分裂（DIANA）：由樹的頂層開始，由上而下將群集逐次分裂
樹狀圖

圖片來源：STHDA

聚合式 Agglomerative（AGNES）

圖片來源：藍鯨的網站分析筆記-層次聚類算法的原理及實現

流程

將每個資料點各自視為一個群集
找出所有群集間，兩個距離最相近的群集
合併兩個群集成一個新的群集
重複步驟 2~3 直至我們所設定的群數

Distance：利用歐式距離判斷資料間的遠近
Centroid：根據所算出之距離，進行資料的聚合，大致有下列方法：
- 單一連結法 Single linkage（最近法）
  - 不同群中最接近之兩點距離
- 完全連結法 Complete linkage（最遠法）
  - 不同群中最遠之兩點距離
- 平均法 Average linkage
  - 不同群中各點間距離總和之平均
- 中心法 Centroid linkage
  - 不同群中之中心點距離
- 華德法 Ward Method
  - 將兩群合併後，各點到合併後的中心點之距離平方和

範例（鳶尾花資料集）

建立距離矩陣

參數	說明
method	euclidean(歐式距離)、manhattan(曼哈頓距離)、canberra、maximum、binary、minkowski

1	dist.res <- dist(iris.x, method = "euclidean") # 歐式距離

Hierarchical clustering

參數	說明
method	single、complete、average、centroid、ward.D2

hc <- hclust(dist.res, method = "complete")

# 分群視覺化
plot(hc)
rect.hclust(hc, k = 3, border = 2:4)

分割式 Divisive（DIANA）

流程

與聚合式（AGNES）相反，不多贅述

範例

請參考上方

最佳分群數

一般有三種方法可供判斷：

Elbow method
Average Silhouette method
Gap statistic method

Elbow method

利用組內變異數 SSE 來衡量分群的好壞
隨著分群數的增加，組內變異數 SSE 會不斷地減少，此方法認為增加分群數，分群效果並不能增強，因此存在一個「拐點」，該點即為最佳分群數
分群數從 1 到 3 下降地最快，之後下降地很慢，因此 3 為最佳分群數
但「拐點」的認定相當主觀，這也是 Elbow method 的最大缺陷
Robert L. Thorndike (December 1953)

k-means

1 2	fviz_nbclust(iris.x, kmeans, method = "wss") + geom_vline(xintercept = 3, linetype = 2)

k-medoid

1 2	fviz_nbclust(iris.x, pam, method = "wss") + geom_vline(xintercept = 3, linetype = 2)

Hierarchical clustering

1 2	fviz_nbclust(iris.x, hcut, method = "wss") + geom_vline(xintercept = 3, linetype = 2)

Average Silhouette method

結合內聚力與分散力，用來衡量群聚效果的好壞
其值越大，代表分群效果越好
Peter J. Rousseeuw (1986)

方法

計算樣本點 i 至同群內其他樣本點的平均距離 ai，ai 越小，代表樣本點 i 越應該被歸類至該群
計算樣本點 i 至其他某群 Cj 的所有樣本之平均距離 bij，bi 越大，代表樣本點 i 越不屬於該群
ai 稱為組內不相似度，bi 稱為組間不相似度，根據兩者定義輪廓係數（Silhouette coefficient）

判斷：
- si 接近 1，代表樣本 i 的分群合理
- si 接近-1，代表樣本 i 更應該被分類至其他群內
- si 接近 0，代表樣本 i 在兩個群的邊界上

k-means

1	fviz_nbclust(iris.x, kmeans, method = "silhouette")

k-medoid

1	fviz_nbclust(iris.x, pam, method = "silhouette")

Hierarchical clustering

1 2	fviz_nbclust(iris.x, hcut, method = "silhouette", hc_method = "complete")

Gap statistic method

將原本 Elbow method 中所提到的缺陷進行改良

當 k 最小，Gap 值最大時，k 即為最佳分群數
R. Tibshirani, G. Walther, and T. Hastie (Standford University, 2001)

k-means

1
2
3

gap_stat <- clusGap(iris.x, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)
fviz_gap_stat(gap_stat)

k-medoid

1 2	gap_stat <- clusGap(iris.x, FUN = pam, K.max = 10, B = 50) fviz_gap_stat(gap_stat)

Hierarchical clustering

1 2	gap_stat <- clusGap(iris.x, FUN = hcut, K.max = 10, B = 50) fviz_gap_stat(gap_stat)

時間序列集群分析

dtwclust 套件使用說明

常見集群方法

階層式分群（Hierarchical Clustering）
- 聚合式 Agglomerative（AGNES）
- 分割式 Divisive（DIANA）
- 其他：Any other suitable time-series centroid method
切割式分群（Partitional Clustering）
- PAM（Patition around medoid）
- DBA+DTW

要件

Distance：
- 動態時間校正 DTW（Dynamic time warping distance）
  - 改寫自歐式距離
  - 找出一條路徑最短，使得兩不同長度序列的距離最短
  - 找出波形之間的相似度
  - 原理
Centroid：
- PAM
- DBA（DTW barycenter averaging）
  - 原理第 15 頁

搭配用法

Prototyping function	Distance	Algorithm
PAM	歐式距離／DTW	Time-series with minimum sum of distances to the other series in the group
DBA	DTW	Average of points grouped according to DTW alignments

切割式分群（Partitional Clustering）

PAM

原理 & 流程

基本上與一般集群分析的 k-medoid 雷同，請參考上方說明

範例（以台灣市值前 50 大股票為例）

資料預處理

# 套件讀取（請自行安裝）
library(tidyverse)
library(tidyquant)
library(dtwclust)
library(cluster)

# 讀取資料
data <- read.table("stockdata.txt", sep = "\t", header = T, stringsAsFactors = F)
# 欄位命名
colnames(data) <- c("code", "name", "date", "open", "high", "low", "close", "volume", "mv")

# 取出2018年1月2日市值前50大的股票
top <- data %>%
  arrange(code, date) %>%
  group_by(code) %>%
  filter(date==20180102) %>%
  arrange(desc(mv)) %>%
  group_by() %>%
  slice(., 1:50) %>%
  pull(name)

# 建立日頻報酬率table
dayLogRet <- data %>%
  arrange(code, date) %>%
  group_by(name) %>%
  filter(name %in% top) %>%
  mutate(dayLogRet=log(close/lag(close, 1))) %>%
  group_by() %>%
  select(code, name, date, dayLogRet) %>%
  na.omit()

##################### 重要 #####################

# 將資料轉成 dtwclust套件的讀取格式
symbol <- unique(dayLogRet$name)
stockRetData <- lapply(symbol, function(iy){
  stockRet <- dayLogRet %>% filter(name==iy) %>% .$dayLogRet
  return(stockRet)
})
names(stockRetData) <- symbol

# 需先將資料進行標準化，以利更好的分群結果
stockRetData_z <- zscore(stockRetData)

pc_k <- tsclust(stockRetData_z,
                k=3,
                distance="dtw_basic",
                centroid = "pam")

# 儲存分群資訊
pCluster <- data.frame(name=names(pc_k@datalist), group=pc_k@cluster)
pCluster <- pCluster %>% arrange(group)

參數說明

可調整參數	參數值	說明
series	stockRetData_z	必須將 data.frame 轉成 list 格式放入
k	3	群數，通常依經驗決定
distance	“dtw_basic”	距離演算法，預設為歐式距離，若為時間序列，建議使用 dtw_basic
centroid	“pam”	決定每群中心點的方法

畫圖

1	plot(pc_k)

DBA+DTW

原理 & 流程

請見上方 DBA 與 DTW 原理

範例（以台灣市值前 50 大股票為例）

DBA+DTW 分群

dba_k <- tsclust(stockRetData_z,
                 k=3,
                 distance = "dtw_basic",
                 centroid = "dba")

# 儲存分群資訊
dbCluster <- data.frame(name=names(dba_k@datalist), group=dba_k@cluster)
dbCluster <- dbCluster %>% arrange(group)

參數說明

可調整參數	參數值	說明
series	stockRetData_z	必須將 data.frame 轉成 list 格式放入
k	3	群數，通常依經驗決定
distance	“dtw_basic”	距離演算法，預設為歐式距離，若為時間序列，建議使用 dtw_basic
centroid	“dba”	決定每群中心點的方法

畫圖

1	plot(dba_k)

階層式分群（Hierarchical Clustering）

原理 & 流程

基本上與一般集群分析方法雷同，請參考上方說明

範例（以台灣市值前 50 大股票為例）

階層式分群（AGNES）

hcaCluster <- tsclust(stockRetData_z,
                      type="h",
                      k=3,
                      distance="dtw_basic",
                      control=hierarchical_control(method=agnes))

參數說明

可調整參數	參數值	說明
series	stockRetData_z	必須將 data.frame 轉成 list 格式放入
type	“h”	h 為階層式集群的代號
k	3	群數，通常依經驗決定
distance	“dtw_basic”	距離演算法，通常用 dtw_basic 即可
control	hierarchical_control(method=agnes))	method 可使用 agnes 或 diana

樹狀圖（AGNES）

1	plot(hcaCluster)

階層式分群（DIANA）

hcaCluster <- tsclust(stockRetData_z,
                      type="h",
                      k=3,
                      distance="dtw_basic",
                      control=hierarchical_control(method=diana))

樹狀圖（DIANA）

1	plot(hcaCluster)

最佳分群數（CVI，Cluster validity indices）

利用 dtwclust 套件中的 cvi 函式
cvi 函式中有多項指標，詳細請見dtwclust 套件使用說明第 24 頁
以下僅示範 cvi 函式的用法

PAM

# 需要設定多群的k
pc_cvi <- tsclust(stockRetData_z,
                  k=2:10,
                  distance = "dtw_basic",
                  centroid = "pam")
# 將結果命名
names(pc_cvi) <- paste0("k_",2:10)

# 使用sapply搭配cvi函式同時計算多群的衡量效果
pc_eva <- sapply(pc_cvi, cvi, type = "internal")

# 使用Average Silhouette method衡量
plot(x=2:10, pc_eva[c("Sil"),], type="l")

DBA+DTW

# 需要設定多群的k
dba_k <- tsclust(stockRetData_z,
                 k=2:10,
                 distance = "dtw_basic",
                 centroid = "dba")
# 將結果命名
names(dba_cvi) <- paste0("k_",2:10)

# 使用sapply搭配cvi函式同時計算多群的衡量效果
dba_eva <- sapply(dba_cvi, cvi, type = "internal")

# 使用Average Silhouette method衡量
plot(x=2:10, dba_eva[c("Sil"),], type="l")

Hierarchical clustering

# 需要設定多群的k
hcaCluster <- tsclust(stockRetData_z,
                      type="h",
                      k=2:10,
                      distance="dtw_basic",
                      control=hierarchical_control(method=diana))

# 將結果命名
names(hcaCluster) <- paste0("k_",2:10)

# 使用sapply搭配cvi函式同時計算多群的衡量效果
hca_eva <- sapply(hcaCluster, cvi, type = "internal")

# 使用Average Silhouette method衡量
plot(x=2:10, hca_eva[c("Sil"),], type="l")

一般集群分析

要件：

大致可分為：

Distance

歐式距離

曼哈頓距離（絕對誤差）

切割式分群（Partitional Clustering）

K-means

原理

流程

範例（鳶尾花資料集）

K-medoid / PAM（Patition around medoid）

原理

流程

範例（鳶尾花資料集）

階層式分群（Hierarchical Clustering）

原理

聚合式 Agglomerative（AGNES）

流程

範例（鳶尾花資料集）

分割式 Divisive（DIANA）

流程

範例

最佳分群數

Elbow method

k-means

k-medoid

Hierarchical clustering

Average Silhouette method

方法

k-means

k-medoid

Hierarchical clustering

Gap statistic method

k-means

k-medoid

Hierarchical clustering

時間序列集群分析

常見集群方法

要件

搭配用法

切割式分群（Partitional Clustering）

PAM

原理 & 流程

範例（以台灣市值前 50 大股票為例）

DBA+DTW

原理 & 流程

範例（以台灣市值前 50 大股票為例）

階層式分群（Hierarchical Clustering）

原理 & 流程

範例（以台灣市值前 50 大股票為例）

最佳分群數（CVI，Cluster validity indices）

PAM

DBA+DTW

Hierarchical clustering

參考資料

Reference