# K-means聚类算法

## 前言

聚类分析是机器学习中无监督学习的重要分支，而K-means算法作为最经典和广泛使用的聚类算法之一，在数据挖掘、图像分割、客户分群等领域有着重要应用。本文将深入探讨K-means算法的原理、实现方法、评估指标以及实际应用案例。

## 1. K-means算法概述

### 1.1 什么是K-means聚类

K-means聚类是一种基于距离的聚类算法，其目标是将n个数据点划分为k个簇，使得每个数据点都属于距离其最近的簇中心所在的簇。算法通过迭代优化来最小化簇内平方和（Within-Cluster Sum of Squares, WCSS）。

### 1.2 算法核心思想

K-means算法的核心思想是：
1. **初始化**：随机选择k个初始聚类中心
2. **分配**：将每个数据点分配到最近的聚类中心
3. **更新**：重新计算每个簇的中心点
4. **迭代**：重复步骤2和3，直到收敛

### 1.3 算法优缺点

**优点：**
- 算法简单，易于理解和实现
- 计算复杂度相对较低
- 适用于球形簇的聚类
- 可扩展性好

**缺点：**
- 需要预先指定聚类数量k
- 对初始聚类中心敏感
- 对异常值敏感
- 假设簇是球形的，对非球形簇效果不佳

## 2. K-means算法实现

### 2.1 基础实现

```python
import os
import matplotlib
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import calinski_harabasz_score

# 设置环境变量
os.environ['OMP_NUM_THREADS'] = '4'
matplotlib.use('TkAgg')

def basic_kmeans_demo():
    """K-means基础演示"""
    # 1. 生成模拟数据
    x, y = make_blobs(
        n_samples=1000,
        n_features=2,
        centers=[[-1, -1], [0, 0], [1, 1], [2, 2]],
        cluster_std=[0.4, 0.2, 0.2, 0.2],
        random_state=23
    )
    
    # 2. 可视化原始数据
    plt.figure(figsize=(12, 5))
    plt.subplot(1, 2, 1)
    plt.scatter(x[:, 0], x[:, 1], c='blue', alpha=0.6)
    plt.title('原始数据分布')
    plt.xlabel('特征1')
    plt.ylabel('特征2')
    
    # 3. K-means聚类
    kmeans = KMeans(n_clusters=4, random_state=23)
    y_pred = kmeans.fit_predict(x)
    
    # 4. 可视化聚类结果
    plt.subplot(1, 2, 2)
    plt.scatter(x[:, 0], x[:, 1], c=y_pred, cmap='viridis', alpha=0.6)
    plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], 
                c='red', marker='x', s=200, linewidths=3)
    plt.title('K-means聚类结果')
    plt.xlabel('特征1')
    plt.ylabel('特征2')
    
    plt.tight_layout()
    plt.show()
    
    # 5. 评估聚类效果
    ch_score = calinski_harabasz_score(x, y_pred)
    print(f"Calinski-Harabasz指数: {ch_score:.2f}")
    
    return kmeans, y_pred
```

### 2.2 最优K值选择

选择合适的K值是K-means算法的关键问题。常用的方法包括肘部法则、轮廓系数和Calinski-Harabasz指数。

```python
from sklearn.metrics import silhouette_score, calinski_harabasz_score

def find_optimal_k():
    """寻找最优K值"""
    # 生成数据
    x, y = make_blobs(
        n_samples=1000,
        n_features=2,
        centers=[[-1, -1], [0, 0], [1, 1], [2, 2]],
        cluster_std=[0.4, 0.2, 0.2, 0.2],
        random_state=23
    )
    
    # 存储不同K值的评估指标
    sse_list = []  # 簇内平方和
    sc_list = []   # 轮廓系数
    ch_list = []   # Calinski-Harabasz指数
    
    k_range = range(2, 20)
    
    for k in k_range:
        kmeans = KMeans(n_clusters=k, max_iter=100, random_state=23)
        kmeans.fit(x)
        y_pred = kmeans.predict(x)
        
        # 计算评估指标
        sse_list.append(kmeans.inertia_)  # 簇内平方和
        sc_list.append(silhouette_score(x, y_pred))  # 轮廓系数
        ch_list.append(calinski_harabasz_score(x, y_pred))  # CH指数
    
    # 可视化评估指标
    plt.figure(figsize=(15, 5))
    
    # 肘部法则
    plt.subplot(1, 3, 1)
    plt.plot(k_range, sse_list, 'bo-')
    plt.xlabel('K值')
    plt.ylabel('簇内平方和(SSE)')
    plt.title('肘部法则')
    plt.grid(True)
    
    # 轮廓系数
    plt.subplot(1, 3, 2)
    plt.plot(k_range, sc_list, 'ro-')
    plt.xlabel('K值')
    plt.ylabel('轮廓系数')
    plt.title('轮廓系数法')
    plt.grid(True)
    
    # Calinski-Harabasz指数
    plt.subplot(1, 3, 3)
    plt.plot(k_range, ch_list, 'go-')
    plt.xlabel('K值')
    plt.ylabel('CH指数')
    plt.title('Calinski-Harabasz指数')
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    # 输出最优K值
    optimal_k_sc = k_range[sc_list.index(max(sc_list))]
    optimal_k_ch = k_range[ch_list.index(max(ch_list))]
    
    print(f"轮廓系数法推荐K值: {optimal_k_sc}")
    print(f"CH指数法推荐K值: {optimal_k_ch}")
    
    return optimal_k_sc, optimal_k_ch
```

## 3. 聚类评估方法详解

### 3.1 内部评估指标

#### 3.1.1 簇内平方和（SSE/Inertia）
```python
def calculate_sse(x, y_pred, centers):
    """计算簇内平方和"""
    sse = 0
    for i in range(len(centers)):
        cluster_points = x[y_pred == i]
        if len(cluster_points) > 0:
            sse += np.sum((cluster_points - centers[i]) ** 2)
    return sse
```

#### 3.1.2 轮廓系数（Silhouette Score）
轮廓系数衡量一个样本与其自身簇的相似度相对于其他簇的相似度。

```python
def silhouette_analysis():
    """轮廓系数分析"""
    x, y = make_blobs(n_samples=1000, n_features=2, centers=4, random_state=23)
    
    # 不同K值的轮廓系数
    k_values = [2, 3, 4, 5, 6]
    silhouette_scores = []
    
    for k in k_values:
        kmeans = KMeans(n_clusters=k, random_state=23)
        y_pred = kmeans.fit_predict(x)
        score = silhouette_score(x, y_pred)
        silhouette_scores.append(score)
        print(f"K={k}, 轮廓系数: {score:.3f}")
    
    return silhouette_scores
```

#### 3.1.3 Calinski-Harabasz指数
CH指数基于簇间和簇内的方差比，值越大表示聚类效果越好。

```python
def ch_analysis():
    """Calinski-Harabasz指数分析"""
    x, y = make_blobs(n_samples=1000, n_features=2, centers=4, random_state=23)
    
    k_values = range(2, 10)
    ch_scores = []
    
    for k in k_values:
        kmeans = KMeans(n_clusters=k, random_state=23)
        y_pred = kmeans.fit_predict(x)
        score = calinski_harabasz_score(x, y_pred)
        ch_scores.append(score)
        print(f"K={k}, CH指数: {score:.2f}")
    
    return ch_scores
```

### 3.2 外部评估指标

当有真实标签时，可以使用外部评估指标：

```python
from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score

def external_evaluation():
    """外部评估指标"""
    x, y_true = make_blobs(n_samples=1000, n_features=2, centers=4, random_state=23)
    
    kmeans = KMeans(n_clusters=4, random_state=23)
    y_pred = kmeans.fit_predict(x)
    
    # 调整兰德指数
    ari = adjusted_rand_score(y_true, y_pred)
    print(f"调整兰德指数: {ari:.3f}")
    
    # 标准化互信息
    nmi = normalized_mutual_info_score(y_true, y_pred)
    print(f"标准化互信息: {nmi:.3f}")
    
    return ari, nmi
```

## 4. 实际应用案例：客户分群

### 4.1 客户分群业务背景

在电商和零售行业中，客户分群是重要的营销策略。通过分析客户的消费行为、收入水平等特征，可以将客户分为不同的群体，从而制定个性化的营销策略。

### 4.2 数据准备与预处理

```python
import pandas as pd

def customer_segmentation():
    """客户分群案例"""
    # 1. 加载客户数据
    df = pd.read_csv("customers.csv")
    print("数据基本信息:")
    print(df.info())
    print("\n数据预览:")
    print(df.head())
    
    # 2. 特征选择
    # 假设使用年收入和消费分数作为聚类特征
    x = df.iloc[:, 3:5]  # 选择第4和第5列作为特征
    print(f"\n特征数据形状: {x.shape}")
    print("特征数据预览:")
    print(x.head())
    
    return x, df
```

### 4.3 寻找最优K值

```python
def find_optimal_k_customers(x):
    """为客户数据寻找最优K值"""
    sse_list = []
    sc_list = []
    
    k_range = range(2, 20)
    
    for k in k_range:
        kmeans = KMeans(n_clusters=k, max_iter=100, random_state=23)
        kmeans.fit(x)
        y_pred = kmeans.predict(x)
        
        sse_list.append(kmeans.inertia_)
        sc_list.append(silhouette_score(x, y_pred))
    
    # 可视化
    plt.figure(figsize=(15, 5))
    
    plt.subplot(1, 2, 1)
    plt.plot(k_range, sse_list, 'bo-')
    plt.xlabel('K值')
    plt.ylabel('簇内平方和(SSE)')
    plt.title('肘部法则 - 客户分群')
    plt.grid(True)
    
    plt.subplot(1, 2, 2)
    plt.plot(k_range, sc_list, 'ro-')
    plt.xlabel('K值')
    plt.ylabel('轮廓系数')
    plt.title('轮廓系数法 - 客户分群')
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    # 选择最优K值
    optimal_k = k_range[sc_list.index(max(sc_list))]
    print(f"推荐K值: {optimal_k}")
    
    return optimal_k
```

### 4.4 客户分群实现

```python
def perform_customer_clustering(x, optimal_k):
    """执行客户分群"""
    # 1. 训练K-means模型
    kmeans = KMeans(n_clusters=optimal_k, max_iter=100, random_state=23)
    kmeans.fit(x)
    
    # 2. 预测客户类别
    y_pred = kmeans.predict(x)
    
    # 3. 可视化分群结果
    plt.figure(figsize=(10, 8))
    colors = ['red', 'blue', 'green', 'orange', 'purple', 'brown', 'pink', 'gray']
    
    for i in range(optimal_k):
        cluster_data = x.values[y_pred == i]
        plt.scatter(cluster_data[:, 0], cluster_data[:, 1], 
                   c=colors[i % len(colors)], label=f'客户群体 {i+1}', alpha=0.7)
    
    # 绘制聚类中心
    plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], 
               c='black', marker='x', s=200, linewidths=3, label='聚类中心')
    
    plt.xlabel('年收入')
    plt.ylabel('消费分数')
    plt.title(f'客户分群结果 (K={optimal_k})')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    # 4. 分析各客户群体特征
    print("\n各客户群体特征分析:")
    for i in range(optimal_k):
        cluster_data = x.values[y_pred == i]
        print(f"\n客户群体 {i+1}:")
        print(f"  客户数量: {len(cluster_data)}")
        print(f"  平均年收入: {cluster_data[:, 0].mean():.2f}")
        print(f"  平均消费分数: {cluster_data[:, 1].mean():.2f}")
    
    return kmeans, y_pred
```

### 4.5 客户群体分析

```python
def analyze_customer_segments(df, y_pred, optimal_k):
    """分析客户群体特征"""
    # 添加聚类标签到原始数据
    df['Cluster'] = y_pred
    
    print("客户群体详细分析:")
    print("=" * 50)
    
    for i in range(optimal_k):
        cluster_data = df[df['Cluster'] == i]
        print(f"\n客户群体 {i+1} (共{len(cluster_data)}人):")
        print("-" * 30)
        
        # 基本统计信息
        if 'Age' in df.columns:
            print(f"平均年龄: {cluster_data['Age'].mean():.1f}岁")
        if 'Annual Income (k$)' in df.columns:
            print(f"平均年收入: {cluster_data['Annual Income (k$)'].mean():.1f}k$")
        if 'Spending Score (1-100)' in df.columns:
            print(f"平均消费分数: {cluster_data['Spending Score (1-100)'].mean():.1f}")
        
        # 性别分布
        if 'Gender' in df.columns:
            gender_dist = cluster_data['Gender'].value_counts()
            print(f"性别分布: {dict(gender_dist)}")
    
    return df
```

## 5. 高级应用：XGBoost与聚类结合

### 5.1 红酒品质多分类案例

```python
import xgboost as xgb
import joblib
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import classification_report, accuracy_score
from sklearn.model_selection import StratifiedKFold
from sklearn.utils import class_weight
from collections import Counter

def wine_quality_classification():
    """红酒品质多分类案例"""
    
    def data_split():
        """数据分割"""
        df = pd.read_csv("红酒品质分类.csv")
        x = df.iloc[:, :-1]
        y = df.iloc[:, -1] - 3  # 将标签从3-8转换为0-5
        
        print(f"标签分布: {Counter(y)}")
        
        # 分层采样分割数据
        x_train, x_test, y_train, y_test = train_test_split(
            x, y, test_size=0.2, random_state=23, stratify=y
        )
        
        # 保存数据
        pd.concat([x_train, y_train], axis=1).to_csv(
            "红酒品质分类_train.csv", index=False
        )
        pd.concat([x_test, y_test], axis=1).to_csv(
            "红酒品质分类_test.csv", index=False
        )
        
        return x_train, x_test, y_train, y_test
    
    def train_model():
        """训练模型"""
        # 加载数据
        train_data = pd.read_csv("红酒品质分类_train.csv")
        test_data = pd.read_csv("红酒品质分类_test.csv")
        
        x_train = train_data.iloc[:, :-1]
        y_train = train_data.iloc[:, -1]
        x_test = test_data.iloc[:, :-1]
        y_test = test_data.iloc[:, -1]
        
        # 创建XGBoost模型
        model = xgb.XGBClassifier(
            max_depth=5,
            n_estimators=100,
            learning_rate=0.1,
            random_state=24,
            objective='multi:softmax'
        )
        
        # 处理样本不平衡
        sample_weights = class_weight.compute_sample_weight("balanced", y_train)
        
        # 训练模型
        model.fit(x_train, y_train, sample_weight=sample_weights)
        
        # 评估模型
        y_pred = model.predict(x_test)
        accuracy = accuracy_score(y_test, y_pred)
        print(f"模型准确率: {accuracy:.4f}")
        
        # 保存模型
        joblib.dump(model, "红酒品质分类.pkl")
        print("模型保存成功")
        
        return model, accuracy
    
    def hyperparameter_tuning():
        """超参数调优"""
        train_data = pd.read_csv("红酒品质分类_train.csv")
        test_data = pd.read_csv("红酒品质分类_test.csv")
        
        x_train = train_data.iloc[:, :-1]
        y_train = train_data.iloc[:, -1]
        x_test = test_data.iloc[:, :-1]
        y_test = test_data.iloc[:, -1]
        
        # 参数网格
        param_grid = {
            'max_depth': [3, 5, 8, 10, 12],
            'n_estimators': [100, 130, 150, 180, 230],
            'learning_rate': [0.1, 0.3, 0.5, 1.0, 1.2]
        }
        
        # 分层K折交叉验证
        skf = StratifiedKFold(n_splits=2, shuffle=True, random_state=23)
        
        # 网格搜索
        model = xgb.XGBClassifier()
        grid_search = GridSearchCV(model, param_grid, cv=skf)
        grid_search.fit(x_train, y_train)
        
        # 评估最佳模型
        y_pred = grid_search.predict(x_test)
        accuracy = accuracy_score(y_test, y_pred)
        
        print(f"网格搜索最佳参数: {grid_search.best_params_}")
        print(f"网格搜索最佳分数: {grid_search.best_score_:.4f}")
        print(f"测试集准确率: {accuracy:.4f}")
        
        return grid_search
    
    # 执行完整流程
    print("1. 数据分割")
    data_split()
    
    print("\n2. 模型训练")
    train_model()
    
    print("\n3. 超参数调优")
    hyperparameter_tuning()
```

## 6. K-means算法优化策略

### 6.1 初始化优化

```python
def kmeans_plus_plus_init(x, k):
    """K-means++初始化"""
    centers = []
    centers.append(x[np.random.randint(x.shape[0])])
    
    for _ in range(k - 1):
        distances = []
        for point in x:
            min_dist = min([np.linalg.norm(point - center) for center in centers])
            distances.append(min_dist ** 2)
        
        probabilities = distances / np.sum(distances)
        cumulative_probabilities = np.cumsum(probabilities)
        r = np.random.random()
        
        for i, prob in enumerate(cumulative_probabilities):
            if r < prob:
                centers.append(x[i])
                break
    
    return np.array(centers)
```

### 6.2 处理不同形状的簇

```python
from sklearn.cluster import DBSCAN, AgglomerativeClustering

def alternative_clustering_methods(x):
    """替代聚类方法"""
    
    # DBSCAN - 基于密度的聚类
    dbscan = DBSCAN(eps=0.5, min_samples=5)
    y_dbscan = dbscan.fit_predict(x)
    
    # 层次聚类
    hierarchical = AgglomerativeClustering(n_clusters=4)
    y_hierarchical = hierarchical.fit_predict(x)
    
    # 可视化比较
    plt.figure(figsize=(15, 5))
    
    plt.subplot(1, 3, 1)
    plt.scatter(x[:, 0], x[:, 1], c=y_dbscan, cmap='viridis')
    plt.title('DBSCAN聚类')
    
    plt.subplot(1, 3, 2)
    plt.scatter(x[:, 0], x[:, 1], c=y_hierarchical, cmap='viridis')
    plt.title('层次聚类')
    
    plt.subplot(1, 3, 3)
    kmeans = KMeans(n_clusters=4, random_state=23)
    y_kmeans = kmeans.fit_predict(x)
    plt.scatter(x[:, 0], x[:, 1], c=y_kmeans, cmap='viridis')
    plt.title('K-means聚类')
    
    plt.tight_layout()
    plt.show()
    
    return y_dbscan, y_hierarchical, y_kmeans
```

## 7. 实际应用建议

### 7.1 数据预处理

1. **标准化处理**：确保不同特征的量纲一致
2. **异常值处理**：K-means对异常值敏感，需要预处理
3. **特征选择**：选择与业务目标相关的特征

### 7.2 参数调优

1. **K值选择**：结合业务需求和统计指标
2. **初始化策略**：使用K-means++减少随机性
3. **迭代次数**：设置合理的最大迭代次数

### 7.3 结果解释

1. **业务含义**：为每个簇赋予业务含义
2. **特征分析**：分析各簇的特征分布
3. **应用策略**：制定针对性的应用策略

## 8. 总结

K-means聚类算法作为最经典的聚类算法之一，具有以下特点：

**优势：**
- 算法简单，易于理解和实现
- 计算效率高，适合大规模数据
- 在球形簇上表现良好
- 可解释性强

**局限性：**
- 需要预先指定聚类数量
- 对初始值敏感
- 假设簇是球形的
- 对异常值敏感

**应用场景：**
- 客户分群和用户画像
- 图像分割
- 数据压缩
- 异常检测

在实际应用中，需要根据数据特点和业务需求选择合适的聚类算法，并通过多种评估指标来验证聚类效果。同时，结合业务知识对聚类结果进行解释和应用，才能真正发挥聚类分析的价值。

## 9. 扩展阅读

- [Scikit-learn聚类算法官方文档](https://scikit-learn.org/stable/modules/clustering.html)
- 《数据挖掘：概念与技术》- Jiawei Han
- 《机器学习》- 周志华

---

*本文更新日期：2025年9月8日

*本文基于实际项目代码编写，所有代码均经过测试验证。如有问题，欢迎交流讨论。*