12  无监督学习(Unsupervised Learning)

在前面的章节中，我们主要讨论了监督学习方法，如回归和分类。在监督学习中，我们通常有\(p\) 个特征\(X_1, X_2, \ldots, X_p\) 以及一个响应变量\(Y\)，目标是利用 \(X\) 来预测\(Y\)。

本章将转向无监督学习，这是一类统计工具，用于我们只有特征 \(X_1, X_2, \ldots, X_p\) 但没有响应变量\(Y\) 的情形。无监督学习的目标是发现关于这些测量的有趣模式：

12.1 学习闭环

先修与可观察目标

完成本章后，学生应能：

1. 从中心化矩阵推导 PCA 的特征向量问题，并用碎石图、载荷与得分解释降维结果。
2. 对给定二维点和初始中心手算至少两轮 K-Means，逐轮报告分配、中心和 WCSS，数值无误。
3. 对新数据选择标准化、距离、聚类算法与 \(K\)，并用轮廓系数和多初值稳定性评价，而非预写簇数。
4. 区分得分、载荷、距离与 linkage 的输入对象，避免把二维图或簇标签解释为因果或真实类型。
5. 完成 M12 的冻结数据聚类证据包，使同伴能从输入快照复算表图与稳定性结论。

入口检查（5 分钟）

点 \((0,0)\) 与 \((2,0)\) 的均值中心是什么？到该中心的平方距离和是多少？

展开入口评分、补救与异形复测

答案为 \((1,0)\) 与 \(2\)。若第二问误用欧氏距离和，先复习平方距离，再对 \((0,0),(0,2)\) 重测；答案仍为 \(2\)。

无提示检索、渐隐链与迁移

闭书写出 K-Means 的“固定中心分配—固定分配更新—计算 WCSS—检查停止”链。第一轮完整对照 小节 12.4.1；第二轮只给初始中心；第三轮把点 \(E\) 改为 \((8,5)\) 并自行完成两轮。陌生迁移任务是对长三角上市公司供应链指标选择 PCA 或聚类方案，明确单位、异常值、稳定性与不能作出的经济结论。

目标—评价映射

目标	评价证据	达标标准
1	练习 1、2、6、9	推导、碎石图和载荷含义一致
2	练习 3	两轮分配、质心、WCSS 与停止判定无误
3	练习 5、6、7	选型与多初值稳定性证据完整
4	练习 4、8、11	距离/linkage 与结论边界正确
5	M12、练习 6—8	输入、输出和审计记录完整

• 是否有一种信息丰富的方式来可视化数据。
• 能否发现变量或观测值中的子群体。

无监督学习的挑战

与监督学习相比，无supervised 学习通常更具挑战性：

1. 没有明确的目标：没有响应变量来指导学习过程
2. 结果难以评估：缺乏像交叉验证这样通用的评估方法
3. 更主观：结果更多依赖研究者的判断

然而，无监督学习在数据探索和模式发现中至关重要。

12.2 无监督学习的主要类型

1. 主成分分析（PCA）：降维和可视化
2. 聚类分析：发现数据中的子群体
3. 关联规则挖掘：发现变量间的关联（如购物篮分析）
4. 异常检测：识别数据中的异常点

本章重点介绍前两种方法。

12.3 主成分分析（PCA）

PCA 的经典统计起点是 Pearson 对“最接近点系”的最小二乘几何构造；现代教材通常用最大方差、最小重构误差与 SVD 三种等价视角表述 (Pearson 1901年)。

12.3.1 PCA 的基本思想

假设我们有\(n\) 个观测，每个观测有\(p\) 个特征。当 \(p\) 很大时：

• 如何可视化数据？
• 是否能找到数据的低维表示，保留大部分信息？

核心思想：找到数据变化最大的方向，用这些方向来表示数据。

PCA 的直观解释

想象你在拍摄一个三维物体：

• 从某个角度拍摄，只能看到物体的部分信息
• 如果你找到了最佳拍摄角度，可以用一张二维照片最大程度地展现三维物体
• PCA 就是在寻找数据的”最佳拍摄角度

在数学上，PCA 寻找的是数据方差最大的方向，这些方向称为主成分。

12.3.2 第一主成分

第一主成分是特征的标准化线性组合：

\[ Z_1 = \phi_{11}X_1 + \phi_{21}X_2 + \cdots + \phi_{p1}X_p \tag{12.1}\]

满足约束： \[ \sum_{j=1}^{p} \phi_{j1}^2 = 1 \tag{12.2}\]

使得 \(Z_1\) 的方差最大化。

优化问题：

\[ \max_{\phi_{11}, \ldots, \phi_{p1}} \left\{ \frac{1}{n} \sum_{i=1}^{n} \left( \sum_{j=1}^{p} \phi_{j1} x_{ij} \right)^2 \right\} \]

约束：\(\sum_{j=1}^{p} \phi_{j1}^2 = 1\)

数学推导：PCA 的奇异值分解 (SVD) 视角

虽然上面是从方差最大化的角度定义主成分，但在实际计算和深层代数理论中，PCA 的奇异值分解（Singular Value Decomposition, SVD） 密不可分。

设中心化（列均值为 0）后的数据矩阵为 \(\mathbf{X}\in\mathbb{R}^{n\times p}\)。它的完整 SVD 为 \[ \mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T \] 其中 \(\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}^{T}\mathbf{X}\) 的特征向量。

关联：PCA 的联系。 样本的协方差矩阵（按 \(1/n\) 缩放）为 \(\mathbf{C} = \frac{1}{n} \mathbf{X}^T \mathbf{X}\)。 将 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. \] 若写紧 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}\mathbf{V} = \mathbf{U}\mathbf{\Sigma}\)。

SVD 视角的优势在于：避免了直接计算高维协方差矩阵 \(\mathbf{X}^T \mathbf{X}\)（这可能导致数值不稳定和巨大内存开销），使得工业界的软件（如 Scikit-Learn）在底层能够利用高度优化为SVD 算法稳健而高效地求解大规模PCA。

12.3.3 中国案例：基于财务指标的上市公司结构分析（非项目演示）

为了理解A股市场的内在结构，我们将对上市公司进行主成分分析。我们的目标是从众多的财务指标中提取出反映公司核心特征（如”规模”、“盈利能力”、“杠杆”）的主成分。

在这个实验环节，我们将丢掉以往那种目标明确的“预测涨跌”，转而像一位在星空中寻找星座的占星师一样，去探索A 股长三角地区（沪、苏、浙、皖）庞大上市公司群体的财务基因重构图谱。 下面的 PCA 流水线使用固定的 2022Q4 横截面，并提取总资产对数、资产负债率、净利率、资产周转率和 ROA。它是帮助读者理解载荷、得分与动态维数的非项目演示，采用全行业长三角样本和顺序分位数删除，不生成 M12 handoff，也不与 小节 12.9 的单行业完整案例清洗合同混称为同一证据。正式 M12 只以练习 6 的 m12-preregistration-v2.json、三个 roster 与排除账本为权威。 请格外注意代码的最后一步（StandardScaler）：若不标准化，大尺度的资产变量会支配基于方差和距离的结果。标准化只统一量纲，不会消除异常值、行业构成或会计口径差异。

import pandas as pd  # 数据分析库
import numpy as np  # 数值计算库
import hashlib  # 计算课程冻结文件的真实字节指纹
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  # 数据标准化工具

# 1. 加载数据
import os  # 文件系统操作
from pathlib import Path  # 使用跨平台路径对象解析显式数据根
def stream_sha256(path, block_size=1024 * 1024):  # 固定块流式哈希，避免复制整个HDF进内存
    digest = hashlib.sha256()  # 初始化摘要对象
    with path.open('rb') as source_file:  # 以二进制流读取真实字节
        for block in iter(lambda: source_file.read(block_size), b''):  # 逐块读取到EOF
            digest.update(block)  # 累加当前块
    return digest.hexdigest()  # 返回完整SHA-256
book_data_root_value = os.environ.get('BOOK_DATA_DIR', '').strip()
if not book_data_root_value:
    raise RuntimeError({'status': 'stopped', 'reason': 'BOOK_DATA_DIR_missing'})
DATA_DIR = Path(book_data_root_value).expanduser().resolve()  # 从必需环境变量取得数据根
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    raise FileNotFoundError(f'BOOK_DATA_DIR 不存在或不是目录: {DATA_DIR}')  # 失败即停止，不回退固定路径
path_fin = DATA_DIR / 'stock' / 'financial_statement.h5'  # 财务报表文件路径
path_basic = DATA_DIR / 'stock' / 'stock_basic_data.h5'  # 股票基本信息文件路径
m12_financial_slice = Path(os.environ['BOOK_M12_FINANCIAL_SLICE']).expanduser().resolve()  # 课程运行器提供的冻结小切片
m12_basic_slice = Path(os.environ['BOOK_M12_BASIC_SLICE']).expanduser().resolve()  # 只含项目列与公司的基本信息切片

# 在任何拟合前验证课程冻结清单；这些值必须由数据交付方提供，不能从文件名猜测
snapshot_contract = {
    'snapshot_id': os.environ.get('BOOK_DATA_SNAPSHOT_ID', '').strip(),
    'freeze_date': os.environ.get('BOOK_DATA_FREEZE_DATE', '').strip(),
    'financial_sha256': os.environ.get('BOOK_FINANCIAL_STATEMENT_SHA256', '').strip().lower(),
    'basic_sha256': os.environ.get('BOOK_STOCK_BASIC_SHA256', '').strip().lower(),
    'financial_slice_sha256': os.environ.get('BOOK_M12_FINANCIAL_SLICE_SHA256', '').strip().lower(),
    'basic_slice_sha256': os.environ.get('BOOK_M12_BASIC_SLICE_SHA256', '').strip().lower(),
    'pca_variance_threshold': os.environ.get('BOOK_M12_PCA_VARIANCE_THRESHOLD', '').strip(),
    'report_period': os.environ.get('BOOK_M12_REPORT_PERIOD', '').strip().lower(),
}
missing_snapshot_fields = [name for name, value in snapshot_contract.items() if not value]
if missing_snapshot_fields:
    raise RuntimeError({'status': 'stopped', 'reason': 'snapshot_contract_incomplete', 'missing': missing_snapshot_fields})
missing_input_files = [str(path) for path in (path_fin, path_basic, m12_financial_slice, m12_basic_slice) if not path.is_file()]
if missing_input_files:
    raise FileNotFoundError({'status': 'stopped', 'reason': 'missing_input_files', 'paths': missing_input_files})
actual_input_hashes = {'financial_sha256': stream_sha256(path_fin), 'basic_sha256': stream_sha256(path_basic), 'financial_slice_sha256': stream_sha256(m12_financial_slice), 'basic_slice_sha256': stream_sha256(m12_basic_slice)}  # 流式核验源与切片
if any(actual_input_hashes[name] != snapshot_contract[name] for name in actual_input_hashes):  # 拒绝任一输入字节漂移
    raise RuntimeError({'status': 'stopped', 'reason': 'input_sha256_mismatch', 'actual': actual_input_hashes})  # 在数据读取前结构化停止

financial_statements = pd.read_hdf(m12_financial_slice)  # 只读取课程冻结的报告期/列切片
stock_basic_data = pd.read_hdf(m12_basic_slice)  # 只读取课程冻结的总体映射切片
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) and required_basic_columns.issubset(stock_basic_data.columns)  # 缺列立即停止
print({'financial_slice_shape': financial_statements.shape, 'basic_slice_shape': stock_basic_data.shape, 'source_hashes': actual_input_hashes})  # 记录实际I/O规模

下面在独立代码块中冻结横截面并合并行业与省份字段；拆分只改善课堂阅读，不改变对象依赖。

# 2. 数据清洗与特征构造
# 固定同一会计报告期，避免把不同年份的公司截面混在一起

# 转换日期列为datetime类型，便于后续时间筛选
financial_statements['info_date'] = pd.to_datetime(financial_statements['info_date'])  # 将列转换为日期时间类型
freeze_date = pd.Timestamp(snapshot_contract['freeze_date'])  # 解析语义冻结日
assert financial_statements['info_date'].notna().all()  # 禁止以缺失披露日绕过时点守卫
financial_statements = financial_statements[financial_statements['info_date'].le(freeze_date)].copy()  # 先删冻结日后修订

target_report_period = snapshot_contract['report_period']  # 使用同一预注册报告期
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))  # 同一报告期内保留最后披露版本
print(f'统一报告期: {target_report_period}; 公司数: {recent_financials.order_book_id.nunique()}')

# 合并基本信息 (获取行业)：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()  # 筛选长三角地区上市公司

在成功加载并筛选出长三角地区上市公司数据后，下一步是构造关键财务指标并进行数据标准化处理。

表 12.1: 上市公司财务指标数据加载与清洗

# 定义原始列名到标准列名的映射关系
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']].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:  # 遍历循环
    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))  # 打印五个财务指标的描述统计

列表 12.1: PCA 输入的训练样本标准化

# 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.3.4 拟合 PCA 模型

我们将分析这 5 个核心财务指标的主成分。

标准化后的 5 维财务数据进入 PCA 后，首先要回答“保留多少维才能达到事前信息保留标准”。代码现场输出单项与累计解释方差；碎石图只帮助检查谱形状，不能预写前两维或前三维一定达到某个比例。用于聚类的维数必须由预注册阈值和现场累计解释率共同决定，二维得分仅服务于展示。

表 12.2: 主成分解释的方差比例

# 拟合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))  # 打印方差解释汇总表

表 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')  # X轴标签
plt.ylabel('解释方差比例', fontsize=12, fontproperties='Source Han Serif SC')  # Y轴标签
plt.title('PCA 碎石图(Scree Plot)', fontsize=14, fontproperties='Source Han Serif SC')  # 图标题
plt.grid(True, alpha=0.3)  # 添加网格线
plt.show()  # 显示图形

图 12.1

解释方差比例由当前统一报告期样本现场计算。选择主成分数量时应结合累计解释方差、碎石图、下游任务和跨报告期稳定性；“肘部”是诊断而非自动成立的结论。

12.3.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}')

图 12.2

图 12.2 把得分与系数载荷箭头叠加在同一二维坐标中。主成分名称必须依据上方当前载荷表给出。载荷整体变号不改变主成分空间，因此“正方向”本身没有经济优劣含义；只有系数载荷的绝对值、相对符号和解释方差可用于描述主成分方向，原始变量相关应另行计算。

12.3.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)  # X轴标签
plt.ylabel('PC2', fontsize=12)  # Y轴标签
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

图 12.3 是探索性投影。行业是否聚集、哪些行业重叠，应从本次图形和组内离散度计算判断；二维分离不能证明行业是主要驱动因素，也不能自动转化为分类或投资结论。

PCA 的实际应用建议

1. 数据标准化至关重要：不同量纲的变量必须先标准化
2. 主成分数量选择：
   • 碎石图法：选择拐点处的主成分数
   • 累积方差法：选择累积解释方差达到70-90%的主成分
   • Kaiser 准则：对标准化变量，常用经验规则是保留特征值大于 1 的主成分
3. 解释主成分：通过载荷矩阵理解每个主成分的实际含义
4. 主成分得分：可用于后续的聚类、回归等分析

12.4 聚类分析

聚类分析的目标是将相似的观测归为一组。在金融中，我们可以用它来发现财务特征相似的公司群体，这对于寻找对标公司（Benchmarking）非常有帮助。

12.4.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.3}\]

算法只交替做两件事：把每个点分给最近中心，再把每个中心更新为簇内均值。取六个点

\[ 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\)。再做一次分配不会改变标签，算法停止。这个例子也暴露初始化风险：较差初值会让第一轮形成极不平衡簇；因此真实分析应使用多初值并报告稳定性。

列表 12.2: 二维 K-Means 的两轮可复算迭代

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.4}\]

其中 \(\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.5}\]

\[ \mu_k^{(t+1)} = \frac{\sum_{i=1}^{n} \gamma_{ik} \cdot x_i}{\sum_{i=1}^{n} \gamma_{ik}} \tag{12.6}\]

\[ \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.7}\]

直观理解：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.4.2 K-Means 聚类

我们按预注册累计解释率阈值选择最小维数 \(d\) 进行聚类；前两维只用于显示。

如果我们不仅想用眼睛去“看”这些星团，还想让机器自动把它们“圈”出来，那就要轮到同为无监督三剑客之一——K-Means 聚类（K-Means Clustering）登场了。 代码使用现场累计解释率确定的 \(d\) 维得分，并在事前登记的 \(K=2,\ldots,8\) 中只保留满足 \(n\ge5K\) 的候选值。快照合同、有限数、样本量和有效簇检查均发生在拟合或绘图之前；若没有合法候选，代码以结构化失败记录停止。

from sklearn.cluster import KMeans  # K-Means聚类算法
from sklearn.metrics import adjusted_rand_score, silhouette_score  # 聚类质量与标签稳定性

variance_threshold = float(snapshot_contract['pca_variance_threshold'])  # 读取事前阈值
if not 0 < variance_threshold <= 1:
    raise ValueError({'status': 'stopped', 'reason': 'invalid_pca_variance_threshold'})
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维聚类输入
if clustering_features.ndim != 2 or clustering_features.shape[1] != selected_component_count:
    raise ValueError({'status': 'stopped', 'reason': 'invalid_cluster_matrix_shape', 'shape': clustering_features.shape})
if not np.isfinite(clustering_features).all():
    raise ValueError({'status': 'stopped', 'reason': 'non_finite_cluster_input'})

# 只遍历事前登记且满足 n>=5K 的候选，不存在合法 K 时不拟合
cluster_count_range = [k for k in range(2, 9) if len(clustering_features) >= 5 * k]
if not cluster_count_range:
    raise ValueError({'status': 'stopped', 'reason': 'no_valid_k', 'n': len(clustering_features), 'rule': '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)  # 对聚类特征数据拟合模型
    if np.unique(candidate_labels).size != k:
        raise ValueError({'status': 'stopped', 'reason': 'empty_or_invalid_cluster', 'k': k})
    sum_squared_errors.append(kmeans_model.inertia_)  # 记录SSE（inertia）
    candidate_silhouette = silhouette_score(clustering_features, candidate_labels)
    if not np.isfinite(candidate_silhouette):
        raise ValueError({'status': 'stopped', 'reason': 'non_finite_silhouette', '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)  # X轴标签
ax1.set_ylabel('Inertia (SSE)', fontsize=12)  # Y轴标签
ax1.set_title('肘部法则', fontsize=14, fontproperties='Source Han Serif SC')  # 图标题
ax1.grid(True, alpha=0.3)  # 添加网格线

ax2.plot(cluster_count_range, silhouette_scores, 'rs-')  # 绘制轮廓系数折线图
ax2.set_xlabel('K', fontsize=12)  # X轴标签
ax2.set_ylabel('轮廓系数 (Silhouette)', fontsize=12)  # Y轴标签
ax2.set_title('轮廓系数', fontsize=14, fontproperties='Source Han Serif SC')  # 图标题
ax2.grid(True, alpha=0.3)  # 添加网格线

plt.tight_layout()  # 自动调整子图间距
plt.show()  # 显示图形

图 12.4

下面将轮廓系数最大的候选值作为本次样本的数据驱动选择。这是一项探索性规则，不意味着存在固定不变的“真实 \(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)  # 拟合并预测每个公司的聚类标签
if np.unique(cluster_assignments).size != optimal_cluster_count:
    raise ValueError({'status': 'stopped', 'reason': 'invalid_final_clusters', 'k': 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)
    if np.unique(repeated_labels).size != optimal_cluster_count:
        raise ValueError({'status': 'stopped', 'reason': 'empty_cluster_in_stability_run', 'seed': initialization_seed})
    initialization_stability.append(adjusted_rand_score(cluster_assignments, repeated_labels))
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)),
}
if stability_summary['min_ari'] < 0.8:
    raise ValueError({'status': 'stopped', 'reason': 'low_initialization_stability', **stability_summary})
cleaned_financial_data['Cluster'] = cluster_assignments  # 将聚类结果存入数据框

print({'selected_k': optimal_cluster_count, 'snapshot': snapshot_contract, 'stability': stability_summary})

图 12.4 同时显示可行候选的 WCSS 与轮廓系数，簇数必须由现场输出和后续稳定性守卫共同决定。

12.4.3 聚类结果解释

聚类编号只是当次运行的数字代号，且簇数由上一节的候选规则动态决定。下面回到原始财务特征，按实际簇编号计算均值，并在前两个主成分上显示分群。颜色代表算法分配，不代表已知的“真实同盟”。

代码还会报告每簇的众数行业、平均资产和负债率。只有看到现场输出后，才可描述某簇是否呈现“重资产”或“高周转”等模式。这些是探索性画像，不是无偏的行业分类，也不足以直接支持并购或投资决策。

表 12.3: 各聚类财务特征画像

# 分析每个簇的特征均值(使用原始特征而非标准化后的值，便于直观理解)
cluster_feature_means = cleaned_financial_data.groupby('Cluster')[financial_features].mean()  # 按聚类分组计算各财务指标均值
print(cluster_feature_means.round(3))  # 打印各簇财务画像

表 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()  # 显示图形

图 12.5

列表 12.3: 按现场簇编号输出代表行业与财务画像

# 逐簇输出代表性行业与核心财务特征
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}')  # 打印该簇画像

表 12.3、图 12.5 与 列表 12.3 是解释聚类的共同起点。读者应根据当次运行的特征均值、载荷和代表行业描述各簇，并检查换随机种子、重抽样或更换期间后结论是否稳定。簇编号没有天然商业含义；未看到运行输出时，不应预设“银行”“软件”或“分区清晰”等结论。

12.4.4 层次聚类（Hierarchical Clustering）

层次聚类可以帮助我们理解公司之间的层级相似性，例如构建”同类公司”。

如果一K-Means 是在做粗暴的“切蛋糕”动作，那么层次聚类（Hierarchical Clustering）则是在编纂一部极其细腻的“公司财务进化谱系图”。 为了展现这幅画卷，并且不让全市场几千家公司把图谱挤得密不透风，我们在下面的代码中只随机抽取 50 家公司作为演示样本。Ward 最小方差法每一步选择使合并后组内平方和增加最小的簇对；这一准则来自 Ward 的层次分组方案 (Ward 1963年)。 最后输出的那张形似倒垂树根的图形，被统计学家们形象地称为“树状图（Dendrogram）”。最底层的每一片树叶都是一家活生生的公司，当两片树叶在极低的高度融合在一起时，意味着它们最近一期财报的结构几乎完全重合（它们可能是最完美的竞争对手或并购标的）；而随着树根向着上方图表的“高空”不断合拢，不同阵营间的壁垒也变得越来越难以跨越。这种无需预设 K 值，全凭距离远近自由生长的层级视觉，是探索小型高价值金融样本时无可替代的艺术品。

from scipy.cluster.hierarchy import dendrogram, linkage  # 层次聚类树状图与链接方法

# 为了可视化清晰，我们只随机选取 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')  # 图标题
plt.ylabel('距离', fontsize=12, fontproperties='Source Han Serif SC')  # Y轴表示合并距离
plt.tight_layout()  # 自动调整布局防止标签溢出
plt.show()  # 显示图形

图 12.6

图 12.6 的叶节点是当次随机抽取的公司，因此不预言银行、电力等具体行业的聚合顺序。读图时可记录低高度合并的公司对，再检查这种结构在更换抽样后是否稳定。与 K-Means 预先指定 \(K\) 不同，树状图可在不同高度切割成不同粒度的探索性分群，但切割高度仍是分析选择。

12.4.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}')  # 打印CH指数评估结果

代码会报告当次数据与所选 \(K\) 对应的 CH 指数。CH 指数没有跨数据集通用的绝对阈值；在同一数据和预处理下，可将它与不同 \(K\) 或聚类算法横向比较，但仍需结合稳定性与业务解释。

聚类算法选择指南

• K-Means: 数据量大，寻找明确分群。
• 层次聚类: 数据量小，寻找层级关系。
• DBSCAN: 寻找任意形状的簇，处理噪声（适合地理位置聚类等）。

12.5 无监督学习的挑战与局限

12.5.1 主成分分析的局限

1. 线性假设：PCA 只能捕获线性关系，无法处理非线性流形。
2. 解释困难：主成分是变量的线性组合，物理意义往往不直观。
3. 方差非信息：方差最大的方向不一定包含最有用的信息（例如信噪比问题）。

12.5.2 聚类分析的局限

1. 结果主观：K值的选择往往没有标准答案。
2. 不仅是数据：聚类结果需要业务解释，数学上分离好的簇在业务上可能没意义。
3. 稳定性：个别数据的变动可能导致聚类结果剧烈变化。

12.6 中国案例：基于股价波动的股票聚类

除了基于财务指标，我们还可以根据股价的历史波动形态对股票进行聚类。这种方法常用于构建投资组合（如分散化投资）或配对交易（Pairs Trading）。

最后，让我们把无监督学习的利刃从“静态的财报剖析”，转移向“动态的市场博弈”。在现代投资组合理论（MPT）中，为了防止黑天鹅事件将你的净值一次性清零，你必须买入那些“平时不怎么同涨同跌”的证券。然而，基于表面行业标签的去相关性往往极其脆弱。 这里从 5 个预先列出的行业各取若干股票，并使用 2023-01-01 之后、截至数据快照末日的收盘价计算日收益率。实际窗口长度由现场日期范围决定，不称为“过去一年”或“龙头股总体”。

这里的核心魔法在于第 3 步： 我们用简单的减法 1 - 相关系数 创造了一种全新的“距离度量”。如果两只股票像连体婴儿一样总是同涨同跌（相关系数接了1），它俩之间的距离就会逼近 0；而如果它们是在做跷跷板般的完美反向对冲（相关系数接近 -1），距离就会被无限拉大到 2。

随后执行的层次聚类和底层那张令人惊叹的“相关性热图（Clustermap）”将以前所未有的清晰度向你展示：在这个真实交易出来的平行宇宙里，资金究竟是如何暗中抱团的？半导体和白酒是否真的存在风格轮动效应？而那些披着不同行业外衣、实际上却有着相同宏观周期软肋的公司，将在这面照妖镜下无所遁形。

表 12.4

# 1. 加载股价数据与基本信息数据
import os  # 文件系统操作
from pathlib import Path  # 使用跨平台路径对象解析显式数据根
exercise_data_root_value = os.environ.get('BOOK_DATA_DIR', '').strip()
if not exercise_data_root_value:
    raise RuntimeError({'status': 'stopped', 'reason': 'BOOK_DATA_DIR_missing'})
DATA_DIR = Path(exercise_data_root_value).expanduser().resolve()  # 从必需环境变量取得数据根
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    raise FileNotFoundError(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'  # 股票基本信息文件路径

stock_price_history = pd.read_hdf(path_price).reset_index()  # 读取后复权股价数据并重置MultiIndex
stock_basic_data = pd.read_hdf(path_basic)  # 读取股票基本信息数据

# 2. 数据准备：构建收益率矩阵（行=日期，列=股票）
analysis_start_date = '2023-01-01'  # 固定分析起点；结束日由数据快照决定
recent_price_data = stock_price_history[stock_price_history['date'] > analysis_start_date].copy()  # 筛选近期数据

# 选取5个代表性行业，每个行业取6只龙头股
target_stock_universe = []  # 初始化目标股票池
target_industry_list = ['银行', '电子', '食品饮料', '非银行金融', '房地产']  # 5个代表性行业（中信一级分类）

for industry_name in target_industry_list:  # 遍历每个目标行业
    industry_top_stocks = stock_basic_data[stock_basic_data['citics_2019_l1_name'] == industry_name]['order_book_id'].head(6).tolist()  # 取该行业前6只股票
    target_stock_universe.extend(industry_top_stocks)  # 加入目标股票池

subset_price_data = recent_price_data[recent_price_data['order_book_id'].isin(target_stock_universe)]  # 筛选目标股票的价格数据
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}')  # 打印矩阵维度（天数 x 股票数）

上方输出动态报告当前收益率矩阵的交易日数和股票数。后续相关距离矩阵的维度由实际列数决定；缺失值筛选或数据版本变化时，不在正文保留旧的固定维度。

接下来基于收益率矩阵计算相关距离 \(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

图 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

图 12.8 与树状图使用同一距离和 linkage，因此排序应一致。是否出现块结构由当前图形决定；相关性会随窗口变化，不能单凭热图认定行业是“主要因子”或某组合已经有效分散。

12.6.1 应用场景

• 组合研究假设：距离较远的股票可能提供不同的短期共动结构，但仍需样本外协方差、权重与交易成本检验。
• 配对研究假设：距离近只表示相关，不等于价差平稳；配对交易还需协整、结构稳定性和样本外成本检验。

12.7 本章小结

本章我们系统学习了两种主要的无监督学习方法：

12.7.1 主成分分析（PCA）

核心思想：找到数据方差最大的方向，实现降维

关键步骤：

1. 数据标准化（零均值，单位方差）
2. 计算主成分
3. 选择主成分数量（碎石图、累积方差）
4. 解释主成分（载荷矩阵）
5. 计算主成分得分

应用场景：

• 数据可视化
• 降维以减少计算量
• 去除多重共线性
• 特征提取

12.7.2 聚类分析

核心思想：将相似观测归为一组

主要算法：

• K-Means：快速、简单，适合大数据
• 层次聚类：可视化层次结构，适合小数据

评估指标：

• 轮廓系数（Silhouette Score）
• Calinski-Harabasz 指数
• Davies-Bouldin 指数

应用场景：

• 客户细分
• 图像分割
• 文档聚类
• 异常检测

无监督学习最佳实践

1. 数据预处理至关重要：标准化、处理缺失值、异常值
2. 选择合适的方法：根据数据特点和业务目标
3. 验证结果：虽然难以验证，但应与业务知识结合
4. 迭代优化：尝试不同参数和方法，比较结果
5. 与监督学习结合：无监督学习可用于特征工程，提升监督学习效果

12.8 理论来源与前沿

无监督学习的“理论母体”主要来自三条线索：

1. 线性代数与谱理论：PCA 将高维数据投影到低维子空间，本质是对协方差矩阵进行特征分解；谱聚类进一步把图拉普拉斯矩阵的特征向量用于发现群落结构。
2. 优化与几何：K-Means 的目标函数是组内平方和（within-cluster sum of squares），其交替最小化过程可解释为一个EM 风格的坐标下降算法；层次聚类则对应一类凝聚式或分裂式的贪心优化。
3. 概率建模：高斯混合模型把聚类看作潜在类别变量的后验推断问题，天然支持软聚类与不确定性量化。

近年来的前沿发展集中在：

• 非线性降维：t-SNE、UMAP 等方法在保持局部邻域结构方面效果突出，但需要谨慎解读全局距离。
• 可扩展与鲁棒聚类：在大规模数据上，MiniBatch K-Means、近似最近邻与基于密度的方法（如 HDBSCAN）更实用。
• 与业务决策联动：在客户细分、风险分层等场景中，无监督结果往往需要与后续监督学习或策略优化结合，并通过稳定性分析与可解释性工具进行验证。

12.9 M12：无监督结构证据包

本里程碑与 小节 6.3 一一对应，且只在本章实施。

输入与依赖：使用课程冻结清单指向的同一财报报告期、公司范围与行业映射版本；数据根必须来自 BOOK_DATA_DIR。输入至少含 order_book_id、报告期、披露日、行业和四个比率的原始分子/分母。任何清洗或拟合前，代码先幂等冻结 m12-preregistration-v2.json，再生成课程级多行业 universe_roster、单行业 analysis_subset_roster、清洗后 clustering_roster 与逐公司排除账本。三个 roster 分别服务于 M13 总体、M12 行业选择和 PCA/K-Means 输入，满足 clustering ⊆ analysis_subset ⊆ universe，不能互换。

冻结规则：在查看或清洗数值前冻结报告期、单行业名称、行业映射版本、原始字段与四项公式、分母/缺失/缩尾规则、StandardScaler 口径、PCA 累计解释率阈值、候选 \(K=2,\ldots,8\)、\(n\ge5K\)、20 个完整种子、轮廓系数、ARI 与最低 ARI 0.8。任一字段改变会改变 preregistration hash；同名制品已有不同字节时必须在拟合前失败。

输出与提交证据：提交预注册 JSON、三个 roster、逐公司排除账本、数据血缘表、碎石图、载荷表、二维得分图、\(K\)—轮廓系数图、簇画像表、逐行 20 种子稳定性表、稳定 handoff 与独立运行遥测。载荷的机器证据是 m12-loadings-v2.csv；图 12.11 明确是得分投影，不得把它误称为载荷图。若另交载荷图，箭头坐标必须逐项等于载荷表的 PC1/PC2 列。稳定 handoff 只记录三个 roster 的用途、包含关系、行业映射版本/哈希，以及预注册、排除账本与六份机器输出的实际 SHA-256；六个固定名称为 loadings、scores、silhouettes、cluster-profile、stability-seeds、lineage。每次运行的时间、资源与命令只写入带唯一 run_id 的追加式终态遥测，不由 handoff 绑定。

验收与失败条件：主成分方差比之和在数值容差 \(10^{-8}\) 内不超过 1；聚类输入有限且已标准化；每个候选 \(K\) 的有效簇数等于 \(K\)。若 \(n<5K\)、出现空簇、轮廓系数不可定义或初值一致性过低，报告失败而非强行画像。结论只描述当前上市公司样本结构。

教师启动器与课堂预检：tests/fixtures/m12/manifest-v1.json 是版本化权威接口，并逐项列出唯一 canonical M12 代码链、全部必需环境值、完全相同的必交/授权文件集和签发后连续 setup 边界。教师不得用全章渲染代替项目启动器。真实授权快照到位后，在空输出目录运行：

conda run -n peter python scripts/m12_classroom_runner.py --manifest tests/fixtures/m12/manifest-v1.json --profile licensed --artifact-dir /absolute/new/m12-run --student-active-seconds 3300

启动器先完成 manifest、全部必需环境值、输入存在性与四个 SHA-256、日期、canonical 标签及源码编译检查。任一签发前检查失败只返回不含 run_id 的 not_started，不得创建制品目录或终态。全部通过后，启动器才同时创建 fresh 目录、签发私有 run_id 并立即注册父进程退出守卫，再启动仅含 canonical M12 链的 fresh 子进程。此后任何验证、抽取、子进程或输出审计失败都必须留下恰好一份规范 stopped 终态；成功时输出名必须恰好等于十三项必交集合。目录中的任何额外文件、第二份终态或临时候选都会触发父守卫关闭，不能保留 completed 冒充成功。终态分别记录 launcher_wall_seconds、project_process_seconds 与教师计入的 student_active_seconds，三者不得混称。

当前没有可分发的授权全量快照，因此“项目进程不超过 30 分钟、学生主动完成不超过 55 分钟”仍是尚未实测的设计目标，不是已验证承诺。教师可先用 CC0 合同夹具执行可再分发替代活动；它走完全相同的项目链、schema、PCA/K-Means 与原子交付代码，只支持合同教学，不得冒充真实 A 股经验结果：

conda run -n peter python scripts/m12_classroom_runner.py --profile contract-success --artifact-dir /absolute/new/m12-contract-success --run-id contract-success
conda run -n peter python scripts/m12_classroom_runner.py --profile contract-stop --artifact-dir /absolute/new/m12-contract-stop --run-id contract-stop
conda run -n peter python scripts/m12_classroom_runner.py --probe-setup-boundaries
conda run -n peter python scripts/m12_classroom_runner.py --probe-failure-contracts
conda run -n peter python scripts/m12_classroom_preflight.py --self-test

contract-success 必须产生完整十三项及一个 completed 终态；contract-stop 用不存在的分析行业确定性触发结构化停止；边界探针在签发后的每个连续 setup 块立即正常退出，并要求每次只留一个原子 stopped 终态。失败合同探针另覆盖缺环境、缺文件、hash 不符、非法日期、非法 canonical 标签、子进程失败和未授权输出：前五项必须是无 ID、无目录、无终态的 not_started，后两项必须由父守卫关闭为唯一规范 stopped。manifest 的 required_outputs 与 authorized_outputs 必须逐项相等。临时候选只允许使用同目录 .<authorized-final-name>.<uuid>.tmp 形式，原子发布后必须删除，且不属于学生交付物。

量规（20 分）：时点与血缘 5 分，PCA 数学/图形 4 分，聚类选择与稳定性 5 分，复现审计 3 分，解释边界 3 分。任一未来信息泄漏或快照不明，整项最高 8 分。评分只认下表中“代码输出—制品文件”的可复算配对；正文声称已经完成不构成证据。

可评分动词	必须由代码打印或断言的输出	必交制品文件	分值
预注册并核验时点/输入	冻结日、报告期、四个输入 hash、preregistration hash	m12-preregistration-v2.json、m12-lineage-v2.csv	3
冻结三层总体与排除	clustering ⊆ analysis_subset ⊆ universe、三层 \(n\)、逐原因计数	三个 m12-*-roster-v2.csv、m12-analysis-exclusions-v2.csv	2
拟合并解释 PCA	解释率和、动态维数、逐 PC 主导载荷	m12-loadings-v2.csv、m12-scores-v2.csv、碎石图与得分图	4
选择 \(K\) 并画像	每个可行 \(K\) 的轮廓系数、最终有效簇数、原尺度画像	m12-silhouettes-v2.csv、m12-cluster-profile-v2.csv	3
复测 20 个种子	恰好 20 个唯一种子、逐种子 ARI、最小/均值/标准差	m12-stability-seeds-v2.csv	2
记录并复现	签发前失败为无 ID/终态的 not_started；签发后任一失败只有一份规范 stopped；成功路径另记录 Python/包版本、命令、起止时刻、耗时、峰值内存、加载行数、输入/输出 hash，且稳定 handoff 可同字节重建	m12-run-<run_id>-v4.json、m12-handoff-v4.json；not_started 不产生文件	3
限定结论边界	当前样本结构、二维投影限制、簇编号无实体含义	决策说明中的对应段落	3

12.10 练习

12.10.1 概念题

1. [核心｜难度:2｜时间:10分钟｜分值:5｜项目:无] 解释主成分分析的几何意义。为什么第一主成分的方向对应于数据方差最大的方向？

2. [核心｜难度:1｜时间:8分钟｜分值:4｜项目:无] 在进行PCA 之前为什么要对数据进行标准化？如果不标准化会怎样？

3. [核心｜难度:2｜时间:18分钟｜分值: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 目标函数。

4. [核心｜难度:2｜时间:10分钟｜分值:5｜项目:无] 比较层次聚类和K-Means 聚类的优缺点。

5. [核心｜难度:2｜时间:10分钟｜分值:5｜项目:无] 解释轮廓系数的含义。它的取值范围是什么？如何解释？

12.10.2 应用题

6. [核心｜难度:3｜时间:55分钟｜分值:20｜项目:M12] 使用本章已登记的本地上市公司财务快照，严格按 小节 12.9 的证据矩阵完成一个行业的结构审计：先在查看/清洗数值前冻结 m12-preregistration-v2.json，再冻结课程级多行业 universe、单行业 analysis subset、清洗后 clustering roster 与逐公司排除账本；用现场累计解释率选择 PCA 维数，以载荷表解释因子，以得分进入聚类；只比较满足 \(n\ge5K\) 的候选，保留逐行 20 个种子的 ARI；最后提交表图、血缘、只绑定稳定制品 hash 的 handoff，以及按 run_id 独立追加的终态遥测。碎石图、载荷表和二维得分图是三个不同证据；不得把得分散点称为载荷图。

7. [拓展｜难度:3｜时间:35分钟｜分值:14｜项目:无] 受控模拟：K-Means 的几何边界。生成三个方差和重叠程度可控的二维高斯簇：

   • 比较标准化前后结果
   • 使用 K-Means，并用轮廓系数选择 \(K\)
   • 改变簇方差与随机种子，报告标签稳定性
   • 只解释算法几何性质，不把模拟簇称为真实客户群

8. [拓展｜难度:3｜时间:40分钟｜分值:15｜项目:无] 以长三角四省市中每家上市公司为分析单位，使用同一冻结报告期的资产规模、收入规模、净利率与 ROA 描述财务画像（不是省域经济水平）：

   • 冻结公司名册、报告期、输入 hash 与逐公司排除原因（3 分）
   • 在同一标准化公司特征矩阵上拟合 Ward 层次聚类并绘制公司级树状图（4 分）
   • 只比较满足 \(n\ge5K\) 的 KMeans 候选，以轮廓系数选择 \(K\)，再按同一个 \(K\) 切割 Ward 树（4 分）
   • 报告 Ward 与 KMeans 的 ARI、每簇公司数和省份构成，并说明省份只用于解释而不参与选型（2 分）
   • 若无可行 \(K\)，提交包含 n、候选规则和停止原因的结构化记录即可获得 KMeans 比较的 4 分，不得强行拟合（2 分）

12.10.3 理论题

9. [核心｜难度:3｜时间:20分钟｜分值:10｜项目:无] 证明 PCA 的主成分方向是数据协方差矩阵的特征向量。

10. [拓展｜难度:3｜时间:25分钟｜分值:12｜项目:无] 推导 K-Means 算法的EM（期望最大化）解释。

11. [拓展｜难度:3｜时间:25分钟｜分值:12｜项目:无] 证明在正态分布假设下，Ward 层次聚类方法等价于最小化合并后的方差增加量。

12. [拓展｜难度:2｜时间:25分钟｜分值:10｜项目:无] 研究并总结以下主题：

• t-SNE 和UMAP（非线性降维方法）
• DBSCAN 和谱聚类（基于密度的聚类）
• 高斯混合模型（基于概率的聚类）

12.11 练习参考解答

展开完整解答、评分点与常见失败模式

评分以题面分值为准；代数与 WCSS 允许 \(10^{-3}\) 绝对误差。应用题若缺碎石图、载荷含义、\(K\) 选择或多初值稳定性，每项扣 20%；用行业标签选择 \(K\)、把簇当真实类型或把省份上市公司画像写成省域经济结论，解释项不得分。

12.11.1 概念题解答

1. PCA 的几何意义

   PCA 的本质是在高维空间中寻找一组新的正交坐标轴（即主成分方向），使得数据在这些新轴上的投影方差依次递减。第一主成分的方向是数据方差最大的方向，因为最大化投影方差等价于最小化数据点到该方向的正交投影距离之和（即最小化重构误差）。从数学上看，这等价于求解数据协方差矩阵 \(\Sigma\) 的最大特征值对应的特征向量。几何上，第一主成分就是数据”椭球”的最长轴方向。

2. 标准化的必要性

   PCA 对变量的尺度敏感。如果不标准化，方差大的变量（如以”元”计的总资产）将主导主成分方向，而方差小的变量（如以百分比计的利润率）几乎没有影响。标准化（减均值、除以标准差）使所有变量在同一尺度上竞争，避免结果被单位差异所扭曲。例如，在分析上市公司财务指标时，营业收入（亿元级）和资产负债率（0-1之间）的量纲完全不同，不标准化的 PCA 结果几乎毫无意义。

3. 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 分。

4. 层次聚类 vs. K-Means 的优缺点

   特征	K-Means	层次聚类
   簇数 \(K\)	需预先指定	无需预先指定，树状图可后验选择
   计算复杂度	\(O(nKT)\)，较快	\(O(n^2 \log n)\) 或 \(O(n^3)\)，较慢
   簇形状	倾向于发现球形簇	可发现不规则形状簇（视联接方式）
   确定性	非确定性（依赖初始化）	确定性
   可解释性	中等	高（树状图直观展示层级关系）
   大数据适用性	适合大规模数据	不适合大规模数据

5. 轮廓系数的含义

   轮廓系数衡量数据点与自身簇的紧密度相对于与最近邻簇的分离度。对于数据点 \(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.11.2 应用题解答

6. 基于本地数据的行业财务因子PCA + 聚类

   我们使用本地的A 股财务数据进行分析，以”医药生物”行业为例。

列表 12.4: M12 fresh-kernel 导入与运行时准备

import pandas as pd  # 数据分析库
import numpy as np  # 数值计算库
import hashlib  # 重新计算冻结输入与制品的字节指纹
import json  # 规范化预注册与handoff合同
import os  # 读取运行合同与失败日志目录
import sys  # 记录解释器版本以支持环境复核
import sklearn  # 记录统计学习库版本以支持环境复核
import time  # 记录端到端墙钟耗时
import tracemalloc  # 记录本进程 Python 内存峰值
import uuid  # 为每次运行生成不会覆盖历史证据的唯一标识
import atexit  # 为已签发运行补写非覆盖异常退出终态
from pathlib import Path  # 在失败前也能定位授权日志目录
from datetime import datetime  # 解析课程运行器提供的固定起始时刻
from zoneinfo import ZoneInfo  # 强制运行日志采用中国标准时间
m12_started_perf = time.perf_counter()  # 必须早于任何输入读取
tracemalloc.start()  # 从 fresh-kernel 设置块开始计量峰值
def m12_timing_fields():  # 将端到端墙钟、子进程墙钟与学生主动用时分开
    try:
        launcher_started_ns = int(os.environ.get('BOOK_LAUNCHER_STARTED_NS', str(time.time_ns())))
        student_active_seconds = float(os.environ.get('BOOK_STUDENT_ACTIVE_SECONDS', '0'))
    except ValueError:
        launcher_started_ns, student_active_seconds = time.time_ns(), 0.0
    return {'launcher_wall_seconds': round(max(0.0, (time.time_ns() - launcher_started_ns) / 1_000_000_000), 6), 'project_process_seconds': round(time.perf_counter() - m12_started_perf, 6), 'student_active_seconds': round(max(0.0, student_active_seconds), 6)}

列表 12.5: M12 终态字节的schema、身份与规范化校验

def m12_validate_terminal_bytes(payload, expected_run_id):  # 校验最终路径上的权威终态字节
    try:  # 将截断或非JSON终态转成可诊断冲突
        terminal_record = json.loads(payload.decode('utf-8'))  # 只解析已经落盘的UTF-8字节
    except (UnicodeDecodeError, json.JSONDecodeError) as error:  # 捕获非法编码与截断JSON
        raise RuntimeError({'status': 'not_started', 'reason': 'terminal_json_invalid', 'run_id': expected_run_id}) from error  # 保留原文件并报告损坏
    required_values = {'schema_version': '4', 'run_id': expected_run_id}  # 冻结终态身份字段
    if not isinstance(terminal_record, dict) or any(terminal_record.get(key) != value for key, value in required_values.items()):  # 核对映射类型与身份
        raise RuntimeError({'status': 'not_started', 'reason': 'terminal_schema_or_run_id_invalid', 'run_id': expected_run_id})  # 禁止错身份终态占位
    if terminal_record.get('status') not in {'completed', 'stopped'}:  # 只允许一个规范终态枚举
        raise RuntimeError({'status': 'not_started', 'reason': 'terminal_status_invalid', 'run_id': expected_run_id})  # 中间状态不能关闭守卫
    canonical_bytes = (json.dumps(terminal_record, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 重建规范字节
    if payload != canonical_bytes:  # 禁止接受尾随内容或非规范序列化
        raise RuntimeError({'status': 'not_started', 'reason': 'terminal_bytes_not_canonical', 'run_id': expected_run_id})  # 拒绝可歧义或尾随字节
    return terminal_record  # 返回已核验且不可从文件外推测的终态

列表 12.6: M12 完整终态的非覆盖原子发布器

def m12_sync_directory(directory):  # 持久化原子名称更新
    directory_fd = os.open(directory, os.O_RDONLY)
    try:
        os.fsync(directory_fd)
    finally:
        os.close(directory_fd)
def m12_publish_terminal(path, record):  # completed、stopped与退出守卫共用唯一发布原语
    terminal_bytes = (json.dumps(record, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 固定候选字节
    m12_validate_terminal_bytes(terminal_bytes, record.get('run_id'))  # 发布前拒绝非法schema或状态
    temporary_path = path.with_name(f'.{path.name}.{uuid.uuid4().hex}.tmp')  # 在相同目录准备不可见完整候选
    with temporary_path.open('xb') as terminal_file:  # 临时候选本身也排他创建
        terminal_file.write(terminal_bytes)  # 最终名称出现前写完全部内容
        terminal_file.flush(); os.fsync(terminal_file.fileno())  # 发布前持久化候选字节
    was_published = False  # 区分本调用发布与先到终态
    try:  # 硬链接保证最终名称只看到完整文件且永不覆盖
        os.link(temporary_path, path)  # 以单次非覆盖操作公布完整终态
        was_published = True  # 记录本调用赢得唯一终态
        m12_sync_directory(path.parent)  # 将最终名称发布同步到磁盘
    except FileExistsError:  # 并发先到记录必须保留原字节
        was_published = False  # 交由调用者决定有效占用是否为冲突
    finally:  # 成功、竞争或异常都清理临时候选
        temporary_path.unlink(missing_ok=True)  # 成功、冲突或异常均不留下候选
    persisted_bytes = path.read_bytes()  # 只从最终路径重读权威字节
    m12_validate_terminal_bytes(persisted_bytes, record['run_id'])  # 显式拒绝截断、错ID或错状态
    if was_published and persisted_bytes != terminal_bytes:  # 发布者必须取回完全相同的字节
        raise RuntimeError({'status': 'stopped', 'reason': 'terminal_verification_failed', 'run_id': record['run_id']})  # 不将不一致字节视为终态
    return persisted_bytes, was_published  # 调用者据此处理成功或先到终态

列表 12.7: M12 签发前制品目录与原子硬链接能力验证

m12_artifact_value = os.environ.get('BOOK_M12_ARTIFACT_DIR', '').strip()
if not m12_artifact_value:
    raise RuntimeError({'status': 'not_started', 'reason': 'BOOK_M12_ARTIFACT_DIR_missing'})
m12_artifact_dir = Path(m12_artifact_value).expanduser().resolve()
m12_artifact_dir.mkdir(parents=True, exist_ok=True)
if not m12_artifact_dir.is_dir():
    raise RuntimeError({'status': 'not_started', 'reason': 'terminal_directory_not_directory', 'path': str(m12_artifact_dir)})
m12_probe_token = uuid.uuid4().hex
m12_probe_source = m12_artifact_dir / f'.m12-run-directory-probe-v4.json.{m12_probe_token}.tmp'
m12_probe_link = m12_artifact_dir / f'.m12-run-directory-probe-v4.json.{m12_probe_token}.link.tmp'
try:  # 实际验证创建、刷盘、同目录硬链接、目录刷盘与清理，而非只依赖os.access
    with m12_probe_source.open('xb') as probe_file:
        probe_file.write(b'm12-directory-probe\n'); probe_file.flush(); os.fsync(probe_file.fileno())
    os.link(m12_probe_source, m12_probe_link)
    if m12_probe_link.read_bytes() != b'm12-directory-probe\n':
        raise OSError('directory_probe_bytes_changed')
    m12_sync_directory(m12_artifact_dir)
except OSError as probe_error:
    raise RuntimeError({'status': 'not_started', 'reason': 'terminal_directory_not_atomically_writable', 'path': str(m12_artifact_dir), 'error_type': type(probe_error).__name__}) from probe_error
finally:
    m12_probe_link.unlink(missing_ok=True); m12_probe_source.unlink(missing_ok=True)

列表 12.8: M12 run_id 签发与立即注册的自包含退出守卫

def m12_issue_run(artifact_dir):  # 返回前保证签发身份已受守卫保护
    candidate = os.environ.get('BOOK_RUN_ID', '').strip() or uuid.uuid4().hex
    candidate_path = artifact_dir / f'm12-run-{candidate}-v4.json'
    if candidate_path.exists():
        existing = m12_validate_terminal_bytes(candidate_path.read_bytes(), candidate)
        raise RuntimeError({'status': 'not_started', 'reason': 'run_id_preexisting', 'run_id': candidate, 'existing_status': existing['status']})
    state = {'written': False, 'run_id': candidate, 'path': candidate_path, 'guard_error': None}
    def exit_guard(_state=state, _publish=m12_publish_terminal, _validate=m12_validate_terminal_bytes, _timing=m12_timing_fields, _datetime=datetime, _zone=ZoneInfo, _dumps=json.dumps, _stderr=sys.stderr):  # 捕获全部依赖，不引用后续代码块或可变全局
        if _state['written']:
            return
        record = {'schema_version': '4', 'run_id': _state['run_id'], 'status': 'stopped', 'reason': 'unexpected_exit', 'finished_at': _datetime.now(_zone('Asia/Shanghai')).isoformat(), **_timing()}
        try:
            persisted, published = _publish(_state['path'], record)
            _validate(persisted, _state['run_id'])
            _state['written'] = published or _state['path'].is_file()
        except BaseException as guard_error:
            _state['guard_error'] = {'reason': 'exit_guard_terminal_conflict', 'error_type': type(guard_error).__name__}
            print(_dumps(_state['guard_error'], sort_keys=True), file=_stderr)
    atexit.register(exit_guard)  # 身份仅在此闭包已登记后才返回给后续代码
    return candidate, candidate_path, state, exit_guard
m12_run_id, m12_run_log_path, m12_terminal_state, m12_exit_guard = m12_issue_run(m12_artifact_dir)

列表 12.9: M12 终态排他写入与落盘校验

def m12_write_terminal(record):  # 直接写终态，不回调停止函数
    candidate_bytes = (json.dumps(record, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 固定显式终态候选
    m12_validate_terminal_bytes(candidate_bytes, m12_run_id)  # 写入前绑定已签发身份、schema与状态
    persisted_bytes, was_published = m12_publish_terminal(m12_run_log_path, record)  # 原子发布并重读最终字节
    if not was_published:  # 已有合法终态也属于明确的身份冲突
        existing_terminal = m12_validate_terminal_bytes(persisted_bytes, m12_run_id)  # 报告先到状态但不改写
        raise RuntimeError({'status': 'stopped', 'reason': 'run_id_reused', 'run_id': m12_run_id, 'existing_status': existing_terminal['status'], 'terminal_path': str(m12_run_log_path)})  # 不递归写遥测
    m12_terminal_state['written'] = True  # 写后校验通过才关闭守卫
    atexit.unregister(m12_exit_guard)  # 避免正常终态后二次写入
    return m12_run_log_path, persisted_bytes  # 只返回落盘字节

列表 12.10: M12 追加式失败证据与必需环境变量入口

def m12_stop(reason, **details):  # 所有结构化停止统一写stopped遥测
    loaded_rows = details.pop('loaded_rows', {})  # 将行数与其他详情分开
    failure_record = {'schema_version': '4', 'run_id': m12_run_id, 'status': 'stopped', 'reason': reason, 'finished_at': datetime.now(ZoneInfo('Asia/Shanghai')).isoformat(), **m12_timing_fields(), 'peak_memory_bytes': int(tracemalloc.get_traced_memory()[1]), 'loaded_rows': loaded_rows, 'details': details}  # 分开记录墙钟、子进程与学生主动用时
    m12_write_terminal(failure_record)  # 在抛错前写唯一终态
    raise RuntimeError(failure_record)  # 阻止失败运行继续产出
def m12_require_env(name):  # 在任何直接索引前统一形成结构化缺参证据
    value = os.environ.get(name, '').strip()
    if not value:
        m12_stop('required_environment_missing', name=name)
    return value

列表 12.11: M12 PCA、KMeans 与绘图库导入

from sklearn.decomposition import PCA  # 主成分分析
from sklearn.preprocessing import StandardScaler  # 数据标准化工具
from sklearn.cluster import KMeans  # K-Means聚类算法
from sklearn.metrics import adjusted_rand_score, silhouette_score  # 聚类选择与稳定性
import matplotlib.pyplot as plt  # 绘图库
import os  # 文件系统操作
from pathlib import Path  # 使用跨平台路径对象解析显式数据根

列表 12.12: M12 输入路径、稳定序列化与幂等写入工具

# 1. 加载数据
def exercise_stream_sha256(path, block_size=1024 * 1024):  # fresh-kernel流式哈希
    digest = hashlib.sha256()  # 初始化摘要
    with path.open('rb') as source_file:  # 打开真实字节流
        for block in iter(lambda: source_file.read(block_size), b''):  # 固定块读取
            digest.update(block)  # 更新摘要
    return digest.hexdigest()  # 返回完整指纹
def exercise_stable_csv(frame):  # 固定roster、排除与输出表字节
    return frame.to_csv(index=False, date_format='%Y-%m-%d', lineterminator='\n').encode('utf-8')  # 固定编码与换行
def exercise_read_slice(path):  # 真实课堂HDF与可再分发CSV合同夹具走同一schema链
    suffix = path.suffix.lower()
    if suffix in {'.h5', '.hdf', '.hdf5'}:
        return pd.read_hdf(path)
    if suffix == '.csv':
        return pd.read_csv(path)
    m12_stop('unsupported_slice_format', path=str(path), suffix=suffix)
def exercise_write_same(path, payload):  # 同名不同字节时fail closed
    temporary_path = path.with_name(f'.{path.name}.{uuid.uuid4().hex}.tmp')  # 在同目录准备完整候选
    with temporary_path.open('xb') as artifact_file:  # 临时名也不允许复用
        artifact_file.write(payload)  # 写入全部候选字节
        artifact_file.flush(); os.fsync(artifact_file.fileno())  # 发布前刷新内容
    try:  # 完整临时文件只用硬链接原子发布一次
        os.link(temporary_path, path)  # 读者不会看到部分稳定制品
    except FileExistsError:
        if path.read_bytes() != payload: m12_stop('artifact_changed', path=str(path))  # 冲突统一停止
    finally:
        temporary_path.unlink(missing_ok=True)  # 只清理未发布的临时名字
    persisted_bytes = path.read_bytes()  # 始终从磁盘重读权威字节
    if persisted_bytes != payload:  # 捕获写后竞态或损坏
        m12_stop('artifact_verification_failed', path=str(path))  # 不将未验证制品交接
    return persisted_bytes  # 后续只哈希落盘字节

列表 12.13: M12 授权数据根与冻结切片路径

exercise_data_root_value = os.environ.get('BOOK_DATA_DIR', '').strip()
if not exercise_data_root_value:
    m12_stop('BOOK_DATA_DIR_missing')
DATA_DIR = Path(exercise_data_root_value).expanduser().resolve()  # 从必需环境变量取得数据根
if not DATA_DIR.is_dir():  # 在读取前验证数据根
    m12_stop('BOOK_DATA_DIR_not_directory', path=str(DATA_DIR))
path_fin = DATA_DIR / 'stock' / 'financial_statement.h5'  # 财务报表文件路径
path_basic = DATA_DIR / 'stock' / 'stock_basic_data.h5'  # 股票基本信息文件路径
exercise_financial_slice = Path(m12_require_env('BOOK_M12_FINANCIAL_SLICE')).expanduser().resolve()  # 冻结财务切片
exercise_basic_slice = Path(m12_require_env('BOOK_M12_BASIC_SLICE')).expanduser().resolve()  # 冻结总体切片

列表 12.14: M12 冻结输入合同与 SHA-256 核验

exercise_snapshot_contract = {
    'snapshot_id': os.environ.get('BOOK_DATA_SNAPSHOT_ID', '').strip(),
    'freeze_date': os.environ.get('BOOK_DATA_FREEZE_DATE', '').strip(),
    'financial_sha256': os.environ.get('BOOK_FINANCIAL_STATEMENT_SHA256', '').strip().lower(),
    'basic_sha256': os.environ.get('BOOK_STOCK_BASIC_SHA256', '').strip().lower(),
    'pca_variance_threshold': os.environ.get('BOOK_M12_PCA_VARIANCE_THRESHOLD', '').strip(),
    'financial_slice_sha256': os.environ.get('BOOK_M12_FINANCIAL_SLICE_SHA256', '').strip().lower(),
    'basic_slice_sha256': os.environ.get('BOOK_M12_BASIC_SLICE_SHA256', '').strip().lower(),
    'report_period': os.environ.get('BOOK_M12_REPORT_PERIOD', '').strip().lower(),
    'industry_name': os.environ.get('BOOK_M12_INDUSTRY_NAME', '').strip(),
    'industry_mapping_version': os.environ.get('BOOK_M12_INDUSTRY_VERSION', '').strip(),
}
exercise_missing_snapshot_fields = [name for name, value in exercise_snapshot_contract.items() if not value]
if exercise_missing_snapshot_fields:
    m12_stop('snapshot_contract_incomplete', missing=exercise_missing_snapshot_fields)
exercise_missing_files = [str(path) for path in (path_fin, path_basic, exercise_financial_slice, exercise_basic_slice) if not path.is_file()]
if exercise_missing_files:
    m12_stop('missing_input_files', paths=exercise_missing_files)
exercise_actual_hashes = {'financial_sha256': exercise_stream_sha256(path_fin), 'basic_sha256': exercise_stream_sha256(path_basic), 'financial_slice_sha256': exercise_stream_sha256(exercise_financial_slice), 'basic_slice_sha256': exercise_stream_sha256(exercise_basic_slice)}  # 流式核验源与切片
if any(exercise_actual_hashes[name] != exercise_snapshot_contract[name] for name in exercise_actual_hashes):  # 对照课程冻结清单
    m12_stop('input_sha256_mismatch', actual=exercise_actual_hashes)  # 字节漂移时停止
if Path(m12_require_env('BOOK_M12_ARTIFACT_DIR')).expanduser().resolve() != m12_artifact_dir:  # 禁止签发后切换制品目录
    m12_stop('artifact_directory_changed_after_run_id')  # 终态仍写入最初授权目录

列表 12.15: M12 清洗前预注册合同冻结

m12_preregistration = {'schema_version': '2', 'snapshot': exercise_snapshot_contract, 'input_hashes': exercise_actual_hashes, 'raw_fields': ['total_equity', 'operating_revenue', 'total_assets', 'total_liabilities', 'net_profit'], 'feature_formulas': {'ROE': 'net_profit/total_equity', 'Net_Margin': 'net_profit/operating_revenue', 'Asset_Turnover': 'operating_revenue/total_assets', 'Debt_Ratio': 'total_liabilities/total_assets'}, 'denominator_rule': 'total_equity!=0; operating_revenue!=0; total_assets>0', 'missing_rule': 'replace +/-inf with NA; complete case on four features', 'winsor_rule': 'none', 'standardization': 'sklearn.StandardScaler with population variance (ddof=0)', 'pca_rule': 'smallest d with cumulative explained variance + 1e-12 >= threshold; cap at all components', 'pca_variance_threshold': float(exercise_snapshot_contract['pca_variance_threshold']), 'k_candidates': list(range(2, 9)), 'sample_size_rule': 'n>=5*K', 'k_metric': 'maximum silhouette; ties choose smallest K', 'candidate_kmeans': {'random_state': 42, 'n_init': 20}, 'stability_seeds': list(range(20)), 'stability_metric': 'adjusted_rand_score versus final labels', 'minimum_ari': 0.8, 'final_kmeans': {'random_state': 42, 'n_init': 50}, 'report_period': exercise_snapshot_contract['report_period'], 'analysis_industry': exercise_snapshot_contract['industry_name'], 'industry_mapping_version': exercise_snapshot_contract['industry_mapping_version']}  # 完整冻结所有事前分析选择
if not 0 < m12_preregistration['pca_variance_threshold'] <= 1:  # 在清洗前拒绝无效阈值
    m12_stop('invalid_pca_variance_threshold')  # 不允许拟合后修阈值
m12_preregistration_bytes = (json.dumps(m12_preregistration, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 规范化完整合同
m12_preregistration_bytes = exercise_write_same(m12_artifact_dir / 'm12-preregistration-v2.json', m12_preregistration_bytes)  # 清洗拟合前原子冻结并重读
m12_preregistration_sha256 = hashlib.sha256(m12_preregistration_bytes).hexdigest()  # 只哈希落盘权威字节

列表 12.16: M12 多行业总体与单行业分析名册冻结

financial_statements = exercise_read_slice(exercise_financial_slice)  # 真实HDF与CC0 CSV夹具共用生产schema
stock_basic_data = exercise_read_slice(exercise_basic_slice)  # 文件格式不改变后续权威项目链
required_financial_fields = {'order_book_id', 'quarter', 'info_date'} | set(m12_preregistration['raw_fields'])  # 从预注册恢复财务schema
required_basic_fields = {'order_book_id', 'citics_2019_l1_name'}  # 冻结公司—行业映射schema
if not required_financial_fields.issubset(financial_statements.columns) or not required_basic_fields.issubset(stock_basic_data.columns):
    m12_stop('slice_schema_missing', financial_missing=sorted(required_financial_fields - set(financial_statements.columns)), basic_missing=sorted(required_basic_fields - set(stock_basic_data.columns)))
if stock_basic_data['order_book_id'].duplicated().any():
    m12_stop('duplicate_company_in_universe')
if stock_basic_data['citics_2019_l1_name'].isna().any():
    m12_stop('industry_mapping_missing')
exercise_freeze_date = pd.Timestamp(exercise_snapshot_contract['freeze_date'])  # 解析事前冻结日期
financial_statements['info_date'] = pd.to_datetime(financial_statements['info_date'])  # 统一披露日类型以审计可得性
financial_statements = financial_statements[financial_statements['info_date'].le(exercise_freeze_date)].copy()  # 在选版本前排除冻结日后修订
universe_roster = stock_basic_data[['order_book_id', 'citics_2019_l1_name']].rename(columns={'citics_2019_l1_name': 'industry_name'}).sort_values('order_book_id').reset_index(drop=True)  # 冻结课程级多行业总体
universe_roster['industry_mapping_version'] = m12_preregistration['industry_mapping_version']  # 绑定行业口径版本
universe_roster['preregistration_sha256'] = m12_preregistration_sha256  # 绑定完整事前规则
if universe_roster['industry_name'].nunique() < 2:
    m12_stop('universe_not_multi_industry', industry_n=int(universe_roster['industry_name'].nunique()))
analysis_subset_roster = universe_roster[universe_roster['industry_name'].eq(m12_preregistration['analysis_industry'])].copy()  # 单行业只服务M12
if analysis_subset_roster.empty or not set(analysis_subset_roster['order_book_id']).issubset(universe_roster['order_book_id']):
    m12_stop('analysis_subset_invalid')
universe_bytes = exercise_write_same(m12_artifact_dir / 'm12-universe-roster-v2.csv', exercise_stable_csv(universe_roster))  # 原子冻结并重读多行业总体
analysis_subset_bytes = exercise_write_same(m12_artifact_dir / 'm12-analysis-subset-roster-v2.csv', exercise_stable_csv(analysis_subset_roster))  # 原子冻结并重读单行业子集

列表 12.17: M12 冻结报告期与披露时点审计

analysis_stock_codes = analysis_subset_roster['order_book_id']  # 从冻结单行业roster恢复分析公司
recent_pharma_financials = financial_statements[financial_statements['order_book_id'].isin(analysis_stock_codes)]  # 只筛M12子集财报
recent_pharma_financials = recent_pharma_financials[recent_pharma_financials['quarter'].astype(str).str.lower().eq(m12_preregistration['report_period'])].sort_values('info_date').groupby('order_book_id').tail(1)  # 选择冻结期最后可得版本
if recent_pharma_financials['info_date'].gt(exercise_freeze_date).any():  # 禁止使用冻结日以后披露的数据
    m12_stop('info_date_after_freeze_date')  # 披露时点越界即停止
print({'financial_slice_shape': financial_statements.shape, 'basic_slice_shape': stock_basic_data.shape, 'universe_industries': universe_roster['industry_name'].nunique(), 'preregistration_sha256': m12_preregistration_sha256})  # 数据血缘首行

取得课程冻结报告期内且披露日不晚于冻结日的医药行业上市公司财报后，接下来构建核心财务指标（ROE、净利润率、资产周转率、负债率），并运用 PCA 与 K-Means 描述当前样本内部结构。

列表 12.18: M12 公司级财务比率构造

# 构建财务指标；先标记无效分母，不用任意常数改变比率含义
valid_denominators = ((recent_pharma_financials['total_equity'] != 0) & (recent_pharma_financials['operating_revenue'] != 0) & (recent_pharma_financials['total_assets'] > 0))  # 执行预注册分母规则
valid_financials = recent_pharma_financials.loc[valid_denominators].copy()  # 仅保留可定义四项比率的报告
pharma_analysis_data = pd.DataFrame({'order_book_id': valid_financials['order_book_id']})  # 初始化清洗前特征表
pharma_analysis_data['ROE'] = valid_financials['net_profit'] / valid_financials['total_equity']  # 按冻结公式计算净资产收益率
pharma_analysis_data['Net_Margin'] = valid_financials['net_profit'] / valid_financials['operating_revenue']  # 按冻结公式计算净利润率
pharma_analysis_data['Asset_Turnover'] = valid_financials['operating_revenue'] / valid_financials['total_assets']  # 按冻结公式计算资产周转率
pharma_analysis_data['Debt_Ratio'] = valid_financials['total_liabilities'] / valid_financials['total_assets']  # 按冻结公式计算资产负债率
pharma_analysis_data = pharma_analysis_data.replace([np.inf, -np.inf], np.nan).dropna().sort_values('order_book_id').reset_index(drop=True)  # 执行完整案例规则
report_ids = set(recent_pharma_financials['order_book_id'])  # 记录冻结报告期有可得版本的公司
denominator_ids = set(valid_financials['order_book_id'])  # 记录分母规则通过的公司
complete_ids = set(pharma_analysis_data['order_book_id'])  # 记录完整特征公司的集合

列表 12.19: M12 清洗后名册与逐公司排除账本

analysis_exclusion_ledger = analysis_subset_roster.copy()  # 从完整单行业子集建立并绑定合同的排除账本
analysis_exclusion_ledger['reason'] = np.select([~analysis_exclusion_ledger['order_book_id'].isin(report_ids), ~analysis_exclusion_ledger['order_book_id'].isin(denominator_ids), ~analysis_exclusion_ledger['order_book_id'].isin(complete_ids)], ['missing_report_version', 'invalid_denominator', 'missing_or_nonfinite_feature'], default='included')  # 按事前顺序给出唯一原因
excluded_analysis_rows = analysis_exclusion_ledger[analysis_exclusion_ledger['reason'].ne('included')].sort_values('order_book_id').reset_index(drop=True)  # 只输出严格排除记录
clustering_roster = analysis_subset_roster[analysis_subset_roster['order_book_id'].isin(complete_ids)].sort_values('order_book_id').reset_index(drop=True)  # 第三层只服务M12拟合
if not set(clustering_roster['order_book_id']).issubset(analysis_subset_roster['order_book_id']):
    m12_stop('clustering_roster_outside_analysis_subset')
if len(clustering_roster) + len(excluded_analysis_rows) != len(analysis_subset_roster):
    m12_stop('analysis_roster_accounting_mismatch')
clustering_bytes = exercise_write_same(m12_artifact_dir / 'm12-clustering-roster-v2.csv', exercise_stable_csv(clustering_roster))  # 原子冻结并重读聚类名册
exclusion_bytes = exercise_write_same(m12_artifact_dir / 'm12-analysis-exclusions-v2.csv', exercise_stable_csv(excluded_analysis_rows))  # 原子冻结并重读排除账本
exercise_roster_sha256 = hashlib.sha256(clustering_bytes).hexdigest()  # 保留后续展示所需的聚类名册指纹
exercise_excluded_count = len(excluded_analysis_rows)  # 从排除账本计算排除数

列表 12.20: M12 有限性守卫与标准化

# 2. PCA 分析
financial_features = ['ROE', 'Net_Margin', 'Asset_Turnover', 'Debt_Ratio']  # 定义分析特征列表
predictor_matrix = pharma_analysis_data[financial_features]  # 提取特征矩阵
predictor_array = predictor_matrix.to_numpy(dtype=float)
if len(predictor_array) < 10:
    m12_stop('no_valid_k', n=len(predictor_array), rule='n>=5K for K>=2', loaded_rows={'clustering_roster': len(predictor_array)})
if not np.isfinite(predictor_array).all():
    m12_stop('non_finite_financial_features', loaded_rows={'clustering_roster': len(predictor_array)})

feature_scaler = StandardScaler()  # 初始化标准化器
scaled_predictor_matrix = feature_scaler.fit_transform(predictor_array)  # 对特征进行标准化

列表 12.21: M12 PCA 拟合与动态维数选择

pca_model = PCA()  # 创建PCA模型
pca_model.fit(scaled_predictor_matrix)  # 拟合PCA

# 打印碎石图信息
print('Explained Variance Ratio:', pca_model.explained_variance_ratio_)  # 输出各主成分解释方差比例

# 3. 按预注册累计解释率选择聚类维数；二维仅用于展示
exercise_variance_threshold = m12_preregistration['pca_variance_threshold']  # 只从已冻结权威合同读取阈值
if not 0 < exercise_variance_threshold <= 1:  # 限制阈值为有效比例
    m12_stop('invalid_pca_variance_threshold')  # 无效阈值时停止
exercise_variance_tolerance = 1e-12  # 固定浮点边界容差
exercise_cumulative_variance = np.cumsum(pca_model.explained_variance_ratio_)  # 计算现场累计解释率
qualifying_components = np.flatnonzero(exercise_cumulative_variance + exercise_variance_tolerance >= exercise_variance_threshold)  # 容差下寻找首个达标边界
selected_component_count = int(qualifying_components[0] + 1) if qualifying_components.size else len(exercise_cumulative_variance)  # 阈值1未命中时使用全部成分
selected_component_count = min(max(selected_component_count, 1), len(exercise_cumulative_variance))  # 硬性限制结果在可用维数内
pca_cluster_features = np.asarray(pca_model.transform(scaled_predictor_matrix)[:, :selected_component_count], dtype=float)  # 用前d维聚类
if pca_cluster_features.ndim != 2 or pca_cluster_features.shape[1] != selected_component_count:  # 核对动态聚类矩阵
    m12_stop('invalid_cluster_matrix_shape', shape=pca_cluster_features.shape)  # 维数异常时停止
if not np.isfinite(pca_cluster_features).all():  # 禁止非有限值进入距离计算
    m12_stop('non_finite_cluster_input')  # 输入异常时停止

完整交付先显示 图 12.9 和载荷，再按轮廓系数选择满足 \(n\ge5K\) 的 \(K\)；不存在固定 \(K\) 的先行拟合路径。

component_numbers = np.arange(1, len(pca_model.explained_variance_ratio_) + 1)  # 建立主成分序号
plt.plot(component_numbers, pca_model.explained_variance_ratio_, marker='o', label='单项')  # 绘制单项解释率
plt.plot(component_numbers, np.cumsum(pca_model.explained_variance_ratio_), marker='s', label='累计')  # 绘制累计解释率
plt.xlabel('主成分')  # 标明横轴为主成分序号
plt.ylabel('解释方差比例')  # 标明纵轴的统计含义
plt.legend()  # 区分单项与累计曲线
plt.show()  # 输出可评分碎石图

图 12.9

列表 12.22: 习题6：载荷、候选 K 与簇画像审计

loading_table = pd.DataFrame(pca_model.components_.T, index=financial_features, columns=[f'PC{i + 1}' for i in range(len(financial_features))])  # 生成载荷表
candidate_cluster_counts = [k for k in m12_preregistration['k_candidates'] if len(pharma_analysis_data) >= 5 * k]  # 落实冻结的 n>=5K
if not candidate_cluster_counts:
    m12_stop('no_valid_k', n=len(pharma_analysis_data), minimum_n=10, loaded_rows={'clustering_roster': len(pharma_analysis_data)})
silhouette_by_k = {}  # 保存每个候选 K 的现场分数
for cluster_count in candidate_cluster_counts:  # 遍历事前允许的候选簇数
    candidate_model = KMeans(n_clusters=cluster_count, **m12_preregistration['candidate_kmeans'])  # 使用冻结的候选模型参数
    candidate_labels = candidate_model.fit_predict(pca_cluster_features)  # 基于达到阈值的前d个主成分聚类
    if np.unique(candidate_labels).size != cluster_count:
        m12_stop('empty_cluster', k=cluster_count, loaded_rows={'clustering_roster': len(pharma_analysis_data)})
    silhouette_by_k[cluster_count] = silhouette_score(pca_cluster_features, candidate_labels)  # 计算轮廓分数
if not np.isfinite(list(silhouette_by_k.values())).all():
    m12_stop('non_finite_silhouette', loaded_rows={'clustering_roster': len(pharma_analysis_data)})
selected_cluster_count = sorted(silhouette_by_k, key=lambda cluster_count: (-silhouette_by_k[cluster_count], cluster_count))[0]  # 轮廓并列时选最小K
final_cluster_labels = KMeans(n_clusters=selected_cluster_count, **m12_preregistration['final_kmeans']).fit_predict(pca_cluster_features)  # 用冻结参数拟合最终模型
if np.unique(final_cluster_labels).size != selected_cluster_count:
    m12_stop('invalid_final_clusters', k=selected_cluster_count, loaded_rows={'clustering_roster': len(pharma_analysis_data)})

silhouette_table = pd.DataFrame({'k': list(silhouette_by_k), 'silhouette': list(silhouette_by_k.values())})  # 形成可评分选择表
plt.plot(silhouette_table['k'], silhouette_table['silhouette'], marker='o')  # 绘制现场选择证据
plt.xlabel('K')  # 标注候选簇数
plt.ylabel('轮廓系数')  # 标注评价准则
plt.show()  # 输出K—轮廓图

图 12.10

图 12.10 是最终 \(K\) 的现场选择证据；未出现在图中的候选因预注册样本量守卫而不可行。

列表 12.23: M12 二十种子稳定性审计

# 在本题真实医药公司样本上以20个单初值种子检验稳定性
exercise_stability_records = []
for initialization_seed in m12_preregistration['stability_seeds']:  # 遍历冻结的完整20种子集合
    repeated_labels = KMeans(n_clusters=selected_cluster_count, n_init=1, random_state=initialization_seed).fit_predict(pca_cluster_features)  # 在同一动态维数复测初始化
    if np.unique(repeated_labels).size != selected_cluster_count:
        m12_stop('empty_cluster_in_stability_run', seed=initialization_seed, loaded_rows={'clustering_roster': len(pharma_analysis_data)})
    exercise_stability_records.append(adjusted_rand_score(final_cluster_labels, repeated_labels))
exercise_stability_summary = {
    'seed_count': len(exercise_stability_records),
    'mean_ari': float(np.mean(exercise_stability_records)),
    'min_ari': float(np.min(exercise_stability_records)),
    'std_ari': float(np.std(exercise_stability_records)),
}
if exercise_stability_summary['min_ari'] < m12_preregistration['minimum_ari']:  # 执行冻结的稳定性门槛
    m12_stop('low_initialization_stability', loaded_rows={'clustering_roster': len(pharma_analysis_data)}, **exercise_stability_summary)
exercise_seed_table = pd.DataFrame({'seed': m12_preregistration['stability_seeds'], 'ari': exercise_stability_records})  # 保留每次种子的完整证据

表 12.5: M12 数据血缘、载荷、稳定性与簇画像现场输出

pharma_cluster_profile = pharma_analysis_data.assign(cluster=final_cluster_labels).groupby('cluster')[financial_features].agg(['mean', 'median', 'count'])  # 输出原尺度画像
exercise_lineage_table = pd.DataFrame([{'snapshot_id': exercise_snapshot_contract['snapshot_id'], 'freeze_date': exercise_snapshot_contract['freeze_date'], 'report_period': exercise_snapshot_contract['report_period'], 'industry': exercise_snapshot_contract['industry_name'], 'n': len(pharma_analysis_data), 'missing_excluded': exercise_excluded_count, 'winsor_rule': 'none; complete-case after valid denominators', 'selected_components': selected_component_count, 'selected_k': selected_cluster_count, 'seed_count': len(exercise_seed_table), 'preregistration_sha256': m12_preregistration_sha256}])  # 机器可读血缘
print(loading_table)  # 用最大绝对载荷现场命名主成分
print(exercise_lineage_table.to_string(index=False))  # 输出可评分数据血缘表
print({'actual_hashes': exercise_actual_hashes, 'roster_sha256': exercise_roster_sha256, 'silhouette_by_k': silhouette_by_k, 'stability': exercise_stability_summary})  # 输出完整可审计合同
print(pharma_cluster_profile)  # 报告每簇均值、中位数与样本量
print(exercise_seed_table.to_string(index=False))  # 输出完整20种子稳定性表

列表 12.24: M12 六份机器表的稳定序列化

score_table = pd.DataFrame(pca_cluster_features, columns=[f'PC{index + 1}' for index in range(selected_component_count)])  # 定型实际聚类得分
score_table.insert(0, 'order_book_id', pharma_analysis_data['order_book_id'].to_numpy())  # 将得分绑定公司键
profile_table = pharma_cluster_profile.reset_index()  # 把簇画像索引还原为字段
profile_table.columns = ['cluster'] + [f'{feature}_{summary}' for feature, summary in profile_table.columns[1:]]  # 固定多层列名
m12_output_frames = {'loadings': loading_table.reset_index(names='feature'), 'scores': score_table, 'silhouettes': silhouette_table, 'cluster-profile': profile_table, 'stability-seeds': exercise_seed_table, 'lineage': exercise_lineage_table}  # 汇总六份正式输出
m12_expected_output_names = {'loadings', 'scores', 'silhouettes', 'cluster-profile', 'stability-seeds', 'lineage'}  # 固定生产者—消费者接口
if set(m12_output_frames) != m12_expected_output_names or len(m12_output_frames) != 6:
    m12_stop('m12_output_schema_mismatch', expected=sorted(m12_expected_output_names), actual=sorted(m12_output_frames))
m12_output_hashes = {}  # 初始化输出指纹表
for output_name, output_frame in m12_output_frames.items():  # 逐一绑定同一预注册合同
    bound_frame = output_frame.copy()  # 不改写现场展示对象
    bound_frame['preregistration_sha256'] = m12_preregistration_sha256  # 把合同指纹写入每份机器输出
    output_bytes = exercise_stable_csv(bound_frame)  # 固定输出表字节
    output_bytes = exercise_write_same(m12_artifact_dir / f'm12-{output_name}-v2.csv', output_bytes)  # 原子冻结并重读正式输出
    m12_output_hashes[output_name] = hashlib.sha256(output_bytes).hexdigest()  # 只登记落盘输出指纹

列表 12.25: M12 稳定 handoff 构建与写后核验

required_run_fields = ['BOOK_RUN_STARTED_AT', 'BOOK_RUN_COMMAND', 'BOOK_GIT_COMMIT']  # 冻结运行身份所需字段
missing_run_fields = [field for field in required_run_fields if not os.environ.get(field, '').strip()]  # 检查运行器是否给出完整身份
if missing_run_fields:  # 缺少运行身份时不能声称可复现
    m12_stop('run_identity_incomplete', missing=missing_run_fields, loaded_rows={'financial_slice': len(financial_statements), 'basic_slice': len(stock_basic_data), 'clustering_roster': len(clustering_roster)})
run_started_at = datetime.fromisoformat(os.environ['BOOK_RUN_STARTED_AT'])  # 解析运行器固定的 ISO 8601 时刻
if run_started_at.tzinfo is None or run_started_at.astimezone(ZoneInfo('Asia/Shanghai')).utcoffset().total_seconds() != 28800:  # 要求可转换为中国标准时间
    m12_stop('run_started_at_not_cst')  # 拒绝无时区或错误时区
industry_mapping_table = universe_roster[['order_book_id', 'industry_name', 'industry_mapping_version']].copy()  # 定型课程级行业映射
industry_mapping_sha256 = hashlib.sha256(exercise_stable_csv(industry_mapping_table)).hexdigest()  # 哈希映射版本对应的实际字节
m12_roster_hashes = {'universe': hashlib.sha256(universe_bytes).hexdigest(), 'analysis_subset': hashlib.sha256(analysis_subset_bytes).hexdigest(), 'clustering': hashlib.sha256(clustering_bytes).hexdigest()}  # 分别登记三层roster
m12_handoff = {'schema_version': '4', 'snapshot_id': exercise_snapshot_contract['snapshot_id'], 'freeze_date': exercise_snapshot_contract['freeze_date'], 'report_period': exercise_snapshot_contract['report_period'], 'analysis_industry': exercise_snapshot_contract['industry_name'], 'preregistration_sha256': m12_preregistration_sha256, 'source_hashes': exercise_actual_hashes, 'industry_mapping': {'version': m12_preregistration['industry_mapping_version'], 'sha256': industry_mapping_sha256}, 'rosters': {'universe': {'file': 'm12-universe-roster-v2.csv', 'sha256': m12_roster_hashes['universe'], 'purpose': 'M13课程级多行业计划总体'}, 'analysis_subset': {'file': 'm12-analysis-subset-roster-v2.csv', 'sha256': m12_roster_hashes['analysis_subset'], 'purpose': 'M12单行业分析边界'}, 'clustering': {'file': 'm12-clustering-roster-v2.csv', 'sha256': m12_roster_hashes['clustering'], 'purpose': 'M12清洗后PCA与聚类样本'}}, 'containment': 'clustering subset of analysis_subset subset of universe', 'exclusions': {'file': 'm12-analysis-exclusions-v2.csv', 'sha256': hashlib.sha256(exclusion_bytes).hexdigest(), 'count': exercise_excluded_count}, 'outputs': m12_output_hashes, 'output_names': sorted(m12_expected_output_names), 'selected_components': selected_component_count, 'selected_k': selected_cluster_count}  # 稳定入口完全排除运行遥测字段
m12_handoff_bytes = (json.dumps(m12_handoff, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 规范化交接合同字节
m12_handoff_path = m12_artifact_dir / 'm12-handoff-v4.json'  # 定位唯一稳定跨章入口
m12_handoff_bytes = exercise_write_same(m12_handoff_path, m12_handoff_bytes)  # 排他写入并取回落盘handoff字节
m12_handoff_sha256 = exercise_stream_sha256(m12_handoff_path)  # 完成遥测只引用磁盘实际指纹

列表 12.26: M12 完整必交输出名称与落盘存在性核对

m12_required_stable_files = ['m12-preregistration-v2.json', 'm12-universe-roster-v2.csv', 'm12-analysis-subset-roster-v2.csv', 'm12-clustering-roster-v2.csv', 'm12-analysis-exclusions-v2.csv', 'm12-loadings-v2.csv', 'm12-scores-v2.csv', 'm12-silhouettes-v2.csv', 'm12-cluster-profile-v2.csv', 'm12-stability-seeds-v2.csv', 'm12-lineage-v2.csv', 'm12-handoff-v4.json']  # 冻结十二份稳定交付名称
m12_missing_stable_files = [name for name in m12_required_stable_files if not (m12_artifact_dir / name).is_file()]  # 在成功终态前核对全部稳定文件
if m12_missing_stable_files:  # 缺任一必交文件都不能宣告完成
    m12_stop('required_output_missing', files=m12_missing_stable_files)  # 追加停止终态并保留已完成制品
m12_required_output_contract = m12_required_stable_files + [m12_run_log_path.name]  # 加入本次唯一终态形成十三项完整合同
if len(m12_required_output_contract) != len(set(m12_required_output_contract)) or len(m12_required_output_contract) != 13:  # 防止重复或接口漏项
    m12_stop('required_output_contract_invalid', outputs=m12_required_output_contract)  # 不允许模糊输出接口进入handoff

列表 12.27: M12 经 handoff 核验后追加唯一 completed 终态

m12_finished_at = datetime.now(ZoneInfo('Asia/Shanghai'))  # 只在稳定制品验证后记录成功时刻
m12_peak_memory_bytes = tracemalloc.get_traced_memory()[1]  # 读取 fresh-kernel Python 峰值
m12_loaded_rows = {'financial_slice': int(len(financial_statements)), 'basic_slice': int(len(stock_basic_data)), 'universe_roster': int(len(universe_roster)), 'analysis_subset_roster': int(len(analysis_subset_roster)), 'clustering_roster': int(len(clustering_roster))}  # 明示每个关键输入与名册的现场行数
m12_run_log = {'schema_version': '4', 'run_id': m12_run_id, 'status': 'completed', 'started_at': run_started_at.astimezone(ZoneInfo('Asia/Shanghai')).isoformat(), 'finished_at': m12_finished_at.isoformat(), **m12_timing_fields(), 'peak_memory_bytes': int(m12_peak_memory_bytes), 'loaded_rows': m12_loaded_rows, 'command': os.environ['BOOK_RUN_COMMAND'], 'git_commit': os.environ['BOOK_GIT_COMMIT'], 'python': sys.version.split()[0], 'packages': {'pandas': pd.__version__, 'numpy': np.__version__, 'scikit_learn': sklearn.__version__}, 'preregistration_sha256': m12_preregistration_sha256, 'input_hashes': exercise_actual_hashes, 'output_hashes': m12_output_hashes, 'handoff_sha256': m12_handoff_sha256}  # 墙钟、项目进程与学生主动用时分列，遥测只引用落盘handoff指纹
m12_run_log_path, m12_run_log_bytes = m12_write_terminal(m12_run_log)  # 为本run_id追加唯一completed终态
m12_run_log_sha256 = hashlib.sha256(m12_run_log_bytes).hexdigest()  # 计算实际终态字节指纹
print({'m12_handoff_sha256': m12_handoff_sha256, 'm12_run_log': m12_run_log_path.name, 'm12_run_log_sha256': m12_run_log_sha256, 'rosters': m12_roster_hashes, 'preregistration_sha256': m12_preregistration_sha256})  # 分别输出稳定入口与本次终态遥测

exercise_display_scores = pca_model.transform(scaled_predictor_matrix)[:, :2]  # 二维得分只用于可视化
plt.scatter(exercise_display_scores[:, 0], exercise_display_scores[:, 1], c=final_cluster_labels, cmap='viridis')  # 展示动态全维聚类的二维投影
plt.title('医药行业上市公司财务聚类', fontproperties='Source Han Serif SC')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()

图 12.11

主成分命名必须引用每列绝对载荷最大的变量及符号；载荷接近时写“混合维度”。若没有候选 \(K\) 满足 \(n\ge5K\)、任一运行出现无效簇，或 20 个初值中的最小 ARI 低于 0.8，本题结构化停止并登记原因。

7. 受控模拟：K-Means 的几何边界

   本题只演示三个人工高斯簇的几何性质，不声称数据来自真实客户或支持商业细分结论。

import pandas as pd  # 数据分析库
import numpy as np  # 数值计算库
from sklearn.cluster import KMeans  # K-Means聚类算法
from sklearn.metrics import silhouette_score  # 轮廓系数评估聚类质量
from sklearn.preprocessing import StandardScaler  # 本题独立比较是否标准化
import matplotlib.pyplot as plt  # 绘图库

# 1. 模拟数据（此处为教学演示用途，模拟三类客户群体）
np.random.seed(42)  # 设置随机数种子保证可复现性
number_of_samples = 500  # 模拟客户总数
# 模拟三个不同消费力的客户群体：低收入低消费、中收入中消费、高收入高消费
simulated_customer_features = np.concatenate([  # 拼接数组
    np.random.normal([20, 20], [5, 5], size=(number_of_samples//3, 2)),  # 低收入低消费群体
    np.random.normal([60, 50], [10, 10], size=(number_of_samples//3, 2)),  # 中收入中消费群体
    np.random.normal([90, 80], [10, 10], size=(number_of_samples//3, 2))  # 高收入高消费群体
])  # 完成构建
customer_data_frame = pd.DataFrame(simulated_customer_features, columns=['Income', 'Score'])  # 构建客户数据框

模拟数据集构建完成后，利用轮廓系数（Silhouette Score）在 $K = 2$ 到 $9$ 的范围内搜索最优聚类数，然后使用最优 $K$ 值执行最终聚类并可视化客户细分结果。

# 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(customer_data_frame)  # 拟合模型
    silhouette_scores.append(silhouette_score(customer_data_frame, kmeans_model.labels_))  # 记录轮廓系数

plt.plot(cluster_count_range, silhouette_scores, 'o-')  # 绘制轮廓系数曲线
plt.title('Silhouette Score')  # 图标题
plt.show()  # 峰值位置以本次运行输出为准

图 12.12: 习题7：受控三簇样本的候选 K 轮廓系数

图 12.12 冻结 \(K\) 的选择证据；图 12.13 只展示依该选择得到的二维标签。

# 3. 使用现场轮廓系数最高的 K 进行最终聚类
selected_k = list(cluster_count_range)[int(np.argmax(silhouette_scores))]
kmeans_model = KMeans(n_clusters=selected_k, n_init=10, random_state=42)  # 初始化最终模型
cluster_assignments = kmeans_model.fit_predict(customer_data_frame)  # 拟合并预测聚类标签

plt.scatter(customer_data_frame['Income'], customer_data_frame['Score'], c=cluster_assignments)  # 按聚类标签着色绘制散点图
plt.xlabel('Annual Income (k$)')  # X轴：年收入
plt.ylabel('Spending Score (1-100)')  # Y轴：消费得分
plt.title('Customer Segments')  # 图标题
plt.show()  # 显示图形

图 12.13: 习题7：现场轮廓最优 K 的二维 KMeans 分群

固定一次种子不能完成“稳定性”要求。下面在两种重叠程度、是否标准化和 10 个随机种子上复算，并用调整兰德指数（ARI）比较各运行与同场景参考标签；标签编号置换不会影响 ARI。

列表 12.28: 习题7：尺度、重叠与随机初值的稳定性审计

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，结论应写“对初始化不稳定”，而不是只展示最好的一次运行。

8. 基于本地数据的长三角上市公司财务画像聚类

分析单位是一家公司；省份只用于解释簇构成，不进入选型。这样 Ward 与 KMeans 在同一公司级矩阵上比较，并且 \(n\ge5K\) 有机会成立。样本仍受上市公司覆盖与行业构成影响，不能代表完整省域经济水平。

列表 12.29: 习题8：公司级 Ward–KMeans 比较的 fresh-kernel 设置

import os
import hashlib  # 核验练习8使用的同一冻结输入文件
import json  # 冻结练习合同与输出清单
import uuid  # 为每次独立尝试生成追加式标识
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
exercise8_attempt_id = os.environ.get('BOOK_EX8_ATTEMPT_ID', '').strip() or uuid.uuid4().hex  # 运行器可注入尝试ID

列表 12.30: 习题8：Ward、KMeans 与稳定性工具导入

from scipy.cluster.hierarchy import dendrogram, fcluster, linkage  # 层次树、按K切割与Ward链接
from sklearn.cluster import KMeans  # 在同一公司矩阵上拟合KMeans
from sklearn.metrics import adjusted_rand_score, silhouette_score  # 选择K并比较两种标签
from sklearn.preprocessing import StandardScaler

列表 12.31: 习题8：稳定序列化与结构化状态工具

def exercise8_stream_sha256(path, block_size=1024 * 1024):  # 本题独立流式哈希
    digest = hashlib.sha256()  # 初始化摘要
    with path.open('rb') as source_file:  # 打开字节流
        for block in iter(lambda: source_file.read(block_size), b''):  # 固定块读取
            digest.update(block)  # 更新摘要
    return digest.hexdigest()  # 返回指纹
def exercise8_stable_csv(frame):  # 固定名册、账本和审计表字节
    return frame.to_csv(index=False, date_format='%Y-%m-%d', lineterminator='\n').encode('utf-8')

列表 12.32: 习题8：按尝试ID追加的唯一终态

exercise8_artifact_value = os.environ.get('BOOK_EX8_ARTIFACT_DIR', '').strip()
def exercise8_terminal_path():  # 每次尝试只对应一个终态文件
    if not exercise8_artifact_value:  # 未授权目录时不自行推测路径
        return None  # 终态仍随异常返回
    status_dir = Path(exercise8_artifact_value).expanduser().resolve()  # 解析授权目录
    status_dir.mkdir(parents=True, exist_ok=True)  # 允许首次尝试建目录
    return status_dir / f'ex8-attempt-{exercise8_attempt_id}-v1.json'  # 停止与成功不共享全局状态名
def exercise8_write_terminal(record):  # 直接追加尝试终态以避免递归
    status_path = exercise8_terminal_path()  # 定位当前尝试的唯一终态
    status_bytes = (json.dumps(record, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 规范化尝试证据
    if status_path is not None:  # 有授权目录时排他落盘
        try:
            with status_path.open('xb') as status_file:  # 原子创建消除并发尝试的竞态窗口
                status_file.write(status_bytes)  # 追加本次尝试的唯一终态
        except FileExistsError: raise RuntimeError({'status': 'stopped', 'reason': 'attempt_id_reused', 'attempt_id': exercise8_attempt_id, 'path': str(status_path)})  # 保留原终态
    return status_path, status_bytes  # 供成功路径输出校验
def exercise8_stop(reason, **details):  # 停止尝试保留可后续成功的历史证据
    record = {'schema_version': '1', 'attempt_id': exercise8_attempt_id, 'status': 'stopped', 'reason': reason, 'details': details}  # 绑定尝试身份与失败现场
    exercise8_write_terminal(record)  # 在停止前追加独立终态
    raise RuntimeError(record)  # 阻止失败尝试继续拟合

列表 12.33: 习题8：稳定制品排他写入与落盘校验

def exercise8_write_same(path, payload):  # 禁止同名稳定制品被结果后改写
    temporary_path = path.with_name(f'.{path.name}.{uuid.uuid4().hex}.tmp')  # 在同目录完整写候选
    with temporary_path.open('xb') as artifact_file:
        artifact_file.write(payload)  # 临时内容尚不对稳定文件读者可见
        artifact_file.flush(); os.fsync(artifact_file.fileno())  # 原子发布前刷新
    try:
        os.link(temporary_path, path)  # 以硬链接排他公布完整字节
    except FileExistsError:
        if path.read_bytes() != payload: exercise8_stop('artifact_changed', path=str(path))
    finally:
        temporary_path.unlink(missing_ok=True)  # 不触碰先到的稳定制品
    persisted_bytes = path.read_bytes()  # 重读权威落盘字节
    if persisted_bytes != payload:  # 捕获写后竞态或损坏
        exercise8_stop('artifact_verification_failed', path=str(path))  # 禁止未验证制品进入清单
    return persisted_bytes  # 后续只哈希落盘字节

列表 12.34: 习题8：必需环境变量入口

def exercise8_require_env(name):  # 所有必需环境变量先走同一失败接口
    value = os.environ.get(name, '').strip()
    if not value:
        exercise8_stop('required_environment_missing', name=name)
    return value

列表 12.35: 习题8：冻结输入路径与授权制品目录

# 1. 独立加载并验证快照，不依赖正文或前题留在内存中的对象
exercise8_data_root_value = exercise8_require_env('BOOK_DATA_DIR')
exercise8_data_dir = Path(exercise8_data_root_value).expanduser().resolve()
exercise8_paths = {
    'financial': exercise8_data_dir / 'stock' / 'financial_statement.h5',
    'basic': exercise8_data_dir / 'stock' / 'stock_basic_data.h5',
    'financial_slice': Path(exercise8_require_env('BOOK_M12_FINANCIAL_SLICE')).expanduser().resolve(),
    'basic_slice': Path(exercise8_require_env('BOOK_M12_BASIC_SLICE')).expanduser().resolve(),
}
exercise8_artifact_dir = Path(exercise8_require_env('BOOK_EX8_ARTIFACT_DIR')).expanduser().resolve()
exercise8_artifact_dir.mkdir(parents=True, exist_ok=True)

列表 12.36: 习题8：冻结输入合同与字节核验

exercise8_snapshot_contract = {
    'snapshot_id': os.environ.get('BOOK_DATA_SNAPSHOT_ID', '').strip(),
    'freeze_date': os.environ.get('BOOK_DATA_FREEZE_DATE', '').strip(),
    'financial_sha256': os.environ.get('BOOK_FINANCIAL_STATEMENT_SHA256', '').strip().lower(),
    'basic_sha256': os.environ.get('BOOK_STOCK_BASIC_SHA256', '').strip().lower(),
    'financial_slice_sha256': os.environ.get('BOOK_M12_FINANCIAL_SLICE_SHA256', '').strip().lower(),
    'basic_slice_sha256': os.environ.get('BOOK_M12_BASIC_SLICE_SHA256', '').strip().lower(),
    'report_period': os.environ.get('BOOK_M12_REPORT_PERIOD', '').strip().lower(),
}
exercise8_missing_contract = [name for name, value in exercise8_snapshot_contract.items() if not value]
exercise8_missing_files = [str(path) for path in exercise8_paths.values() if not path.is_file()]
if exercise8_missing_contract or exercise8_missing_files:
    exercise8_stop('input_contract_incomplete', missing_contract=exercise8_missing_contract, missing_files=exercise8_missing_files)
exercise8_actual_hashes = {name + '_sha256': exercise8_stream_sha256(path) for name, path in exercise8_paths.items()}  # 流式重算源与切片指纹
if any(exercise8_actual_hashes[name] != exercise8_snapshot_contract[name] for name in exercise8_actual_hashes):  # 对照课程清单
    exercise8_stop('input_sha256_mismatch', actual=exercise8_actual_hashes)  # 字节漂移时停止

列表 12.37: 习题8：冻结长三角公司级分析名册

exercise8_financials = pd.read_hdf(exercise8_paths['financial_slice'])
exercise8_basic = pd.read_hdf(exercise8_paths['basic_slice'])
required_financial_columns = {'order_book_id', 'quarter', 'info_date', 'operating_revenue', 'net_profit', 'total_assets'}
required_basic_columns = {'order_book_id', 'province'}
missing_columns = sorted((required_financial_columns - set(exercise8_financials)) | (required_basic_columns - set(exercise8_basic)))
if missing_columns:
    exercise8_stop('missing_required_columns', columns=missing_columns)

exercise8_financials['info_date'] = pd.to_datetime(exercise8_financials['info_date'])
exercise8_financials = exercise8_financials[exercise8_financials['info_date'].le(pd.Timestamp(exercise8_snapshot_contract['freeze_date']))].copy()  # 先过滤冻结日后修订
exercise8_period = exercise8_snapshot_contract['report_period']
exercise8_cross_section = (exercise8_financials.loc[exercise8_financials['quarter'].astype(str).str.lower() == exercise8_period]
                           .sort_values('info_date').groupby('order_book_id').tail(1))
required_provinces = {'上海市', '江苏省', '浙江省', '安徽省'}
exercise8_basic_roster = exercise8_basic.loc[exercise8_basic['province'].isin(required_provinces), ['order_book_id', 'province']].drop_duplicates().sort_values('order_book_id')
if exercise8_basic_roster['order_book_id'].duplicated().any():
    exercise8_stop('duplicate_company_in_basic_roster')
exercise8_raw_roster = exercise8_basic_roster.merge(exercise8_cross_section, on='order_book_id', how='left', validate='one_to_one')  # 左联结保留每个原始候选及其排除证据
observed_provinces = set(exercise8_raw_roster['province'].dropna())
if observed_provinces != required_provinces:
    exercise8_stop('incomplete_province_coverage', observed=sorted(observed_provinces))

列表 12.38: 习题8：公司级财务特征与可行 K 守卫

exercise8_feature_frame = exercise8_raw_roster[['order_book_id', 'province', 'quarter', 'info_date', 'operating_revenue', 'net_profit', 'total_assets']].copy()
exercise8_feature_frame['log_revenue'] = np.log(exercise8_feature_frame['operating_revenue'].where(exercise8_feature_frame['operating_revenue'].gt(0)))
exercise8_feature_frame['log_assets'] = np.log(exercise8_feature_frame['total_assets'].where(exercise8_feature_frame['total_assets'].gt(0)))
exercise8_feature_frame['net_margin'] = exercise8_feature_frame['net_profit'] / exercise8_feature_frame['operating_revenue']
exercise8_feature_frame['roa'] = exercise8_feature_frame['net_profit'] / exercise8_feature_frame['total_assets']
exercise8_feature_names = ['log_revenue', 'log_assets', 'net_margin', 'roa']  # 冻结两项规模与两项盈利特征
exercise8_feature_frame = exercise8_feature_frame.replace([np.inf, -np.inf], np.nan)
exercise8_exclusion_ledger = exercise8_feature_frame[['order_book_id', 'province']].copy()
exercise8_exclusion_ledger['reason'] = np.select(
    [exercise8_feature_frame['quarter'].isna(), exercise8_feature_frame['total_assets'].le(0), exercise8_feature_frame['operating_revenue'].le(0), exercise8_feature_frame[exercise8_feature_names].isna().any(axis=1)],
    ['missing_report_version', 'nonpositive_total_assets', 'nonpositive_operating_revenue', 'missing_or_nonfinite_feature'], default='included')
exercise8_company_profiles = exercise8_feature_frame.loc[exercise8_exclusion_ledger['reason'].eq('included'), ['order_book_id', 'province'] + exercise8_feature_names].sort_values('order_book_id').reset_index(drop=True)
exercise8_analysis_roster = exercise8_company_profiles[['order_book_id', 'province']].copy()
exercise8_contract = {'schema_version': '1', 'snapshot': exercise8_snapshot_contract, 'input_hashes': exercise8_actual_hashes, 'required_provinces': sorted(required_provinces), 'features': exercise8_feature_names, 'feature_rules': {'log_revenue': 'log(operating_revenue), operating_revenue>0', 'log_assets': 'log(total_assets), total_assets>0', 'net_margin': 'net_profit/operating_revenue', 'roa': 'net_profit/total_assets'}, 'k_candidates': list(range(2, 9)), 'sample_size_rule': 'n>=5*K', 'selection': 'maximum KMeans silhouette; ties choose smallest K', 'comparison': 'Ward and KMeans on identical StandardScaler matrix and selected K'}
exercise8_contract_bytes = (json.dumps(exercise8_contract, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')
exercise8_contract_bytes = exercise8_write_same(exercise8_artifact_dir / 'ex8-contract-v1.json', exercise8_contract_bytes)  # 原子冻结并重读合同
exercise8_contract_sha256 = hashlib.sha256(exercise8_contract_bytes).hexdigest()  # 只哈希落盘合同字节
for exercise8_name, exercise8_frame in {'raw-roster': exercise8_raw_roster[['order_book_id', 'province']], 'analysis-roster': exercise8_analysis_roster, 'exclusions': exercise8_exclusion_ledger.loc[exercise8_exclusion_ledger['reason'].ne('included')]}.items():
    exercise8_write_same(exercise8_artifact_dir / f'ex8-{exercise8_name}-v1.csv', exercise8_stable_csv(exercise8_frame))
exercise8_feature_matrix = exercise8_company_profiles[exercise8_feature_names].to_numpy(dtype=float)  # 建立公司×特征矩阵
exercise8_valid_k = [cluster_count for cluster_count in range(2, 9) if len(exercise8_company_profiles) >= 5 * cluster_count]  # 只保留满足样本量守卫的候选
if not exercise8_valid_k:  # 无可行K时给出可评分停止记录
    exercise8_stop('no_valid_k_for_comparison', n=len(exercise8_company_profiles), rule='n>=5K', candidates=list(range(2, 9)), contract_sha256=exercise8_contract_sha256)  # 序列化可评分停止证据
feature_scaler = StandardScaler()  # 初始化标准化器
scaled_company_profiles = feature_scaler.fit_transform(exercise8_feature_matrix)  # 对公司财务特征进行标准化

hierarchical_linkage_matrix = linkage(scaled_company_profiles, method='ward')  # 计算公司级Ward链接矩阵
plt.figure(figsize=(10, 6))  # 创建画布
dendrogram(hierarchical_linkage_matrix, no_labels=True, color_threshold=None)  # 公司较多时隐藏叶标签但保留全部观测
plt.title('长三角上市公司财务画像 Ward 树', fontproperties='Source Han Serif SC')  # 图标题
plt.ylabel('Ward 合并距离', fontproperties='Source Han Serif SC')  # 明确纵轴含义
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图形

图 12.14

表 12.6: 习题8：公司级 Ward 与 KMeans 同 K 比较

exercise8_silhouettes = {}  # 保存每个可行K的KMeans轮廓系数
for cluster_count in exercise8_valid_k:  # 遍历事前候选
    candidate_labels = KMeans(n_clusters=cluster_count, n_init=20, random_state=42).fit_predict(scaled_company_profiles)  # 在同一矩阵拟合候选
    exercise8_silhouettes[cluster_count] = silhouette_score(scaled_company_profiles, candidate_labels)  # 记录选择证据
exercise8_selected_k = sorted(exercise8_silhouettes, key=lambda cluster_count: (-exercise8_silhouettes[cluster_count], cluster_count))[0]  # 并列时选择更小K
exercise8_kmeans_labels = KMeans(n_clusters=exercise8_selected_k, n_init=50, random_state=42).fit_predict(scaled_company_profiles)  # 拟合冻结K的最终模型
exercise8_ward_labels = fcluster(hierarchical_linkage_matrix, t=exercise8_selected_k, criterion='maxclust')  # 把Ward树切为同样簇数
exercise8_comparison_ari = adjusted_rand_score(exercise8_ward_labels, exercise8_kmeans_labels)  # 用标签置换不变指标比较
exercise8_cluster_audit = exercise8_company_profiles[['order_book_id', 'province']].assign(ward_cluster=exercise8_ward_labels, kmeans_cluster=exercise8_kmeans_labels)  # 绑定公司键
exercise8_ward_province = exercise8_cluster_audit.groupby(['ward_cluster', 'province']).size().rename('province_n').reset_index()
exercise8_ward_province['cluster_total_n'] = exercise8_ward_province.groupby('ward_cluster')['province_n'].transform('sum')
exercise8_kmeans_province = exercise8_cluster_audit.groupby(['kmeans_cluster', 'province']).size().rename('province_n').reset_index()
exercise8_kmeans_province['cluster_total_n'] = exercise8_kmeans_province.groupby('kmeans_cluster')['province_n'].transform('sum')
exercise8_summary = pd.DataFrame([{'n_raw': len(exercise8_raw_roster), 'n_analysis': len(exercise8_company_profiles), 'n_excluded': int(exercise8_exclusion_ledger['reason'].ne('included').sum()), 'selected_k': exercise8_selected_k, 'ward_kmeans_ari': exercise8_comparison_ari, 'contract_sha256': exercise8_contract_sha256}])
for exercise8_name, exercise8_frame in {'assignments': exercise8_cluster_audit, 'ward-province': exercise8_ward_province, 'kmeans-province': exercise8_kmeans_province, 'summary': exercise8_summary}.items():
    exercise8_write_same(exercise8_artifact_dir / f'ex8-{exercise8_name}-v1.csv', exercise8_stable_csv(exercise8_frame))
print({'n': len(exercise8_company_profiles), 'valid_k': exercise8_valid_k, 'silhouette_by_k': exercise8_silhouettes, 'selected_k': exercise8_selected_k, 'ward_kmeans_ari': exercise8_comparison_ari, 'contract_sha256': exercise8_contract_sha256})
print(exercise8_ward_province.to_string(index=False))
print(exercise8_kmeans_province.to_string(index=False))

列表 12.39: 习题8：稳定完成清单与追加式成功终态

exercise8_required_artifacts = ['ex8-contract-v1.json', 'ex8-raw-roster-v1.csv', 'ex8-analysis-roster-v1.csv', 'ex8-exclusions-v1.csv', 'ex8-assignments-v1.csv', 'ex8-ward-province-v1.csv', 'ex8-kmeans-province-v1.csv', 'ex8-summary-v1.csv']  # 冻结完成验收的全部机器制品
exercise8_missing_artifacts = [name for name in exercise8_required_artifacts if not (exercise8_artifact_dir / name).is_file()]  # 核对稳定清单前的完整性
if exercise8_missing_artifacts:  # 任一必需CSV或合同缺失都不能称为完成
    exercise8_stop('completed_artifact_missing', files=exercise8_missing_artifacts)  # 保留本次失败尝试
exercise8_artifact_hashes = {name: exercise8_stream_sha256(exercise8_artifact_dir / name) for name in exercise8_required_artifacts}  # 重读实际字节计算SHA-256
exercise8_completed_manifest = {'schema_version': '1', 'status': 'completed', 'artifacts': exercise8_artifact_hashes}  # 排除尝试ID与现场指标以保持稳定
exercise8_manifest_bytes = (json.dumps(exercise8_completed_manifest, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode('utf-8')  # 固定完成清单字节
exercise8_manifest_path = exercise8_artifact_dir / 'ex8-completed-manifest-v1.json'  # 与追加式尝试状态分离
exercise8_manifest_bytes = exercise8_write_same(exercise8_manifest_path, exercise8_manifest_bytes)  # 原子冻结并重读完成清单
exercise8_manifest_sha256 = exercise8_stream_sha256(exercise8_manifest_path)  # 只哈希落盘清单字节
exercise8_success_record = {'schema_version': '1', 'attempt_id': exercise8_attempt_id, 'status': 'completed', 'manifest_sha256': exercise8_manifest_sha256, 'contract_sha256': exercise8_contract_sha256, 'n': len(exercise8_company_profiles), 'valid_k': exercise8_valid_k, 'selected_k': exercise8_selected_k, 'ward_kmeans_ari': exercise8_comparison_ari}  # 尝试遥测只引用落盘稳定清单
exercise8_success_path, exercise8_success_bytes = exercise8_write_terminal(exercise8_success_record)  # 追加本次成功而不改写历史停止尝试
print({'attempt': exercise8_success_path.name if exercise8_success_path else None, 'attempt_sha256': hashlib.sha256(exercise8_success_bytes).hexdigest(), 'manifest': exercise8_completed_manifest})  # 分别输出尝试与稳定清单

本题的可验收输出是原始公司名册、分析名册、逐公司排除账本、规则合同及其输入 hash，另含 图 12.14、表 12.6、公司标签表，以及 Ward 与 KMeans 各自的“簇规模—省份构成”表。ex8-completed-manifest-v1.json 列出 ex8-contract-v1.json、全部七份必需 CSV 及从落盘字节重算的 SHA-256；ex8-attempt-<attempt_id>-v1.json 只追加本次尝试的 stopped 或 completed 终态。因此旧的停止尝试会被保留，换用新 attempt_id 后仍可成功，不会与一个全局状态文件冲突。Ward 与 KMeans 使用同一个标准化公司矩阵和同一个现场选择的 \(K\)；ARI 只衡量两种算法标签的一致程度。省份构成是事后描述，不参与选 \(K\)，也不能把上市公司财务画像推广为完整省域经济或“真实”省域类型。

12.11.3 理论题解答

9. 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\) 为对应特征值。

10. 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。

11. 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 不满足此推导。

12. 主题总结（写作建议）

    • t-SNE/UMAP：强调‘局部结构保持’与‘可视化解释限制’，并讨论超参数敏感性。
    • DBSCAN/谱聚类：强调密度可发现任意形状簇、谱方法适合图结构与非凸簇。
    • GMM：强调软聚类与模型选择（BIC/AIC），以及协方差结构对簇形状的影响。

12.12 章末闭环

逐项目标自检：说明标准化为何改变距离几何；从解释方差规则选择 PCA 维数并解释载荷；区分二维展示坐标与实际聚类输入；按 \(n\ge5K\)、有效簇与有限轮廓筛选候选；用多初值稳定性区分算法波动与结构证据。另设不可抵消的 M12 门槛：必须提交课程级多行业 universe、单行业 analysis subset、清洗后 clustering roster 的包含关系，以及清洗前冻结合同、排除账本、六份输出和 handoff 的逐字节 hash；缺少任一项均不得进入 M13。

禁用情境：输入快照、总体或逐文件指纹不完整时，不应报告实样本画像；样本量不足、有效簇少于 \(K\) 或稳定性失败时，不应强行命名簇。常见误区是把 PC1/PC2 图上的分离当作全维结构，或把簇编号、行业标签和因果机制相互替代。

无提示检索：1）为何 PCA 前是否标准化会改变答案？2）二维图只解释了哪些信息？3）轮廓系数最高是否足以证明存在真实类别？

展开检索反馈与学习决策

1）标准化改变变量相对尺度，从而改变协方差/相关矩阵的特征方向和距离。2）只展示所画两个主成分承载的投影，不代表未画维度。3）不能；它是候选几何与样本内分离证据，还需样本量、稳定性、外部解释和失败审计。概念题满分 3 分；任一错项先回到本章“距离与标准化”“主成分分析”或“基于本地数据的行业财务因子 PCA + 聚类”的相应公式和代码锚点补修，再用改变尺度、旋转投影或更换候选 \(K\) 的异形题复测，答对才重入。M12 门槛单独满分 4 分（合同、三层名册、排除账本、证据包各 1 分）且必须 4 分；任何缺项均回到 小节 12.9 重建并从 fresh kernel 核验全部 hash，不得用概念题得分抵消。

陌生迁移：换一个冻结公司横截面，先登记逐文件 SHA-256、总体和排除，再预注册 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.