import pandas as pd # 用表格对象组织公司横截面与聚类画像
import numpy as np # 提供有限值筛选、矩阵运算与可复现随机数
import matplotlib.pyplot as plt # 绘制降维和聚类的诊断图形
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC', 'Arial Unicode MS'] # 优先使用项目中文字体并保留兼容后备
plt.rcParams['axes.unicode_minus'] = False # 确保主成分负载荷的负号正常显示
from sklearn.decomposition import PCA # 用正交方向概括标准化财务变异
from sklearn.preprocessing import StandardScaler # 消除财务指标量纲对距离的直接支配12 无监督学习(Unsupervised Learning)
12.1 导读
监督学习利用特征 \(X\) 预测响应 \(Y\)；无监督学习则在没有响应变量的条件下描述观测或变量之间的结构。本章聚焦两类常用工具：PCA 用较少维度概括主要变异，聚类用距离或相似度探索样本分组。两者都是探索性方法，结果依赖尺度、样本范围和分析选择，不能自动证明因果关系或“真实类型”。
12.2 学习目标
完成本章后，学生应能：
- 从中心化矩阵推导 PCA 的特征向量问题，并用碎石图、载荷与得分解释降维结果。
- 对给定二维点和初始中心手算至少两轮 K-Means，逐轮报告分配、中心和 WCSS，数值无误。
- 对新数据选择标准化、距离、聚类算法与 \(K\)，并用轮廓系数和多初值稳定性评价，而非预写簇数。
- 区分得分、载荷、距离与 linkage 的输入对象，避免把二维图或簇标签解释为因果或真实类型。
- 在同一财务横截面上完成 PCA 与聚类分析，并据实际输出解释载荷、得分和稳定性。
先修自检
点 \((0,0)\) 与 \((2,0)\) 的均值中心是 \((1,0)\)，到该中心的平方距离和是 \(2\)。若第二问误用欧氏距离和，可先复习平方距离。
无监督学习的挑战
与监督学习相比，无监督学习通常更具挑战性：
- 没有明确的目标：没有响应变量来指导学习过程
- 结果难以评估：缺乏像交叉验证这样通用的评估方法
- 更主观：结果更多依赖研究者的判断
然而，无监督学习在数据探索和模式发现中至关重要。
12.3 无监督学习的主要类型
- 主成分分析（PCA）：降维和可视化
- 聚类分析：发现数据中的子群体
- 关联规则挖掘：发现变量间的关联（如购物篮分析）
- 异常检测：识别数据中的异常点
本章重点介绍前两种方法。
12.4 主成分分析（PCA）
PCA 的经典统计起点是 Pearson 对“最接近点系”的最小二乘几何构造；现代教材通常用最大方差、最小重构误差与 SVD 三种等价视角表述 (Pearson 1901年)。
12.4.1 PCA 的基本思想
假设我们有\(n\) 个观测，每个观测有\(p\) 个特征。当 \(p\) 很大时：
- 如何可视化数据？
- 是否能找到数据的低维表示，保留大部分信息？
核心思想：找到数据变化最大的方向，用这些方向来表示数据。
PCA 的直观解释
想象你在拍摄一个三维物体：
- 从某个角度拍摄，只能看到物体的部分信息
- 如果你找到了最佳拍摄角度，可以用一张二维照片最大程度地展现三维物体
- PCA 就是在寻找数据的“最佳拍摄角度”。
在数学上，PCA 先把每个变量减去其样本均值，再在中心化数据中寻找投影方差最大的方向；这些彼此正交的方向称为主成分。中心化是协方差 PCA 的定义性步骤，是否再除以标准差则是分析者对变量尺度的选择。
12.4.2 第一主成分：从中心化到方差最大化
令原始观测为 \(x_{ij}\)，第 \(j\) 个变量的样本均值为 \(\bar x_j\)，中心化观测为 \(x_{ij}^{c}=x_{ij}-\bar x_j\)。给定单位方向 \(\boldsymbol\phi_1=(\phi_{11},\ldots,\phi_{p1})^\top\)，第 \(i\) 个观测的第一主成分得分为
\[ z_{i1}=\sum_{j=1}^{p}\phi_{j1}(x_{ij}-\bar x_j)=\boldsymbol\phi_1^\top\mathbf{x}_i^c. \tag{12.1}\]
式 12.1 明确把均值从每个变量中移除；这一步使得后续平方投影的平均值能够解释为得分方差。
载荷方向满足单位长度约束
\[ \boldsymbol\phi_1^\top\boldsymbol\phi_1=\sum_{j=1}^{p}\phi_{j1}^2=1. \tag{12.2}\]
式 12.2 排除任意放大载荷即可放大投影方差的退化解，并把问题限定为比较方向。
由于中心化得分的均值为零，最大化得分方差可直接写成
\[ \max_{\boldsymbol\phi_1^\top\boldsymbol\phi_1=1} \frac{1}{n}\sum_{i=1}^{n} \left[\sum_{j=1}^{p}\phi_{j1}(x_{ij}-\bar x_j)\right]^2 =\max_{\boldsymbol\phi_1^\top\boldsymbol\phi_1=1} \boldsymbol\phi_1^\top\mathbf S_n\boldsymbol\phi_1, \tag{12.3}\]
其中 \(\mathbf S_n=(\mathbf X^c)^\top\mathbf X^c/n\)。若采用样本协方差的 \(1/(n-1)\) 约定，载荷方向不变，只有特征值整体缩放。对 式 12.3 构造拉格朗日函数并求一阶条件，可得 \(\mathbf S_n\boldsymbol\phi_1=\lambda_1\boldsymbol\phi_1\)；因此第一载荷是最大特征值对应的特征向量，后续主成分依次在与已有方向正交的子空间中最大化剩余投影方差。
易混淆概念辨析：协方差 PCA 与相关矩阵 PCA
协方差 PCA 只对变量中心化，直接比较原单位下的方差；它适合变量量纲一致且绝对波动大小本身有意义的场景。相关矩阵 PCA 先构造 \(x_{ij}^{s}=(x_{ij}-\bar x_j)/s_j\)，等价于对样本相关矩阵做特征分解；它适合量纲不同或不希望大尺度变量机械支配方向的场景。标准化会改变各变量的相对权重，因而可能改变载荷、得分与下游聚类；它不是 PCA 的普遍必要条件。若某变量样本标准差为零，则不能做相关矩阵 PCA，应先移除或重新定义该变量。
[拓展] 数学推导：PCA 的奇异值分解 (SVD) 视角
虽然上面是从方差最大化的角度定义主成分，但在实际计算和深层代数理论中，PCA 的奇异值分解（Singular Value Decomposition, SVD） 密不可分。
设中心化（列均值为 0）后的数据矩阵为 \(\mathbf{X}^c\in\mathbb{R}^{n\times p}\)。它的完整 SVD 为 \[ \mathbf{X}^c = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T \tag{12.4}\] 式 12.4 使用的对象与 式 12.3 完全相同，区别只在计算视角。 其中 \(\mathbf{U}\in\mathbb{R}^{n\times n}\) 和 \(\mathbf{V}\in\mathbb{R}^{p\times p}\) 为正交矩阵，\(\mathbf{\Sigma}\in\mathbb{R}^{n\times p}\) 为“矩形对角”矩阵，非负对角元 \(\sigma_1\geq\sigma_2\geq\cdots\) 为奇异值。\(\mathbf{V}\) 的列是 \((\mathbf{X}^{c})^{T}\mathbf{X}^{c}\) 的特征向量。
关联：PCA 的联系。 样本的协方差矩阵（按 \(1/n\) 缩放）为 \(\mathbf{C} = \frac{1}{n} (\mathbf{X}^c)^T \mathbf{X}^c\)。 将 SVD 代入 \(\mathbf{C}\) ： \[ \mathbf{C} = \frac{1}{n} (\mathbf{U} \mathbf{\Sigma} \mathbf{V}^T)^T (\mathbf{U} \mathbf{\Sigma} \mathbf{V}^T) = \mathbf{V} \left( \frac{\mathbf{\Sigma}^T\mathbf{\Sigma}}{n} \right) \mathbf{V}^T. \tag{12.5}\] 式 12.5 因而把最大方差问题与右奇异向量连接起来。 若写紧 SVD，中间的对角阵可简记为 \(\operatorname{diag}(\sigma_1^2,\ldots,\sigma_r^2)/n\)，而不是对一个非方阵的 \(\mathbf{\Sigma}\) 直接写 \(\mathbf{\Sigma}^2\)。
这恰好是协方差矩阵的特征值分解！这意味着。 1. 右奇异矩的\(\mathbf{V}\) 的列向量就是 PCA 的主成分载荷（Principal Component Loadings），即 \(\phi\) 向量。 2. 特征的\(\lambda_i = \sigma_i^2 / n\) 代表了主成分轴上的方差。这就是为什么我们可以通过奇异值的平方比来计算各个主成分的方差贡献率。 3. 数据在主成分空间上的投影（主成分得分, PC Scores）为 \(\mathbf{Z} = \mathbf{X}^c\mathbf{V} = \mathbf{U}\mathbf{\Sigma}\)。若先单位方差标准化，则这里的 \(\mathbf X^c\) 相应替换为标准化矩阵，SVD 逻辑不变，但求解对象变为相关矩阵 PCA。
SVD 视角的优势在于：避免了直接计算高维协方差矩阵 \(\mathbf{X}^T \mathbf{X}\)（这可能导致数值不稳定和巨大内存开销），使得工业界的软件（如 Scikit-Learn）在底层能够利用高度优化为SVD 算法稳健而高效地求解大规模PCA。
12.4.3 中国案例：基于财务指标的上市公司结构分析（非项目演示）
为了理解A股市场的内在结构，我们将对上市公司进行主成分分析。我们的目标是从众多的财务指标中提取出反映公司核心特征（如”规模”、“盈利能力”、“杠杆”）的主成分。
下面的 PCA 流程使用 2022Q4 横截面，提取总资产对数、资产负债率、净利率、资产周转率和 ROA，用于说明载荷、得分与维数选择。该横截面按本地数据当前版本构造，因此属于回顾性教学样本，不是历史时点投资证据。
本案例明确选择相关矩阵 PCA：五项指标包含对数金额、比率和收益率，单位与自然波动范围不同；若只中心化，绝对方差较大的指标会机械主导方向，而这里的教学目标是比较相对财务画像。StandardScaler 因此是本案例的尺度决策，不是 PCA 的定义。它只统一量纲，不会消除异常值、行业构成或会计口径差异；若研究问题恰好关心原单位的绝对波动，则应改用协方差 PCA 并解释尺度含义。
# 1. 加载数据
import os # 将在线教材的固定数据根同步给本章后续独立代码块
from pathlib import Path # 使用跨平台路径对象解析显式数据根
BOOK_DATA_DIR = Path('/home/ubuntu/r2_data_mount/data').resolve() # 明文定义在线教材的BOOK_DATA_DIR绝对路径
DATA_DIR = BOOK_DATA_DIR # 保留本章后续代码使用的数据根名称
os.environ['BOOK_DATA_DIR'] = str(BOOK_DATA_DIR) # 为本章后续独立代码块登记同一数据根
assert DATA_DIR.is_dir(), f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}' # 检查数据根
path_fin = DATA_DIR / 'stock' / 'financial_statement.h5' # 财务报表文件路径
path_basic = DATA_DIR / 'stock' / 'stock_basic_data.h5' # 股票基本信息文件路径
missing_pca_paths = [str(required_path) for required_path in (path_fin, path_basic) if not required_path.is_file()] # 收集当前实际缺失的绝对路径
assert not missing_pca_paths, f'BOOK_DATA_DIR={DATA_DIR}；缺少必需文件: {missing_pca_paths}；请修正 BOOK_DATA_DIR 或挂载这些 HDF 文件' # 读取前同时给出数据根与精确缺失项
target_report_period = '2022q4' # 在读取前锁定教学横截面以便下推HDF筛选
financial_statements = pd.read_hdf( # 在存储层只读取目标季度与分析字段
path_fin, where=f'quarter={target_report_period!r}', # 利用quarter数据列下推2022Q4筛选
columns=['order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities', 'operating_revenue', 'net_profit'] # 仅读取PCA所需财务字段
).copy() # 隔离目标横截面供版本筛选
stock_basic_data = pd.read_hdf( # 公司表规模较小但仍只读取合并所需字段
path_basic, columns=['order_book_id', 'citics_2019_l1_name', 'province'] # 仅读取公司键、行业与省份
).copy() # 隔离公司属性映射
required_financial_columns = {'order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities', 'operating_revenue', 'net_profit'} # 定义财务字段
required_basic_columns = {'order_book_id', 'citics_2019_l1_name', 'province'} # 定义公司字段
assert required_financial_columns.issubset(financial_statements.columns), '财务报表缺少 PCA 所需字段' # 检查财务字段要求
assert required_basic_columns.issubset(stock_basic_data.columns), '公司表缺少行业或省份字段' # 检查公司字段要求下面在独立代码块中预先固定横截面并合并行业与省份字段；拆分只改善课堂阅读，不改变对象依赖。
# 2. 数据清洗与特征构造
# 固定同一会计报告期，避免把不同年份的公司截面混在一起
# 转换日期并固定同一会计报告期
financial_statements['info_date'] = pd.to_datetime(financial_statements['info_date'], errors='coerce') # 转换信息可得日
assert financial_statements['info_date'].notna().all(), 'info_date 含缺失或无效日期' # 检查日期完整性
period_mask = financial_statements['quarter'].astype(str).str.lower() == target_report_period # 只保留事先固定的共同会计期
recent_financials = (financial_statements.loc[period_mask] # 从共同报告期开始构造横截面
.sort_values('info_date') # 按披露日排序以识别最后版本
.groupby('order_book_id').tail(1)) # 同一报告期内保留最后披露版本
finite_columns = ['total_assets', 'total_liabilities', 'operating_revenue', 'net_profit'] # 指定数值字段
recent_financials = recent_financials[np.isfinite(recent_financials[finite_columns]).all(axis=1)].copy() # 保留有限值
print(f'统一报告期: {target_report_period}; 公司数: {recent_financials.order_book_id.nunique()}; 最晚信息日: {recent_financials.info_date.max().date()}') # 输出样本身份
# 合并基本信息 (获取行业)：stock_basic_data 中行业列名为 citics_2019_l1_name，此处重命名为 industry_name 便于后续引用
stock_basic_with_industry = stock_basic_data[['order_book_id', 'citics_2019_l1_name', 'province']].rename(columns={'citics_2019_l1_name': 'industry_name'}) # 提取并重命名行业分类列
merged_financial_data = pd.merge(recent_financials, stock_basic_with_industry, on='order_book_id', how='inner') # 合并财务数据与基本信息
# 仅保留长三角地区（上海、江苏、浙江、安徽）的上市公司进行分析
yangtze_river_delta_provinces = ['上海市', '江苏省', '浙江省', '安徽省'] # 定义长三角省份列表（省份名含行政后缀）
merged_financial_data = merged_financial_data[merged_financial_data['province'].isin(yangtze_river_delta_provinces)].copy() # 筛选长三角地区上市公司统一报告期: 2022q4; 公司数: 5370; 最晚信息日: 2026-04-30在成功加载并筛选出长三角地区上市公司数据后，下一步是构造关键财务指标并进行数据标准化处理。
# 定义原始列名到标准列名的映射关系
column_name_mapping = { # 统一原始财务字段与教学特征命名
'total_assets': 'Total_Assets', # 将资产总额映射为规模分母
'total_liabilities': 'Total_Liabilities', # 将负债总额映射为杠杆分子
'operating_revenue': 'Revenue', # 将营业收入映射为周转率和利润率分母
'net_profit': 'Net_Income' # 将净利润映射为盈利能力分子
} # 完成四个原始金额字段的语义映射
# 检查数据中是否包含所有必需的财务列
available_financial_columns = [c for c in column_name_mapping.keys() if c in merged_financial_data.columns] # 筛选存在的列
# 提取需要的列并重命名为标准英文名
financial_analysis_data = merged_financial_data[available_financial_columns + ['order_book_id', 'industry_name', 'province']].copy() # 提取财务数据与画像所需省份
financial_analysis_data.rename(columns=column_name_mapping, inplace=True) # 将中文/原始列名统一为英文标准名
# 计算关键财务比率指标
valid_denominators = (financial_analysis_data['Total_Assets'] > 0) & (financial_analysis_data['Revenue'] != 0) # 排除无法形成规模和利润率的分母
financial_analysis_data = financial_analysis_data.loc[valid_denominators].copy() # 仅在有效分母公司上构造财务比率
financial_analysis_data['Log_Assets'] = np.log(financial_analysis_data['Total_Assets']) # 1. 规模指标
financial_analysis_data['Debt_Ratio'] = financial_analysis_data['Total_Liabilities'] / financial_analysis_data['Total_Assets'] # 2. 杠杆指标
financial_analysis_data['Net_Margin'] = financial_analysis_data['Net_Income'] / financial_analysis_data['Revenue'] # 3. 盈利指标
financial_analysis_data['Asset_Turnover'] = financial_analysis_data['Revenue'] / financial_analysis_data['Total_Assets'] # 4. 效率指标
financial_analysis_data['ROA'] = financial_analysis_data['Net_Income'] / financial_analysis_data['Total_Assets'] # 5. 回报指标
# 清洗：去除inf, NaN, 极端值
financial_features = ['Log_Assets', 'Debt_Ratio', 'Net_Margin', 'Asset_Turnover', 'ROA'] # 定义用于PCA的特征列表
cleaned_financial_data = financial_analysis_data.replace([np.inf, -np.inf], np.nan).dropna(subset=financial_features) # 替换无穷值为NaN后删除缺失行
# 按1%和99%分位数去除离群值（极端的利润率等）
for col in financial_features: # 对每个PCA输入独立执行尾部截断
lower_bound = cleaned_financial_data[col].quantile(0.01) # 计算下1%分位数
upper_bound = cleaned_financial_data[col].quantile(0.99) # 计算上99%分位数
cleaned_financial_data = cleaned_financial_data[(cleaned_financial_data[col] >= lower_bound) & (cleaned_financial_data[col] <= upper_bound)] # 保留分位数范围内的数据
print(f'最终样本量: {len(cleaned_financial_data)} 家公司') # 报告分母与尾部筛选后的分析规模
print(cleaned_financial_data[financial_features].describe().round(4)) # 核对五项输入的中心、尺度和支持范围最终样本量: 1814 家公司
Log_Assets Debt_Ratio Net_Margin Asset_Turnover ROA
count 1814.0000 1814.0000 1814.0000 1814.0000 1814.0000
mean 22.0787 0.3928 0.0780 0.5973 0.0443
std 1.1886 0.1889 0.1304 0.3303 0.0545
min 19.9477 0.0514 -0.8610 0.0595 -0.1501
25% 21.2092 0.2442 0.0302 0.3685 0.0155
50% 21.8503 0.3768 0.0783 0.5431 0.0430
75% 22.7197 0.5340 0.1422 0.7530 0.0756
max 26.4492 0.9172 0.4659 2.2059 0.2127# 3. 标准化：将所有特征缩放为均值0、方差1，消除量纲差异
feature_scaler = StandardScaler() # 创建标准化器实例
scaled_features_matrix = feature_scaler.fit_transform(cleaned_financial_data[financial_features]) # 对财务特征进行标准化变换
scaled_features_df = pd.DataFrame(scaled_features_matrix, columns=financial_features) # 将标准化结果转为DataFrame便于查看样本量与描述统计以代码输出为准，并应与打印的统一报告期一起记录。标准化使不同量纲的指标可以共同进入 PCA，但不会消除异常值处理、行业构成或会计口径带来的结构差异。
12.4.4 拟合 PCA 模型
我们将分析这 5 个核心财务指标的主成分。
标准化后的 5 维财务数据进入 PCA 后，首先要回答“保留多少维才能达到事前信息保留标准”。代码现场输出单项与累计解释方差；碎石图只帮助检查谱形状，不能预写前两维或前三维一定达到某个比例。用于聚类的维数必须由预注册阈值和现场累计解释率共同决定，二维得分仅服务于展示。
# 拟合PCA模型，保留所有主成分
pca_model = PCA() # 创建PCA实例，默认保留全部主成分
pca_scores = pca_model.fit_transform(scaled_features_matrix) # 拟合并保存样本得分
# 计算每个主成分的解释方差比例和累计解释方差
explained_variance_summary = pd.DataFrame({ # 构建DataFrame数据表
'PC': [f'PC{i+1}' for i in range(len(pca_model.explained_variance_ratio_))], # 主成分编号
'Explained Variance': pca_model.explained_variance_ratio_, # 各主成分解释方差比例
'Cumulative Variance': np.cumsum(pca_model.explained_variance_ratio_) # 累计解释方差比例
}) # 完成单项与累计解释率的同表对照
print(explained_variance_summary.round(4)) # 输出选择保留维数所需的数值证据 PC Explained Variance Cumulative Variance
0 PC1 0.4141 0.4141
1 PC2 0.2740 0.6881
2 PC3 0.1962 0.8843
3 PC4 0.0853 0.9695
4 PC5 0.0305 1.0000表 12.2 给出现场计算的单项与累计解释率；图 12.1 只把同一组数值画成碎石图，不承担表格输出。
# 绘制碎石图（Scree Plot），用于可视化判断主成分数量的选择
plt.figure(figsize=(10, 5)) # 为五个主成分的碎石诊断设置横向画布
plt.plot(range(1, len(pca_model.explained_variance_ratio_)+1), pca_model.explained_variance_ratio_, 'bo-') # 绘制解释方差比例折线图
plt.xlabel('主成分', fontsize=12, fontproperties='Source Han Serif SC') # 横轴按方差从大到小排列成分编号
plt.ylabel('解释方差比例', fontsize=12, fontproperties='Source Han Serif SC') # 纵轴给出各成分的样本方差占比
plt.title('PCA 碎石图(Scree Plot)', fontsize=14, fontproperties='Source Han Serif SC') # 标明图形用于检查特征值谱的衰减
plt.grid(True, alpha=0.3) # 辅助比较相邻成分的解释率落差
plt.show() # 渲染各主成分的解释方差比例以辅助选择维数
解释方差比例由当前统一报告期样本现场计算。选择主成分数量时应结合累计解释方差、碎石图、下游任务和跨报告期稳定性；“肘部”是诊断而非自动成立的结论。
12.4.5 主成分载荷与可视化
解释 PCA 需要同时查看得分与载荷。双标图把样本投影和变量方向放在同一平面，但它只是高维结构的二维摘要。本图的箭头坐标是 PCA.components_ 给出的特征向量系数载荷 \((v_{j1},v_{j2})\)，而不是与主成分的相关载荷。因此，两根箭头在图上接近只表示它们的前两个系数载荷投影方向接近，不能单凭夹角断言原始变量“极强正相关”。若要用箭头内积近似变量相关，应改用 \(√{\lambda_k}v_{jk}\) 的相关载荷，并且只有在所画主成分已承载几乎全部方差时，二维夹角才是良好近似；精确相关需使用全部主成分。 主成分的商业名称必须由当次载荷决定。若资产规模或杠杆在 PC1 上绝对载荷最大，才可将 PC1 暂述为“规模/杠杆”轴；若利润率或周转率主导 PC2，才可进一步讨论“盈利/效率”。载荷符号可整体翻转，因此解释应关注相对大小和方向关系，而不是预设固定的“商业十字架”。
# 构建载荷矩阵：行为原始特征，列为主成分
pca_loadings_df = pd.DataFrame(pca_model.components_.T, columns=[f'PC{i+1}' for i in range(len(financial_features))], index=financial_features) # 转置载荷矩阵便于阅读
print('\n主成分载荷(Loadings):') # 标识随后输出为前两维的系数载荷
print(pca_loadings_df.iloc[:, :2].round(3)) # 显示前两个主成分的载荷值
# 双标图同时展示样本得分与系数载荷；夹角只表示二维投影几何
fig, ax = plt.subplots(figsize=(10, 8)) # 为得分与载荷共用坐标系设置近方形画布
display_count = min(500, len(pca_scores)) # 限制散点数量以保持双标图可读
display_indices = np.random.default_rng(42).choice(len(pca_scores), display_count, replace=False) # 固定抽样使展示可复现
score_scale = np.max(np.abs(pca_scores[display_indices, :2]), axis=0) # 按展示样本统一缩放前两维得分
normalized_scores = pca_scores[display_indices, :2] / score_scale # 将得分压到与载荷箭头相近的显示尺度
ax.scatter(normalized_scores[:, 0], normalized_scores[:, 1], s=12, alpha=0.2, # 用半透明点显示公司得分密度
color='steelblue', label='公司得分') # 固定颜色避免与红色载荷箭头混淆
# 遍历每个财务指标，绘制从原点出发的红色载荷箭头
for i, feature in enumerate(financial_features): # 遍历并枚举每个元素
ax.arrow(0, 0, pca_model.components_[0, i], pca_model.components_[1, i], # 从原点显示指标在前两维的系数载荷方向
head_width=0.05, head_length=0.08, fc='red', ec='red') # 绘制载荷箭头
ax.text(pca_model.components_[0, i]*1.15, pca_model.components_[1, i]*1.15, # 添加文本标签
feature, color='red', fontsize=12, ha='center', va='center') # 在箭头末端标注特征名
ax.set_xlabel(f'PC1 ({pca_model.explained_variance_ratio_[0]:.1%})', fontsize=12) # X轴：第一主成分及解释方差比例
ax.set_ylabel(f'PC2 ({pca_model.explained_variance_ratio_[1]:.1%})', fontsize=12) # Y轴：第二主成分及解释方差比例
ax.set_title('PCA 二维投影：公司得分与系数载荷', fontsize=14, fontproperties='Source Han Serif SC') # 明示箭头的统计约定
ax.grid(True, alpha=0.3) # 辅助读取公司得分与载荷箭头的相对位置
ax.axhline(0, color='black', lw=1) # 绘制水平零线
ax.axvline(0, color='black', lw=1) # 绘制垂直零线
plt.axis('equal') # 设置等比例坐标轴以避免视觉失真
plt.show() # 渲染公司得分与载荷箭头的共同二维投影
# 输出当前数据的主导载荷，避免用固定模板命名主成分
for component_index in range(2): # 分别核对展示平面的两条主成分轴
component_loadings = pd.Series(pca_model.components_[component_index], index=financial_features) # 对齐载荷与原始指标名称
leading_features = component_loadings.abs().nlargest(2).index # 选择绝对载荷最大的两个解释线索
signed_terms = ', '.join(f'{name}={component_loadings[name]:+.3f}' for name in leading_features) # 保留方向供条件式命名
print(f'PC{component_index + 1} 主导载荷: {signed_terms}') # 输出当次数据决定的主导指标
主成分载荷(Loadings):
PC1 PC2
Log_Assets -0.293 0.519
Debt_Ratio -0.465 0.478
Net_Margin 0.592 0.249
Asset_Turnover 0.001 0.541
ROA 0.589 0.384
PC1 主导载荷: Net_Margin=+0.592, ROA=+0.589
PC2 主导载荷: Asset_Turnover=+0.541, Log_Assets=+0.519图 12.2 把得分与系数载荷箭头叠加在同一二维坐标中。主成分名称必须依据上方当前载荷表给出。载荷整体变号不改变主成分空间，因此“正方向”本身没有经济优劣含义；只有系数载荷的绝对值、相对符号和解释方差可用于描述主成分方向，原始变量相关应另行计算。
12.4.6 主成分得分与行业分布
在依据当前载荷确定 PC1、PC2 的描述后，可以检查行业标签在该二维投影中的重叠与离散。行业位置、方向和分离程度必须从本次载荷、组内离散度与图形读取；在看到输出前不预设任何行业落在哪个象限。
pca_transformed_features = pca_model.transform(scaled_features_matrix) # 将标准化数据投影到主成分空间
cleaned_financial_data['PC1'] = pca_transformed_features[:, 0] # 将第一主成分得分存入数据框
cleaned_financial_data['PC2'] = pca_transformed_features[:, 1] # 将第二主成分得分存入数据框
# 选取公司数量最多的前6个行业进行可视化
top_industry_indices = cleaned_financial_data['industry_name'].value_counts().head(6).index # 获取最大的六个行业
plot_subset_data = cleaned_financial_data[cleaned_financial_data['industry_name'].isin(top_industry_indices)] # 筛选这六个行业的数据
plt.figure(figsize=(12, 8)) # 为六个行业的二维重叠保留足够图例空间
# 按行业分颜色绘制散点图
for industry_name in top_industry_indices: # 逐行业叠加散点以观察组间重叠
industry_subset = plot_subset_data[plot_subset_data['industry_name'] == industry_name] # 筛选当前行业子集
plt.scatter(industry_subset['PC1'], industry_subset['PC2'], label=industry_name, alpha=0.6, s=30) # 绘制该行业散点
plt.xlabel('PC1', fontsize=12) # 横轴为第一主成分公司得分
plt.ylabel('PC2', fontsize=12) # 纵轴为第二主成分公司得分
plt.title('上市公司行业分布 (PCA空间)', fontsize=14, fontproperties='Source Han Serif SC') # 标明颜色仅用于观察行业重叠
plt.legend(prop={'family': 'Source Han Serif SC'}) # 添加中文图例
plt.grid(True, alpha=0.3) # 辅助读取行业组内离散与组间重叠
plt.show() # 渲染行业在前两主成分上的重叠与离散程度
图 12.3 是探索性投影。行业是否聚集、哪些行业重叠，应从本次图形和组内离散度计算判断；二维分离不能证明行业是主要驱动因素，也不能自动转化为分类或投资结论。
PCA 的实际应用建议
- 按分析目标决定是否标准化：本例希望各项不同量纲的财务指标具有可比权重，因此先标准化；若原始方差本身具有明确含义，则应说明保留原尺度的理由
- 主成分数量选择：
- 碎石图法：选择拐点处的主成分数
- 累积方差法：选择累积解释方差达到70-90%的主成分
- Kaiser 准则：对标准化变量，常用经验规则是保留特征值大于 1 的主成分
- 解释主成分：通过载荷矩阵理解每个主成分的实际含义
- 主成分得分：可用于后续的聚类、回归等分析
12.5 聚类分析
聚类分析的目标是将相似的观测归为一组。在金融中，我们可以用它来发现财务特征相似的公司群体，这对于寻找对标公司（Benchmarking）非常有帮助。
12.5.1 二维 K-Means 两轮手算
K-Means 在固定 \(K\) 时最小化组内平方和
\[ J=\sum_{k=1}^{K}\sum_{i:z_i=k}\lVert x_i-\mu_k\rVert_2^2. \tag{12.6}\]
算法只交替做两件事：把每个点分给最近中心，再把每个中心更新为簇内均值。取六个点
\[ A=(0,0),\ B=(0,2),\ C=(2,0),\ D=(8,8),\ E=(8,6),\ F=(6,8), \]
令 \(K=2\)，初始中心为 \(\mu_1^{(0)}=A=(0,0)\) 与 \(\mu_2^{(0)}=C=(2,0)\)。平方距离足以比较，避免无意义地开平方。
第 1 轮。 分配结果为 \(C_1^{(1)}=\{A,B\}\)、\(C_2^{(1)}=\{C,D,E,F\}\)。更新中心：
\[ \mu_1^{(1)}=(0,1),\qquad \mu_2^{(1)}=(6,5.5). \]
用更新后的中心计算目标函数，第一簇贡献 \(1+1=2\)，第二簇贡献 \(46.25+10.25+4.25+6.25=67\)，所以 \(J^{(1)}=69\)。
第 2 轮。 重新比较距离后，\(C\) 转入第一簇：\(C_1^{(2)}=\{A,B,C\}\)、\(C_2^{(2)}=\{D,E,F\}\)。中心更新为
\[ \mu_1^{(2)}=(2/3,2/3),\qquad \mu_2^{(2)}=(22/3,22/3). \]
两簇各贡献 \(16/3\)，故 \(J^{(2)}=32/3\approx10.667<69\)。再做一次分配不会改变标签，算法停止。这个例子也暴露初始化风险：较差初值会让第一轮形成极不平衡簇；因此真实分析应使用多初值并报告稳定性。
import numpy as np # 用纯内存数组复核两轮手算
two_dimensional_points = np.array([[0, 0], [0, 2], [2, 0], [8, 8], [8, 6], [6, 8]], dtype=float) # 固定六点
cluster_centers = np.array([[0, 0], [2, 0]], dtype=float) # 使用题设的刻意较差初值
iteration_audit = [] # 保存每轮标签、中心与目标函数
for iteration_number in range(1, 3): # 严格执行两轮以对应手算
squared_distances = ((two_dimensional_points[:, None, :] - cluster_centers[None, :, :]) ** 2).sum(axis=2) # 算平方距离
cluster_labels = squared_distances.argmin(axis=1) # 把每点分到最近中心
cluster_centers = np.vstack([two_dimensional_points[cluster_labels == k].mean(axis=0) for k in range(2)]) # 更新均值中心
within_cluster_sum = ((two_dimensional_points - cluster_centers[cluster_labels]) ** 2).sum() # 计算更新后 WCSS
iteration_audit.append((iteration_number, cluster_labels.tolist(), cluster_centers.copy(), within_cluster_sum)) # 留存核对
print(iteration_audit) # 输出两轮证据供逐项比对[(1, [0, 0, 1, 1, 1, 1], array([[0. , 1. ],
[6. , 5.5]]), 69.0), (2, [0, 0, 0, 1, 1, 1], array([[0.66666667, 0.66666667],
[7.33333333, 7.33333333]]), 10.666666666666666)]当数据未标准化时，大尺度变量会支配平方距离；当簇并非近似球形、含强异常点或需要概率归属时，K-Means 可能不是合适工具。学生现在已完成 K-Means 核心路线，才进入下面的概率拓展。
拓展：EM 与高斯混合模型（完成 K-Means 核心路线后再展开）
选读：期望最大化 (EM) 算法的详细阐释 {#sec-em-algorithm}
第一次阅读请先跳到后面的 K-Means 算法与目标函数；掌握“分配—更新”迭代后再回到本选读。EM 是理解高斯混合模型的统一框架，但不是学习 K-Means 基本操作的先修门槛。
EM 算法要解决什么问题？
EM 算法用于解决含有隐变量(latent variables)的概率模型的最大似然估计(Maximum Likelihood Estimation, MLE)问题。所谓”隐变量”，是指我们知道它存在、但无法直接观测的变量。
在聚类问题中，经典的隐变量就是每个数据点所属的簇标签。我们能观测到数据点的特征（如公司的财务指标），但不知道它属于哪个类群。给定簇标签可以按组估计参数；给定参数则可以计算各簇的后验概率。EM 算法通过交替执行这两个条件步骤，处理标签与参数相互依赖的估计问题。
EM 算法的核心思想：交替迭代
EM 算法通过在两个步骤之间交替迭代来逐步逼近最优解：
E 步 (Expectation Step，期望步)： 固定当前的模型参数估计值，计算每个数据点属于各个簇的后验概率（也称”责任”或”软分配”）。具体来说，对于数据点\(x_i\)属于第\(k\)个簇的概率为：
\[ \gamma_{ik} = P(Z_i = k | x_i, \theta^{(t)}) = \frac{\pi_k^{(t)} \cdot f(x_i | \mu_k^{(t)}, \Sigma_k^{(t)})}{\sum_{l=1}^{K} \pi_l^{(t)} \cdot f(x_i | \mu_l^{(t)}, \Sigma_l^{(t)})} \tag{12.7}\]
其中 \(\theta^{(t)} = \{\pi_k^{(t)}, \mu_k^{(t)}, \Sigma_k^{(t)}\}\) 是第 \(t\) 次迭代的参数估计，\(Z_i\) 是数据点 \(i\) 的隐含簇标签，\(f(\cdot)\) 是高斯密度函数。
直观理解：E步就像是一位金融分析师，在已知各个”公司群体”（簇）的特征画像后，评估每家上市公司最可能属于哪个群体的概率。
M 步 (Maximization Step，最大化步)： 固定 E 步算出的后验概率，重新估计模型参数，使得数据的期望对数似然最大化。更新公式为：
\[ \pi_k^{(t+1)} = \frac{1}{n} \sum_{i=1}^{n} \gamma_{ik} \tag{12.8}\]
\[ \mu_k^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik} \cdot x_i}{\sum_{i=1}^{n} \gamma_{ik}} \tag{12.9}\]
\[ \Sigma_k^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik} \cdot (x_i - \mu_k^{(t+1)})(x_i - \mu_k^{(t+1)})^T}{\sum_{i=1}^{n} \gamma_{ik}} \tag{12.10}\]
直观理解：M步就像是根据各公司的”归属概率”，重新计算每个群体的中心特征和分散程度。属于某个群体概率越高的公司，在计算该群体中心时的权重就越大。
EM 算法的收敛保证
EM 算法的一个重要理论性质是：每次迭代都保证数据的对数似然函数单调递增（或至少不减少）。数学上可以证明：
\[ \log L(\theta^{(t+1)}) \geq \log L(\theta^{(t)}) \]
在目标有上界、协方差受非退化约束等正则条件下，这支持目标值趋于极限，极限点通常是驻点；它不保证参数序列总收敛，也不保证驻点是局部最优。未约束 GMM 还可能因协方差塌缩使似然无界。实践中应限制协方差、使用多初值，并同时报告失败与退化运行。
K-Means 是 EM 算法的硬分配极端特例
理解了完整的 EM 算法后，我们就能清晰地看到 K-Means 与 EM 的关系：
| 特征 | EM 算法 (GMM) | K-Means |
|---|---|---|
| E 步（分配方式） | 软分配：\(\gamma_{ik} \in [0, 1]\)，每个点按概率”部分地”属于各个簇 | 硬分配：\(\gamma_{ik} \in \{0, 1\}\)，每个点完全属于距离最近的簇 |
| M 步（参数更新） | 按概率加权计算均值和协方差 | 仅计算已分配点的算术平均值 |
| 簇的形状 | 可以是椭圆形（取决于协方差参数化） | 平方欧氏距离偏好凸、近球形结构；实际簇样本量不必相等 |
| 目标函数 | 最大化数据的对数似然 | 最小化组内平方和 (WCSS) |
当我们将 GMM 中每个高斯分量的协方差矩阵退化为极小的各向同性矩阵 \(\epsilon I\)（\(\epsilon \to 0\)），E 步中的软概率就会坍缩为 0 或 1 的硬分配——每个点被确定性地分配给距离最近的均值中心。这正是 K-Means 的分配规则。
这一联系说明 K-Means 对尺度和几何形状敏感，但不意味着算法强制输出等样本量簇。当群体呈现椭圆、不同密度或强烈不均衡时，应比较标准化、其他距离、GMM 或其他聚类模型，并用现场稳定性证据判断。
12.5.2 K-Means 聚类
我们按事先设定的 80% 累计解释率选择最小维数 \(d\) 进行聚类；前两维只用于显示。代码在 \(K=2,\ldots,8\) 中只保留满足 \(n\ge5K\) 的候选值，并检查输入有限性与实际簇数。
from sklearn.cluster import KMeans # 最小化标准化财务得分的簇内平方和
from sklearn.metrics import adjusted_rand_score, silhouette_score # 聚类质量与标签稳定性
variance_threshold = 0.80 # 事先设定累计解释率阈值
selected_component_count = int(np.searchsorted(np.cumsum(pca_model.explained_variance_ratio_), variance_threshold) + 1) # 最小达标d
clustering_features = np.asarray(pca_transformed_features[:, :selected_component_count], dtype=float) # 动态d维聚类输入
assert clustering_features.ndim == 2 and clustering_features.shape[1] == selected_component_count, '聚类矩阵维度不符' # 检查矩阵形状
assert np.isfinite(clustering_features).all(), '聚类输入含非有限值' # 检查数值有限性
# 只遍历事前登记且满足 n>=5K 的候选，不存在合法 K 时不拟合
cluster_count_range = [k for k in range(2, 9) if len(clustering_features) >= 5 * k] # 按样本量筛选候选K
assert cluster_count_range, f'样本量 {len(clustering_features)} 下没有满足 n>=5K 的候选' # 检查候选非空
sum_squared_errors = [] # 存储各K值对应的组内误差平方和
silhouette_scores = [] # 存储各K值对应的轮廓系数
for k in cluster_count_range: # 遍历每个候选K值
kmeans_model = KMeans(n_clusters=k, random_state=42, n_init=10) # 初始化K-Means模型
candidate_labels = kmeans_model.fit_predict(clustering_features) # 对聚类特征数据拟合模型
assert np.unique(candidate_labels).size == k, f'K={k} 未形成足够的非空簇' # 检查实际簇数
sum_squared_errors.append(kmeans_model.inertia_) # 记录SSE（inertia）
candidate_silhouette = silhouette_score(clustering_features, candidate_labels) # 衡量当次K下样本与本簇相对邻簇的分离程度
assert np.isfinite(candidate_silhouette), f'K={k} 的轮廓系数不是有限数' # 检查评价值
silhouette_scores.append(candidate_silhouette) # 记录轮廓系数
# 绘制双子图：左侧为肘部法则图，右侧为轮廓系数图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # 创建1行2列子图
ax1.plot(cluster_count_range, sum_squared_errors, 'bo-') # 绘制SSE折线图
ax1.set_xlabel('K', fontsize=12) # 横轴为预先给定的候选簇数
ax1.set_ylabel('Inertia (SSE)', fontsize=12) # 纵轴为各候选K的组内平方和
ax1.set_title('肘部法则', fontsize=14, fontproperties='Source Han Serif SC') # 标明左图用于检查组内平方和降幅
ax1.grid(True, alpha=0.3) # 辅助比较增加一个簇带来的SSE变化
ax2.plot(cluster_count_range, silhouette_scores, 'rs-') # 绘制轮廓系数折线图
ax2.set_xlabel('K', fontsize=12) # 与左图共享同一候选簇数集合
ax2.set_ylabel('轮廓系数 (Silhouette)', fontsize=12) # 纵轴比较簇内紧密与簇间分离
ax2.set_title('轮廓系数', fontsize=14, fontproperties='Source Han Serif SC') # 标明右图用于选择现场分离度较高的K
ax2.grid(True, alpha=0.3) # 辅助读取候选K之间的轮廓系数差异
plt.tight_layout() # 自动调整子图间距
plt.show() # 并排渲染候选K的组内平方和与轮廓系数诊断
下面将轮廓系数最大的候选值作为本次样本的数据驱动选择。这是一项探索性规则，不意味着存在固定不变的“真实 \(K\)”；应同时检查肘部图、稳定性与业务可解释性。
optimal_cluster_count = list(cluster_count_range)[int(np.argmax(silhouette_scores))] # 候选集中轮廓系数最大的K
optimal_kmeans_model = KMeans(n_clusters=optimal_cluster_count, random_state=42, n_init=10) # 初始化最终K-Means模型
cluster_assignments = optimal_kmeans_model.fit_predict(clustering_features) # 拟合并预测每个公司的聚类标签
assert np.unique(cluster_assignments).size == optimal_cluster_count, '最终模型未形成指定簇数' # 检查最终簇数
# 在同一真实样本上以至少20个单初值种子检验标签稳定性
initialization_stability = []
for initialization_seed in range(20): # 用二十个单初值复算评估初始化敏感性
repeated_labels = KMeans(n_clusters=optimal_cluster_count, random_state=initialization_seed, n_init=1).fit_predict(clustering_features) # 在相同样本上只改变初始中心
assert np.unique(repeated_labels).size == optimal_cluster_count, f'种子{initialization_seed}出现空簇' # 检查重复拟合
initialization_stability.append(adjusted_rand_score(cluster_assignments, repeated_labels)) # 用置换不变的ARI比较参考标签
stability_summary = { # 汇总而不只报告最有利的一次初始化
'seed_count': len(initialization_stability), # 记录实际比较的随机种子数
'mean_ari': float(np.mean(initialization_stability)), # 描述平均标签一致度
'min_ari': float(np.min(initialization_stability)), # 保留最差运行暴露不稳定风险
'std_ari': float(np.std(initialization_stability)), # 描述初始化结果的离散程度
} # 完成初始化稳定性摘要
cleaned_financial_data['Cluster'] = cluster_assignments # 将聚类结果存入数据框
print({'selected_k': optimal_cluster_count, 'report_period': target_report_period, 'stability': stability_summary}) # 输出选择与稳定性{'selected_k': 3, 'report_period': '2022q4', 'stability': {'seed_count': 20, 'mean_ari': 0.8720045005489385, 'min_ari': 0.47656499129212265, 'std_ari': 0.2194989416977441}}图 12.4 同时显示可行候选的 WCSS 与轮廓系数。若最小 ARI 低于 0.8，应把结果解释为对初始化敏感，而不是中止或只报告最好的一次。
12.5.3 聚类结果解释
聚类编号只是当次运行的数字代号，且簇数由上一节的候选规则动态决定。下面回到原始财务特征，按实际簇编号计算均值，并在前两个主成分上显示分群。颜色代表算法分配，不代表已知的“真实同盟”。
代码还会报告每簇的众数行业、平均资产和负债率。只有看到现场输出后，才可描述某簇是否呈现“重资产”或“高周转”等模式。这些是探索性画像，不是无偏的行业分类，也不足以直接支持并购或投资决策。
# 分析每个簇的特征均值(使用原始特征而非标准化后的值，便于直观理解)
cluster_feature_means = cleaned_financial_data.groupby('Cluster')[financial_features].mean() # 按聚类分组计算各财务指标均值
print(cluster_feature_means.round(3)) # 在原尺度输出各簇的条件式财务画像 Log_Assets Debt_Ratio Net_Margin Asset_Turnover ROA
Cluster
0 21.727 0.433 -0.163 0.426 -0.052
1 22.967 0.551 0.060 0.731 0.037
2 21.514 0.269 0.144 0.538 0.071表 12.3 报告原始财务尺度上的簇均值；图 12.5 是同一簇标签在前两个主成分上的展示投影。
# 在PCA空间中可视化聚类结果
plt.figure(figsize=(10, 8)) # 为簇标签的二维展示设置近方形画布
scatter = plt.scatter(cleaned_financial_data['PC1'], cleaned_financial_data['PC2'], c=cleaned_financial_data['Cluster'], cmap='viridis', s=50, alpha=0.7) # 按聚类标签着色
plt.colorbar(scatter, label='Cluster') # 添加颜色条标注聚类编号
plt.xlabel('PC1', fontsize=12, fontproperties='Source Han Serif SC') # 载荷需由当次运行结果解释
plt.ylabel('PC2', fontsize=12, fontproperties='Source Han Serif SC') # 不预先把轴命名为特定商业含义
plt.title('上市公司聚类分布', fontsize=14, fontproperties='Source Han Serif SC') # 标明颜色是算法簇编号而非真实类别
plt.grid(True, alpha=0.3) # 辅助读取二维投影中的簇内离散与重叠
plt.show() # 渲染按现场K-Means簇标签着色的PCA投影
# 逐簇输出代表性行业与核心财务特征
for cluster_id in range(optimal_cluster_count): # 遍历每个聚类
dominant_industry = cleaned_financial_data[cleaned_financial_data['Cluster'] == cluster_id]['industry_name'].mode()[0] # 获取该簇中出现频率最高的行业
print(f'Cluster {cluster_id}: 代表行业 [{dominant_industry}], 平均资产 {cluster_feature_means.loc[cluster_id, "Log_Assets"]:.2f}, 负债率 {cluster_feature_means.loc[cluster_id, "Debt_Ratio"]:.2f}') # 报告当次簇的行业众数和两项核心尺度Cluster 0: 代表行业 [机械], 平均资产 21.73, 负债率 0.43
Cluster 1: 代表行业 [机械], 平均资产 22.97, 负债率 0.55
Cluster 2: 代表行业 [机械], 平均资产 21.51, 负债率 0.27表 12.3、图 12.5 与 列表 12.3 是解释聚类的共同起点。读者应根据当次运行的特征均值、载荷和代表行业描述各簇，并检查换随机种子、重抽样或更换期间后结论是否稳定。簇编号没有天然商业含义；未看到运行输出时，不应预设“银行”“软件”或“分区清晰”等结论。
12.5.4 层次聚类（Hierarchical Clustering）
层次聚类可以帮助我们理解公司之间的层级相似性，例如构建”同类公司”。
K-Means 直接给出指定数量的互斥簇；层次聚类（Hierarchical Clustering）则记录样本从单个对象逐步合并的嵌套关系，并用树状图展示不同切割高度下的分组。 为了让树状图的叶标签保持可读，我们在下面的代码中随机抽取至多 50 家公司作为演示样本。Ward 最小方差法每一步选择使合并后组内平方和增加最小的簇对；这一准则来自 Ward 的层次分组方案 (Ward 1963年)。 树状图（Dendrogram）的叶节点代表公司，较低高度发生的合并表示公司在所选标准化财务特征上更相似。相似性不等同于竞争关系、并购价值或因果联系；切割高度也仍需研究者结合任务选择。
from scipy.cluster.hierarchy import dendrogram, linkage # 计算Ward链接并展示逐次合并高度
# 为了可视化清晰，我们只随机选取 50 家公司
sample_size = min(50, len(cleaned_financial_data)) # 确保样本量不超过总数据量
sample_financial_data = cleaned_financial_data.sample(sample_size, random_state=42) # 随机抽取样本公司
sample_features_matrix = sample_financial_data[financial_features].values # 提取特征矩阵
# 标准化：如果 feature_scaler 已经 fit 过则用 transform，否则重新 fit
if not hasattr(feature_scaler, 'mean_'): # 检查 scaler 是否已被拟合
feature_scaler.fit(cleaned_financial_data[financial_features]) # 对全量数据拟合标准化器
scaled_sample_features = feature_scaler.transform(sample_features_matrix) # 对样本数据进行标准化变换
# 使用Ward最小方差法计算层次聚类的链接矩阵
hierarchical_linkage_matrix = linkage(scaled_sample_features, method='ward') # Ward法使组内方差增加最小
plt.figure(figsize=(12, 6)) # 为五十个叶标签保留横向显示空间
dendrogram(hierarchical_linkage_matrix, labels=sample_financial_data['industry_name'].values, leaf_rotation=90, leaf_font_size=8) # 绘制树状图，叶节点标注行业名
plt.title('上市公司财务相似性树状图 (部分样本)', fontsize=14, fontproperties='Source Han Serif SC') # 标明树高对应Ward合并代价
plt.ylabel('距离', fontsize=12, fontproperties='Source Han Serif SC') # Y轴表示合并距离
plt.tight_layout() # 自动调整布局防止标签溢出
plt.show() # 渲染财务特征Ward链接下的公司合并顺序/home/ubuntu/miniconda3/envs/peter/lib/python3.10/site-packages/sklearn/utils/validation.py:2749: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names
warnings.warn(
图 12.6 的叶节点是当次随机抽取的公司，因此不预言银行、电力等具体行业的聚合顺序。读图时可记录低高度合并的公司对，再检查这种结构在更换抽样后是否稳定。与 K-Means 预先指定 \(K\) 不同，树状图可在不同高度切割成不同粒度的探索性分群，但切割高度仍是分析选择。
12.5.5 聚类方法比较
K-Means 适合大型数据集，计算效率高，但要求簇呈球状且大小相近。层次聚类适合小型数据集及探索层级结构，但计算量大。在A股全市场分析中，K-Means 更为常用；而在寻找特定几家公司的对标时，层次聚类更有直观价值。
除了轮廓系数外，Calinski-Harabasz 指数（CH 指数）也是评估聚类质量的重要指标。该指数衡量的是簇间离散度与簇内离散度的比值，数值越大表示聚类效果越好。下面使用 CH 指数对前面的 K-Means 聚类结果进行量化评估。
# 使用Calinski-Harabasz指数量化评估K-Means聚类效果
from sklearn.metrics import calinski_harabasz_score # CH指数：簇间方差/簇内方差，越大越好
ch_score_value = calinski_harabasz_score(clustering_features, cluster_assignments) # 计算CH得分
print(f'K-Means CH Score: {ch_score_value:.2f}') # 报告组间离散相对组内离散的辅助指标K-Means CH Score: 706.51代码会报告当次数据与所选 \(K\) 对应的 CH 指数。CH 指数没有跨数据集通用的绝对阈值；在同一数据和预处理下，可将它与不同 \(K\) 或聚类算法横向比较，但仍需结合稳定性与业务解释。
聚类算法选择指南
- K-Means: 数据量大，寻找明确分群。
- 层次聚类: 数据量小，寻找层级关系。
- DBSCAN: 寻找任意形状的簇，处理噪声（适合地理位置聚类等）。
12.6 无监督学习的挑战与局限
12.6.1 主成分分析的局限
- 线性假设：PCA 只能捕获线性关系，无法处理非线性流形。
- 解释困难：主成分是变量的线性组合，物理意义往往不直观。
- 方差非信息：方差最大的方向不一定包含最有用的信息（例如信噪比问题）。
12.6.2 聚类分析的局限
- 结果主观：K值的选择往往没有标准答案。
- 不仅是数据：聚类结果需要业务解释，数学上分离好的簇在业务上可能没意义。
- 稳定性：个别数据的变动可能导致聚类结果剧烈变化。
12.7 中国案例：基于股价波动的股票聚类
除了基于财务指标，我们还可以根据股价的历史波动形态对股票进行聚类。这种方法常用于构建投资组合（如分散化投资）或配对交易（Pairs Trading）。
最后把无监督学习从静态财务横截面扩展到动态市场数据。在投资组合分析中，收益相关性可以描述证券在当前样本窗口内的共同波动，但它会随市场状态变化，不能仅凭历史低相关就保证未来分散化效果。 这里从 5 个预先列出的行业各取若干股票，并使用 2023-01-01 之后、截至数据快照末日的收盘价计算日收益率。实际窗口长度由现场日期范围决定，不称为“过去一年”或“龙头股总体”。
第 3 步用 1 - 相关系数 定义相关距离。相关系数接近 1 时距离接近 0，接近 -1 时距离接近 2。随后用层次聚类和相关热图展示当前样本窗口中的共同波动结构；它们不能单独证明资金流向、风格轮动或未来分散化效果。
# 1. 加载股价数据与基本信息数据
import os # 读取跨平台数据根环境变量
from pathlib import Path # 使用跨平台路径对象解析显式数据根
book_data_dir_value = os.environ.get('BOOK_DATA_DIR') # 安全读取本案例独立入口的数据根
assert book_data_dir_value, '请先设置 BOOK_DATA_DIR，使其指向包含 stock/ 子目录的数据根' # 缺失时说明修复方法
DATA_DIR = Path(book_data_dir_value).expanduser().resolve() # 将配置解析为统一绝对路径
assert DATA_DIR.is_dir(), f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}' # 在读取前验证数据根
path_price = DATA_DIR / 'stock' / 'stock_price_post_adjusted.h5' # 后复权股价文件路径
path_basic = DATA_DIR / 'stock' / 'stock_basic_data.h5' # 股票基本信息文件路径
assert path_price.is_file(), f'缺少后复权行情文件: {path_price}' # 空目录时指出价格输入文件
assert path_basic.is_file(), f'缺少公司基本信息文件: {path_basic}' # 空目录时指出行业输入文件
# 选取5个预先指定行业，每个行业按股票代码排序后取6只教学样本
stock_basic_data = pd.read_hdf( # 基本信息表仅读取股票池构造字段
path_basic, columns=['order_book_id', 'citics_2019_l1_name'] # 排除本案例不使用的公司属性列
).copy() # 隔离行业股票池映射
target_stock_universe = [] # 初始化目标股票池
target_industry_list = ['银行', '电子', '食品饮料', '非银行金融', '房地产'] # 5个代表性行业（中信一级分类）
for industry_name in target_industry_list: # 遍历每个目标行业
industry_top_stocks = stock_basic_data.loc[stock_basic_data['citics_2019_l1_name'].eq(industry_name), 'order_book_id'].sort_values().head(6).tolist() # 按代码排序后取前6只教学样本
target_stock_universe.extend(industry_top_stocks) # 加入目标股票池
assert target_stock_universe, '预设行业在公司基本信息中未形成可读股票池' # 防止空股票池进入行情读取
price_frames = [pd.read_hdf(path_price, where=f'order_book_id={stock_id!r}', columns=['close']).reset_index() for stock_id in target_stock_universe] # 在HDF存储层逐只读取同一预设股票池
stock_price_history = pd.concat(price_frames, ignore_index=True) # 合并选择性读取的目标证券行情# 2. 数据准备：构建收益率矩阵（行=日期，列=股票）
analysis_start_date = '2023-01-01' # 固定分析起点；结束日由数据快照决定
subset_price_data = stock_price_history[stock_price_history['date'] > analysis_start_date].copy() # 在已下推公司筛选的行情中限定日期窗口
price_pivot_table = subset_price_data.pivot(index='date', columns='order_book_id', values='close') # 转为透视表：行=日期，列=股票代码
daily_returns_matrix = price_pivot_table.pct_change().dropna() # 计算日收益率并删除首行缺失值
print(f'收益率矩阵维度: {daily_returns_matrix.shape}') # 报告相关距离实际使用的交易日数与证券数收益率矩阵维度: (726, 30)上方输出动态报告当前收益率矩阵的交易日数和股票数。后续相关距离矩阵的维度由实际列数决定；缺失值筛选或数据版本变化时，不在正文保留旧的固定维度。
接下来基于收益率矩阵计算相关距离 \(1-\rho\)，并使用平均链接绘制树状图。Ward 准则依赖欧氏平方距离，不适用于这里预先计算的相关距离。
# 3. 计算相关性距离矩阵：1-相关系数，相关性越高距离越近
nonconstant_columns = daily_returns_matrix.std().loc[lambda values: values.gt(0)].index
return_correlation_matrix = daily_returns_matrix[nonconstant_columns].corr(min_periods=60) # 只在至少六十个共同交易日上估计相关性
return_correlation_matrix = return_correlation_matrix.dropna(axis=0, how='any').dropna(axis=1, how='any') # 排除无法形成完整距离的证券
return_correlation_matrix = return_correlation_matrix.loc[ # 强制行列证券顺序完全一致
return_correlation_matrix.index, return_correlation_matrix.index # 用同一索引排列方阵两轴
] # 完成对称相关矩阵的顺序校准
correlation_distance_matrix = (1 - return_correlation_matrix).clip(0, 2) # 把相关性映射到合法的非负差异尺度
correlation_distance_matrix = (correlation_distance_matrix + correlation_distance_matrix.T) / 2 # 消除浮点误差造成的轻微非对称
np.fill_diagonal(correlation_distance_matrix.values, 0.0) # 同一证券到自身的距离固定为零
assert np.isfinite(correlation_distance_matrix.to_numpy()).all() # 在压缩距离前排除非有限元素
# 4. 对相关性距离使用平均链接；Ward 只适用于欧氏特征空间
from scipy.cluster.hierarchy import dendrogram, linkage # 对收益相关距离执行平均链接并绘制树状图
from scipy.spatial.distance import squareform # 将对称距离方阵压缩为层次聚类所需向量
condensed_correlation_distance = squareform(correlation_distance_matrix, checks=True) # 同时核验对称性与零对角线
hierarchical_linkage_matrix = linkage(condensed_correlation_distance, method='average') # 对非欧氏相关距离采用平均链接
# 绘制层次聚类树状图
plt.figure(figsize=(12, 6)) # 为证券代码叶标签预留横向空间
dendrogram(hierarchical_linkage_matrix, labels=return_correlation_matrix.columns, leaf_rotation=90, leaf_font_size=8) # 用证券代码标记共同链接矩阵的叶节点
plt.title('股票收益相关性聚类（平均链接）', fontsize=14, fontproperties='Source Han Serif SC') # 标明树高对应平均相关距离
plt.ylabel('距离', fontsize=12, fontproperties='Source Han Serif SC') # Y轴表示合并距离
plt.tight_layout() # 防止标签溢出
plt.show() # 渲染收益相关距离平均链接下的证券合并顺序
图 12.7 只描述所选窗口中的收益共动结构。可对照行业标签检查低距离合并是否更常发生在行业内，但树状图本身既不识别因果因子，也不证明分散效果或可交易关系。
最后，我们用热图（Clustermap）显示经过聚类排序后的股票相关性矩阵，以更直观地揭示行业内外的相关性模式。
import seaborn as sns # 复用既定链接矩阵排列收益相关热图
# 显式复用树状图的 linkage，避免 seaborn 默认重算另一套聚类
clustermap_figure = sns.clustermap(
return_correlation_matrix, cmap='coolwarm', center=0, figsize=(10, 10), # 用零为中心显示正负相关
row_linkage=hierarchical_linkage_matrix, col_linkage=hierarchical_linkage_matrix # 行列复用同一链接以保持排序一致
) # 完成相关矩阵聚类热图
clustermap_figure.fig.suptitle('股票相关性热图(聚类排序)', fontsize=16, fontproperties='Source Han Serif SC', y=1.02) # 设置总标题
plt.show() # 渲染与树状图复用同一链接矩阵的相关性排序
print('\n观察结果：') # 输出分析提示
print('聚类是否成功将同一行业的股票（如所有银行股）归为一类？') # 提示学生验证假设
print('行业内聚集若出现，只是该窗口中的描述性共动证据。') # 限定热图不能直接支持因果或交易结论
观察结果：
聚类是否成功将同一行业的股票（如所有银行股）归为一类？
行业内聚集若出现，只是该窗口中的描述性共动证据。图 12.8 与树状图使用同一距离和 linkage，因此排序应一致。是否出现块结构由当前图形决定；相关性会随窗口变化，不能单凭热图认定行业是“主要因子”或某组合已经有效分散。
12.7.1 应用场景
- 组合研究假设：距离较远的股票可能提供不同的短期共动结构，但仍需样本外协方差、权重与交易成本检验。
- 配对研究假设：距离近只表示相关，不等于价差平稳；配对交易还需协整、结构稳定性和样本外成本检验。
12.8 本章小结
本章我们系统学习了两种主要的无监督学习方法：
12.8.1 主成分分析（PCA）
核心思想：找到数据方差最大的方向，实现降维
关键步骤：
- 所有 PCA 都先中心化；仅在尺度问题需要时再做单位方差标准化
- 计算主成分
- 选择主成分数量（碎石图、累积方差）
- 解释主成分（载荷矩阵）
- 计算主成分得分
应用场景：
- 数据可视化
- 降维以减少计算量
- 去除多重共线性
- 特征提取
12.8.2 聚类分析
核心思想：将相似观测归为一组
主要算法：
- K-Means：快速、简单，适合大数据
- 层次聚类：可视化层次结构，适合小数据
评估指标：
- 轮廓系数（Silhouette Score）
- Calinski-Harabasz 指数
- Davies-Bouldin 指数
应用场景：
- 客户细分
- 图像分割
- 文档聚类
- 异常检测
无监督学习最佳实践
- 数据预处理至关重要：标准化、处理缺失值、异常值
- 选择合适的方法：根据数据特点和业务目标
- 验证结果：虽然难以验证，但应与业务知识结合
- 迭代优化：尝试不同参数和方法，比较结果
- 与监督学习结合：无监督学习可用于特征工程，提升监督学习效果
12.9 理论来源与前沿
无监督学习的“理论母体”主要来自三条线索：
- 线性代数与谱理论：PCA 将高维数据投影到低维子空间，本质是对协方差矩阵进行特征分解；谱聚类进一步把图拉普拉斯矩阵的特征向量用于发现群落结构。
- 优化与几何：K-Means 的目标函数是组内平方和（within-cluster sum of squares），其交替最小化过程可解释为一个EM 风格的坐标下降算法；层次聚类则对应一类凝聚式或分裂式的贪心优化。
- 概率建模：高斯混合模型把聚类看作潜在类别变量的后验推断问题，天然支持软聚类与不确定性量化。
近年来的前沿发展集中在：
- 非线性降维：t-SNE、UMAP 等方法在保持局部邻域结构方面效果突出，但需要谨慎解读全局距离。
- 可扩展与鲁棒聚类：在大规模数据上，MiniBatch K-Means、近似最近邻与基于密度的方法（如 HDBSCAN）更实用。
- 与业务决策联动：在客户细分、风险分层等场景中，无监督结果往往需要与后续监督学习或策略优化结合，并通过稳定性分析与可解释性工具进行验证。
12.10 M12：无监督结构证据
本里程碑直接复用本章已经执行的真实数据对象、表格和图形。学生在同一个 QMD/IPYNB 内按顺序运行代码，即可检查数据、修改参数并立即看到结果。
学生沿四步统计路径完成任务：
- 从 表 12.1 核对报告期、公司范围与清洗后的分析样本；
- 根据 表 12.2 与 图 12.1 选择 PCA 维数，并结合载荷解释主成分；
- 比较满足 \(n\ge5K\) 的候选 \(K\)，用轮廓系数和不同随机种子的 ARI 检查稳定性；
- 结合 表 12.3 与 图 12.5 解释原尺度簇画像，同时说明二维投影、簇编号和样本外推边界。
selected_components=d 只定义聚类使用 PC1 至 PC\(d\) 的输入维数，不定义展示维数。若 \(d=1\)，聚类只能使用 PC1；PC2 即使用于二维展示，也不得反向进入聚类。若没有候选 \(K\) 满足样本量条件、出现空簇、轮廓系数不可定义或稳定性不足，应直接在 notebook 输出中报告停止原因，不能放宽规则或补造结果。
学生量规（20 分）：总体与清洗核对 4 分；标准化理由、碎石图和 PCA 维数选择 5 分；载荷与得分解释 4 分；\(K\) 选择及稳定性 4 分；原尺度画像和结论边界 3 分。
12.11 练习
12.11.1 概念题
[核心｜难度:2｜分值:5｜任务:独立] 解释主成分分析的几何意义。为什么第一主成分的方向对应于数据方差最大的方向？
[核心｜难度:1｜分值:4｜任务:独立] 区分中心化、协方差 PCA 与相关矩阵 PCA。什么情形应再做单位方差标准化？尺度选择会怎样改变载荷？
[核心｜难度:2｜分值:10｜任务:独立] 对 \(A=(0,0)\)、\(B=(0,2)\)、\(C=(4,0)\)、\(D=(6,0)\) 做 \(K=2\) 的 K-Means，初始质心为 \(\mu_1^{(0)}=(0,0)\)、\(\mu_2^{(0)}=(6,0)\)；距离相同时分到编号较小的簇。完整写出两轮的点—簇分配、更新后质心和 WCSS，判定是否停止，并用一句话解释 K-Means 目标函数。
[核心｜难度:2｜分值:5｜任务:独立] 比较层次聚类和K-Means 聚类的优缺点。
[核心｜难度:2｜分值:5｜任务:独立] 解释轮廓系数的含义。它的取值范围是什么？如何解释？
12.11.2 应用题
[核心｜难度:3｜分值:20｜项目:M12] 在同一 QMD/IPYNB 中重新运行 小节 12.4.3 至 小节 12.10 的真实数据代码。根据现场输出核对样本与清洗规则，说明标准化口径；用累计解释率选择 PCA 维数并区分载荷与得分；比较满足 \(n\ge5K\) 的候选 \(K\)，用不同随机种子的 ARI 检查稳定性；最后结合原尺度画像和二维得分投影写一页以内的结论。不得把得分图称为载荷图，也不得把当前样本中的簇解释为真实企业类型或因果关系。
[拓展｜难度:3｜分值:14｜任务:独立] 受控模拟：K-Means 的几何边界。生成三个方差和重叠程度可控的二维高斯簇：
- 比较标准化前后结果
- 使用 K-Means，并用轮廓系数选择 \(K\)
- 改变簇方差与随机种子，报告标签稳定性
- 只解释算法几何性质，不把人工簇称为真实群体或业务分层
[拓展｜难度:3｜分值:15｜任务:独立] 以长三角四省市中每家上市公司为分析单位，使用同一预先固定报告期的资产规模、收入规模、净利率与 ROA 描述财务画像（不是省域经济水平）：
- 说明公司范围、报告期与逐公司排除原因（3 分）
- 在同一标准化公司特征矩阵上拟合 Ward 层次聚类并绘制公司级树状图（4 分）
- 只比较满足 \(n\ge5K\) 的 KMeans 候选，以轮廓系数选择 \(K\)，再按同一个 \(K\) 切割 Ward 树（4 分）
- 报告 Ward 与 KMeans 的 ARI、每簇公司数和省份构成，并说明省份只用于解释而不参与选型（2 分）
- 若无可行 \(K\)，报告样本量 \(n\)、候选规则和停止原因即可获得 KMeans 比较的 4 分，不得强行拟合（2 分）
12.11.3 理论题
[核心｜难度:3｜分值:10｜任务:独立] 证明 PCA 的主成分方向是数据协方差矩阵的特征向量。
[拓展｜难度:3｜分值:12｜任务:独立] 推导 K-Means 算法的EM（期望最大化）解释。
[拓展｜难度:3｜分值:12｜任务:独立] 证明在正态分布假设下，Ward 层次聚类方法等价于最小化合并后的方差增加量。
[拓展｜难度:2｜分值:10｜任务:独立] 研究并总结以下主题：
- t-SNE 和UMAP（非线性降维方法）
- DBSCAN 和谱聚类（基于密度的聚类）
- 高斯混合模型（基于概率的聚类）
12.12 练习参考解答
展开完整解答、评分点与常见失败模式
评分以题面分值为准；代数与 WCSS 允许 \(10^{-3}\) 绝对误差。应用题若缺碎石图、载荷含义、\(K\) 选择或多初值稳定性，每项扣 20%；用行业标签选择 \(K\)、把簇当真实类型或把省份上市公司画像写成省域经济结论，解释项不得分。
12.12.1 概念题解答
PCA 的几何意义
PCA 的本质是在高维空间中寻找一组新的正交坐标轴（即主成分方向），使得数据在这些新轴上的投影方差依次递减。第一主成分的方向是数据方差最大的方向，因为最大化投影方差等价于最小化数据点到该方向的正交投影距离之和（即最小化重构误差）。从数学上看，这等价于求解数据协方差矩阵 \(\Sigma\) 的最大特征值对应的特征向量。几何上，第一主成分就是数据”椭球”的最长轴方向。
中心化与尺度选择
中心化，即减去各变量均值，是协方差 PCA 的必要步骤；否则投影平方和混入均值位置，不能直接解释为方差。是否再除以标准差不是定义要求。只中心化等价于协方差矩阵 PCA，保留原单位下的方差权重；中心化并除以标准差等价于相关矩阵 PCA，使每个非恒定变量以单位方差进入分析。若变量量纲不同且研究问题关注相对画像，例如同时使用对数资产和财务比率，应选择相关矩阵 PCA；若变量同量纲且绝对波动大小正是研究对象，则协方差 PCA 可能更合适。标准化改变相对方差权重，因此载荷方向、得分和下游聚类都可能改变，不能把两种答案视为机械等价。
K-Means 两轮手算
K-Means 的目标函数是组内平方和（Within-Cluster Sum of Squares, WCSS）： \[ \min_{C_1,...,C_K} \sum_{k=1}^{K} \sum_{x_i \in C_k} \|x_i - \mu_k\|^2 \] 第 1 轮分配为 \(C_1^{(1)}=\{A,B\}\)、\(C_2^{(1)}=\{C,D\}\)，更新得 \(\mu_1^{(1)}=(0,1)\)、\(\mu_2^{(1)}=(5,0)\)。因此 \(\mathrm{WCSS}^{(1)}=(1+1)+(1+1)=4\)。第 2 轮的分配仍为 \(\{A,B\}\) 和 \(\{C,D\}\)，质心仍为 \((0,1)\) 和 \((5,0)\)，\(\mathrm{WCSS}^{(2)}=4\)；分配与质心均不再变化，故停止。目标函数最小化每个点到所属簇均值中心的平方距离和，表示簇内紧凑度。
评分量规（10 分）：目标函数及含义 2 分；第 1 轮分配 2 分、两个更新质心 2 分、WCSS 1 分；第 2 轮分配、质心与 WCSS 共 2 分；停止理由 1 分。
层次聚类 vs. K-Means 的优缺点
特征 K-Means 层次聚类 簇数 \(K\) 需预先指定 无需预先指定，树状图可后验选择 计算复杂度 \(O(nKT)\)，较快 \(O(n^2 \log n)\) 或 \(O(n^3)\)，较慢 簇形状 倾向于发现球形簇 可发现不规则形状簇（视联接方式） 确定性 非确定性（依赖初始化） 确定性 可解释性 中等 高（树状图直观展示层级关系） 大数据适用性 适合大规模数据 不适合大规模数据 轮廓系数的含义
轮廓系数衡量数据点与自身簇的紧密度相对于与最近邻簇的分离度。对于数据点 \(i\)： \[ s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))} \] 其中 \(a(i)\) 是点 \(i\) 到同簇其他点的平均距离（紧密度），\(b(i)\) 是点 \(i\) 到最近邻簇所有点的平均距离（分离度）。取值范围为 \([-1, 1]\)：\(s(i) \approx 1\) 表示聚类效果好（簇内紧密、簇间分离），\(s(i) \approx 0\) 表示点在两个簇的边界上，\(s(i) < 0\) 表示可能被错误分配。
12.12.2 应用题解答
同一文档中的财务 PCA 与聚类
先读取 表 12.1 的报告期、候选公司数和清洗后样本数，再用 表 12.2 与 图 12.1 报告累计解释率达到阈值时的最小维数。主成分命名必须引用载荷绝对值较大的变量及符号；载荷接近时写“混合维度”。
随后从 notebook 的候选表读取所有满足 \(n\ge5K\) 的 \(K\) 及轮廓系数，并报告不同随机种子所得 ARI 的最小值、均值和离散程度。图 12.5 只是在 PC1/PC2 平面展示同一次拟合的簇标签；聚类实际使用的维数仍由 \(d\) 决定。最后结合 表 12.3 概括当前报告期内的原尺度财务画像，并明确簇编号可置换、二维距离会丢失高维信息、结论不能外推为稳定企业类型或因果关系。
受控模拟：K-Means 的几何边界
本题只演示三个人工高斯簇的几何性质，不声称这些坐标来自真实业务对象或支持实际分层结论。
import pandas as pd # 组织受控三簇机制样本与重复实验结果
import numpy as np # 生成固定高斯样本并汇总ARI分布
from sklearn.cluster import KMeans # 在候选K与不同初值下生成算法簇标签
from sklearn.metrics import silhouette_score # 在预设候选范围比较样本内分离度
from sklearn.preprocessing import StandardScaler # 本题独立比较是否标准化
import matplotlib.pyplot as plt # 展示候选K证据与受控机制样本的簇标签
# 1. 构造仅用于观察 K-Means 几何边界的二维机制样本
np.random.seed(42) # 设置随机数种子保证可复现性
number_of_samples = 500 # 设定机制样本总数
base_group_size = number_of_samples // 3 # 前两组使用相同整数样本量
# 三组坐标只编码不同中心和离散程度，不对应任何真实业务群体
simulated_cluster_points = np.concatenate([ # 合并三组受控二维高斯点
np.random.normal([20, 20], [5, 5], size=(base_group_size, 2)), # 生成离散度较小的机制组A
np.random.normal([60, 50], [10, 10], size=(base_group_size, 2)), # 生成中等离散的机制组B
np.random.normal([90, 80], [10, 10], size=(number_of_samples - 2 * base_group_size, 2)) # 让机制组C吸收整数除法余数
]) # 形成只用于比较聚类几何的二维数组
mechanism_data_frame = pd.DataFrame(simulated_cluster_points, columns=['x1', 'x2']) # 使用中性坐标名避免业务实体暗示
assert len(mechanism_data_frame) == number_of_samples # 核对生成行数与预先声明的样本总数一致机制样本构建后，利用轮廓系数在 $K=2$ 到 $K=9$ 的候选中比较样本内分离度，并可视化算法分组。该结果不代表任何真实业务分层。# 2. 寻找最优K：遍历K=2到9计算轮廓系数
silhouette_scores = [] # 存储各K对应的轮廓系数
cluster_count_range = range(2, 10) # 候选聚类数范围
for k in cluster_count_range: # 遍历每个候选K值
kmeans_model = KMeans(n_clusters=k, n_init=10, random_state=42) # 初始化K-Means
kmeans_model.fit(mechanism_data_frame) # 在同一机制坐标上拟合候选模型
silhouette_scores.append(silhouette_score(mechanism_data_frame, kmeans_model.labels_)) # 记录候选的样本内分离度
plt.plot(cluster_count_range, silhouette_scores, 'o-') # 绘制轮廓系数曲线
plt.title('Silhouette Score') # 标明曲线用于在预设范围内选择机制实验的K
plt.show() # 峰值位置以本次运行输出为准
图 12.9 预先固定 \(K\) 的选择证据；图 12.10 只展示依该选择得到的二维标签。
# 3. 使用现场轮廓系数最高的 K 进行最终聚类
selected_k = list(cluster_count_range)[int(np.argmax(silhouette_scores))] # 在预设候选内选择轮廓系数最高的K
kmeans_model = KMeans(n_clusters=selected_k, n_init=10, random_state=42) # 初始化最终模型
cluster_assignments = kmeans_model.fit_predict(mechanism_data_frame) # 对机制样本生成算法标签
plt.scatter(mechanism_data_frame['x1'], mechanism_data_frame['x2'], c=cluster_assignments) # 按算法标签显示二维机制点
plt.xlabel('机制坐标 x1') # 横轴仅表示受控模拟的第一坐标
plt.ylabel('机制坐标 x2') # 纵轴仅表示受控模拟的第二坐标
plt.title('受控高斯簇的 K-Means 标签') # 明示图形不代表真实业务细分
plt.show() # 渲染受控高斯机制样本的现场K-Means标签
固定一次种子不能完成“稳定性”要求。下面在两种重叠程度、是否标准化和 10 个随机种子上复算，并用调整兰德指数（ARI）比较各运行与同场景参考标签；标签编号置换不会影响 ARI。
from sklearn.metrics import adjusted_rand_score # 使用不受标签编号影响的稳定性指标
stability_records = [] # 保存每个受控场景的比较结果
for overlap_scale in [0.6, 1.4]: # 改变簇方差以控制重叠程度
scenario_rng = np.random.default_rng(20260812) # 固定数据生成随机源
scenario_points = np.vstack([scenario_rng.normal(center, overlap_scale, size=(120, 2)) for center in [[0, 0], [4, 0], [2, 3.5]]]) # 生成三簇
for should_standardize in [False, True]: # 比较原尺度与标准化尺度
scenario_input = StandardScaler().fit_transform(scenario_points) if should_standardize else scenario_points # 按场景处理尺度
reference_labels = KMeans(n_clusters=3, n_init=50, random_state=0).fit_predict(scenario_input) # 建立多初值参考
for initialization_seed in range(10): # 改变初始化随机种子
repeated_labels = KMeans(n_clusters=3, n_init=1, random_state=initialization_seed).fit_predict(scenario_input) # 刻意单初值测试敏感性
stability_records.append({'overlap': overlap_scale, 'standardized': should_standardize, 'seed': initialization_seed, 'ari': adjusted_rand_score(reference_labels, repeated_labels)}) # 保存 ARI
stability_table = pd.DataFrame(stability_records) # 汇总全部受控运行
print(stability_table.groupby(['overlap', 'standardized'])['ari'].agg(['mean', 'min', 'std'])) # 报告平均、最差与离散度 mean min std
overlap standardized
0.6 False 1.000000 1.000000 0.000000
True 1.000000 1.000000 0.000000
1.4 False 0.967808 0.925961 0.035083
True 0.977926 0.957632 0.020060交付必须包含四个场景的均值、最小值与标准差；若最小 ARI 低于 0.8，结论应写“对初始化不稳定”，而不是只展示最好的一次运行。
- 基于本地数据的长三角上市公司财务画像聚类
本题的数据身份、2022Q4 版本选择、长三角公司范围、比率构造与顺序缩尾以 小节 12.4.3 为唯一规范实现；具体清洗输出见 表 12.1，尺度选择见 列表 12.1。解答不另建第二条 HDF 读取路径，而是复用其中的 merged_financial_data 与 cleaned_financial_data，再补充题目特有的正收入要求、逐公司排除账本和 Ward/K-Means 比较。
m12_ex8_audit = merged_financial_data[['order_book_id', 'province', 'total_assets', 'operating_revenue', 'net_profit']].copy() # 从规范横截面建立逐公司账本
m12_ex8_amount_columns = ['total_assets', 'operating_revenue', 'net_profit'] # 声明画像构造所需原始金额
m12_ex8_is_finite = np.isfinite(m12_ex8_audit[m12_ex8_amount_columns]).all(axis=1) # 识别缺失或非有限金额
m12_ex8_has_denominators = m12_ex8_audit['total_assets'].gt(0) & m12_ex8_audit['operating_revenue'].gt(0) # 对数和比率要求正资产与正收入
m12_ex8_sample = cleaned_financial_data.loc[cleaned_financial_data['Revenue'].gt(0)].copy() # 复用规范清洗并施加题目正收入约束
m12_ex8_sample['Log_Revenue'] = np.log(m12_ex8_sample['Revenue']) # 构造收入规模指标
m12_ex8_included_ids = set(m12_ex8_sample['order_book_id']) # 锁定最终公司级分析样本
m12_ex8_audit['exclusion_reason'] = np.select( # 为每家公司登记首要排除原因
[~m12_ex8_is_finite, ~m12_ex8_has_denominators, ~m12_ex8_audit['order_book_id'].isin(m12_ex8_included_ids)], # 按数据有效性、分母和尾部规则排序
['missing_or_nonfinite_amount', 'nonpositive_asset_or_revenue', 'outside_sequential_tail_filter'], # 使用可复核原因代码
default='included' # 其余公司进入四指标聚类
) # 完成互斥的公司级审计状态
m12_ex8_exclusion_summary = m12_ex8_audit['exclusion_reason'].value_counts().rename_axis('reason').rename('companies') # 汇总各原因公司数
print({'report_period': target_report_period, 'candidate_companies': len(m12_ex8_audit), 'included_companies': len(m12_ex8_sample)}) # 报告题目样本身份
print(m12_ex8_exclusion_summary) # 输出排除原因计数供评分
m12_ex8_audit[['order_book_id', 'province', 'exclusion_reason']] # 保留可逐公司检查的完整账本对象{'report_period': '2022q4', 'candidate_companies': 2014, 'included_companies': 1814}
reason
included 1814
outside_sequential_tail_filter 198
nonpositive_asset_or_revenue 2
Name: companies, dtype: int64| order_book_id | province | exclusion_reason | |
|---|---|---|---|
| 2 | 300356.XSHE | 江苏省 | outside_sequential_tail_filter |
| 3 | 300336.XSHE | 上海市 | outside_sequential_tail_filter |
| 4 | 300330.XSHE | 上海市 | included |
| 7 | 600532.XSHG | 上海市 | included |
| 22 | 000918.XSHE | 浙江省 | outside_sequential_tail_filter |
| ... | ... | ... | ... |
| 5334 | 002605.XSHE | 上海市 | included |
| 5337 | 002883.XSHE | 江苏省 | included |
| 5341 | 688076.XSHG | 江苏省 | included |
| 5342 | 688022.XSHG | 江苏省 | included |
| 5345 | 002514.XSHE | 江苏省 | included |
2014 rows × 3 columns
下面在同一四指标标准化矩阵上拟合两种算法。候选 \(K\) 只由 \(n\ge5K\) 与轮廓系数决定，省份不进入标准化、选型或拟合。
from scipy.cluster.hierarchy import fcluster # 按选定K切割同一Ward树
m12_ex8_features = ['Log_Assets', 'Log_Revenue', 'Net_Margin', 'ROA'] # 锁定题面四项公司特征
m12_ex8_matrix = StandardScaler().fit_transform(m12_ex8_sample[m12_ex8_features]) # 在同一公司样本统一四项尺度
m12_ex8_candidates = [cluster_count for cluster_count in range(2, 9) if len(m12_ex8_sample) >= 5 * cluster_count] # 只保留满足样本量规则的K
assert m12_ex8_candidates, f'n={len(m12_ex8_sample)} 时没有满足 n>=5K 的候选；停止算法比较' # 无可行候选时诚实停止
m12_ex8_silhouettes = {} # 保存每个候选的样本内分离度
for cluster_count in m12_ex8_candidates: # 对共享输入逐一比较候选K
candidate_labels = KMeans(n_clusters=cluster_count, n_init=50, random_state=0).fit_predict(m12_ex8_matrix) # 使用相同多初值预算
m12_ex8_silhouettes[cluster_count] = silhouette_score(m12_ex8_matrix, candidate_labels) # 计算候选轮廓系数
m12_ex8_selected_k = max(m12_ex8_silhouettes, key=m12_ex8_silhouettes.get) # 依事前规则选择轮廓系数最高的K
m12_ex8_kmeans_labels = KMeans(n_clusters=m12_ex8_selected_k, n_init=50, random_state=0).fit_predict(m12_ex8_matrix) # 拟合选定K的K-Means
m12_ex8_ward_linkage = linkage(m12_ex8_matrix, method='ward') # 在同一欧氏标准化矩阵计算Ward链接
m12_ex8_ward_labels = fcluster(m12_ex8_ward_linkage, t=m12_ex8_selected_k, criterion='maxclust') - 1 # 把Ward树切成相同K组
m12_ex8_algorithm_ari = adjusted_rand_score(m12_ex8_kmeans_labels, m12_ex8_ward_labels) # 用置换不变ARI比较算法标签
m12_ex8_sample['kmeans_cluster'] = m12_ex8_kmeans_labels # 将K-Means标签对齐回公司
m12_ex8_sample['ward_cluster'] = m12_ex8_ward_labels # 将Ward标签对齐回同一公司列表 12.5 保持样本、特征、尺度和 \(K\) 完全一致，因而 ARI 的差异只反映两种聚类准则在当前矩阵上的分组差别。
plt.figure(figsize=(14, 7)) # 为公司级Ward树提供横向画布
dendrogram(m12_ex8_ward_linkage, no_labels=True, color_threshold=None) # 展示全部公司叶节点并隐藏拥挤代码文本
plt.xlabel('公司叶节点（股票代码保存在公司级审计对象中）', fontproperties='Source Han Serif SC') # 说明叶节点单位与标签位置
plt.ylabel('Ward 合并距离', fontproperties='Source Han Serif SC') # 标明纵轴为组内平方和增量对应距离
plt.title('长三角上市公司财务画像的 Ward 层次结构', fontproperties='Source Han Serif SC') # 限定图形为当期公司画像
plt.tight_layout() # 避免轴标题超出画布
plt.show() # 渲染同一标准化矩阵的公司级树
图 12.11 展示全部公司参与的合并结构；隐藏叶标签只是为了避免数百个代码重叠，m12_ex8_sample 保留逐公司标签。题面要求的候选、算法一致性、簇规模与省份构成由 表 12.6 统一输出。
m12_ex8_candidate_table = pd.Series(m12_ex8_silhouettes, name='silhouette').rename_axis('k').to_frame() # 整理所有可行候选的轮廓证据
m12_ex8_cluster_sizes = m12_ex8_sample.groupby(['kmeans_cluster', 'ward_cluster']).size().rename('companies') # 报告两种标签组合的公司数
m12_ex8_province_composition = pd.crosstab(m12_ex8_sample['kmeans_cluster'], m12_ex8_sample['province'], normalize='index') # 描述K-Means簇内省份比例
print(m12_ex8_candidate_table.round(3)) # 输出全部候选而非只报最佳结果
print({'selected_k': m12_ex8_selected_k, 'ward_kmeans_ari': round(m12_ex8_algorithm_ari, 3)}) # 输出共享K与算法一致性
print(m12_ex8_cluster_sizes) # 输出两种算法标签的联合簇规模
print(m12_ex8_province_composition.round(3)) # 输出仅用于解释的省份构成 silhouette
k
2 0.331
3 0.367
4 0.302
5 0.287
6 0.289
7 0.291
8 0.274
{'selected_k': 3, 'ward_kmeans_ari': 0.359}
kmeans_cluster ward_cluster
0 0 5
1 82
2 131
1 0 1
2 1024
2 0 295
2 276
Name: companies, dtype: int64
province 上海市 安徽省 江苏省 浙江省
kmeans_cluster
0 0.284 0.092 0.339 0.284
1 0.193 0.073 0.372 0.362
2 0.222 0.119 0.291 0.368诊断顺序是：先确认 表 12.5 中每家公司只有一个状态，再确认候选非空、各标签实际包含 \(K\) 个非空簇，最后解释 ARI。轮廓系数是样本内几何指标，ARI 只衡量 Ward 与 K-Means 的一致程度；二者都不能证明存在真实企业类型。省份只作事后描述，不能据省份构成反选 \(K\)，也不能把上市公司样本推广为省域经济。若 ARI 很低，应报告算法假设敏感；若某簇过小、树高没有清晰断点或更换报告期后标签剧变，应把画像结论降格为不稳定探索结果。
12.12.3 理论题解答
第 12 题详细比较量规
- t-SNE 侧重保留局部邻域，UMAP 兼顾局部与部分全局结构；二者的二维轴和簇间距离通常没有直接业务尺度，超参数与随机种子都会改变图形。
- DBSCAN 以密度定义簇并可识别噪声点，适合非球形结构，但对距离尺度、
eps和维度敏感；谱聚类借助相似度图的特征向量处理非凸结构，但相似度构造和样本规模决定计算成本。 - 高斯混合模型给出软聚类概率，允许不同协方差结构；其分布假设、成分数选择和局部最优需要用信息准则、稳定性与业务可解释性共同检查。
- 满分答案必须按“目标函数或结构假设、关键超参数、优势、失效情形、计算成本、一个中国商业金融应用及结论边界”七列比较，而不能只罗列方法名称。
PCA 主成分方向是协方差矩阵特征向量
最大化 \(v^\top\Sigma v\) 的\(v^\top v=1\) 的拉格朗日函数为 \(\mathcal{L}(v,\lambda)=v^\top\Sigma v-\lambda(v^\top v-1)\)，一阶条件给。 \[ \Sigma v = \lambda v, \] 的\(v\) 为特征向量，\(\lambda\) 为对应特征值。
K-Means 作为球形 GMM 的硬分配极限
设混合模型 \(p(x_i)=\sum_{k=1}^K\pi_k\phi(x_i;\mu_k,\sigma^2I)\)。E 步责任度为
\[ \gamma_{ik}=\frac{\pi_k\exp[-\lVert x_i-\mu_k\rVert^2/(2\sigma^2)]} {\sum_{\ell}\pi_\ell\exp[-\lVert x_i-\mu_\ell\rVert^2/(2\sigma^2)]}. \]
当 \(\sigma^2\to0\) 且最近中心唯一时，距离最小项的指数衰减最慢，故其责任度趋于 1，其余趋于 0。这给出硬分配
\[ z_{ik}=I\!\left(k=\arg\min_\ell\lVert x_i-\mu_\ell\rVert^2\right). \]
固定 \(z_{ik}\) 后，最小化完全数据负对数似然中依赖 \(\mu_k\) 的部分
\[ \sum_i z_{ik}\lVert x_i-\mu_k\rVert^2 \]
对 \(\mu_k\) 求导并令零，得到
\[ \mu_k=\frac{\sum_i z_{ik}x_i}{\sum_i z_{ik}}, \]
即簇内均值更新。两步交替不会增加 WCSS，但只保证到局部最优；有限 \(\sigma^2\)、不同协方差或非等混合权重时，GMM 不等同于 K-Means。
Ward 合并的方差增加量推导
对簇 \(A\) 定义 \(W(A)=\sum_{i\in A}\lVert x_i-\bar x_A\rVert^2\)。合并 \(A,B\) 后，利用平方和分解
\[ \sum_{i\in A}\lVert x_i-\bar x_{A\cup B}\rVert^2 =W(A)+n_A\lVert\bar x_A-\bar x_{A\cup B}\rVert^2, \]
对 \(B\) 同理。又因为
\[ \bar x_{A\cup B}=\frac{n_A\bar x_A+n_B\bar x_B}{n_A+n_B}, \]
可得
\[ \Delta(A,B)=W(A\cup B)-W(A)-W(B) =\frac{n_An_B}{n_A+n_B}\lVert\bar x_A-\bar x_B\rVert^2. \]
Ward 每步选择 \(\Delta(A,B)\) 最小的簇对，因此正是最小化该次合并造成的组内平方和增长。结论依赖欧氏平方距离与均值中心；把相关距离直接传给 Ward 不满足此推导。
主题总结（写作建议）
- t-SNE/UMAP：强调‘局部结构保持’与‘可视化解释限制’，并讨论超参数敏感性。
- DBSCAN/谱聚类：强调密度可发现任意形状簇、谱方法适合图结构与非凸簇。
- GMM：强调软聚类与模型选择（BIC/AIC），以及协方差结构对簇形状的影响。
12.13 章末回顾
学习自检：说明标准化为何改变距离几何；从解释方差规则选择 PCA 维数并解释载荷；区分二维展示坐标与实际聚类输入；按 \(n\ge5K\)、有效簇与有限轮廓筛选候选；用多初值稳定性区分算法波动与结构证据。还应说明多行业总体、单行业分析样本与清洗后聚类样本的包含关系，并列出逐公司排除原因和核心结果。
禁用情境：样本边界或排除理由不清时，不应报告实样本画像；样本量不足、有效簇少于 \(K\) 或稳定性失败时，不应强行命名簇。常见误区是把 PC1/PC2 图上的分离当作全维结构，或把簇编号、行业标签和因果机制相互替代。
无提示检索：1）为何 PCA 前是否标准化会改变答案？2）二维图只解释了哪些信息？3）轮廓系数最高是否足以证明存在真实类别？
展开检索反馈与学习决策
1）标准化改变变量相对尺度，从而改变协方差或相关矩阵的特征方向和距离。2）二维图只展示所画两个主成分承载的投影，不代表未画维度。3）不能；轮廓系数只是候选几何与样本内分离证据，还需样本量、稳定性、外部解释和失败分析。需要复习时，可分别改变变量尺度、旋转展示投影或更换候选 \(K\)，观察结论如何变化。
陌生迁移：换一个公司横截面，先声明总体、排除规则、PCA 维数与 \(K\) 选择规则；若结论改变，报告敏感性而非选择更好看的图。下一章处理同时检验许多探索发现时的错误累积。
Pearson, Karl. 1901年. 《On Lines and Planes of Closest Fit to Systems of Points in Space》. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2 (11): 559～72. https://doi.org/10.1080/14786440109462720.
Ward, Joe H. 1963年. 《Hierarchical Grouping to Optimize an Objective Function》. Journal of the American Statistical Association 58 (301): 236～44. https://doi.org/10.1080/01621459.1963.10500845.