【ChatGPT】クラスター分析のクラスター数を自動で決定するコードを書いてもらった

クラスター分析を自動で行うアプリを作成したいため、クラスター分析のクラスター数決定を自動で行うコードをChatGPTに書いてもらいました。

エルボー法を利用した自動決定

エルボー法を用いたクラスタ数の自動決定は主観的な部分が含まれるため、完全な自動化は難しいですが、WCSSの変化率を計算することで、一定の基準に基づいて最適と思われるクラスタ数を自動で選ぶことは可能です。

以下のコードは、WCSSの変化率を計算し、その変化率が最も大きい位置（エルボーの位置に相当）をトップ3として選択、それぞれクラスタリングを実施するものです。また、WCSSの変化率のプロットを表示し、変化率が最も大きいトップ3のクラスタ数を赤い破線で示します。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# サンプルデータを生成（適宜、実際のデータを読み込む）
data = {
    'ユーザーID': np.arange(1, 101),
    '特徴量1': np.random.randint(1000, 10000, 100),
    '特徴量2': np.random.randint(1, 10, 100),
    '特徴量3': np.random.rand(100) * 10,
    '特徴量4': np.random.randint(1, 5, 100),
    '特徴量5': np.random.randint(1000, 5000, 100),
    '特徴量6': np.random.randint(500, 5000, 100)
}

df = pd.DataFrame(data)

# 'ユーザーID'はクラスタリングに使用しないため、除外
df_cluster = df.drop(columns=['ユーザーID'])

# データを標準化
scaler = StandardScaler()
df_std = scaler.fit_transform(df_cluster)

# クラスタ数の最適化（エルボー法）
wcss = []
max_clusters = 10
for i in range(1, max_clusters + 1):
    kmeans = KMeans(n_clusters=i, init='k-means++', random_state=42)
    kmeans.fit(df_std)
    wcss.append(kmeans.inertia_)

# WCSSの変化率を計算
wcss_diff_percent = np.diff(wcss) / wcss[:-1] * 100
top_n = 3
top_clusters = np.argsort(wcss_diff_percent)[:top_n] + 2  # "+2" は0始まりのindexとクラスタ数のoffsetの調整

# エルボーのプロットと最適と思われるクラスタ数の表示
plt.figure(figsize=(10, 5))
plt.plot(range(2, max_clusters + 1), wcss_diff_percent, marker='o', linestyle='--')
for cluster in top_clusters:
    plt.axvline(x=cluster, color='r', linestyle='--')
plt.title('WCSS Percent Change')
plt.xlabel('Number of clusters')
plt.ylabel('WCSS % Change')
plt.show()

for idx, cluster_num in enumerate(top_clusters, 1):
    kmeans = KMeans(n_clusters=cluster_num, init='k-means++', random_state=42)
    df[f'cluster_top{idx}'] = kmeans.fit_predict(df_std)

print(df.head())

シルエット法を利用した自動決定

シルエット法を用いてクラスタ数を最適化し、その結果をもとにクラスタリングを行うPythonのコードを以下に示します。

シルエットスコアの変化を示すグラフを表示し、スコアが最も高いトップ1, 2, 3のクラスタ数を赤い破線で示します。そして、それぞれのクラスタ数でクラスタリングを実施し、その結果をデータフレームに追加して表示します。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

# サンプルデータを生成（適宜、実際のデータを読み込む）
data = {
    'ユーザーID': np.arange(1, 101),
    '特徴量1': np.random.randint(1000, 10000, 100),
    '特徴量2': np.random.randint(1, 10, 100),
    '特徴量3': np.random.rand(100) * 10,
    '特徴量4': np.random.randint(1, 5, 100),
    '特徴量5': np.random.randint(1000, 5000, 100),
    '特徴量6': np.random.randint(500, 5000, 100)
}

df = pd.DataFrame(data)

# 'ユーザーID'はクラスタリングに使用しないため、除外
df_cluster = df.drop(columns=['ユーザーID'])

# データを標準化
scaler = StandardScaler()
df_std = scaler.fit_transform(df_cluster)

# シルエット法を用いたクラスタ数の最適化
silhouette_scores = []
max_clusters = 10
for i in range(2, max_clusters + 1):  # シルエットスコアはクラスタ数が1の場合は計算できないので、2から始める
    kmeans = KMeans(n_clusters=i, init='k-means++', random_state=42)
    kmeans_labels = kmeans.fit_predict(df_std)
    score = silhouette_score(df_std, kmeans_labels)
    silhouette_scores.append(score)

# 最適と思われるクラスタ数をトップ3で選択
top_n = 3
top_clusters = np.argsort(silhouette_scores)[-top_n:] + 2  # "+2" は0始まりのindexとクラスタ数のoffsetの調整

# シルエットスコアのプロットと最適と思われるクラスタ数の表示
plt.figure(figsize=(10, 5))
plt.plot(range(2, max_clusters + 1), silhouette_scores, marker='o', linestyle='--')
for cluster in top_clusters:
    plt.axvline(x=cluster, color='r', linestyle='--')
plt.title('Silhouette Scores')
plt.xlabel('Number of clusters')
plt.ylabel('Silhouette Score')
plt.show()

for idx, cluster_num in enumerate(top_clusters, 1):
    kmeans = KMeans(n_clusters=cluster_num, init='k-means++', random_state=42)
    df[f'cluster_top{idx}'] = kmeans.fit_predict(df_std)

print(df.head())

ギャップ統計量を利用した自動決定

ギャップ統計量（Gap Statistic）を計算するためには、実データのクラスタ内誤差平方和と、そのデータをランダムにサンプリングした場合の誤差平方和との差（ギャップ）を求めます。この差が最も大きくなるクラスタ数を最適と判断します。

以下は、ギャップ統計量を計算してクラスタ数を最適化し、その結果をもとにクラスタリングを行うPythonのコード例です。

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from gap_statistic import OptimalK
from sklearn.cluster import KMeans

# サンプルデータを生成（適宜、実際のデータを読み込む）
data = {
    'ユーザーID': np.arange(1, 101),
    '特徴量1': np.random.randint(1000, 10000, 100),
    '特徴量2': np.random.randint(1, 10, 100),
    '特徴量3': np.random.rand(100) * 10,
    '特徴量4': np.random.randint(1, 5, 100),
    '特徴量5': np.random.randint(1000, 5000, 100),
    '特徴量6': np.random.randint(500, 5000, 100)
}

df = pd.DataFrame(data)

# 'ユーザーID'はクラスタリングに使用しないため、除外
df_cluster = df.drop(columns=['ユーザーID'])

# データを標準化
scaler = StandardScaler()
df_std = scaler.fit_transform(df_cluster)

# ギャップ統計量を用いて最適なクラスタ数を求める
optimalK = OptimalK(parallel_backend='rust')
n_clusters = optimalK(df_std, cluster_array=np.arange(1, 11))
print(f"Optimal number of clusters: {n_clusters}")

# 最適なクラスタ数でクラスタリングを実施
kmeans = KMeans(n_clusters=n_clusters, init='k-means++', random_state=42)
df['cluster'] = kmeans.fit_predict(df_std)

print(df.head())

まとめ

決定が難しい面もあるクラスター数ですが、今回の3つの方法でも結果が異なるため、複数の方法で試してみるのがよいかと思いました。