什么是无监督学习?
在前面的章节中,我们研究的都是监督学习——为每个观测 \(i\) 拥有特征向量 \(x_i\) 和响应变量 \(y_i\)。
无监督学习面对的是完全不同的局面:
- 只有特征 \(X_1, X_2, \ldots, X_p\),没有响应变量 \(Y\)
- 目标不是预测,而是发现数据中隐藏的结构和模式
- 两大核心工具:主成分分析 (PCA) 和 聚类分析 (Clustering)
无监督学习的三大核心任务
| 降维 |
减少特征数,保留信息 |
PCA、t-SNE |
因子构建、数据压缩 |
| 聚类 |
发现自然分组 |
K-Means、层次聚类 |
客户分群、行业分类 |
| 异常检测 |
找出离群点 |
孤立森林、LOF |
欺诈检测、异常交易 |
与监督学习的关键区别:
- 监督学习:有”标准答案”(标签),可以算准确率
- 无监督学习:没有标准答案,只能靠”内部一致性”评估
挑战:如何评价结果的好坏?
- 聚类成3组合适还是5组?没有唯一正确答案
- 需要结合业务知识判断
无监督学习的独特挑战
与监督学习相比,无监督学习面临更深层次的困难:
| 目标 |
明确(最小化预测误差) |
模糊(探索结构) |
| 评估 |
交叉验证、MSE、AUC |
没有统一标准 |
| 主观性 |
较低 |
较高(结果解读因人而异) |
| 验证 |
测试集验证 |
需要领域知识判断 |
为什么这里不能直接谈准确率、召回率和 AUC?
这几个指标在监督学习里非常重要,但在无监督学习里不能直接照搬,原因很简单:
- 准确率需要真实标签,才能知道模型“猜对了多少”。
- 召回率 / 精确率需要先知道谁是真正的正类、谁是负类。
- AUC需要把正负类对象做排序比较。
而无监督学习恰恰一开始就没有这些标签。
所以无监督学习真正能做的是:
- 看结构是否稳定。
- 看组内是否更相似、组间是否更分离。
- 看结果能否被业务知识解释。
这也是为什么监督学习更像“判卷”,而无监督学习更像“看地图找结构”。两者不是同一套评分体系。
为什么没有标签仍然能学到东西?
无监督学习不是在回答’谁对谁错’,而是在回答:
- 数据里有没有高相关的变量组合?
- 样本之间有没有自然形成的群体?
- 有没有一些观测在结构上明显离群?
| 能不能预测 \(Y\)? |
数据内部有什么结构? |
| 哪些特征最能解释标签? |
哪些方向最能解释差异? |
| 模型准确率是多少? |
分组和低维表示是否有业务意义? |
所以:无监督学习的价值不在于替代预测,而在于为后续预测、解释和策略设计提供结构化线索。
PCA 的数学定义
给定 \(n\) 个观测、\(p\) 个特征的数据矩阵 \(\mathbf{X}\)(已中心化),第一主成分定义为:
\[Z_1 = \phi_{11}X_1 + \phi_{21}X_2 + \cdots + \phi_{p1}X_p\]
其中载荷向量 \(\boldsymbol{\phi}_1 = (\phi_{11}, \phi_{21}, \ldots, \phi_{p1})^T\) 满足:
\[\max_{\boldsymbol{\phi}_1} \text{Var}(Z_1) \quad \text{s.t.} \quad \|\boldsymbol{\phi}_1\|^2 = 1\]
- 目标:使投影后的方差最大化
- 约束:载荷向量为单位向量(防止无穷大解)
- 解:\(\boldsymbol{\phi}_1\) 是协方差矩阵 \(\mathbf{S}\) 的最大特征值对应的特征向量
主成分、载荷与得分分别是什么?
PCA 里最容易混淆的三个词是:主成分、载荷、得分。
| 主成分 |
一个新的综合方向 |
新坐标轴 |
| 载荷 |
原变量在主成分上的权重 |
这个新坐标轴由哪些原变量拼出来 |
| 得分 |
每个样本投影到主成分上的坐标 |
每家公司在新坐标轴上的位置 |
举例来说:
- 如果 PC1 主要由 ROA 和净利润率组成,那么它可以被解释成’盈利能力方向’
- 某家公司在 PC1 上得分很高,说明它在这个’盈利能力方向’上位置靠前
一句话记忆:载荷决定主成分’长什么样’,得分决定样本’站在哪里’。
PCA的优化目标与求解
PCA的两种等价表述:
表述1:最大方差
\[
\max_{\phi} \text{Var}(Z_1) = \max_{\phi} \frac{1}{n} \sum_{i=1}^n \left(\sum_{j=1}^p \phi_{j1} x_{ij}\right)^2
\]
约束:\(\|\phi_1\| = 1\)(单位向量)
表述2:最小重构误差
\[
\min_{\phi} \sum_{i=1}^n \|x_i - \hat{x}_i\|^2
\]
其中 \(\hat{x}_i\) 是 \(x_i\) 在主成分方向上的投影
这两个目标是等价的!
- 最大方差方向 = 最小重构误差方向
- 求解:对协方差矩阵做特征值分解
SVD 视角:PCA 的高效计算
PCA 的计算可以通过奇异值分解 (SVD) 高效实现。
对中心化数据矩阵 \(\mathbf{X}\) 进行 SVD:
\[\mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^T\]
- \(\mathbf{V}\) 的列 = 主成分载荷向量(即特征向量)
- \(\boldsymbol{\Sigma}\) 的对角元素 \(\sigma_i\) 与特征值的关系:\(\lambda_i = \sigma_i^2 / n\)
- 方差解释比:\(\frac{\lambda_k}{\sum_{j=1}^p \lambda_j} = \frac{\sigma_k^2}{\sum_{j=1}^p \sigma_j^2}\)
数据准备:长三角上市公司财务指标
import pandas as pd # 导入pandas库用于数据框操作
import numpy as np # 导入numpy库用于数值计算
import os # 导入os模块用于跨平台路径处理
# 根据操作系统选择数据目录
DATA_DIR = 'C:/qiufei/data' if os.name == 'nt' else '/home/ubuntu/r2_data_mount/data' # 跨平台数据路径
path_financial = os.path.join(DATA_DIR, 'stock/financial_statement.h5') # 财务报表数据路径
path_basic = os.path.join(DATA_DIR, 'stock/stock_basic_data.h5') # 上市公司基本信息路径
financial_statement_raw = pd.read_hdf(path_financial) # 读取财务报表数据
stock_basic_data = pd.read_hdf(path_basic) # 读取上市公司基本信息
# 定义长三角省份列表(上海、江苏、浙江、安徽)
yangtze_river_delta_provinces = ['上海市', '江苏省', '浙江省', '安徽省'] # 长三角四省市(与数据中省份全称匹配)
# 筛选长三角地区上市公司
yrd_stock_codes = stock_basic_data[
stock_basic_data['province'].isin(yangtze_river_delta_provinces) # 按省份筛选
]['order_book_id'].tolist() # 提取股票代码列表
构建五大财务指标
# 筛选长三角公司、最近期的财务数据
yrd_financial_data = financial_statement_raw[
financial_statement_raw['order_book_id'].isin(yrd_stock_codes) # 筛选长三角公司
].copy() # 创建副本防止修改原数据
# 只保留2023年之后的财务数据
yrd_financial_data = yrd_financial_data[
yrd_financial_data['quarter'] >= '2023q1' # 筛选最近期数据
] # 确保数据时效性
# 计算五个核心财务指标
yrd_financial_data['Log_Assets'] = np.log(yrd_financial_data['total_assets'].clip(lower=1)) # 对数总资产(取对数消除量纲差异)
yrd_financial_data['Debt_Ratio'] = ( # 资产负债率
yrd_financial_data['current_liabilities'].fillna(0) + yrd_financial_data['non_current_liabilities'].fillna(0)
) / yrd_financial_data['total_assets'] # 负债总额除以总资产
yrd_financial_data['Net_Margin'] = ( # 净利润率
yrd_financial_data['net_profit'] / yrd_financial_data['revenue'].replace(0, np.nan) # 净利润除以营业收入
) # 衡量盈利能力
yrd_financial_data['Asset_Turnover'] = ( # 总资产周转率
yrd_financial_data['revenue'] / yrd_financial_data['total_assets'] # 营业收入除以总资产
) # 衡量资产使用效率
yrd_financial_data['ROA'] = ( # 总资产收益率
yrd_financial_data['net_profit'] / yrd_financial_data['total_assets'] # 净利润除以总资产
) # 衡量资产盈利能力
# 每家公司取最近一期财报
pca_feature_columns = ['Log_Assets', 'Debt_Ratio', 'Net_Margin', 'Asset_Turnover', 'ROA'] # PCA分析的五个特征
yrd_latest_financial = yrd_financial_data.sort_values('quarter').groupby('order_book_id').last().reset_index() # 取最新一期
yrd_pca_data = yrd_latest_financial[['order_book_id'] + pca_feature_columns].dropna() # 提取特征列并去除缺失值
为什么 PCA 前几乎总要先标准化?
PCA 会优先寻找方差最大的方向,因此如果变量量纲差异过大,量纲大的变量会自动’抢走话语权’。
一个典型对比:
- 总资产可能以十亿、百亿计
- ROA 往往只是几个百分点
如果不标准化:
- PCA 的第一主成分很可能几乎只是在描述’企业规模’
- 盈利能力、运营效率等变量会被压缩得几乎看不见
如果先做标准化:
- 每个变量都先变成均值为 0、标准差为 1 的可比尺度
- PCA 才是在比较’结构差异’,而不是比较’谁的数值单位更大’
数据标准化与 Winsorization
from sklearn.preprocessing import StandardScaler # 导入标准化工具
# 1% 和 99% 分位数 Winsorization(缩尾处理)
for feature_col in pca_feature_columns: # 遍历每个特征列
lower_bound = yrd_pca_data[feature_col].quantile(0.01) # 计算1%分位数作为下界
upper_bound = yrd_pca_data[feature_col].quantile(0.99) # 计算99%分位数作为上界
yrd_pca_data[feature_col] = yrd_pca_data[feature_col].clip(lower_bound, upper_bound) # 将极端值截断到上下界
# 标准化(均值为0,标准差为1)
feature_scaler = StandardScaler() # 初始化标准化器
standardized_features = feature_scaler.fit_transform(yrd_pca_data[pca_feature_columns]) # 拟合并转换数据
print(f'标准化后的数据形状:{standardized_features.shape[0]} 家公司 × {standardized_features.shape[1]} 个特征') # 输出数据维度
标准化后的数据形状:2015 家公司 × 5 个特征
为什么要标准化?
- 不同指标量纲差异巨大(总资产数十亿,ROA不到1%)
- 不标准化时,PCA 会被量纲最大的变量主导
- 标准化确保每个变量在 PCA 中获得公平权重
PCA 拟合与碎石图
Code
from sklearn.decomposition import PCA # 导入PCA模型
import matplotlib.pyplot as plt # 导入绑图库
plt.rcParams['font.family'] = ['Source Han Serif SC', 'Noto Serif CJK SC', 'SimSun'] # 设置中文字体优先级
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
pca_model = PCA() # 初始化PCA模型(保留所有主成分)
pca_scores = pca_model.fit_transform(standardized_features) # 拟合PCA并获取主成分得分
explained_variance_ratio = pca_model.explained_variance_ratio_ # 提取各主成分的方差解释比
cumulative_variance_ratio = np.cumsum(explained_variance_ratio) # 计算累积方差解释比
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) # 创建1行2列子图
# 左图:碎石图
axes[0].bar(range(1, len(explained_variance_ratio) + 1), explained_variance_ratio, color='#008080', alpha=0.8) # 柱状图
axes[0].set_xlabel('主成分编号') # 设置x轴标签
axes[0].set_ylabel('方差解释比') # 设置y轴标签
axes[0].set_title('碎石图 (Scree Plot)') # 设置标题
# 右图:累积方差解释比
axes[1].plot(range(1, len(cumulative_variance_ratio) + 1), cumulative_variance_ratio, 'o-', color='#E3120B') # 折线图
axes[1].axhline(y=0.80, color='gray', linestyle='--', alpha=0.5, label='80% 阈值') # 80%阈值线
axes[1].set_xlabel('主成分个数') # 设置x轴标签
axes[1].set_ylabel('累积方差解释比') # 设置y轴标签
axes[1].set_title('累积方差解释比') # 设置标题
axes[1].legend() # 显示图例
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
print(f'PC1 方差解释比: {explained_variance_ratio[0]:.2%}') # 输出PC1方差解释比
print(f'PC1+PC2 累积: {cumulative_variance_ratio[1]:.2%}') # 输出前两个PC累积方差
print(f'PC1+PC2+PC3 累积: {cumulative_variance_ratio[2]:.2%}') # 输出前三个PC累积方差
PC1 方差解释比: 39.07%
PC1+PC2 累积: 66.18%
PC1+PC2+PC3 累积: 85.63%
方差解释比到底在解释什么?
方差解释比回答的问题不是’解释了多少经济意义’,而是:
这个主成分保留了原始数据总波动中的多少比例?
例如:
- 如果 PC1 的方差解释比是 40%
- 意味着把数据投影到 PC1 这一条轴上时,原始五维数据中约 40% 的总体变异仍被保留下来
需要特别注意:
- 解释比高,说明信息保留多
- 但不必然说明业务解释更重要
- 有时解释比不高的主成分,反而对应着很关键的行业结构或异常模式
因此:方差解释比是降维质量指标,不是业务价值的自动排名。
载荷图 (Biplot):解读主成分含义
Code
fig, ax = plt.subplots(figsize=(8, 6)) # 创建单张图表
# 提取前两个主成分的载荷
loadings_pc1 = pca_model.components_[0] # PC1的载荷向量
loadings_pc2 = pca_model.components_[1] # PC2的载荷向量
# 绘制载荷箭头
for i, feature_name in enumerate(pca_feature_columns): # 遍历每个特征
ax.annotate('', xy=(loadings_pc1[i], loadings_pc2[i]), xytext=(0, 0), # 从原点到载荷点画箭头
arrowprops=dict(arrowstyle='->', color='#E3120B', lw=2)) # 红色箭头
ax.text(loadings_pc1[i] * 1.15, loadings_pc2[i] * 1.15, feature_name, # 在箭头端点标注特征名
fontsize=11, ha='center', color='#2C3E50', fontweight='bold') # 加粗标注
ax.axhline(y=0, color='gray', linestyle='--', alpha=0.3) # 水平参考线
ax.axvline(x=0, color='gray', linestyle='--', alpha=0.3) # 垂直参考线
ax.set_xlabel(f'PC1 ({explained_variance_ratio[0]:.1%})') # x轴标签(含方差解释比)
ax.set_ylabel(f'PC2 ({explained_variance_ratio[1]:.1%})') # y轴标签(含方差解释比)
ax.set_title('PCA 载荷图') # 设置标题
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
载荷解读:
- PC1(盈利质量维度):Net_Margin 和 ROA 载荷较高 → 反映企业盈利能力
- PC2(规模与杠杆维度):Log_Assets 和 Debt_Ratio 载荷较高 → 反映企业规模和负债水平
主成分得分按行业分布
Code
# 合并行业信息
yrd_pca_data_with_industry = yrd_pca_data.merge( # 合并行业分类
stock_basic_data[['order_book_id', 'industry_name']], on='order_book_id', how='left' # 左连接
) # 获取每家公司的行业归属
# 统计各行业频次,选择前6大行业
top_six_industries = yrd_pca_data_with_industry['industry_name'].value_counts().head(6).index.tolist() # 前6大行业名称
plot_mask = yrd_pca_data_with_industry['industry_name'].isin(top_six_industries) # 筛选属于前6行业的公司
fig, ax = plt.subplots(figsize=(10, 6)) # 创建图表
industry_color_palette = ['#E3120B', '#008080', '#F0A700', '#2C3E50', '#8E9EAA', '#6A5ACD'] # 6种行业颜色
for i, industry_name in enumerate(top_six_industries): # 遍历前6大行业
industry_mask = yrd_pca_data_with_industry['industry_name'] == industry_name # 筛选当前行业
valid_indices = yrd_pca_data_with_industry[industry_mask].index # 获取行业内公司的索引
# 获取这些公司在PCA得分矩阵中的位置
pca_row_indices = [yrd_pca_data.index.get_loc(idx) for idx in valid_indices if idx in yrd_pca_data.index] # 映射索引
ax.scatter(pca_scores[pca_row_indices, 0], pca_scores[pca_row_indices, 1], # 绘制散点
c=industry_color_palette[i], label=industry_name, alpha=0.6, s=30) # 按行业着色
ax.set_xlabel(f'PC1 ({explained_variance_ratio[0]:.1%})') # x轴标签
ax.set_ylabel(f'PC2 ({explained_variance_ratio[1]:.1%})') # y轴标签
ax.set_title('主成分得分的行业分布') # 设置标题
ax.legend(fontsize=9, ncol=2) # 显示两列图例
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
PCA 实践要点
如何选择主成分个数?
| 碎石图 (Scree Plot) |
寻找”肘部”(方差急剧下降处) |
直观判断 |
| 累积方差比 |
选择累积达到80%~90%的点 |
通用标准 |
| Kaiser 准则 |
保留特征值 > 1 的成分 |
标准化数据 |
PCA 的局限性:
- 线性假设:只能捕获线性关系,非线性结构需要 t-SNE/UMAP
- 解释困难:主成分是原始变量的线性组合,业务含义不直接
- 方差 ≠ 信息:最大方差方向未必是最有业务价值的方向
PCA在量化投资中的应用
应用1:因子降维
- 上百个基本面指标 → PCA提取几个综合因子
- 类似于构建”超级因子”
应用2:去噪
- 保留前几个主成分(信号),丢弃后面的(噪声)
- 协方差矩阵去噪 → 更稳定的投资组合优化
应用3:异常检测
- 在PCA空间中距离原点很远的样本可能是异常值
- 利用重构误差检测异常
注意事项:
- PCA是线性方法:只能捕捉线性关系
- 对量纲敏感:必须先标准化
- 主成分不一定有业务含义:需要通过载荷矩阵解读
从 PCA 到聚类:EM 算法的视角
EM 算法 是理解聚类(以及许多其他统计方法)的统一框架:
- E 步(期望步):给定当前参数估计,计算每个观测属于各类别的后验概率
\[\gamma_{ik} = \frac{\pi_k \mathcal{N}(x_i | \mu_k, \Sigma_k)}{\sum_{j=1}^K \pi_j \mathcal{N}(x_i | \mu_j, \Sigma_j)}\]
- M 步(最大化步):利用后验概率更新参数
- 更新混合比例:\(\hat{\pi}_k = \frac{1}{n}\sum_{i=1}^n \gamma_{ik}\)
- 更新均值:\(\hat{\mu}_k = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}\)
- 更新协方差:\(\hat{\Sigma}_k = \frac{\sum_i \gamma_{ik}(x_i - \hat{\mu}_k)(x_i - \hat{\mu}_k)^T}{\sum_i \gamma_{ik}}\)
读法建议:先把 \(\gamma_{ik}\) 看成“第 \(i\) 个样本对第 \(k\) 类的责任权重”,再看均值和协方差更新,你就会发现 M 步本质上只是“按责任权重重新做加权平均”。
EM算法的数学框架
EM算法的通用流程:
假设数据来自 \(K\) 个”隐含”组,每个数据点的组归属是未观测的隐变量
E步(Expectation):
\[
\gamma_{ik} = P(z_i = k | x_i, \theta^{(t)}) = \frac{\pi_k f(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^K \pi_l f(x_i | \mu_l, \Sigma_l)}
\]
- \(\gamma_{ik}\):样本 \(i\) 属于第 \(k\) 组的后验概率
M步(Maximization):
\[
\mu_k^{(t+1)} = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}, \quad \pi_k^{(t+1)} = \frac{\sum_i \gamma_{ik}}{n}
\]
迭代直到收敛
K-Means是EM的”硬分配”特例:\(\gamma_{ik} \in \{0, 1\}\)(不是概率而是确定性分配)
实务提醒:EM 的优势是能表达“边界模糊”的样本,代价是更依赖初始化,也更容易停在局部最优;所以解释结果时,不能只看一次运行的答案。
K-Means 是 EM 的特例
K-Means 可以被视为高斯混合模型 EM 算法的硬分配版本:
| 分配方式 |
软分配(概率 \(\gamma_{ik} \in [0,1]\)) |
硬分配(\(\gamma_{ik} \in \{0,1\}\)) |
| 均值更新 |
加权平均 |
算术平均 |
| 簇形状 |
椭圆形(不同协方差) |
球形(等协方差) |
| 目标函数 |
对数似然 |
组内平方和 (WCSS) |
K-Means 的目标函数:
\[\min_{C_1,\ldots,C_K} \sum_{k=1}^K \sum_{i \in C_k} \| x_i - \mu_k \|^2\]
聚类结果为什么没有唯一标准答案?
聚类和分类最大的不同在于:
- 分类通常有外部标签可核对
- 聚类没有一个预先写好的’正确分组’
因此,聚类结果的好坏取决于三个层面:
| 统计层面 |
组内够不够紧、组间够不够分开? |
| 业务层面 |
分出来的组有没有经济含义? |
| 决策层面 |
这样的分组能不能支持后续策略? |
所以:
- 把公司分成 3 组,可能更利于教学解释
- 分成 6 组,可能更利于投资细分
- 二者不一定谁对谁错,而是服务于不同目标
K-Means 实战:最优 K 值选择
Code
from sklearn.cluster import KMeans # 导入K-Means聚类模型
from sklearn.metrics import silhouette_score # 导入轮廓系数
# 使用前3个主成分作为聚类特征
pca_features_for_clustering = pca_scores[:, :3] # 提取前3个PC得分
k_range = range(2, 10) # 候选K值范围:2到9
inertia_values = [] # 存储各K值的组内平方和 (SSE/Inertia)
silhouette_values = [] # 存储各K值的轮廓系数
for k in k_range: # 遍历每个候选K值
kmeans_model = KMeans(n_clusters=k, random_state=42, n_init=10) # 初始化K-Means模型
cluster_labels = kmeans_model.fit_predict(pca_features_for_clustering) # 拟合并预测聚类标签
inertia_values.append(kmeans_model.inertia_) # 记录组内平方和
silhouette_values.append(silhouette_score(pca_features_for_clustering, cluster_labels)) # 记录轮廓系数
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) # 创建1行2列子图
axes[0].plot(list(k_range), inertia_values, 'o-', color='#E3120B', linewidth=2) # 肘部法折线图
axes[0].set_xlabel('聚类数 K') # x轴标签
axes[0].set_ylabel('组内平方和 (SSE)') # y轴标签
axes[0].set_title('肘部法 (Elbow Method)') # 标题
axes[1].plot(list(k_range), silhouette_values, 's-', color='#008080', linewidth=2) # 轮廓系数折线图
axes[1].set_xlabel('聚类数 K') # x轴标签
axes[1].set_ylabel('轮廓系数') # y轴标签
axes[1].set_title('轮廓系数法 (Silhouette Method)') # 标题
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
聚类结果解读:K=4 的业务含义
# 以K=4执行最终聚类
optimal_kmeans = KMeans(n_clusters=4, random_state=42, n_init=10) # K=4的KMeans模型
final_cluster_labels = optimal_kmeans.fit_predict(pca_features_for_clustering) # 获取聚类标签
yrd_pca_data_copy = yrd_pca_data.copy() # 创建副本
yrd_pca_data_copy['Cluster'] = final_cluster_labels # 将聚类标签添加到数据框
# 计算各聚类在原始特征上的均值画像
cluster_profile = yrd_pca_data_copy.groupby('Cluster')[pca_feature_columns].mean().round(4) # 各聚类的特征均值
cluster_profile.index = [f'聚类 {i}' for i in cluster_profile.index] # 重命名索引
cluster_profile # 显示聚类画像表
聚类在 PCA 空间中的可视化
Code
fig, ax = plt.subplots(figsize=(9, 6)) # 创建图表
cluster_colors = ['#E3120B', '#008080', '#F0A700', '#2C3E50'] # 4个聚类的颜色
for cluster_id in range(4): # 遍历每个聚类
cluster_mask = final_cluster_labels == cluster_id # 筛选当前聚类的样本
ax.scatter(pca_scores[cluster_mask, 0], pca_scores[cluster_mask, 1], # 在PC1-PC2空间绑制散点
c=cluster_colors[cluster_id], label=f'聚类 {cluster_id}', alpha=0.6, s=30) # 按聚类着色
ax.set_xlabel(f'PC1 ({explained_variance_ratio[0]:.1%})') # x轴标签
ax.set_ylabel(f'PC2 ({explained_variance_ratio[1]:.1%})') # y轴标签
ax.set_title('K-Means 聚类结果(PCA空间投影)') # 标题
ax.legend() # 显示图例
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
层次聚类与树状图
Code
from scipy.cluster.hierarchy import dendrogram, linkage # 导入层次聚类工具
np.random.seed(42) # 设置随机种子
sample_size = min(50, len(standardized_features)) # 取样本量(最多50家)
sample_indices = np.random.choice(len(standardized_features), sample_size, replace=False) # 随机抽样索引
sample_features = standardized_features[sample_indices] # 提取样本特征
# 合并行业名称用于标签
sample_stock_codes = yrd_pca_data.iloc[sample_indices]['order_book_id'].values # 获取样本公司代码
sample_industry_names = stock_basic_data.set_index('order_book_id').loc[sample_stock_codes, 'industry_name'].values # 获取行业名
# Ward 法层次聚类
linkage_matrix = linkage(sample_features, method='ward') # 使用Ward方法计算链接矩阵
fig, ax = plt.subplots(figsize=(14, 6)) # 创建图表
dendrogram(linkage_matrix, labels=sample_industry_names, leaf_rotation=90, leaf_font_size=8, ax=ax) # 绘制树状图
ax.set_title('层次聚类树状图(Ward法)') # 标题
ax.set_ylabel('距离') # y轴标签
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
层次聚类的”家谱树”隐喻:
- 树的底部:每家公司是一个叶子节点
- 向上合并:最相似的公司先合并,形成小组
- 树的高度:反映组间差异程度——高度越大,差异越显著
距离度量为什么会改变聚类结果?
聚类并不是在问’谁看起来像谁’,而是在问:
按照我们定义的距离,谁离谁更近?
不同距离,代表不同的’相似’含义:
- 欧氏距离:数值绝对水平接近
- 相关距离:走势形态接近
- 曼哈顿距离:逐维偏差的累计更稳健
以股票为例:
- 两只股票价格水平一个 10 元、一个 100 元,但如果涨跌节奏同步,相关距离会认为它们很像
- 若只看欧氏距离,它们可能被认为差得很远
因此:距离不是技术细节,而是在用数学语言定义’什么叫相似’。
层次聚类的距离度量选择
层次聚类的三要素:距离度量 + 连接方式 + 截断高度
常用连接方式(Linkage):
| 单连接 |
两簇最近点的距离 |
易产生”链状”聚类 |
| 全连接 |
两簇最远点的距离 |
产生紧凑球形聚类 |
| 平均连接 |
两簇所有点对平均距离 |
均衡,推荐 |
| Ward |
合并后方差增量最小 |
最常用 |
距离度量的选择:
- 欧氏距离:适用于数值型特征,对量纲敏感
- 相关距离 \(1 - |r|\):适用于时间序列聚类
- Manhattan距离:对异常值更鲁棒
实操建议:结构化数据用Ward连接 + 欧氏距离,时间序列用平均连接 + 相关距离
聚类方法对比与选择
| 簇形状 |
球形 |
任意 |
任意 |
| 需指定 K |
是 |
否(后期切割) |
否 |
| 大规模数据 |
适合 |
不适合 |
适合 |
| 噪声处理 |
无 |
无 |
能识别噪声点 |
| 可解释性 |
中等 |
高(树状图) |
中等 |
from sklearn.metrics import calinski_harabasz_score # 导入CH指数(方差比准则)
# 使用 Calinski-Harabasz 指数评估 K=4 的聚类质量
ch_index_value = calinski_harabasz_score(pca_features_for_clustering, final_cluster_labels) # 计算CH指数
print(f'K=4 的 Calinski-Harabasz 指数: {ch_index_value:.2f}') # 输出CH指数(越大越好)
print('(CH指数越大,表示聚类越紧凑且分离度越高)') # 解读说明
K=4 的 Calinski-Harabasz 指数: 727.57
(CH指数越大,表示聚类越紧凑且分离度越高)
基于收益率相关性的股票聚类
# 选取5大行业各6只代表性股票
industry_stock_mapping = { # 行业-股票映射字典
'银行': ['601398.XSHG', '601288.XSHG', '600036.XSHG', '601818.XSHG', '600000.XSHG', '601166.XSHG'], # 工商银行、农业银行、招商银行等
'电子': ['002415.XSHE', '603986.XSHG', '600703.XSHG', '002049.XSHE', '300750.XSHE', '688008.XSHG'], # 海康威视等电子企业
'食品饮料': ['600519.XSHG', '000858.XSHE', '603369.XSHG', '002304.XSHE', '000568.XSHE', '600887.XSHG'], # 贵州茅台等食品企业
'非银金融': ['601318.XSHG', '600030.XSHG', '601688.XSHG', '600837.XSHG', '000776.XSHE', '601211.XSHG'], # 中国平安等金融企业
'房地产': ['000002.XSHE', '001979.XSHE', '600048.XSHG', '000069.XSHE', '600383.XSHG', '600340.XSHG'] # 万科A等房地产企业
} # 5行业 × 6只 = 30只股票
all_selected_stock_codes = [] # 存放所有选中股票代码
for stock_list in industry_stock_mapping.values(): # 展开所有行业的股票
all_selected_stock_codes.extend(stock_list) # 追加到总列表
# 读取后复权股价数据
path_price = os.path.join(DATA_DIR, 'stock/stock_price_post_adjusted.h5') # 股价文件路径
stock_price_data = pd.read_hdf( # 按日期和股票代码选择性读取所需子集
path_price, # 指定后复权股价文件
key='data', # 指定HDF5中的数据键
where=["date>=Timestamp('2023-01-01')", f'order_book_id={all_selected_stock_codes}'], # 只读取目标股票的2023年以来数据
columns=['close'] # 只保留收盘价这一必要列
).reset_index() # 将MultiIndex恢复为普通列
计算收益率相关距离矩阵
# 筛选选中股票的2023年以来数据
selected_stock_data = stock_price_data[
(stock_price_data['order_book_id'].isin(all_selected_stock_codes)) & # 筛选目标股票
(stock_price_data['date'] >= '2023-01-01') # 筛选2023年后数据
].copy() # 创建副本
# 计算日收益率
selected_stock_data['daily_return'] = selected_stock_data.groupby('order_book_id')['close'].pct_change() # 按股票分组计算收益率
# 构建收益率宽表(行=日期,列=股票代码)
daily_returns_wide = selected_stock_data.pivot_table( # 透视为宽表
index='date', columns='order_book_id', values='daily_return' # 行=日期,列=股票
).dropna(axis=1, thresh=100).dropna() # 去除数据不足的股票和缺失行
# 构建相关性距离矩阵: distance = 1 - correlation
correlation_matrix = daily_returns_wide.corr() # 计算股票间的Pearson相关系数矩阵
distance_matrix = 1 - correlation_matrix # 将相关性转换为距离(相关性越高,距离越近)
print(f'有效股票数: {distance_matrix.shape[0]}') # 输出矩阵维度
股票树状图:行业结构自动浮现
Code
from scipy.spatial.distance import squareform # 导入距离矩阵压缩工具
# 将方阵转换为压缩距离向量
condensed_distance = squareform(distance_matrix.values, checks=False) # 压缩距离矩阵
stock_linkage = linkage(condensed_distance, method='ward') # Ward法层次聚类
# 构建行业标签
industry_labels = [] # 存放行业简称标签
for stock_code in distance_matrix.columns: # 遍历每只股票
label_found = False # 标记是否找到行业
for ind_name, ind_stocks in industry_stock_mapping.items(): # 遍历行业映射
if stock_code in ind_stocks: # 匹配股票代码
industry_labels.append(f'{ind_name[:2]}') # 取行业名前两个字作为简称
label_found = True # 标记已找到
break # 跳出内层循环
if not label_found: # 未匹配时
industry_labels.append(stock_code[:6]) # 使用股票代码前6位作为标签
fig, ax = plt.subplots(figsize=(14, 6)) # 创建图表
dendrogram(stock_linkage, labels=industry_labels, leaf_rotation=0, leaf_font_size=10, ax=ax) # 绘制树状图
ax.set_title('基于收益率相关性的股票聚类树状图') # 标题
ax.set_ylabel('距离(1 - 相关系数)') # y轴标签
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
相关性热力图:对角线块状结构
Code
import seaborn as sns # 导入seaborn高级绑图库
# 使用seaborn的clustermap自动排列行列
cluster_heatmap = sns.clustermap( # 创建聚类热力图
correlation_matrix, # 输入相关系数矩阵
method='ward', # 使用Ward法聚类
cmap='coolwarm', # 冷暖色调(蓝=低相关,红=高相关)
vmin=-0.5, vmax=1, # 颜色范围
figsize=(10, 8), # 图表大小
linewidths=0.5 # 网格线宽度
) # 生成聚类热力图
cluster_heatmap.ax_heatmap.set_title('股票收益率相关性聚类热力图', fontsize=14, pad=60) # 设置标题
plt.show() # 显示图表
对角线上的明显”块状”结构表明:
- 同行业股票之间高度正相关(红色块)
- 不同行业之间相关性较弱(蓝色/白色区域)
- 应用价值:投资组合分散化、配对交易、风险管理