核心概念

  1. 高斯混合模型：

    • ：混合系数（第 k 个高斯成分的权重）
    • ：均值向量
    • ：协方差矩阵

  2. EM 算法：

    • E 步骤：计算后验概率（责任度）
    • M 步骤：更新模型参数

解题思路

  1. 数据读取：获取数据点数量 和二维坐标
  2. 参数初始化：
    • 混合系数 均等初始化
    • 均值 随机选择数据点
    • 协方差 初始化为单位矩阵
  3. EM 迭代：
    • E 步：计算责任度
    • M 步：更新
  4. 标签分配：根据最大后验概率确定聚类标签
  5. 结果输出：输出每个数据点的聚类标签

代码解析

import sys
import numpy as np
from scipy.stats import multivariate_normal

def load_data():
    """读取输入数据"""
    N = int(sys.stdin.readline().strip())  # 数据点数量
    data_points = []
    for _ in range(N):
        x, y = map(float, sys.stdin.readline().split())
        data_points.append([x, y])
    K, T = map(int, sys.stdin.readline().split())  # 聚类数和迭代次数
    return np.array(data_points), K, T

def initialize_parameters(data, K):
    """初始化GMM参数"""
    N, D = data.shape
    pi = np.ones(K) / K  # 均等混合系数
    indices = np.random.choice(N, K, replace=False)  # 随机选择初始均值
    mu = data[indices]
    sigma = np.array([np.eye(D) for _ in range(K)])  # 单位协方差矩阵
    return pi, mu, sigma

def e_step(data, pi, mu, sigma, K):
    """E步骤：计算责任度"""
    N = data.shape[0]
    gamma = np.zeros((N, K))
    
    # 计算每个高斯成分的概率密度
    for k in range(K):
        rv = multivariate_normal(mean=mu[k], cov=sigma[k])
        gamma[:, k] = pi[k] * rv.pdf(data)
    
    # 归一化责任度（加1e-6避免除零）
    gamma_sum = np.sum(gamma, axis=1, keepdims=True)
    gamma = gamma / gamma_sum + 1e-6
    return gamma

def m_step(data, gamma, K):
    """M步骤：更新参数"""
    N, D = data.shape
    N_k = np.sum(gamma, axis=0)  # 有效样本数
    
    # 更新混合系数
    pi = N_k / N
    
    # 更新均值
    mu = np.zeros((K, D))
    for k in range(K):
        mu[k] = np.sum(gamma[:, k][:, np.newaxis] * data, axis=0) / N_k[k]
    
    # 更新协方差
    sigma = np.zeros((K, D, D))
    for k in range(K):
        diff = data - mu[k]  # 中心化数据
        weighted_diff = gamma[:, k][:, np.newaxis] * diff
        sigma[k] = np.dot(weighted_diff.T, diff) / N_k[k]
    
    return pi, mu, sigma

def assign_labels(data, pi, mu, sigma, K):
    """分配聚类标签"""
    N = data.shape[0]
    gamma = np.zeros((N, K))
    for k in range(K):
        rv = multivariate_normal(mean=mu[k], cov=sigma[k])
        gamma[:, k] = pi[k] * rv.pdf(data)
    return np.argmax(gamma, axis=1)  # 选择最大后验概率

if __name__ == "__main__":
    np.random.seed(0)  # 固定随机种子
    data, K, T = load_data()
    
    # 初始化参数
    pi, mu, sigma = initialize_parameters(data, K)
    
    # EM迭代
    for _ in range(T):
        gamma = e_step(data, pi, mu, sigma, K)  # E步
        pi, mu, sigma = m_step(data, gamma, K)  # M步
    
    # 分配并输出标签
    labels = assign_labels(data, pi, mu, sigma, K)
    for label in labels:
        print(label)

关键函数说明

  1. load_data()：

    • 读取数据点数量
    • 读取 个二维数据点
    • 读取聚类数 和迭代次数

  2. initialize_parameters(data, K)：

    • ：初始化为
    • ：随机选择 个数据点
    • ：初始化为 单位矩阵

  3. e_step(data, pi, mu, sigma, K)：

    • 使用多元高斯分布计算概率密度
    • 计算责任度 并归一化
    • 添加 避免数值下溢

  4. m_step(data, gamma, K)：

    • 计算有效样本数
    • 更新混合系数
    • 更新均值 （责任度加权平均）
    • 更新协方差 （责任度加权协方差）

  5. assign_labels(data, pi, mu, sigma, K)：

    • 计算最终责任度
    • 选择最大后验概率对应的聚类标签

总结

本题实现了基于 EM 算法的高斯混合模型聚类：

  1. 通过 E 步计算责任度（数据点属于各聚类的概率）
  2. 通过 M 步更新模型参数
  3. 迭代优化直至收敛
  4. 根据最大后验概率分配最终标签

该算法能有效识别客户群体，为个性化营销提供数据支持，帮助电商平台优化产品策略。