12 主成分分析

# 为“数据点的分布揭示了主要变化方向”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“数据点的分布揭示了主要变化方向”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt

# 固定随机种子，使相关估值指标的 PCA 几何机制示例可重复。
np.random.seed(42)
# 创建两个相关的变量
X1 = np.random.randn(20)
# 令第二维等于 0.8 倍第一维再加噪声，形成 PCA 可识别的斜向点云。
X2 = X1 * 0.8 + np.random.randn(20) * 0.5
# 提取 `X` 作为“数据点的分布揭示了主要变化方向”的特征矩阵，明确模型可见的输入列。
X = np.vstack((X1, X2)).T

# 对二维指标逐列去均值并除以标准差，使 PCA 不受原量纲支配。
X_std = (X - X.mean(axis=0)) / X.std(axis=0)

# 为“数据点的分布揭示了主要变化方向”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(11, 5))
# 将标准化估值指标画成点云，检查两维是否存在共同变化方向。
plt.scatter(X_std[:, 0], X_std[:, 1], c='#0d6efd', s=60, alpha=0.8, edgecolors='w', label='原始数据点 (估值指标)')
# 将横轴标为“标准化市销率 (x1)”，明确横向编码的变量。
plt.xlabel('标准化市销率 (x1)', fontsize=16)
# 将纵轴标为“标准化市盈率 (x2)”，明确纵向编码的变量。
plt.ylabel('标准化市盈率 (x2)', fontsize=16)
# 将图题设为“数据分布揭示了主要变化方向”，直接说明当前图形的比较目的。
plt.title('数据分布揭示了主要变化方向', fontsize=14)
# 在 `0` 处添加水平参考线，标出“两个相关变量的散点图”的基准或阈值。
plt.axhline(0, color='grey', lw=0.5)
# 在 `0` 处添加垂直参考线，标出“两个相关变量的散点图”的基准或阈值。
plt.axvline(0, color='grey', lw=0.5)
# 为“数据点的分布揭示了主要变化方向”添加辅助网格，便于比较位置、斜率或组间差异。
plt.grid(True)
# 取得当前坐标轴对象，以便继续设置图形比例与样式。
plt.gca().set_aspect('equal', adjustable='box')
# 显示“数据点的分布揭示了主要变化方向”图例，使颜色或线型与比较对象一一对应。
plt.legend()
# 显示二维点云及主变化方向，检查样本沿正斜率方向的离散程度明显大于正交方向。
plt.show()

# 为“PCA的目标：找到捕捉最大方差的新坐标轴”，从 `sklearn.decomposition` 导入`PCA` 用于估计主成分方向并投影样本。
from sklearn.decomposition import PCA

# 建立保留两个方向的 PCA，用于从标准化点云估计完整正交基。
pca = PCA(n_components=2)
# 从 `X_std` 的协方差结构估计 `pca` 的主成分方向。
pca.fit(X_std)
# 从 PCA 结果提取 `components`（主成分方向），用于解释投影方向与载荷。
components = pca.components_
# 提取 `explained_variance_` 的特征值尺度方差，用于按投影标准差缩放主成分方向；它不是归一化解释比例。
explained_variance = pca.explained_variance_

# 为“PCA的目标：找到捕捉最大方差的新坐标轴”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(10, 3.5))
# 叠加主成分前先画出原始点云，作为方向比较的参照。
plt.scatter(X_std[:, 0], X_std[:, 1], c='#0d6efd', s=60, alpha=0.8, edgecolors='w', label='原始数据点')
# 将横轴标为“标准化市销率 (x1)”，明确横向编码的变量。
plt.xlabel('标准化市销率 (x1)', fontsize=12)
# 将纵轴标为“标准化市盈率 (x2)”，明确纵向编码的变量。
plt.ylabel('标准化市盈率 (x2)', fontsize=12)
# 将图题设为“PCA找到了捕捉最大方差的新坐标轴”，直接说明当前图形的比较目的。
plt.title('PCA找到了捕捉最大方差的新坐标轴', fontsize=18)
plt.tick_params(labelsize=16)  # 放大刻度以保证投影视距下可读
# 在 `0` 处添加水平参考线，标出“第一主成分捕捉了最大的方差”的基准或阈值。
plt.axhline(0, color='grey', lw=0.5)
# 在 `0` 处添加垂直参考线，标出“第一主成分捕捉了最大的方差”的基准或阈值。
plt.axvline(0, color='grey', lw=0.5)
# 为“PCA的目标：找到捕捉最大方差的新坐标轴”添加辅助网格，便于比较位置、斜率或组间差异。
plt.grid(True)
# 取得当前坐标轴对象，以便继续设置图形比例与样式。
plt.gca().set_aspect('equal', adjustable='box')

# 按第一主成分方差的平方根设定箭头长度，使方向与解释尺度同时可见。
arrow_len_pc1 = np.sqrt(explained_variance[0]) * 3
# 取得第一主成分单位向量，绘制最大方差方向。
pc1_vec = components[0, :]
# 从 PCA 载荷矩阵取出第一主成分单位向量，作为最大样本方差方向的箭头坐标。
plt.arrow(0, 0, pc1_vec[0] * arrow_len_pc1, pc1_vec[1] * arrow_len_pc1,
          # 指定 `head_width` 为箭头头部宽度，细化“第一主成分捕捉了最大的方差”的输出。
          head_width=0.2, head_length=0.3, fc='#dc3545', ec='#dc3545', lw=2,
          # 将 `label` 设为“第一主成分 (PC1)”，明确“第一主成分捕捉了最大的方差”的显示语义。
          label='第一主成分 (PC1)', zorder=3)

# 显示“PCA的目标：找到捕捉最大方差的新坐标轴”图例，使颜色或线型与比较对象一一对应。
plt.legend()
# 显示点云与 PC1 箭头，核对第一主成分沿数据最分散方向而非原坐标轴。
plt.show()

# 为“...然后找到与PC1正交的次要方向”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(10, 6))
# 以同一批标准化样本为底图，核对两个主成分的正交方向。
plt.scatter(X_std[:, 0], X_std[:, 1], c='#0d6efd', s=60, alpha=0.8, edgecolors='w', label='原始数据点')
# 将横轴标为“标准化市销率 (x1)”，明确横向编码的变量。
plt.xlabel('标准化市销率 (x1)', fontsize=12)
# 将纵轴标为“标准化市盈率 (x2)”，明确纵向编码的变量。
plt.ylabel('标准化市盈率 (x2)', fontsize=12)
# 将图题设为“PCA找到了捕捉最大方差的新坐标轴”，直接说明当前图形的比较目的。
plt.title('PCA找到了捕捉最大方差的新坐标轴', fontsize=14)
# 在 `0` 处添加水平参考线，标出“第二主成分与第一主成分正交”的基准或阈值。
plt.axhline(0, color='grey', lw=0.5)
# 在 `0` 处添加垂直参考线，标出“第二主成分与第一主成分正交”的基准或阈值。
plt.axvline(0, color='grey', lw=0.5)
# 为“...然后找到与PC1正交的次要方向”添加辅助网格，便于比较位置、斜率或组间差异。
plt.grid(True)
# 取得当前坐标轴对象，以便继续设置图形比例与样式。
plt.gca().set_aspect('equal', adjustable='box')

# 按 PC1 解释方差设定第一支箭头长度，供与 PC2 的正交方向比较。
arrow_len_pc1 = np.sqrt(explained_variance[0]) * 3
# 取得第一主成分单位向量，绘制最大方差方向。
pc1_vec = components[0, :]
# 取出第一主成分载荷向量，绘制最大方差方向作为正交比较基准。
plt.arrow(0, 0, pc1_vec[0] * arrow_len_pc1, pc1_vec[1] * arrow_len_pc1,
          # 指定 `head_width` 为箭头头部宽度，细化“第二主成分与第一主成分正交”的输出。
          head_width=0.2, head_length=0.3, fc='#dc3545', ec='#dc3545', lw=2,
          # 将 `label` 设为“第一主成分 (PC1)”，明确“第二主成分与第一主成分正交”的显示语义。
          label='第一主成分 (PC1)', zorder=3)

# 用第二主成分标准差的两倍确定次轴箭头长度。
arrow_len_pc2 = np.sqrt(explained_variance[1]) * 3
# 取得与第一主成分正交的第二主成分单位向量。
pc2_vec = components[1, :]
# 取出第二主成分载荷向量，绘制与 PC1 正交的剩余方差方向。
plt.arrow(0, 0, pc2_vec[0] * arrow_len_pc2, pc2_vec[1] * arrow_len_pc2,
          # 指定 `head_width` 为箭头头部宽度，细化“第二主成分与第一主成分正交”的输出。
          head_width=0.2, head_length=0.3, fc='#fd7e14', ec='#fd7e14', lw=2,
          # 将 `label` 设为“第二主成分 (PC2)”，明确“第二主成分与第一主成分正交”的显示语义。
          label='第二主成分 (PC2)', zorder=3)

# 显示“...然后找到与PC1正交的次要方向”图例，使颜色或线型与比较对象一一对应。
plt.legend()
# 显示 PC1 与 PC2 两支箭头，核对方向正交且第一方向解释的离散程度更大。
plt.show()

# 将标准化样本转换为主成分得分，取得每个点在 PC1 轴上的一维坐标。
X_projected = pca.transform(X_std)
# 将第一主成分得分逆变换为 `X_reconstructed_pc1`，取得原坐标系中的投影点。
X_reconstructed_pc1 = pca.inverse_transform(np.c_[X_projected[:,0], np.zeros_like(X_projected[:,0])])

# 绘制标准化点云、PC1 轴及每个样本到轴上的正交投影线。

# 为“降维的实现：将数据投影到PC1上”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(11, 2.8))
# 以 `X_std[:, 0]` 为横轴、`X_std[:, 1]` 为纵轴绘制散点，展示“降维就是将数据点投影到最重要的主成分上”。
plt.scatter(X_std[:, 0], X_std[:, 1], alpha=0.5, label='原始数据点', c='#0d6efd', s=60, edgecolors='w')
# 将横轴标为“标准化市销率 (x1)”，明确横向编码的变量。
plt.xlabel('标准化市销率 (x1)', fontsize=16)
# 将纵轴标为“标准化市盈率 (x2)”，明确纵向编码的变量。
plt.ylabel('标准化市盈率 (x2)', fontsize=16)
# 将图题设为“降维的实现：将数据点投影到第一主成分”，直接说明当前图形的比较目的。
plt.title('降维的实现：将数据点投影到第一主成分', fontsize=18)
plt.tick_params(labelsize=16)  # 统一主成分图的刻度字号
# 在 `0` 处添加水平参考线，标出“降维就是将数据点投影到最重要的主成分上”的基准或阈值。
plt.axhline(0, color='grey', lw=0.5)
# 在 `0` 处添加垂直参考线，标出“降维就是将数据点投影到最重要的主成分上”的基准或阈值。
plt.axvline(0, color='grey', lw=0.5)
# 为“降维的实现：将数据投影到PC1上”添加辅助网格，便于比较位置、斜率或组间差异。
plt.grid(True)
# 取得当前坐标轴对象，以便继续设置图形比例与样式。
plt.gca().set_aspect('equal', adjustable='box')

# 绘制PC1轴
line_range = np.array([-4, 4])
# 生成穿过样本均值的第一主成分直线，承载正交投影示意。
pc1_line = line_range[:, np.newaxis] * components[0, :]
# 绘制第一主成分轴，用作原始点与投影点的共同参照。
plt.plot(pc1_line[:, 0], pc1_line[:, 1], color='#dc3545', lw=2.5, label='第一主成分轴 (PC1)')

# 绘制投影点和投影线
plt.scatter(X_reconstructed_pc1[:, 0], X_reconstructed_pc1[:, 1],
            # 用绿色方形标记第一主成分上的投影点，与原始蓝色散点区分。
            marker='s', c='green', s=50, label='投影后的数据点 (z1)', zorder=3)
# 逐个样本连接原始坐标与 PC1 投影坐标，显式画出几何投影距离。
for i in range(X_std.shape[0]):
    # 以原始点和投影点的第一维为横坐标，绘制两点之间的连接线。
    plt.plot([X_std[i, 0], X_reconstructed_pc1[i, 0]],
             # 以同一对点的第二维为纵坐标，保证线段端点一一对应。
             [X_std[i, 1], X_reconstructed_pc1[i, 1]],
             # 为“降维的实现：将数据投影到PC1上”，把投影连接线设为灰色虚线，并用 0.8 线宽弱化其视觉层级。
             '--', color='gray', lw=0.8)

# 显示“降维的实现：将数据投影到PC1上”图例，使颜色或线型与比较对象一一对应。
plt.legend()
# 显示原点、PC1 轴与每个样本的正交投影，检查二维信息如何压缩为一维坐标。
plt.show()

# 为“步骤 0：机制模拟（不是公开经验验证）”导入 `pandas` 并绑定 `pd`，用于整理当前任务的表格、字段与时间索引。
import pandas as pd
# 为“步骤 0：机制模拟（不是公开经验验证）”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“步骤 0：机制模拟（不是公开经验验证）”，从 `sklearn.model_selection` 导入`train_test_split` 用于把样本随机拆成两个互斥子集。
from sklearn.model_selection import train_test_split
# 为“步骤 0：机制模拟（不是公开经验验证）”，导入 `StandardScaler`，供 Pipeline 在每个训练子集内按特征列估计均值与尺度。
from sklearn.preprocessing import StandardScaler
# 为“步骤 0：机制模拟（不是公开经验验证）”，从 `sklearn.linear_model` 导入`LinearRegression` 用于拟合普通最小二乘线性回归。
from sklearn.linear_model import LinearRegression
# 为“步骤 0：机制模拟（不是公开经验验证）”，从 `sklearn.metrics` 导入`mean_squared_error` 用于计算真实值与预测值的均方误差。
from sklearn.metrics import mean_squared_error
# 为“步骤 0：机制模拟（不是公开经验验证）”，从 `sklearn.decomposition` 导入`PCA` 用于估计主成分方向并投影样本。
from sklearn.decomposition import PCA
# 为“步骤 0：机制模拟（不是公开经验验证）”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt
# 为“步骤 0：机制模拟（不是公开经验验证）”导入 `seaborn` 并绑定 `sns`，用于编码当前任务的统计分布或分组关系。
import seaborn as sns

# 固定随机状态并设定机制样本规模，构造 PCA 不保证预测改进的比较案例。
# 确保结果可复现
np.random.seed(42)
# 固定样本数与特征数，界定“步骤 0：机制模拟（不是公开经验验证）”的机制演示规模。
n_samples = 500

# 生成共同潜在因子，构造三个相关财务特征的共享变化来源。
latent_factor = np.random.randn(n_samples) * 2

# 由共同潜在因子生成 ROA、B/M 与价格代理三项相关特征，并分别叠加噪声。
# roa: Return on Assets (与因子正相关)
roa = 0.6 * latent_factor + np.random.randn(n_samples) * 0.5
# bm: Book-to-Market (价值因子，也与因子正相关)
bm = 0.5 * latent_factor + np.random.randn(n_samples) * 0.6
# pps: Price per Share（仅表示当前模拟机制中的每股价格水平，与潜在因子弱相关）
pps = 0.2 * latent_factor + np.random.randn(n_samples) * 0.8

# 由潜在因子、ROA 和噪声生成类 EPS 响应，保留已知结构供 PCA 机制比较。
simulated_eps_like = 1.2 * latent_factor + 0.8 * roa + np.random.randn(n_samples) * 0.4

# 将三项机制特征与类 EPS 响应按行组成表格，固定后续切分字段顺序。
data = pd.DataFrame({
    # 保存机制生成的 ROA 特征。
    'roa': roa,
    # 保存机制生成的账面市值比特征。
    'bm': bm,
    # 保存机制生成的每股股价特征。
    'pps': pps,
    # 保存由潜在因子、ROA 与噪声生成的模拟类 EPS 响应。
    'simulated_eps_like': simulated_eps_like
# 完成三项模拟特征与类 EPS 响应表，仅用于演示 PCA 的相关结构机制。
})

# 用固定随机状态随机划分机制训练集和测试集，保证两种回归规格使用同一留出样本。
training_data, testing_data = train_test_split(data, test_size=0.3, random_state=42)

# 报告机制数据已按固定随机状态划分为训练集与测试集。
print("模拟数据已生成并分割。")

# 显示训练数据的基本统计信息
training_data.describe().round(2)

# 为“EDA: 检查特征间的相关性”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(7, 4.25))
# 把 `training_data.corr()` 编码为热图色块，比较“特征相关性热力图”中的成对关系。
sns.heatmap(training_data.corr(), annot=True, cmap='vlag', fmt='.2f', linewidths=.5)
# 将图题设为“特征与目标变量的相关性矩阵”，直接说明当前图形的比较目的。
plt.title('特征与目标变量的相关性矩阵')
# 显示特征与目标的相关矩阵，检查三项特征间冗余及其与响应变量的线性关联。
plt.show()

# 分离训练样本中的三项模拟特征，避免类 EPS 响应进入输入矩阵。
X_train = training_data.drop('simulated_eps_like', axis=1)
# 提取训练样本的类 EPS 响应，供原始特征与 PCA 回归共同拟合。
y_train = training_data['simulated_eps_like']
# 提取 `X_test` 作为“步骤 2: 数据准备与标准化”的特征矩阵，明确模型可见的输入列。
X_test = testing_data.drop('simulated_eps_like', axis=1)
# 提取随机留出样本的类 EPS 响应，按同一观测顺序评价两个模型。
y_test = testing_data['simulated_eps_like']

# 创建只由训练样本逐特征拟合的缩放器。
scaler_X = StandardScaler()
# 为“步骤 2: 数据准备与标准化”，逐特征估计训练样本均值与标准差，并转换训练特征。
X_train_scaled = scaler_X.fit_transform(X_train)
# 用确定的训练样本逐特征统计量转换测试特征。
X_test_scaled = scaler_X.transform(X_test)

# 报告“数据标准化完成。训练集shape:”中的 `X_train_scaled`，核对“步骤 2: 数据准备与标准化”的样本形状。
print("数据标准化完成。训练集shape:", X_train_scaled.shape)

# 创建并拟合PCA对象
pca_all = PCA(n_components=3)
# 从 `X_train_scaled` 的协方差结构估计 `pca_all` 的主成分方向。
pca_all.fit(X_train_scaled)

# 获取每个主成分解释的方差比率
variance_explained = pca_all.explained_variance_ratio_
# 报告每个主成分的解释方差比例，比较单轴信息贡献。
print(f'各主成分解释的方差比率: {np.round(variance_explained, 3)}')
# 报告累计解释方差，判断保留多少主成分可覆盖主要信息。
print(f'累计解释方差: {np.round(np.cumsum(variance_explained), 3)}')

# 为“可视化解释方差 (碎石图)”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(8, 5))
# 绘制三个主成分各自的解释方差比柱形，柱高表示单个方向保留的信息比例。
plt.bar(range(1, 4), variance_explained, color='skyblue', alpha=0.9, label='单个PC方差')
# 叠加累计解释方差折线，展示保留多个主成分后的总信息比例。
plt.plot(range(1, 4), np.cumsum(variance_explained), 'ro-', label='累计方差')
# 将横轴标为“主成分序号”，明确横向编码的变量。
plt.xlabel('主成分序号')
# 将纵轴标为“解释的方差比率”，明确纵向编码的变量。
plt.ylabel('解释的方差比率')
# 将图题设为“各主成分解释的方差比率（碎石图）”，直接说明当前图形的比较目的。
plt.title('各主成分解释的方差比率（碎石图）')
# 查询当前自动生成的横轴刻度位置与标签；本句不修改刻度设置。
plt.xticks()
# 为“可视化解释方差 (碎石图)”，限定纵轴显示范围，避免比例差异误导比较。
plt.ylim(0, 1.1)
# 显示“可视化解释方差 (碎石图)”图例，使颜色或线型与比较对象一一对应。
plt.legend()
# 为“可视化解释方差 (碎石图)”添加辅助网格，便于比较位置、斜率或组间差异。
plt.grid(axis='x', linestyle='')
# 显示单项与累计解释方差，检查 PC1 是否已保留主要变异以及追加成分的边际信息。
plt.show()

# 创建一个新的PCA对象，目标维度为1
pca_k1 = PCA(n_components=1)

# 在训练数据上拟合并转换
X_train_pca = pca_k1.fit_transform(X_train_scaled)

# 在测试数据上进行转换
X_test_pca = pca_k1.transform(X_test_scaled)

# 报告“降维前shape:”中的 `X_train_scaled`，核对“步骤 4: 执行降维”的样本形状。
print("降维前shape:", X_train_scaled.shape)
# 报告“降维后shape:”中的 `X_train_pca`，核对“步骤 4: 执行降维”的样本形状。
print("降维后shape:", X_train_pca.shape)

# 模型1: 基于原始3个特征
ols_model = LinearRegression()
# 用训练期三个标准化特征估计 OLS 截距与系数，作为未降维预测基准。
ols_model.fit(X_train_scaled, y_train)

# 模型2: 基于PCA降维后的1个特征
ols_pca_model = LinearRegression()
# 用训练期 PC1 得分估计一元线性回归的截距与斜率；主成分方向已由前一步 PCA 确定。
ols_pca_model.fit(X_train_pca, y_train)

# 报告原始特征回归与 PCA 回归均已在同一训练期拟合。
print('两个模型训练完成。')

# 用三特征 OLS 对标准化测试特征生成预测，作为未降维基准。
test_predictions = ols_model.predict(X_test_scaled)
# 计算 `test_mse`（均方误差），量化“步骤 6: 模型评估”的模型表现。
test_mse = mean_squared_error(y_test, test_predictions)

# 用已拟合的 PC1 回归模型生成测试集预测数组，供计算降维模型的 MSE。
test_pca_predictions = ols_pca_model.predict(X_test_pca)
# 计算 `test_pca_mse`（均方误差），量化“步骤 6: 模型评估”的模型表现。
test_pca_mse = mean_squared_error(y_test, test_pca_predictions)

# 打印三特征 OLS 的测试 MSE，作为同一切分下的未降维基准。
print(f'原始模型 (3个特征) - 测试集MSE: {test_mse:.4f}')
# 打印仅用 PC1 的回归测试 MSE，与三特征基准按相同样本比较。
print(f'PCA模型 (1个特征) -  测试集MSE: {test_pca_mse:.4f}')

# 汇总当前固定种子下的两种测试 MSE 与 PC1 解释方差比，保证展示值和模型输出同源。
pca_comparison_table = pd.DataFrame({'评估口径': ['固定种子、随机留出'], '三特征 OLS 测试 MSE': [test_mse], '仅 PC1 OLS 测试 MSE': [test_pca_mse], 'PC1 解释方差比': [pca_k1.explained_variance_ratio_[0]]})
# 将同源计算结果统一显示到六位小数，便于核对模型比较与解释方差。
pca_comparison_table.round(6)

# 创建 `fig, axes` 画布，承载“两个模型在测试集上的预测表现”的并排视觉比较。
fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharey=True, sharex=True)
# 汇总真实值与两组预测的共同边界 `max_val`，让两幅散点图使用同一尺度。
max_val = max(y_test.max(), test_predictions.max(), test_pca_predictions.max())
# 汇总真实值与两组预测的共同边界 `min_val`，让两幅散点图使用同一尺度。
min_val = min(y_test.min(), test_predictions.min(), test_pca_predictions.min())

# 原始模型
axes[0].scatter(y_test, test_predictions, alpha=0.6, edgecolors='w')
# 绘制曲线元素，编码 `axes` 与 `min_val` 与 `max_val` 的关系，服务于“两个模型在测试集上的预测表现”的视觉比较。
axes[0].plot([min_val, max_val], [min_val, max_val], 'r--', lw=2)
# 用图题概括“可视化模型表现：预测值 vs. 真实值”的比较对象与当前计算结果。
axes[0].set_title(f'原始模型 (3特征)\nTest MSE: {test_mse:.3f}')
# 将横轴标为“真实值 (Actual)”，明确横向编码的变量。
axes[0].set_xlabel('真实值 (Actual)')
# 将纵轴标为“预测值 (Predicted)”，明确纵向编码的变量。
axes[0].set_ylabel('预测值 (Predicted)')
# 为“可视化模型表现：预测值 vs. 真实值”添加辅助网格，便于比较位置、斜率或组间差异。
axes[0].grid(True)

# PCA模型
axes[1].scatter(y_test, test_pca_predictions, alpha=0.6, edgecolors='w', c='green')
# 绘制曲线元素，编码 `axes` 与 `min_val` 与 `max_val` 的关系，服务于“两个模型在测试集上的预测表现”的视觉比较。
axes[1].plot([min_val, max_val], [min_val, max_val], 'r--', lw=2)
# 用图题概括“可视化模型表现：预测值 vs. 真实值”的比较对象与当前计算结果。
axes[1].set_title(f'PCA模型 (1特征)\nTest MSE: {test_pca_mse:.3f}')
# 将横轴标为“真实值 (Actual)”，明确横向编码的变量。
axes[1].set_xlabel('真实值 (Actual)')
# 为“可视化模型表现：预测值 vs. 真实值”添加辅助网格，便于比较位置、斜率或组间差异。
axes[1].grid(True)

# 将图题设为“预测值 vs. 真实值 (测试集)”，直接说明当前图形的比较目的。
plt.suptitle('预测值 vs. 真实值 (测试集)', fontsize=16)
# 为“可视化模型表现：预测值 vs. 真实值”，压缩子图留白，避免标题和坐标标签相互遮挡。
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
# 显示原特征模型与 PC1 模型的测试预测对真实值散点，比较两者偏离 45 度线的程度。
plt.show()

from pathlib import Path
from urllib.request import urlretrieve  # 复用本章隐藏设置单元安装的浏览器标识下载器
from sklearn.impute import SimpleImputer  # 在每个训练窗估计缺失中位数
from sklearn.linear_model import Ridge  # 建立有 PCA 与无 PCA 的同级模型
from sklearn.pipeline import Pipeline  # 把插补、缩放、PCA 和回归锁进折内流程
# 按 Linux 共享数据、Windows 共享数据、项目缓存的顺序选择财务报表文件。
financial_path = next((candidate_path for candidate_path in [Path('/home/ubuntu/r2_data_mount/data/stock/financial_statement.h5'), Path('C:/qiufei/data/stock/financial_statement.h5'), Path('data/course/financial_statement.h5')] if candidate_path.exists()), Path('data/course/financial_statement.h5'))
if not financial_path.exists():
    financial_path.parent.mkdir(parents=True, exist_ok=True)
    urlretrieve('https://assets.qiufei.site/data/stock/financial_statement.h5', financial_path)
pca_companies = ['600276.XSHG', '600009.XSHG', '600585.XSHG', '603259.XSHG', '002415.XSHE', '600104.XSHG', '600019.XSHG', '600406.XSHG']  # 确定八家公司
pca_quarters = [f'{year}q4' for year in range(2005, 2025)]  # 覆盖延迟披露标签所需的历史年度 Q4 版本
pca_columns = ['order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities', 'operating_revenue', 'net_profit_parent_company', 'basic_earnings_per_share']  # 确定输入字段
pca_company_frames = [pd.read_hdf(financial_path, key='financial_data', where=f'order_book_id == "{company}"', columns=pca_columns) for company in pca_companies]  # 逐公司选择性读取必要列
pca_reports = pd.concat(pca_company_frames, ignore_index=True).query('quarter in @pca_quarters').copy()  # 只保留确定范围内的年度 Q4 报告版本
pca_reports['info_date'] = pd.to_datetime(pca_reports['info_date'])  # 统一披露日类型
pca_reports['fiscal_year'] = pca_reports['quarter'].str[:4].astype(int)  # 提取财年用于时间对齐
pca_target_versions = pca_reports.rename(columns={'fiscal_year': 'target_year', 'info_date': 'target_info_date', 'basic_earnings_per_share': 'target_eps'}).copy()  # 保留公司财年的全部标签披露版本
pca_targets = pca_target_versions.sort_values(['order_book_id', 'target_year', 'target_info_date']).drop_duplicates(['order_book_id', 'target_year'], keep='last')  # 最终版本只作为确定后的事后评分目标

pca_grid = pca_targets[['order_book_id', 'target_year']].query('2010 <= target_year <= 2024').copy()  # 为延迟披露的历史验证折保留足够训练年份
pca_grid['decision_date'] = pd.to_datetime(pca_grid['target_year'].astype(str) + '-06-30')  # 确定每年六月末预测时点
pca_candidates = pca_grid.merge(pca_reports, on='order_book_id', how='left').query('fiscal_year < target_year and info_date <= decision_date').copy()  # 只列举当时已披露的历史年报
pca_features_at_decision = pca_candidates.sort_values(['order_book_id', 'target_year', 'info_date', 'fiscal_year']).drop_duplicates(['order_book_id', 'target_year'], keep='last')  # 选择预测时点最新可得版本
pca_target_values = pca_targets[['order_book_id', 'target_year', 'target_info_date', 'target_eps']].copy()  # 建立仅供确定后评分的最终目标表
pca_panel = pca_features_at_decision.merge(pca_target_values, on=['order_book_id', 'target_year'], how='inner')  # 合并点时特征与最终评分目标
pca_panel['log_assets'] = np.log(pca_panel['total_assets'])  # 构造对数规模
pca_panel['debt_ratio'] = pca_panel['total_liabilities'] / pca_panel['total_assets']  # 构造负债率
pca_panel['log_revenue'] = np.log(pca_panel['operating_revenue'].clip(lower=1))  # 构造对数收入
pca_panel['roa'] = pca_panel['net_profit_parent_company'] / pca_panel['total_assets']  # 构造资产收益率
pca_panel['log_profit_scale'] = np.sign(pca_panel['net_profit_parent_company']) * np.log1p(np.abs(pca_panel['net_profit_parent_company']))  # 构造带符号盈利规模
pca_feature_names = ['log_assets', 'debt_ratio', 'log_revenue', 'roa', 'log_profit_scale']  # 确定五项 PCA 输入
pca_panel = pca_panel.replace([np.inf, -np.inf], np.nan).dropna(subset=['target_eps']).sort_values(['target_year', 'order_book_id'])  # 清除无效目标并保持时间顺序
assert pca_panel['info_date'].le(pca_panel['decision_date']).all(), 'PCA 特征晚于预测时点'  # 阻止点时特征泄漏
print({'文件': financial_path.name, '样本期': (pca_panel.target_year.min(), pca_panel.target_year.max()), '公司年': len(pca_panel), '特征': pca_feature_names})  # 输出输入说明

def build_pca_label_snapshot(feature_panel, cutoff_date):  # 为指定历史截止日重建当时可得的标签版本
    cutoff = pd.Timestamp(cutoff_date)  # 统一拟合或选择截止日类型
    label_candidates = pca_target_versions.query('target_info_date <= @cutoff').copy()  # 排除截止日后披露或修订的标签
    available_labels = label_candidates.sort_values(['order_book_id', 'target_year', 'target_info_date']).drop_duplicates(['order_book_id', 'target_year'], keep='last')  # 选择截止日前最后可得版本
    feature_columns = [column for column in feature_panel.columns if column not in ['target_info_date', 'target_eps']]  # 防止最终评分标签进入历史数据版本
    snapshot = feature_panel[feature_columns].merge(available_labels[['order_book_id', 'target_year', 'target_info_date', 'target_eps']], on=['order_book_id', 'target_year'], how='inner')  # 合并时点一致的特征与标签
    snapshot = snapshot.dropna(subset=['target_eps']).sort_values(['target_year', 'order_book_id'])  # 只保留当时已披露的完整标签
    assert snapshot['target_info_date'].le(cutoff).all(), 'PCA 标签晚于指定截止日'  # 阻止未来标签或修订进入历史评估
    assert not snapshot.duplicated(['order_book_id', 'target_year']).any(), 'PCA 标签数据版本公司年不唯一'  # 强制每个公司年仅一个当时版本
    snapshot.attrs['label_flow'] = {'截止日': cutoff.date(), '候选版本': len(label_candidates), '截止日前唯一公司年': len(available_labels), '完整行': len(snapshot)}  # 保存可复查的样本数量变化
    return snapshot  # 返回指定截止日可实现的标签面板

pca_validation_years = [2018, 2019, 2020]  # 只使用能在 2023 最终确定前取得标签的验证年
ridge_alphas = [0.1, 1.0, 10.0]  # 事先确定 Ridge 惩罚候选
pca_component_counts = [2, 3, 4]  # 事先确定 PCA 维数候选
pca_candidate_configs = [{'family': 'Ridge', 'name': f'Ridge(a={alpha})', 'alpha': alpha, 'components': None} for alpha in ridge_alphas]  # 建立无 PCA 候选
pca_candidate_configs += [{'family': 'PCA-Ridge', 'name': f'PCA(k={components},a={alpha})', 'alpha': alpha, 'components': components} for components in pca_component_counts for alpha in ridge_alphas]  # 建立完整 PCA 候选
pca_validation_records = []  # 收集逐窗候选和训练均值基线
pca_validation_flows = []  # 收集每折拟合与选择时点的标签样本数量变化
pca_validation_inputs = []  # 缓存各折已按时点确定的训练与验证行
for validation_year in pca_validation_years:  # 逐年执行扩展窗口验证
    pca_fit_cutoff = pd.Timestamp(f'{validation_year}-06-30')  # 在当前目标年预测日确定候选拟合信息
    pca_selection_cutoff = pd.Timestamp(f'{validation_year + 3}-06-30')  # 按公开数据版本约三年披露滞后定义评分截止日
    pca_training_snapshot = build_pca_label_snapshot(pca_panel, pca_fit_cutoff)  # 重建拟合截止日前标签版本
    pca_selection_snapshot = build_pca_label_snapshot(pca_panel, pca_selection_cutoff)  # 重建选择截止日前标签版本
    pca_training_rows = pca_training_snapshot.query('target_year < @validation_year').copy()  # 只使用更早且拟合时已知的公司年
    pca_validation_rows = pca_selection_snapshot.query('target_year == @validation_year').copy()  # 只用选择截止日前已披露的验证标签
    assert len(pca_training_rows) > 0 and len(pca_validation_rows) > 0, 'PCA 训练折或验证折为空'  # 先阻止空表让后续时点断言真空通过
    assert pca_training_rows['target_info_date'].le(pca_fit_cutoff).all(), 'PCA 训练折含晚披露标签'  # 执行逐折拟合时点断言
    assert pca_validation_rows['target_info_date'].le(pca_selection_cutoff).all(), 'PCA 验证折含晚披露标签'  # 执行逐折选择时点断言
    pca_training_flow = pca_training_snapshot.attrs['label_flow']  # 读取拟合截止日的标签样本数量变化
    pca_selection_flow = pca_selection_snapshot.attrs['label_flow']  # 读取选择截止日的标签样本数量变化
    pca_validation_total = pca_panel['target_year'].eq(validation_year).sum()  # 统计该验证年最终评分池规模
    pca_validation_flows.append({'验证年': validation_year, '拟合截止日': pca_fit_cutoff.date(), '拟合候选版本': pca_training_flow['候选版本'], '拟合唯一公司年': pca_training_flow['截止日前唯一公司年'], '完整训练行': len(pca_training_rows), '选择截止日': pca_selection_cutoff.date(), '选择候选版本': pca_selection_flow['候选版本'], '验证可用行': len(pca_validation_rows), '验证排除行': pca_validation_total - len(pca_validation_rows), '验证标签最晚披露': pca_validation_rows['target_info_date'].max().date()})  # 记录逐折标签版本、排除与可得日
    mean_predictions = np.repeat(pca_training_rows['target_eps'].mean(), len(pca_validation_rows))  # 构造训练均值预测
    pca_validation_records.append({'model': '训练均值', 'family': 'Baseline', 'fold': validation_year, 'mse': mean_squared_error(pca_validation_rows['target_eps'], mean_predictions)})  # 保存基线比较结果
    pca_validation_inputs.append({'fold': validation_year, 'training_rows': pca_training_rows, 'validation_rows': pca_validation_rows})  # 交给下一代码块复用同一确定样本

for fold_input in pca_validation_inputs:  # 在已确定的各折样本上拟合完整候选集
    validation_year = fold_input['fold']  # 读取当前验证年
    pca_training_rows = fold_input['training_rows']  # 读取拟合截止日确定的训练行
    pca_validation_rows = fold_input['validation_rows']  # 读取选择截止日确定的验证行
    for config in pca_candidate_configs:  # 对完整事先确定集合使用同一时间窗
        steps = [('imputer', SimpleImputer(strategy='median')), ('scaler', StandardScaler())]  # 只在当前训练折逐特征估计插补与缩放统计量
        steps += [('pca', PCA(n_components=config['components'], random_state=42))] if config['components'] else []  # 仅 PCA 候选增加折内降维
        candidate_pipeline = Pipeline(steps + [('ridge', Ridge(alpha=config['alpha']))]).fit(pca_training_rows[pca_feature_names], pca_training_rows['target_eps'])  # 拟合当前候选
        fold_predictions = candidate_pipeline.predict(pca_validation_rows[pca_feature_names])  # 用候选流程对当前验证窗的具名财务特征生成预测
        pca_validation_records.append({'model': config['name'], 'family': config['family'], 'fold': validation_year, 'mse': mean_squared_error(pca_validation_rows['target_eps'], fold_predictions)})  # 保存逐窗 MSE

pca_validation_table = pd.DataFrame(pca_validation_records)  # 形成完整候选比较结果
display(pca_validation_table.query("family == 'PCA-Ridge' and model.str.startswith('PCA(k=2')", engine='python').pivot(index='model', columns='fold', values='mse').round(4))  # 输出 k=2 的三组惩罚候选

display(pca_validation_table.query("family == 'PCA-Ridge' and model.str.startswith('PCA(k=3')", engine='python').pivot(index='model', columns='fold', values='mse').round(4))  # 展示三主成分候选的逐折误差

display(pca_validation_table.query("family == 'PCA-Ridge' and model.str.startswith('PCA(k=4')", engine='python').pivot(index='model', columns='fold', values='mse').round(4))  # 展示四主成分候选的逐折误差

display(pca_validation_table.query("family in ['Ridge', 'Baseline']").pivot(index='model', columns='fold', values='mse').round(4))  # 输出全部无 PCA 候选与训练均值基线

display(pd.DataFrame(pca_validation_flows))

pca_validation_means = pca_validation_table.query("family in ['Ridge', 'PCA-Ridge']").groupby(['family', 'model'])['mse'].mean()  # 只用验证窗汇总两类候选
best_ridge_name = pca_validation_means.loc['Ridge'].idxmin()  # 确定无 PCA 规格
best_pca_name = pca_validation_means.loc['PCA-Ridge'].idxmin()  # 确定 PCA-Ridge 规格
best_ridge_config = next(config for config in pca_candidate_configs if config['name'] == best_ridge_name)  # 读取确定 Ridge 参数
best_pca_config = next(config for config in pca_candidate_configs if config['name'] == best_pca_name)  # 读取确定 PCA 参数
pca_loadings_by_fold = []  # 收集确定 PCA 规格的逐折载荷
for validation_year in pca_validation_years:  # 复算每个训练窗的 PCA 载荷
    pca_fit_cutoff = pd.Timestamp(f'{validation_year}-06-30')  # 恢复该折拟合时点
    pca_training_rows = build_pca_label_snapshot(pca_panel, pca_fit_cutoff).query('target_year < @validation_year').copy()  # 使用拟合日当时可得的训练标签
    loading_pipeline = Pipeline([('imputer', SimpleImputer(strategy='median')), ('scaler', StandardScaler()), ('pca', PCA(n_components=best_pca_config['components'], random_state=42)), ('ridge', Ridge(alpha=best_pca_config['alpha']))]).fit(pca_training_rows[pca_feature_names], pca_training_rows['target_eps'])  # 只由当前训练折逐特征拟合缩放统计量与确定规格
    pca_loadings_by_fold.append(loading_pipeline.named_steps['pca'].components_)  # 保存主成分载荷
loading_stability = [float(np.abs(np.sum(pca_loadings_by_fold[index] * pca_loadings_by_fold[index - 1], axis=1)).mean()) for index in range(1, len(pca_loadings_by_fold))]  # 计算相邻折同序主成分绝对余弦
print({'确定Ridge': best_ridge_name, '确定PCA-Ridge': best_pca_name, '相邻折平均载荷稳定性': np.round(loading_stability, 4).tolist()})  # 输出确定与稳定性证据

pca_final_fit_cutoff = pd.Timestamp('2023-06-30')  # 在 2024 最终测试预测日前一年确定候选与训练信息
assert max(pd.Timestamp(f'{year + 3}-06-30') for year in pca_validation_years) <= pca_final_fit_cutoff, 'PCA 验证标签尚未到选择截止日'  # 保证延迟披露的选模证据在模型选择日前可得
pca_development_snapshot = build_pca_label_snapshot(pca_panel, pca_final_fit_cutoff)  # 重建最终拟合截止日可得的标签版本
pca_development_rows = pca_development_snapshot.query('target_year < 2024').copy()  # 只使用模型选择日已披露且早于最终测试年的开发标签
pca_testing_rows = pca_panel.query('target_year == 2024').copy()  # 仅在确定后的下一年打开最终评分标签
assert len(pca_development_rows) > 0 and len(pca_testing_rows) > 0, 'PCA 最终开发期或最终测试期为空'  # 在拟合与评分前明确验证两窗有样本
ridge_pipeline = Pipeline([('imputer', SimpleImputer(strategy='median')), ('scaler', StandardScaler()), ('ridge', Ridge(alpha=best_ridge_config['alpha']))]).fit(pca_development_rows[pca_feature_names], pca_development_rows['target_eps'])  # 只由最终开发样本逐特征拟合插补、缩放与确定 Ridge
final_pca_pipeline = Pipeline([('imputer', SimpleImputer(strategy='median')), ('scaler', StandardScaler()), ('pca', PCA(n_components=best_pca_config['components'], random_state=42)), ('ridge', Ridge(alpha=best_pca_config['alpha']))]).fit(pca_development_rows[pca_feature_names], pca_development_rows['target_eps'])  # 只由最终开发样本逐特征拟合缩放、PCA 与确定 Ridge
sealed_baseline_predictions = np.repeat(pca_development_rows['target_eps'].mean(), len(pca_testing_rows))  # 构造同窗训练均值基线
sealed_ridge_predictions = ridge_pipeline.predict(pca_testing_rows[pca_feature_names])  # 生成无 PCA 最终测试预测
sealed_pca_predictions = final_pca_pipeline.predict(pca_testing_rows[pca_feature_names])  # 生成 PCA-Ridge 最终测试预测
sealed_pca_mse = mean_squared_error(pca_testing_rows['target_eps'], sealed_pca_predictions)  # 计算 PCA-Ridge 最终测试 MSE
sealed_results = pd.DataFrame({'模型': ['训练均值', best_ridge_name, best_pca_name], '最终测试MSE': [mean_squared_error(pca_testing_rows['target_eps'], sealed_baseline_predictions), mean_squared_error(pca_testing_rows['target_eps'], sealed_ridge_predictions), sealed_pca_mse]})  # 汇总同窗结果
sealed_results['相对训练均值改善'] = 1 - sealed_results['最终测试MSE'] / sealed_results.loc[sealed_results['模型'].eq('训练均值'), '最终测试MSE'].iloc[0]  # 计算相对基线改善
fitted_pca = final_pca_pipeline.named_steps['pca']  # 读取开发期 PCA 而不重拟合测试期
print({'训练期': (int(pca_development_rows.target_year.min()), int(pca_development_rows.target_year.max())), '模型选择日': pca_final_fit_cutoff.date(), **pca_development_snapshot.attrs['label_flow'], '最终开发行': len(pca_development_rows), '最终测试期': 2024, '确定后最终标签用途': '仅一次性评分', '累计解释方差': round(float(fitted_pca.explained_variance_ratio_.sum()), 4)})  # 报告真实样本数量变化、窗口与解释方差
display(sealed_results.round(4))  # 输出三类同窗最终测试证据