(标签-R语言|关键词-时间序列)

Question

请对如上数据使用R语言，做时间序列聚类（基于DTW的K-means算法），包括肘部法则（轮廓系数）确定最佳聚类数目。

Answer 1

引自免费微信小程序：皆我百晓生

为了进行基于DTW（动态时序 warping）的K-means聚类，并利用肘部法则和轮廓系数确定最佳聚类数，你需要安装一些必要的R包，例如tsclust, ggplot2, 和 dbscan. 首先确保你已经安装了这些包，如果没有，可以使用以下命令安装：

install.packages(c("tsclust", "ggplot2", "dbscan"))

接下来，我们将导入所需的库并处理数据。由于提供的数据是文本格式，我们首先将它转换为一个矩阵，然后将其转换为时间序列对象。

# 导入所需库
library(tsclust)
library(ggplot2)

# 将数据转换为矩阵
data_matrix <- matrix(as.numeric(read.table(text = paste(data, collapse = "\n"), header = FALSE)), ncol = 14, byrow = TRUE)

# 定义时间序列对象
time_series <- ts(data_matrix[, -1], start = c(1962, 1), frequency = 1)

# 按照行名创建索引
country_index <- data_matrix[, 1]

# 对于每个国家的时间序列进行标准化
time_series_normalized <- lapply(time_series, scale)

现在我们可以开始执行基于DTW的K-means聚类。由于DTW的计算复杂度较高，这里我们选择较小的聚类数量（例如5到15个）进行尝试。

# 聚类函数
dtw_kmeans <- function(n_clusters) {
  dtw_clust <- kmeans(dtwDist(time_series_normalized), centers = n_clusters, iter.max = 1000)
  return(list(clusters = dtw_clust$cluster, centroids = time_series_normalized[dtw_clust$centers, , drop = FALSE]))
}

# 确定最佳聚类数
clusters_range <- 5:15
results <- sapply(clusters_range, dtw_kmeans)

接下来，我们需要评估聚类效果。首先绘制肘部法则图，然后计算轮廓系数。

# 肘部法则图
elbow_plot <- ggplot(data.frame(n_clusters = clusters_range, withinss = sapply(results, function(x) x$withinss[[1]]))) +
  geom_line(aes(x = n_clusters, y = withinss)) +
  labs(title = "Elasticity Rule for DTW K-Means Clustering",
       x = "Number of Clusters",
       y = "Within Sum of Squares") +
  theme_minimal()

print(elbow_plot)

# 计算轮廓系数
silhouette_scores <- sapply(results, function(x) silhouette::silhouette(x$clusters, distmat = dtwDist(time_series_normalized)))
best_n_clusters <- which.max(silhouette_scores$n)

根据肘部法则图和轮廓系数，你可以决定最佳的聚类数。假设最佳聚类数为best_n_clusters，则可以输出最终的聚类结果。

# 输出最佳聚类结果
final_clusters <- results[[best_n_clusters]]
final_centroids <- final_clusters$centroids

# 打印轮廓系数
cat("Best number of clusters based on the silhouette method:", best_n_clusters, "\n")
cat("Silhouette score:", silhouette_scores[names(silhouette_scores) == best_n_clusters], "\n")

# 打印聚类结果
cat("\nCluster assignments:\n")
print(final_clusters$clusters)

# 打印聚类中心
cat("\nCentroids of each cluster:\n")
for (i in seq_along(final_centroids)) {
  cat(paste("Cluster", i, ":", country_index[final_clusters$centers[i]], "\n"))
  print(final_centroids[[i]])
}

这将给出最佳聚类数以及对应的聚类结果和聚类中心。注意，实际应用时可能需要根据具体的数据调整聚类范围、迭代次数等参数以获得更优的结果。