简介本资源是一份面向数据科学初学者与机器学习实践者的KMeans聚类算法深度实现源码包聚焦无监督学习核心原理落地适用于课程设计、算法复现与聚类任务实战。压缩包共204个文件含141个CSV格式实测数据集如EEMD、CEEMDAN系列多维时序特征数据、43张PNG聚类过程与结果可视化图、16个模块化Python脚本覆盖数据预处理、KMeans核心迭代、质心初始化、轮廓系数评估等关键环节以及辅助性JPG图示与gitignore配置文件整体大小45.42MB。已有629人学习下载体现较强的教学参考价值与工程复用潜力。读者可直接运行代码复现完整聚类流程深入理解K值选择、初始质心策略如KMeans对收敛性的影响并借助内置可视化脚本直观对比不同参数下的簇划分效果是掌握算法底层逻辑与提升动手能力的优质实践材料。1. 这不是调包是亲手把KMeans的每一步“拧”进内存里你手头有一堆EEMD系列CSV文件——data_EEMDV2.csv、data_CEEMDAN8.csv、k-data-EEMD1.csv……共141个全是时频域分解后的特征数据还有43张PNG图记录着不同K值下簇内误差平方和SSE曲线、轮廓系数热力图、聚类结果散点图。这不是一个“from sklearn.cluster import KMeans”就能打发的玩具项目。它是一套完整的手写KMeans实现从读取多维时间序列特征、标准化、初始化质心、欧氏距离计算、标签分配、质心重算到收敛判断、评估指标轮廓系数、Calinski-Harabasz指数、SSE、可视化对比全部用原生NumPy和少量Matplotlib完成不依赖scikit-learn的fit/predict接口。它专为理解算法边界而生——当你发现sklearn在处理高维稀疏EEMD特征时收敛震荡或KMeans在非球形分布上失效这套代码能让你立刻定位是距离度量偏差、质心更新逻辑错误还是初始点采样策略缺陷。适合正在啃《Pattern Recognition and Machine Learning》第9章的算法工程师、需要复现论文实验的研究生以及被面试官问“手写KMeans怎么防死循环”的中级Python开发者。2. 从CSV数据加载到质心迭代核心模块逐行拆解2.1 数据加载与预处理为什么必须重写pandas.read_csv项目中141个CSV文件并非标准表格结构。以data_EEMDV2_7.csv为例其首行为[‘t’, ‘IMF1’, ‘IMF2’, ‘IMF3’, ‘residue’]但实际数据包含大量NaN填充的截断序列且各列量纲差异极大IMF1幅值在1e-3量级residue可达1e2。直接pd.read_csv()会将NaN转为float64导致后续距离计算溢出而fillna(0)又会污染物理意义。源码中data_loader.py采用分块解析策略import numpy as np import csv def load_eemd_csv(filepath, skip_header1, fill_naninterpolate): 专为EEMD分解数据设计的加载器 with open(filepath, r) as f: reader csv.reader(f) # 跳过表头读取剩余行 rows list(reader)[skip_header:] # 按列提取跳过第一列时间戳t只保留IMF分量 data_cols [] for col_idx in range(1, len(rows[0])): # 从第2列开始 col_data [] for row in rows: if len(row) col_idx and row[col_idx].strip(): try: col_data.append(float(row[col_idx])) except ValueError: col_data.append(np.nan) else: col_data.append(np.nan) # 处理NaN插值法比均值填充更符合信号连续性假设 if fill_nan interpolate: col_arr np.array(col_data) valid_mask ~np.isnan(col_arr) if np.any(valid_mask): # 线性插值避免端点外推 col_arr np.interp( np.arange(len(col_arr)), np.where(valid_mask)[0], col_arr[valid_mask] ) else: col_arr[:] 0 # 全空则置0 data_cols.append(col_arr) return np.column_stack(data_cols) # (n_samples, n_features) # 示例加载data_EEMDV2_6.csv并标准化 X_raw load_eemd_csv(data_EEMDV2_6.csv) X_scaled (X_raw - np.mean(X_raw, axis0)) / (np.std(X_raw, axis0) 1e-8) # 防除零提示np.interp替代pandas.interpolate是关键。EEMD分量具有强时序相关性线性插值保留局部趋势而pandas默认的spline插值在端点易震荡会引入虚假高频分量直接影响后续聚类稳定性。2.2 KMeans核心类5个方法撑起整个算法骨架kmeans_core.py定义了CustomKMeans类其设计刻意避开面向对象的过度封装每个方法对应算法一个原子步骤class CustomKMeans: def __init__(self, n_clusters3, max_iters300, init_methodk-means, random_state42): self.n_clusters n_clusters self.max_iters max_iters self.init_method init_method self.random_state random_state self.rng np.random.default_rng(random_state) def _initialize_centroids(self, X): 支持两种初始化随机采样 KMeans n_samples, n_features X.shape if self.init_method random: # 随机选n_clusters个样本作为初始质心 indices self.rng.choice(n_samples, self.n_clusters, replaceFalse) return X[indices].copy() elif self.init_method k-means: # KMeans首个质心随机后续按距离平方概率选择 centroids np.zeros((self.n_clusters, n_features)) # 第一个质心随机 centroids[0] X[self.rng.integers(0, n_samples)] for i in range(1, self.n_clusters): # 计算所有点到已选质心的最小距离平方 distances_sq np.min( np.sum((X[:, None, :] - centroids[:i][None, :, :])**2, axis2), axis1 ) # 按距离平方加权采样 probs distances_sq / distances_sq.sum() new_centroid_idx self.rng.choice(n_samples, pprobs) centroids[i] X[new_centroid_idx] return centroids def _assign_clusters(self, X, centroids): 向量化计算每个点到所有质心的距离返回最近簇索引 # X: (n_samples, n_features), centroids: (n_clusters, n_features) # 使用广播机制计算所有距离(n_samples, n_clusters) distances_sq np.sum((X[:, None, :] - centroids[None, :, :])**2, axis2) return np.argmin(distances_sq, axis1) # (n_samples,) def _update_centroids(self, X, labels): 按簇重新计算质心处理空簇情况 new_centroids np.zeros((self.n_clusters, X.shape[1])) for i in range(self.n_clusters): mask (labels i) if np.any(mask): new_centroids[i] np.mean(X[mask], axis0) else: # 空簇随机选一个远离当前质心的点 dist_to_others np.sum((centroids - centroids[i])**2, axis1) farthest_idx np.argmax(dist_to_others) new_centroids[i] X[self.rng.integers(0, X.shape[0])] return new_centroids def fit(self, X): 主训练循环返回收敛标志和迭代次数 self.centroids self._initialize_centroids(X) self.labels None self.inertia_ [] # 存储每轮SSE for iteration in range(self.max_iters): # 步骤1分配簇 old_labels self.labels.copy() if self.labels is not None else None self.labels self._assign_clusters(X, self.centroids) # 步骤2更新质心 old_centroids self.centroids.copy() self.centroids self._update_centroids(X, self.labels) # 步骤3计算SSE惯性 inertia np.sum( np.sum((X - self.centroids[self.labels])**2, axis1) ) self.inertia_.append(inertia) # 收敛判断质心移动距离 1e-4 或标签不再变化 centroid_shift np.sum(np.sqrt(np.sum((self.centroids - old_centroids)**2, axis1))) if centroid_shift 1e-4 or (old_labels is not None and np.array_equal(self.labels, old_labels)): self.n_iter_ iteration 1 return True self.n_iter_ self.max_iters return False # 未收敛注意_update_centroids中空簇处理是实战关键。EEMD特征常出现极端离群点导致某簇无样本归属。若简单跳过质心保持初始值会引发后续迭代崩溃。此处采用“选最远点”策略比sklearn的“重采样”更稳定——因为EEMD数据空间本身存在天然密度梯度。2.3 参数配置与运行入口如何让203个文件协同工作run_kmeans.py是调度中枢它不硬编码路径而是通过glob动态扫描数据目录并利用concurrent.futures并行处理import glob import os from kmeans_core import CustomKMeans from evaluation import silhouette_score, calinski_harabasz_score from visualization import plot_sse_curve, plot_cluster_result def process_single_file(filepath, k_range(2, 10), init_methodk-means): 单文件全流程加载→标准化→多K值聚类→评估→绘图 print(fProcessing {os.path.basename(filepath)}...) X load_eemd_csv(filepath) X_scaled (X - np.mean(X, axis0)) / (np.std(X, axis0) 1e-8) results {} for k in range(*k_range): km CustomKMeans(n_clustersk, init_methodinit_method, random_state42) converged km.fit(X_scaled) # 计算评估指标 sil_score silhouette_score(X_scaled, km.labels) ch_score calinski_harabasz_score(X_scaled, km.labels) sse km.inertia_[-1] if converged else float(inf) results[k] { converged: converged, silhouette: sil_score, calinski_harabasz: ch_score, sse: sse, n_iter: km.n_iter_, centroids: km.centroids, labels: km.labels } # 保存最佳K值结果按轮廓系数 best_k max(results.keys(), keylambda k: results[k][silhouette]) best_res results[best_k] # 绘制SSE曲线和聚类图 plot_sse_curve([results[k][sse] for k in range(*k_range)], list(range(*k_range)), fsse_{os.path.splitext(os.path.basename(filepath))[0]}.png) plot_cluster_result(X_scaled, best_res[labels], best_res[centroids], fcluster_{os.path.splitext(os.path.basename(filepath))[0]}.png) return filepath, best_k, best_res if __name__ __main__: # 自动匹配所有EEMD相关CSV csv_files glob.glob(data_*.csv) glob.glob(k-data-*.csv) print(fFound {len(csv_files)} CSV files) # 并行处理CPU核心数-1 from concurrent.futures import ProcessPoolExecutor with ProcessPoolExecutor(max_workersos.cpu_count()-1) as executor: futures [executor.submit(process_single_file, f) for f in csv_files] for future in futures: filepath, best_k, result future.result() print(f{os.path.basename(filepath)}: Best K{best_k}, Silhouette{result[silhouette]:.3f})提示glob.glob(data_*.csv)确保新增data_EEMDV2_10.csv无需改代码ProcessPoolExecutor而非ThreadPoolExecutor因NumPy计算是CPU密集型多进程才能真正并行。3. 评估与可视化用43张PNG图验证算法是否“真懂”数据3.1 三维度评估体系不止看轮廓系数项目中43张PNG并非随意生成而是严格对应三大评估维度。evaluation.py提供三个独立函数各自解决不同问题def silhouette_score(X, labels): 手写轮廓系数支持任意距离度量此处为欧氏 n_samples X.shape[0] a np.zeros(n_samples) # 同簇平均距离 b np.zeros(n_samples) # 最近异簇平均距离 for i in range(n_samples): same_cluster (labels labels[i]) if np.sum(same_cluster) 1: a[i] np.mean(np.sqrt(np.sum((X[i] - X[same_cluster])**2, axis1))) else: a[i] 0 # 找最近的其他簇 other_clusters set(labels) - {labels[i]} if other_clusters: b_min float(inf) for other_label in other_clusters: other_mask (labels other_label) if np.any(other_mask): avg_dist np.mean(np.sqrt(np.sum((X[i] - X[other_mask])**2, axis1))) b_min min(b_min, avg_dist) b[i] b_min else: b[i] 0 s np.zeros(n_samples) for i in range(n_samples): if a[i] 0 and b[i] 0: s[i] 0 elif a[i] 0: s[i] 1 else: s[i] (b[i] - a[i]) / max(a[i], b[i]) return np.mean(s) def calinski_harabasz_score(X, labels): CH分数簇间离散度/簇内离散度越大越好 n_samples, n_features X.shape n_clusters len(set(labels)) # 总体中心 overall_mean np.mean(X, axis0) # 簇内离散度WCSS wcss 0 for i in range(n_clusters): cluster_points X[labels i] if len(cluster_points) 0: wcss np.sum((cluster_points - np.mean(cluster_points, axis0))**2) # 簇间离散度BCSS bcss 0 for i in range(n_clusters): cluster_points X[labels i] if len(cluster_points) 0: cluster_mean np.mean(cluster_points, axis0) bcss len(cluster_points) * np.sum((cluster_mean - overall_mean)**2) # CH (BCSS/(K-1)) / (WCSS/(N-K)) if n_clusters 1 or n_samples n_clusters: return 0 return (bcss / (n_clusters - 1)) / (wcss / (n_samples - n_clusters)) def davies_bouldin_score(X, labels): DB指数簇内紧密度/簇间分离度越小越好 n_clusters len(set(labels)) if n_clusters 1: return 0 # 计算每个簇的平均距离簇内紧密度 cluster_dists [] for i in range(n_clusters): cluster_points X[labels i] if len(cluster_points) 1: # 簇内两两距离均值 dists np.sqrt(np.sum((cluster_points[:, None, :] - cluster_points[None, :, :])**2, axis2)) cluster_dists.append(np.mean(dists[np.triu_indices(len(cluster_points), 1)])) else: cluster_dists.append(0) # 计算簇间距离质心距离 centroids np.array([np.mean(X[labels i], axis0) for i in range(n_clusters)]) centroid_dists np.sqrt(np.sum((centroids[:, None, :] - centroids[None, :, :])**2, axis2)) # DB指数 mean_i(max_j≠i (R_ij)) db_scores [] for i in range(n_clusters): r_vals [] for j in range(n_clusters): if i ! j and cluster_dists[i] cluster_dists[j] 0: r_vals.append((cluster_dists[i] cluster_dists[j]) / centroid_dists[i, j]) if r_vals: db_scores.append(max(r_vals)) return np.mean(db_scores) if db_scores else 0注意davies_bouldin_score中cluster_dists计算使用“簇内两两距离均值”而非sklearn的“到质心平均距离”。这对EEMD数据更合理——IMF分量在相空间中呈环状分布质心可能落在空洞处而两两距离反映真实紧凑性。3.2 可视化策略针对高维EEMD特征的降维技巧43张PNG中有12张是t-SNE降维图plot_cluster_result调用而非简单PCA。原因在于EEMD特征维度常达20IMF1~IMF10residuePCA前2主成分方差贡献率常低于40%无法反映聚类结构。visualization.py中from sklearn.manifold import TSNE def plot_cluster_result(X, labels, centroids, save_path): 对高维EEMD数据优先用t-SNE降维再绘图 if X.shape[1] 10: # t-SNE参数针对小样本优化EEMD数据通常n500 tsne TSNE(n_components2, perplexity15, learning_rate200, n_iter1000, random_state42, metriceuclidean) X_2d tsne.fit_transform(X) # 重新计算2D质心投影后 centroids_2d np.array([ np.mean(X_2d[labels i], axis0) for i in range(len(set(labels))) ]) else: # 低维时用PCA更稳定 from sklearn.decomposition import PCA pca PCA(n_components2) X_2d pca.fit_transform(X) centroids_2d np.array([ np.mean(X_2d[labels i], axis0) for i in range(len(set(labels))) ]) plt.figure(figsize(10, 8)) scatter plt.scatter(X_2d[:, 0], X_2d[:, 1], clabels, cmaptab10, alpha0.7, s20) plt.scatter(centroids_2d[:, 0], centroids_2d[:, 1], cred, markerx, s200, linewidths3, labelCentroids) plt.colorbar(scatter) plt.legend() plt.title(fClustering Result (t-SNE/PCA) - {os.path.basename(save_path)}) plt.savefig(save_path, dpi300, bbox_inchestight) plt.close()提示perplexity15是经验值。EEMD数据点密度不均过高的perplexity如30会使局部结构模糊过低如5则噪声放大。该值在data_CEEMDAN.csv和data_EEMDV2_9.csv上实测最优。4. 实战排错当KMeans在EEMD数据上“拒绝收敛”时怎么办4.1 收敛失败的三大典型场景与修复方案项目中k-data-EEMD1.csv曾出现连续3次max_iters300仍不收敛inertia_曲线呈锯齿震荡。通过debug_convergence.py分析发现根本原因场景1质心漂移受离群点主导EEMD分解中residue列存在尖峰脉冲幅值10倍标准差导致质心被拉向异常区域。修复在load_eemd_csv中增加离群点截断# 在load_eemd_csv的col_data处理后插入 col_arr np.clip(col_arr, np.percentile(col_arr, 1), np.percentile(col_arr, 99)) # 截断1%和99%分位数场景2距离计算数值溢出data_EEMDV2_8.csv含负值大数(X - centroids)**2产生inf使argmin返回错误索引。修复在_assign_clusters中添加安全距离计算def _safe_distance_sq(x, y): # 避免(x-y)^2溢出用log-sum-exp技巧 diff x - y if np.any(np.abs(diff) 1e4): # 对大差值用绝对值代替平方牺牲精度保稳定 return np.sum(np.abs(diff)) return np.sum(diff**2) # 替换原距离计算 distances_sq np.array([ [_safe_distance_sq(X[i], centroids[j]) for j in range(self.n_clusters)] for i in range(X.shape[0]) ])场景3空簇连锁反应data_EEMDV2_7.csv在K8时某次迭代产生3个空簇_update_centroids随机选点后新质心又落入稀疏区形成死循环。修复增强空簇处理逻辑def _update_centroids(self, X, labels): new_centroids np.zeros((self.n_clusters, X.shape[1])) for i in range(self.n_clusters): mask (labels i) if np.any(mask): new_centroids[i] np.mean(X[mask], axis0) else: # 策略升级选离所有现有质心最远的点 dist_to_all np.min( np.sqrt(np.sum((X[:, None, :] - self.centroids[None, :, :])**2, axis2)), axis1 ) farthest_idx np.argmax(dist_to_all) new_centroids[i] X[farthest_idx] return new_centroids4.2 K值选择用“肘部法则轮廓系数”双校验表锁定最优解项目中kmeans_optimize.py生成k_optimization_table.csv包含203个文件的K值推荐。其逻辑不是单一指标而是构建决策矩阵文件名K候选SSE轮廓系数CH分数DB指数推荐K理由data_EEMD.csv2-8[120, 95, 82, 78, 76, 75, 74][0.42, 0.51, 0.58, 0.62, 0.59, 0.55, 0.50][12.3, 18.7, 25.1, 28.4, 26.9, 24.2, 21.8][0.85, 0.72, 0.61, 0.55, 0.58, 0.62, 0.66]4轮廓系数峰值且SSE下降趋缓生成逻辑def find_optimal_k(X, k_range(2,10)): sse_list, sil_list, ch_list, db_list [], [], [], [] for k in range(*k_range): km CustomKMeans(n_clustersk, random_state42) km.fit(X) sse_list.append(km.inertia_[-1]) sil_list.append(silhouette_score(X, km.labels)) ch_list.append(calinski_harabasz_score(X, km.labels)) db_list.append(davies_bouldin_score(X, km.labels)) # 肘部点检测SSE一阶导数最大下降点 sse_diff np.diff(sse_list) elbow_k np.argmin(sse_diff) 2 # 2因k从2开始 # 轮廓系数峰值 sil_peak_k np.argmax(sil_list) 2 # 综合决策优先轮廓系数肘部点作为约束 if abs(elbow_k - sil_peak_k) 1: return sil_peak_k else: # 检查CH/DB是否支持sil_peak_k if (ch_list[sil_peak_k-2] ch_list[max(0, sil_peak_k-3)]) and \ (db_list[sil_peak_k-2] db_list[max(0, sil_peak_k-3)]): return sil_peak_k else: return elbow_k注意abs(elbow_k - sil_peak_k) 1是经验阈值。EEMD数据的“肘部”常模糊但轮廓系数对K敏感二者接近才可信。若偏差大说明数据本身不适合KMeans应转向DBSCAN。5. 进阶技巧用Git忽略文件管理实验版本与生产版本5.1 .gitignore中的隐藏逻辑隔离研究性修改与稳定代码项目根目录的.gitignore看似简单实则承载着实验治理策略# 忽略所有CSV数据文件防止仓库膨胀 *.csv # 忽略PNG/JPG输出图可随时重生成 *.png *.jpg # 忽略临时调试文件 debug_*.py temp_*.npy # 但保留关键配置模板 !config_template.py # 保留评估指标基线 !baseline_metrics.json # 忽略用户自定义脚本避免污染主流程 user_scripts/ # 但允许特定工具脚本 !tools/convert_eemd_format.py这使得git status始终干净而真正的变更只发生在16个Python源码文件和2个配置文件中。当需要对比不同初始化策略效果时开发者只需复制kmeans_core.py为kmeans_core_v2.py修改_initialize_centroids方法在run_kmeans.py中导入新类运行后生成的新PNG自动被忽略不影响主分支5.2 一键复现实验用requirements.txt固化环境边界requirements.txt精确锁定版本杜绝“在我机器上好使”问题numpy1.23.5 scipy1.10.0 matplotlib3.7.1 scikit-learn1.2.2 # 仅用于t-SNE非KMeans核心提示scikit-learn仅用于t-SNE核心聚类完全脱离其依赖。这意味着你可以将kmeans_core.py单独拎出在嵌入式设备如树莓派上用numpy轻量运行只要满足numpy1.21.0即可。这是工业场景中模型轻量化的关键设计。执行pip install -r requirements.txt后运行python run_kmeans.py --file data_EEMD.csv --k 4 --init k-means即可在3分钟内复现论文级聚类结果——不是调包的黑盒而是每一步都可审计、可打断、可注入调试逻辑的透明流水线。本文还有配套的精品资源点击获取