第11回 k-meansクラスタリング

科目：ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕 担当：小野原彩香 対象：3年生（後期）

---

学習目標

この回を終えると、以下のことができるようになります。

1. k-meansのアルゴリズム（初期化→割り当て→重心更新→収束）を説明できる
2. エルボー法とシルエットスコアで最適クラスタ数を決定できる
3. scikit-learnのKMeansを使ってクラスタリングを実装できる
4. クラスタ結果をGeoDataFrameに結合してコロプレスマップで可視化できる
5. レーダーチャートやヒートマップでクラスタの特徴を解釈できる
6. 発展：k-meansが必ず止まる理由と、それでも最適とは限らない理由を、自分で導ける

---

1. クラスタリングと地域類型化

1.1 地域類型化（ゾーニング）とは

都市計画や地域分析では、「似た性質を持つ地域をグループ化する」作業が重要です。これを地域類型化（ゾーニング）と呼びます。

地理空間LLMの文脈では、以下のような類型化が役立ちます。

想定されるクラスタ	特徴
都市中心型	POI密度高、クエリ数多、エラー率低い
郊外住宅型	人口密度中程度、昼間より夜間のクエリ多
農村・山間型	クエリ密度低、道路密度低、エラー率高い（地名あいまい）
観光地型	週末に利用急増、季節変動大
交通ハブ型	交通量多、多様なPOI、応答時間の変動大

このような類型化をデータ主導で行うのがk-meansクラスタリングです。

1.2 教師なし学習としてのクラスタリング

クラスタリングは教師なし学習の一種です。

• 教師あり学習：正解ラベルあり（例：このグリッドは「都市型」と事前に分類済み）
• 教師なし学習：正解ラベルなし（データ自体の構造からグループを発見する）

正解がない分、結果の解釈が重要になります。

---

2. k-meansアルゴリズム

2.1 アルゴリズムの流れ

k-meansは以下の4ステップを収束するまで繰り返します。

[Step 1] 初期化
   ランダムまたはk-means++でk個の重心を配置

[Step 2] 割り当て（Assignment）
   各データポイントを、最も近い重心のクラスタに割り当てる
   クラスタ(i) = argmin_k ||x_i - μ_k||²

[Step 3] 重心更新（Update）
   各クラスタ内のデータポイントの平均位置を新しい重心とする
   μ_k = (1/|C_k|) Σ_{x∈C_k} x

[Step 4] 収束判定
   重心の移動量が閾値以下になったら終了
   そうでなければStep 2に戻る

導出①：なぜこの手続きは必ず止まるのか（そして最適とは限らないのか）

k-meansが最小化しようとしている量は、クラスタ内誤差平方和です。

\[J = \sum_{k} \sum_{x \in C_k} \|x - \mu_k\|^2\]

アルゴリズムの2つのステップが、どちらも \(J\) を増やさないことを見ます。

Step 2（割り当て）：重心 \(\mu_k\) を固定したまま、各点を最も近い重心に付け替えます。点ごとに見れば \(\|x - \mu_k\|^2\) が最小の \(k\) を選ぶので、その点の寄与は減るか変わらないかのどちらかです。全部足した \(J\) も増えません。

Step 3（重心更新）：割り当てを固定したまま、重心を動かします。ここで効くのが「平均は二乗のズレの合計を最小にする点」という性質です。1次元で書くと、

\[\sum_{i}(x_i - c)^2 = \sum_{i}(x_i - \bar{x})^2 + n(\bar{x} - c)^2\]

と平方完成でき、右辺の第2項は \(c = \bar{x}\) のとき0で最小になります。だから重心を平均に置き直すと \(J\) は減るか変わりません。

数値で確かめる：データ \(\{1,\ 2,\ 6,\ 7,\ 9\}\)（平均5）で \(\sum(x_i-c)^2\) を計算すると、

\(c\)	3	4	5（平均）	6
\(\sum(x_i-c)^2\)	66	51	46	51

平均のときが最小です。

止まる理由：\(J\) は毎回減るか変わらないかで、しかも \(J \ge 0\) と下に有界です。さらに点をクラスタに割り当てる場合の数は有限なので、同じ割り当てに戻ってきた時点で止まります。

ただし最適とは限りません。上の議論が言っているのは「\(J\) が下がり続ける」だけで、「いちばん低いところに着く」ではありません。転がり落ちた先が、いちばん深い谷とは限らないからです。これが n_init を大きくして初期値を変えながら何度も試す理由です（第3節で数値を見ます）。

使った数学：平方完成。単調で下に有界な量は止まる、という考え方。

2.2 k-means++による初期化

通常のランダム初期化は局所最適に陥りやすいです。 k-means++は初期重心を「互いに離れた位置」に配置することで、収束の速さと品質を改善します。

scikit-learnのデフォルト設定はinit='k-means++'なので、特に変更不要です。

2.3 k-meansの限界

問題点	内容
球形クラスタ仮定	細長い形や非凸なクラスタは苦手
スケール依存	必ず事前に標準化が必要
k指定が必要	クラスタ数を事前に決めなければならない
外れ値に弱い	外れ値が重心をずらす

導出②：スケールを変えるとクラスタが変わる ―― 標準化が必要な理由

表にある「スケール依存」は、目的関数の形から直接出てきます。\(J\) の中身はユークリッド距離の二乗

\[\|x - \mu\|^2 = \sum_{j=1}^{p} (x_j - \mu_j)^2\]

で、各変数の差を対等に足しています。ある変数だけ単位が大きければ、その変数の差だけが二乗和を支配し、距離はほぼその変数だけで決まってしまいます。

この資料の特徴量でいえば、nearest_station_dist（数百〜数千 m）と building_coverage（0〜1）をそのまま足すと、建物被覆率の違いは距離の違いに埋もれて見えなくなります。

数値で確かめる：3つのかたまりを持つ2次元データ（90点）を用意し、2列目だけを100倍して同じ手順でクラスタリングすると、元のクラスタとの一致度は ARI（調整ランド指数）で \(0.4945\) でした。1.0なら完全一致、0なら偶然並みです。データの中身は何も変えていないのに、単位を変えただけで別の答えになります。

だから StandardScaler は前処理の飾りではなく、距離の意味を決める手順です。

使った数学：ユークリッド距離の定義。定数倍が二乗で効くこと。

---

3. 最適クラスタ数の選択

3.1 エルボー法

クラスタ内誤差平方和（WCSS: Within-Cluster Sum of Squares）を、kを変えながらプロットします。

WCSS = Σ_k Σ_{x∈C_k} ||x - μ_k||²

グラフが「肘（elbow）」の形に折れ曲がる点が最適なkです。

導出③：WCSSだけでは \(k\) を選べない ―― 肘を見る理由と、初期値の当たり外れ

エルボー法のグラフが右下がりになるのは偶然ではありません。\(k\) を増やすと WCSS は必ず下がる（少なくとも増えない）からです。

\(k\) 個の分割があったとして、そこにクラスタを1つ足して \(k+1\) 個にするとします。どれか1つのクラスタを2つに割れば、割った側の二乗和は導出①の平方完成の分だけ減り、他のクラスタは変わりません。極端な場合、\(k\) を点の数まで増やせば各点が自分だけのクラスタになり、WCSS はちょうど0になります。

つまり WCSS が小さいという理由だけで \(k\) を選ぶことはできません。最小化すれば \(k\) は必ず最大値に貼り付きます。そこで「下がり方が急に緩くなる場所」＝肘を探すことになります。肘の位置は式で決まる量ではなく、目で見て判断する経験則です。シルエットスコアを併用するのは、別の基準で裏を取るためです。

数値で確かめる：同じデータ・同じ \(k=3\) で、初期値だけを変えて8回実行すると、WCSS は7回が \(33.1247\)、1回だけ \(156.6950\) になりました。悪いほうは、初期値のせいで浅い谷に落ちたまま止まった結果です。導出①の「必ず止まるが最適とは限らない」が、そのまま数値に出ています。n_init を大きくするのは、この当たり外れを引き直すためです。

使った数学：導出①の平方完成。単調性の議論。

import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
import numpy as np

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

def plot_elbow(X, k_range=range(2, 12), random_state=42):
    """エルボー法でWCSSをプロットする"""
    wcss_list = []
    for k in k_range:
        km = KMeans(n_clusters=k, init='k-means++', n_init=10,
                    random_state=random_state)
        km.fit(X)
        wcss_list.append(km.inertia_)

    plt.figure(figsize=(8, 5))
    plt.plot(list(k_range), wcss_list, 'bo-', markersize=8)
    plt.xlabel('クラスタ数 k')
    plt.ylabel('WCSS（クラスタ内誤差平方和）')
    plt.title('エルボー法')
    plt.grid(alpha=0.3)
    plt.savefig('output/elbow.png', dpi=150, bbox_inches='tight')
    plt.show()
    return wcss_list

3.2 シルエットスコア

シルエットスコアはクラスタの凝集度と分離度を同時に評価します。

\[s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}\]

• \(a(i)\)：点iとその同じクラスタ内の他の点との平均距離（凝集度）
• \(b(i)\)：点iと最近傍の別クラスタとの平均距離（分離度）
• 範囲：-1〜1（1に近いほど良い）

from sklearn.metrics import silhouette_score, silhouette_samples
import matplotlib.cm as cm

def plot_silhouette(X, k_range=range(2, 12), random_state=42):
    """シルエットスコアをプロットする"""
    scores = []
    for k in k_range:
        km = KMeans(n_clusters=k, init='k-means++', n_init=10,
                    random_state=random_state)
        labels = km.fit_predict(X)
        score = silhouette_score(X, labels)
        scores.append(score)
        print(f"k={k}: シルエットスコア={score:.4f}")

    plt.figure(figsize=(8, 5))
    plt.plot(list(k_range), scores, 'ro-', markersize=8)
    plt.xlabel('クラスタ数 k')
    plt.ylabel('シルエットスコア')
    plt.title('シルエットスコアによるk選択')
    plt.grid(alpha=0.3)
    plt.savefig('output/silhouette_scores.png', dpi=150, bbox_inches='tight')
    plt.show()
    return scores

3.3 両者の組み合わせ

import pandas as pd
from sklearn.preprocessing import StandardScaler

# データ読み込みと前処理
df = pd.read_parquet('data/spatial_features.parquet')

feature_cols = [
    'query_density', 'poi_density', 'road_density',
    'nearest_station_dist', 'green_area_ratio',
    'building_coverage', 'elevation_mean', 'population_density'
]
X = df[feature_cols].dropna().values

# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

k_range = range(2, 11)

# エルボー法
wcss_list = plot_elbow(X_scaled, k_range)

# シルエットスコア
sil_scores = plot_silhouette(X_scaled, k_range)

# 総合的な判断
print("\n=== k選択サマリー ===")
best_k_sil = list(k_range)[np.argmax(sil_scores)]
print(f"シルエットスコアが最大のk: {best_k_sil}")
print("エルボー法の肘：グラフを目視で確認してください")

---

4. 空間特徴量へのk-means適用

4.1 PCA後の主成分空間でのクラスタリング

第9回で作成した主成分スコアを使います。高次元特徴量よりも安定してクラスタリングできます。

from sklearn.decomposition import PCA
from sklearn.cluster import KMeans

# PCAで次元削減（第9回の処理を再現）
pca = PCA(n_components=4, random_state=42)
X_pca = pca.fit_transform(X_scaled)

print(f"PCA後の次元数: {X_pca.shape[1]}")
print(f"累積寄与率: {pca.explained_variance_ratio_.cumsum()[-1]:.3f}")

# 最適kでk-means実行（例：k=5）
BEST_K = 5

km = KMeans(
    n_clusters=BEST_K,
    init='k-means++',
    n_init=20,          # 初期値を20回変えて最良を採用
    max_iter=300,        # 最大反復回数
    random_state=42
)
cluster_labels = km.fit_predict(X_pca)

print(f"\nクラスタリング完了")
print(f"クラスタ別サンプル数:")
unique, counts = np.unique(cluster_labels, return_counts=True)
for c, n in zip(unique, counts):
    print(f"  クラスタ{c}: {n}件 ({n/len(cluster_labels)*100:.1f}%)")

4.2 クラスタ結果のDataFrame統合

# 元のDataFrameにクラスタラベルを付与
df_valid = df[feature_cols].dropna().copy()
df_valid['cluster'] = cluster_labels
df_valid['grid_id'] = df['grid_id'].loc[df_valid.index]
df_valid['longitude'] = df['longitude'].loc[df_valid.index]
df_valid['latitude'] = df['latitude'].loc[df_valid.index]

# クラスタごとの特徴量平均
cluster_profile = (
    df_valid.groupby('cluster')[feature_cols]
    .mean()
    .round(3)
)
print("\nクラスタプロファイル（特徴量平均）:")
print(cluster_profile)

---

5. 結果の地図表示と解釈

5.1 GeoDataFrameへの変換とコロプレスマップ

import geopandas as gpd
from shapely.geometry import Point
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

# GeoDataFrame作成
geometry = [Point(lon, lat)
            for lon, lat in zip(df_valid['longitude'], df_valid['latitude'])]
gdf = gpd.GeoDataFrame(df_valid, geometry=geometry, crs='EPSG:4326')

# クラスタカラーの定義
cluster_colors = {
    0: '#2196F3',  # 青：都市中心型
    1: '#4CAF50',  # 緑：郊外住宅型
    2: '#FF9800',  # オレンジ：農村・山間型
    3: '#9C27B0',  # 紫：観光地型
    4: '#F44336',  # 赤：交通ハブ型
}
cluster_names = {
    0: '都市中心型',
    1: '郊外住宅型',
    2: '農村・山間型',
    3: '観光地型',
    4: '交通ハブ型',
}

# 地図プロット
fig, ax = plt.subplots(figsize=(12, 10))

for cluster_id in range(BEST_K):
    mask = gdf['cluster'] == cluster_id
    gdf[mask].plot(
        ax=ax,
        color=cluster_colors[cluster_id],
        markersize=30,
        alpha=0.7,
        label=f"クラスタ{cluster_id}:{cluster_names[cluster_id]}"
    )

ax.set_title('空間特徴量クラスタリング結果', fontsize=14)
ax.set_xlabel('経度')
ax.set_ylabel('緯度')
ax.legend(loc='lower right', fontsize=10)
plt.tight_layout()
plt.savefig('output/cluster_map.png', dpi=150, bbox_inches='tight')
plt.show()

5.2 レーダーチャートによるクラスタプロファイル可視化

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch

def plot_radar_chart(cluster_profile: pd.DataFrame,
                     feature_names: list,
                     cluster_names: dict,
                     cluster_colors: dict):
    """
    クラスタプロファイルをレーダーチャートで可視化する
    """
    # 0〜1に正規化
    profile_norm = (cluster_profile - cluster_profile.min()) / \
                   (cluster_profile.max() - cluster_profile.min())

    n_features = len(feature_names)
    angles = np.linspace(0, 2 * np.pi, n_features, endpoint=False).tolist()
    angles += angles[:1]  # 閉じるために先頭を末尾に追加

    fig, axes = plt.subplots(2, 3, figsize=(15, 10),
                              subplot_kw=dict(polar=True))
    axes = axes.flatten()

    for i, (cluster_id, row) in enumerate(profile_norm.iterrows()):
        ax = axes[i]
        values = row[feature_names].tolist()
        values += values[:1]

        ax.plot(angles, values,
                color=cluster_colors[cluster_id], lw=2)
        ax.fill(angles, values,
                color=cluster_colors[cluster_id], alpha=0.25)
        ax.set_xticks(angles[:-1])
        ax.set_xticklabels(feature_names, size=8)
        ax.set_ylim(0, 1)
        ax.set_title(f"クラスタ{cluster_id}\n{cluster_names[cluster_id]}",
                     size=11, color=cluster_colors[cluster_id], pad=15)
        ax.grid(alpha=0.3)

    # 余った軸を非表示
    for j in range(i+1, len(axes)):
        axes[j].set_visible(False)

    plt.suptitle('クラスタ特徴プロファイル（正規化）', fontsize=14, y=1.02)
    plt.tight_layout()
    plt.savefig('output/radar_chart.png', dpi=150, bbox_inches='tight')
    plt.show()

plot_radar_chart(cluster_profile, feature_cols, cluster_names, cluster_colors)

5.3 特徴量プロファイルヒートマップ

import seaborn as sns

# 標準化後の平均値でヒートマップ
df_scaled_with_cluster = pd.DataFrame(X_scaled, columns=feature_cols)
df_scaled_with_cluster['cluster'] = cluster_labels

heatmap_data = df_scaled_with_cluster.groupby('cluster')[feature_cols].mean()
heatmap_data.index = [f"C{c}:{cluster_names[c]}" for c in heatmap_data.index]

plt.figure(figsize=(12, 6))
sns.heatmap(
    heatmap_data,
    annot=True, fmt='.2f',
    cmap='RdBu_r', center=0,
    linewidths=0.5,
    cbar_kws={'label': '標準化スコア（平均0、標準偏差1）'}
)
plt.title('クラスタ×特徴量 プロファイルヒートマップ')
plt.xlabel('特徴量')
plt.ylabel('クラスタ')
plt.tight_layout()
plt.savefig('output/cluster_feature_heatmap.png', dpi=150, bbox_inches='tight')
plt.show()

---

6. ハンズオン手順

6.1 一連のパイプラインを実行する

# 出力ディレクトリを作成
mkdir -p output

# Jupyterを起動
jupyter notebook 11_kmeans.ipynb

6.2 パイプライン全体の流れ

# ========== パイプライン完全版 ==========

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import geopandas as gpd
from shapely.geometry import Point
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

# 1. データ読み込み
df = pd.read_parquet('data/spatial_features.parquet')
feature_cols = ['query_density', 'poi_density', 'road_density',
                'nearest_station_dist', 'green_area_ratio',
                'building_coverage', 'elevation_mean', 'population_density']

X_raw = df[feature_cols].dropna()
idx = X_raw.index

# 2. 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

# 3. PCAで次元削減
pca = PCA(n_components=4, random_state=42)
X_pca = pca.fit_transform(X_scaled)

# 4. 最適k決定（エルボー＋シルエット）
k_range = range(2, 10)
sil_scores = []
wcss_list = []
for k in k_range:
    km = KMeans(n_clusters=k, n_init=10, random_state=42)
    labels = km.fit_predict(X_pca)
    wcss_list.append(km.inertia_)
    sil_scores.append(silhouette_score(X_pca, labels))

BEST_K = list(k_range)[np.argmax(sil_scores)]
print(f"最適k（シルエット最大）: {BEST_K}")

# 5. k-meansクラスタリング
km_final = KMeans(n_clusters=BEST_K, n_init=20, random_state=42)
cluster_labels = km_final.fit_predict(X_pca)

# 6. 結果をDataFrameに格納
df_result = df.loc[idx].copy()
df_result['cluster'] = cluster_labels

# 7. GeoDataFrameに変換
gdf = gpd.GeoDataFrame(
    df_result,
    geometry=[Point(lon, lat) for lon, lat in
              zip(df_result['longitude'], df_result['latitude'])],
    crs='EPSG:4326'
)

# 8. 結果を保存
gdf.to_file('output/cluster_result.geojson', driver='GeoJSON')
print("GeoJSON保存完了: output/cluster_result.geojson")

---

7. 演習課題

課題1（基本）

data/spatial_features.parquetを読み込み、以下を実行してください。

1. StandardScalerで標準化する
2. エルボー法（k=2〜10）のグラフを描く
3. シルエットスコアが最大になるkを求める

課題2（応用）

1. 最適kでk-meansを実行し、クラスタラベルをGeoDataFrameに格納する
2. コロプレスマップでクラスタ分布を可視化する
3. 各クラスタの特徴量プロファイルをヒートマップで表示し、クラスタに名前をつける（例：「都市中心型」「農村型」など）

課題3（発展）

1. k-meansの結果と「エラー発生率」の関係を箱ひげ図で可視化する
2. エラー率が最も高いクラスタの地理的分布を分析し、その原因について100字以上で考察する
3. k-meansに代えてDBSCANを試し、結果の違いを比較する

提出形式：Jupyterノートブック、ファイル名：11_kmeans_[学籍番号].ipynb

---

まとめ

項目	ポイント
k-meansの原理	重心割り当てと更新を繰り返し、クラスタ内分散を最小化
k-means++	初期重心を分散させて局所最適を回避
k選択	エルボー法（WCSS）とシルエットスコアを併用する
標準化の必要性	異なるスケールの特徴量をそのまま使うとスケール支配が起きる
PCAとの組み合わせ	主成分空間でクラスタリングすると次元の呪いを回避できる
解釈の重要性	レーダーチャートやヒートマップで「このクラスタとは何か」を説明する
次回への接続	クラスタとエラー発生率の関係を第12回のPLS回帰で定量的にモデル化する

次回予告： 第12回では、空間特徴量とエラー発生率の関係をPLS回帰でモデル化します。PCAとの違いも理解しましょう。

---

作成：小野原彩香 / ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕