Rで機械学習:k-means

library(ggplot2)
library(cluster)
library(factoextra)
library(gridExtra)

# ---------- Elbow法の推奨k検出（二次多項式フィット） ----------
detect_elbow <- function(sse, k_range) {
  fit <- lm(sse ~ poly(k_range, 2, raw = TRUE))
  pred <- predict(fit, newdata = data.frame(k_range = k_range))
  residuals <- sse - pred
  elbow_k <- k_range[which.min(residuals)]
  return(elbow_k)
}

# ---------- 推奨kを求めてプロットする関数 ----------
find_optimal_k <- function(data, k_min = 2, k_max = 10, B = 50, seed) {
  set.seed(seed)
  k_range <- k_min:k_max

  # 1. Elbow法（SSE）
  sse <- numeric(length(k_range))
  for (i in seq_along(k_range)) {
    km <- kmeans(data, centers = k_range[i], nstart = 20)
    sse[i] <- km$tot.withinss
  }
  elbow_k <- detect_elbow(sse, k_range)
  df_sse <- data.frame(k = k_range, value = sse, method = "SSE (Elbow)")

  # 2. シルエット係数
  sil <- numeric(length(k_range))
  for (i in seq_along(k_range)) {
    km <- kmeans(data, centers = k_range[i], nstart = 20)
    ss <- silhouette(km$cluster, dist(data))
    sil[i] <- mean(ss[, 3])
  }
  sil_k <- k_range[which.max(sil)]
  df_sil <- data.frame(k = k_range, value = sil, method = "Silhouette")

  # 3. ギャップ統計量
  gap_stat <- clusGap(data, FUN = kmeans, K.max = k_max, B = B, nstart = 20)
  gap_values <- gap_stat$Tab[, "gap"]
  gap_se <- gap_stat$Tab[, "SE.sim"]
  gap_k <- maxSE(f = gap_values, SE.f = gap_se, method = "firstSEmax")
  df_gap <- data.frame(k = 1:k_max, value = gap_values, method = "Gap Statistic")
  df_gap <- df_gap[df_gap$k %in% k_range, ]

  # 4. 評価指標比較グラフ
  df_all <- rbind(df_sse, df_sil, df_gap)
  plot_metrics <- ggplot(df_all, aes(x = k, y = value, color = method)) +
    geom_line() +
    geom_point() +
    scale_x_continuous(breaks = k_range) +
    facet_wrap(~method, scales = "free_y") +
    labs(
      title = "クラスタ数選択の3手法比較",
      x = "クラスタ数 k",
      y = "指標値"
    ) +
    theme_minimal()

  # 5. クラスタリング前のデータ
  plot_original <- ggplot(data, aes(x = x, y = y)) +
    geom_point(alpha = 0.7) +
    labs(
      title = "クラスタリング前のデータ",
      x = "X", y = "Y"
    ) +
    theme_minimal()

  # 6. 各推奨kでクラスタリング＆散布図
  make_cluster_plot <- function(k_value, method_name) {
    km <- kmeans(data, centers = k_value, nstart = 20)
    cluster_df <- data.frame(data, cluster = factor(km$cluster))
    ggplot(cluster_df, aes(x = x, y = y, color = cluster)) +
      geom_point(size = 2, alpha = 0.7) +
      labs(
        title = paste(method_name, " 推奨k =", k_value),
        x = "X", y = "Y"
      ) +
      theme_minimal() +
      scale_color_brewer(palette = "Set1")
  }

  plot_elbow <- make_cluster_plot(elbow_k, "Elbow法")
  plot_sil <- make_cluster_plot(sil_k, "Silhouette法")
  plot_gap <- make_cluster_plot(gap_k, "Gap統計量")

  # 7. 3つのクラスタリング結果をまとめる
  # plot_clusters_all <- grid.arrange(plot_elbow, plot_sil, plot_gap, ncol = 3) では出力される
  # プロットをリストにまとめる
  plot_clusters_all <- list(plot_elbow, plot_sil, plot_gap)

  return(list(
    recommended_k = list(
      Elbow = elbow_k,
      Silhouette = sil_k,
      Gap_Statistic = gap_k
    ),
    plot_original = plot_original,
    metrics_plot = plot_metrics,
    cluster_plots = plot_clusters_all
  ))
}

# ---------- データ生成 ----------

seed <- 20250814
set.seed(seed)
n <- 100
group1 <- data.frame(
  x = rnorm(n, mean = 0, sd = 0.5),
  y = rnorm(n, mean = 0, sd = 0.5)
)
group2 <- data.frame(
  x = rnorm(n, mean = 3, sd = 0.5),
  y = rnorm(n, mean = 0, sd = 0.5)
)
group3 <- data.frame(
  x = rnorm(n, mean = 1.5, sd = 0.5),
  y = rnorm(n, mean = 2.5, sd = 0.5)
)
data <- rbind(group1, group2, group3)

# ---------- 実行 ----------
result <- find_optimal_k(data, k_min = 2, k_max = 10, B = 50, seed = seed)

# クラスタリング前のデータ
print(result$plot_original)

# 指標比較グラフ
print(result$metrics_plot)

# プロットを出力
result$cluster_plots