11 超越“黑箱”——监督学习总结与模型可解释性

# 导入 Matplotlib，绘制模拟观测、已知生成函数和三种候选预测曲线。
import matplotlib.pyplot as plt
# 创建 `fig, axes` 画布，承载“模型复杂度与拟合效果的关系”的并排视觉比较。
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# 遍历三种候选阶数，把同源计算指标直接写入对应面板标题。
for plot_index, polynomial_degree in enumerate(candidate_degrees):
    # 选取当前子图轴，分别承载欠拟合、适度拟合与过拟合曲线。
    axis = axes[plot_index]
    # 读取当前候选的共享执行结果，避免图题与曲线分别维护数值。
    candidate_evidence = model_evidence_by_degree[polynomial_degree]
    # 绘制已知无噪声生成函数，为独立网格误差提供可见参照。
    axis.plot(independent_feature_grid, independent_true_responses, color='#005f73', linestyle='--', linewidth=1.8, label='已知生成函数')
    # 绘制当前阶数在独立网格上的预测曲线，展示训练点之间的形态。
    axis.plot(independent_feature_grid, candidate_evidence['grid_predictions'], label='模型拟合', color='crimson', linewidth=2)
    # 绘制固定种子生成的含噪观测，明确散点并非真实金融数据。
    axis.scatter(simulated_feature_observations, simulated_noisy_responses, edgecolor='#174F87', s=30, label='模拟含噪观测', facecolors='none')
    # 标注模拟特征横轴，使生成机制的输入范围清晰可见。
    axis.set_xlabel('模拟特征')
    # 标注模拟响应纵轴，使其不被误读为实际金融指标。
    axis.set_ylabel('模拟响应')
    # 为“可视化理解：过拟合的陷阱”，限定横轴范围，使比较对象使用一致尺度。
    axis.set_xlim((0, 1))
    # 为“可视化理解：过拟合的陷阱”，限定纵轴范围，使比较对象使用一致尺度。
    axis.set_ylim((-2, 2))
    # 用同一执行结果动态报告训练与网格 MSE，避免静态结论和代码输出漂移。
    axis.set_title(f"{candidate_labels[polynomial_degree]}（阶数={polynomial_degree}）\n训练 MSE={candidate_evidence['train_mse']:.4f}｜网格 MSE={candidate_evidence['grid_mse']:.4f}", fontsize=12)
    # 显示“可视化理解：过拟合的陷阱”图例，使颜色或线型与比较对象一一对应。
    axis.legend(loc='upper right', fontsize=9)

# 用总标题再次标明模拟属性与证据口径，防止教学机制被误读为真实经验结果。
fig.suptitle('模拟余弦机制：训练拟合不等于独立网格表现', fontsize=18)
# 为“可视化理解：过拟合的陷阱”，压缩子图留白，避免标题和坐标标签相互遮挡。
plt.tight_layout(rect=[0, 0, 1, 0.96])
# 显示一阶、四阶与十五阶的模拟机制证据，比较训练误差与独立网格误差。
plt.show()

# 为“步骤 1：导入必要的工具库”导入 `pandas` 并绑定 `pd`，用于整理当前任务的表格、字段与时间索引。
import pandas as pd
# 为“步骤 1：导入必要的工具库”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 导入 SHAP 并绑定为 `shap`，用于构造树模型解释器及全局、局部贡献图。
import shap
# 为“步骤 1：导入必要的工具库”，从 `sklearn.ensemble` 导入`RandomForestRegressor` 用于拟合随机森林回归器。
from sklearn.ensemble import RandomForestRegressor
# 为“步骤 1：导入必要的工具库”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt
# 为“步骤 1：导入必要的工具库”，从 `IPython.display` 导入 `display`，在幻灯片中渲染表格结果。
from IPython.display import display

# 为了让matplotlib的图表更好看

from pathlib import Path
from urllib.request import urlretrieve  # 复用本章隐藏设置单元安装的浏览器标识下载器
# 按 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)
# 确定年度 Q4 行条件，使大型 HDF5 在读取前先按样本期过滤。
financial_quarters = [f'{year}q4' for year in range(2015, 2025)]
# 只读取公司、季度、披露日、五项解释字段与 EPS，控制大表内存占用。
financial_columns = ['order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities', 'operating_revenue', 'net_profit_parent_company', 'cash_from_operating_activities', 'basic_earnings_per_share']
# 从表式 HDF5 键同时选择年度行与必要列，数据获取与后续解释逻辑保持解耦。
financial_raw = pd.read_hdf(financial_path, key='financial_data', where='quarter in financial_quarters', columns=financial_columns)
# 统一披露日期类型，为版本选择与时间切分建立可复算时点。
financial_raw['info_date'] = pd.to_datetime(financial_raw['info_date'])
# 检查 HDF 行过滤只返回事先确定年度，防止季度口径静默漂移。
assert financial_raw['quarter'].isin(financial_quarters).all(), 'HDF 行筛选返回了事先确定范围外季度'
# 从季度字段提取财年，作为严格时间留出的分组变量。
financial_raw['fiscal_year'] = financial_raw['quarter'].str[:4].astype(int)
# 对公司—财年保留最后披露版本，避免同一报告修订重复进入样本。
financial_panel = financial_raw.sort_values(['order_book_id', 'fiscal_year', 'info_date']).drop_duplicates(['order_book_id', 'fiscal_year'], keep='last')
# 固定五项财务解释字段，避免把目标或披露后信息放进模型输入。
financial_features = ['total_assets', 'total_liabilities', 'operating_revenue', 'net_profit_parent_company', 'cash_from_operating_activities']
# 删除目标或输入缺失的记录，保持模型和 SHAP 使用完全相同的行。
financial_panel = financial_panel.dropna(subset=financial_features + ['basic_earnings_per_share']).copy()
# 提取公开财务特征矩阵，按亿元缩放只改善数值可读性而不改变树的排序。
X = financial_panel[financial_features].div(1e8)
# 提取同一已披露年度报告的 EPS 响应，明确本节是同期模型解释。
y = financial_panel['basic_earnings_per_share'].astype(float)
# 用 2015—2021 财年拟合候选模型，留出后续年份作解释检查。
X_train, y_train = X.loc[financial_panel['fiscal_year'].le(2021)], y.loc[financial_panel['fiscal_year'].le(2021)]
# 用 2022—2024 财年形成留出样本，不随机打乱公司—年份记录。
X_test, y_test = X.loc[financial_panel['fiscal_year'].ge(2022)], y.loc[financial_panel['fiscal_year'].ge(2022)]

# 快速浏览一下数据
print("特征数据 (X_train) 概览:")
# 将五项财务特征改用投影友好的短中文列名，只影响展示而不改变模型输入。
training_preview = X_train.head().rename(columns={'total_assets': '资产', 'total_liabilities': '负债', 'operating_revenue': '营收', 'net_profit_parent_company': '归母净利', 'cash_from_operating_activities': '经营现金流'})
# 展示短列名训练样本，核对五项指标及亿元口径。
display(training_preview)

# 为“步骤 3：训练候选集成模型”导入 `time` 并绑定 `time`，用于记录当前任务的实际运行耗时。
import time

# 初始化随机森林回归器
# 为了在低配置服务器上更快运行，我们减少了树的数量和深度
rf_model = RandomForestRegressor(
    n_estimators=50,        # 从100减少到50棵树
    max_depth=10,           # 限制树的最大深度
    min_samples_split=10,   # 保持原设置
    min_samples_leaf=5,     # 增加叶子节点的最小样本数
    # 固定随机森林的抽样与特征选择过程，确保结果可复算。
    random_state=42,
    n_jobs=-1              # 使用所有可用CPU核心
# 完成 `rf_model` 配置，以 50 棵限深树控制课堂运行时间与过拟合风险。
)

# 在计时前提示随机森林拟合即将开始。
print("正在训练随机森林模型...")
# 记录 `start_time` 的时间戳，用于计算该步骤的实际运行耗时。
start_time = time.time()

# 使用训练数据拟合模型
rf_model.fit(X_train, y_train)

# 记录 `end_time` 的时间戳，用于计算该步骤的实际运行耗时。
end_time = time.time()
# 报告随机森林训练耗时，说明后续解释计算的额外成本基准。
print(f"模型训练完成！用时: {end_time - start_time:.2f} 秒")

from sklearn.metrics import mean_squared_error, r2_score  # 用同一测试集计算候选与基线的可比指标
rf_predictions = rf_model.predict(X_test)  # 生成确定随机森林在测试行上的预测
training_mean_predictions = np.full(len(y_test), float(y_train.mean()))  # 以训练期目标均值生成事先确定常数基线
y_test_array = y_test.to_numpy(dtype=float)  # 固定测试响应的数组顺序供成对重抽样
rf_r2 = r2_score(y_test_array, rf_predictions)  # 计算随机森林测试集 R²
baseline_r2 = r2_score(y_test_array, training_mean_predictions)  # 计算训练均值基线测试集 R²
rf_mse = mean_squared_error(y_test_array, rf_predictions)  # 计算随机森林测试集 MSE
baseline_mse = mean_squared_error(y_test_array, training_mean_predictions)  # 计算训练均值基线测试集 MSE
bootstrap_rng = np.random.default_rng(2026)  # 固定重抽样随机种子以支持复算
bootstrap_rf_r2, bootstrap_delta_mse = [], []  # 保存有效重抽样的 R² 与成对 MSE 差
for bootstrap_iteration in range(1000):  # 重抽样测试行以近似当前切片的抽样不确定性
    bootstrap_index = bootstrap_rng.integers(0, len(y_test_array), len(y_test_array))  # 成对抽取响应和两组预测的同一行
    if np.ptp(y_test_array[bootstrap_index]) == 0:  # 跳过响应无变异、R² 无定义的重抽样
        continue  # 不把无定义的 R² 纳入区间
    bootstrap_rf_r2.append(r2_score(y_test_array[bootstrap_index], rf_predictions[bootstrap_index]))  # 保存当前重抽样的随机森林 R²
    bootstrap_delta_mse.append(mean_squared_error(y_test_array[bootstrap_index], rf_predictions[bootstrap_index]) - mean_squared_error(y_test_array[bootstrap_index], training_mean_predictions[bootstrap_index]))  # 保存 RF 相对基线的成对 MSE 差
assert len(bootstrap_rf_r2) >= 950 and np.isfinite(bootstrap_rf_r2).all() and np.isfinite(bootstrap_delta_mse).all()  # 拒绝有效重复不足或非有限结果
rf_r2_ci = np.quantile(bootstrap_rf_r2, [0.025, 0.975])  # 计算随机森林 R² 的百分位区间
delta_mse_ci = np.quantile(bootstrap_delta_mse, [0.025, 0.975])  # 计算 RF 减基线 MSE 的百分位区间
print({'RF测试R²': round(rf_r2, 4), '训练均值基线R²': round(baseline_r2, 4), 'RF测试MSE': round(rf_mse, 4), '训练均值基线MSE': round(baseline_mse, 4), 'RF R² paired-bootstrap 95%区间': np.round(rf_r2_ci, 4).tolist(), 'ΔMSE=RF−基线 95%区间': np.round(delta_mse_ci, 4).tolist()})  # 输出本次干净执行生成的比较证据

# 读取已拟合随机森林的 MDI 特征重要性向量，作为全局解释的模型内指标。
importances = rf_model.feature_importances_
# 为“步骤 4.1：全局视角 - 内置特征重要性”，固定 `features` 的特征名称，使图表标签与模型输入列一一对应。
feature_display_names = {
    'total_assets': '总资产',
    'total_liabilities': '总负债',
    'operating_revenue': '营业收入',
    'net_profit_parent_company': '归母净利润',
    'cash_from_operating_activities': '经营活动现金流'
}
features = pd.Index([feature_display_names.get(name, name) for name in X_train.columns])

# 按 MDI 数值升序取得特征索引，使水平条形图从低到高排列。
indices = np.argsort(importances)

# 创建 MDI 水平条形图画布，为排序后的特征标签和条长预留空间。
plt.figure(figsize=(12, 8))
# 将图题设为“随机森林模型计算的特征重要性”，直接说明当前图形的比较目的。
plt.title('随机森林模型计算的特征重要性', fontsize=18)
# 按重要性降序绘制 MDI 水平条，条长对应每个特征的随机森林不纯度贡献。
plt.barh(range(len(indices)), importances[indices], color='#003366', align='center')
# 将纵轴刻度标签设为与排序后 MDI 数值一一对应的特征名。
plt.yticks(range(len(indices)), [features[i] for i in indices])
# 将横轴标为“相对重要性”，明确横向编码的变量。
plt.xlabel('相对重要性', fontsize=14)
# 显示排序后的 MDI 重要性条形图，识别模型依赖特征并提醒该排序不表示因果贡献。
plt.show()

# 加载JS库以在notebook中美观地显示SHAP图
shap.initjs()

# 对于随机森林，使用TreeExplainer更加高效
# 同时设置较小的背景数据集以加速计算
background_sample = X_train.sample(n=min(100, len(X_train)), random_state=42)  # 从训练期抽取不超过100行作为背景
# 基于已拟合树模型建立 `explainer`，把预测拆解为基准值与特征贡献。
explainer = shap.TreeExplainer(rf_model, background_sample)

# 打印 TreeExplainer 背景样本的行列形状，核对解释基准只取自训练数据。
print(f"使用TreeExplainer，背景数据集大小: {background_sample.shape}")

# 为“步骤 5.2：计算SHAP值”导入 `time` 并绑定 `time`，用于记录当前任务的实际运行耗时。
import time

# 为了加快计算速度，我们只使用测试集的前500个样本
# 在配置较低的服务器上，这可以显著减少计算时间
X_test_sample = X_test.iloc[:500]

# 报告“原始测试集大小: {X_test.shape}”中的 `X_test`，核对“步骤 5.2：计算SHAP值”的样本形状。
print(f"原始测试集大小: {X_test.shape}")
# 报告“采样后测试集大小: {X_test_sample.shape}”中的 `X_test_sample`，核对“步骤 5.2：计算SHAP值”的样本形状。
print(f"采样后测试集大小: {X_test_sample.shape}")

# 记录SHAP计算时间
start_time = time.time()
# 在计时后提示采样测试集的 SHAP 计算即将开始。
print("正在计算SHAP值，请稍候...")

# 使用解释器计算采样测试集的SHAP值
shap_values = explainer(X_test_sample)
shap_values.feature_names = [feature_display_names.get(name, name) for name in X_test_sample.columns]
X_test_sample_display = X_test_sample.rename(columns=feature_display_names)

# 记录 `end_time` 的时间戳，用于计算该步骤的实际运行耗时。
end_time = time.time()
# 报告 SHAP 计算耗时，与模型训练耗时作量级比较。
print(f"SHAP计算完成！用时: {end_time - start_time:.2f} 秒")

# 打印 SHAP 解释结果的容器类型，核对当前版本返回的是解释对象而非预测数组。
print("SHAP值对象类型:", type(shap_values))
# 报告“SHAP值维度:”中的 `shap_values`，核对“步骤 5.2：计算SHAP值”的样本形状。
print("SHAP值维度:", shap_values.shape)

# 可视化测试集中第一个样本的预测解释
# shap_values[0] 表示我们只看第一个样本
shap.plots.waterfall(shap_values[0], show=False)
plt.gca().set_xlabel('SHAP 贡献值（模型输出尺度）')
plt.tight_layout()
plt.show()

# 绘制蜂群图来展示采样样本的SHAP值
shap.summary_plot(shap_values.values, X_test_sample_display, show=False)
summary_figure = plt.gcf()
summary_figure.axes[0].set_xlabel('SHAP 贡献值（模型输出尺度）')
if len(summary_figure.axes) > 1:
    summary_figure.axes[-1].set_ylabel('特征值：低 → 高')
plt.tight_layout()
plt.show()

# 绘制营业收入字段的依赖图，检查模型归因随输入变化的形状。
# 关闭自动交互选择，避免把着色变量误称为已验证交互
shap.dependence_plot('营业收入', shap_values.values, X_test_sample_display, interaction_index=None, show=False)
plt.gca().set_xlabel('营业收入（亿元）')
plt.gca().set_ylabel('营业收入的 SHAP 贡献值')
plt.tight_layout()
plt.show()

from pathlib import Path  # 管理公开数据路径
from urllib.request import urlretrieve  # 复用本章已安装的浏览器标识下载器
import numpy as np  # 处理无穷比率
import pandas as pd  # 合并审计意见与财务报表
# 按 Linux 共享数据、Windows 共享数据、项目缓存的顺序选择股票数据根目录。
stock_data_root = next((candidate_root for candidate_root in [Path('/home/ubuntu/r2_data_mount/data/stock'), Path('C:/qiufei/data/stock'), Path('data/course')] if candidate_root.exists()), Path('data/course'))
# 将审计意见文件定位到已选择的股票数据根目录。
audit_path = stock_data_root / 'audit_opinion.h5'
# 将财务报表文件定位到已选择的股票数据根目录。
financial_path = stock_data_root / 'financial_statement.h5'
# 在读者环境缺少数据根目录时创建项目缓存目录。
stock_data_root.mkdir(parents=True, exist_ok=True)
if not audit_path.exists():
    urlretrieve('https://assets.qiufei.site/data/stock/audit_opinion.h5', audit_path)
if not financial_path.exists():
    urlretrieve('https://assets.qiufei.site/data/stock/financial_statement.h5', financial_path)
audit_opinions = pd.read_hdf(audit_path, key='audit_opinion')
financial_columns = ['order_book_id', 'quarter', 'info_date', 'total_assets', 'total_liabilities', 'operating_revenue', 'net_profit_parent_company']  # 限定财务字段
audit_opinions['info_date'] = pd.to_datetime(audit_opinions['info_date'])  # 统一审计意见披露日类型
audit_opinions = audit_opinions.query("type == 'financial_statements' and quarter.str.endswith('q4')", engine='python').copy()  # 仅保留年度财报审计意见
audit_opinions['modified_opinion'] = audit_opinions['opinion_type'].ne('unqualified').astype(int)  # 构造风险标签
audit_opinions['target_year'] = audit_opinions['quarter'].str[:4].astype(int)  # 从目标季度提取财年
audit_opinions['decision_date'] = pd.to_datetime(audit_opinions['target_year'].astype(str) + '-06-30')  # 固定财年年中的预测时点
audit_opinions['label_deadline'] = pd.to_datetime((audit_opinions['target_year'] + 1).astype(str) + '-06-30')  # 确定次年六月末标签可得窗口
late_label_count = int((audit_opinions['info_date'] > audit_opinions['label_deadline']).sum())  # 统计晚于标签窗口的披露记录
audit_opinions = audit_opinions.query('info_date <= label_deadline').copy()  # 排除不能按题设及时观察的审计标签
duplicate_label_count = int(audit_opinions.duplicated(['order_book_id', 'target_year'], keep=False).sum())  # 统计公司财年重复标签记录
audit_opinions = audit_opinions.sort_values(['order_book_id', 'target_year', 'info_date']).drop_duplicates(['order_book_id', 'target_year'], keep='last')  # 保留窗口内最后披露版本
assert not audit_opinions.duplicated(['order_book_id', 'target_year']).any(), '审计标签键仍不唯一'  # 强制公司财年唯一性

annual_quarters = sorted((audit_opinions['target_year'] - 1).astype(str).add('q4').unique().tolist())  # 只读取可能进入特征的年度季度
financial_reports = pd.read_hdf(financial_path, key='financial_data', where='quarter in annual_quarters', columns=financial_columns)  # 按年度季度选择性读取财务列
financial_reports['info_date'] = pd.to_datetime(financial_reports['info_date'])  # 统一财报披露日类型
financial_reports['target_year'] = financial_reports['quarter'].str[:4].astype(int) + 1  # 将上一财年对齐到下一目标财年
audit_financial_panel = audit_opinions.merge(financial_reports, on=['order_book_id', 'target_year'], suffixes=('_audit', '_financial'))  # 按公司与相邻财年合并
audit_financial_panel = audit_financial_panel.query('info_date_financial <= decision_date').copy()  # 只保留决策日前已披露信息
audit_financial_panel['debt_ratio'] = audit_financial_panel['total_liabilities'] / audit_financial_panel['total_assets']  # 构造负债率
audit_financial_panel['net_margin'] = audit_financial_panel['net_profit_parent_company'] / audit_financial_panel['operating_revenue']  # 构造净利率
audit_financial_panel = audit_financial_panel.replace([np.inf, -np.inf], np.nan).dropna(subset=['debt_ratio', 'net_margin', 'modified_opinion'])  # 清除无效比率
audit_financial_panel = audit_financial_panel.sort_values(['order_book_id', 'target_year', 'info_date_financial']).drop_duplicates(['order_book_id', 'target_year'], keep='last')  # 同一公司财年仅保留最后可得特征版本
assert not audit_financial_panel.duplicated(['order_book_id', 'target_year']).any(), '建模面板键仍不唯一'  # 防止重复公司财年泄漏权重
print({'文件': ['audit_opinion.h5', 'financial_statement.h5'], '键': ['audit_opinion', 'financial_data'], '晚披露排除': late_label_count, '重复记录': duplicate_label_count})  # 输出数据说明与排除数量
print(audit_financial_panel[['target_year', 'info_date_financial', 'decision_date', 'info_date_audit']].head())  # 检查三类时点顺序
print({'样本期': (int(audit_financial_panel['target_year'].min()), int(audit_financial_panel['target_year'].max())), '公司数': int(audit_financial_panel['order_book_id'].nunique()), '样本数': len(audit_financial_panel)})  # 输出样本覆盖
print(audit_financial_panel['modified_opinion'].value_counts().sort_index())  # 报告类别计数

import shap  # 计算树模型的局部贡献
from sklearn.ensemble import GradientBoostingClassifier  # 建立梯度提升分类基准
from sklearn.metrics import average_precision_score, balanced_accuracy_score, confusion_matrix  # 报告不平衡分类指标
from sklearn.utils.class_weight import compute_sample_weight  # 平衡训练期少数类权重
feature_names = ['debt_ratio', 'net_margin']  # 固定可解释特征顺序
audit_feature_display_names = {'debt_ratio': '负债率', 'net_margin': '净利率'}
split_year = audit_financial_panel['target_year'].quantile(0.8, interpolation='lower')  # 按完整财年确定最终测试边界
training_panel = audit_financial_panel.query('target_year <= @split_year')  # 较早完整财年用于训练
testing_panel = audit_financial_panel.query('target_year > @split_year')  # 较晚完整财年作为最终测试期
assert len(training_panel) > 0 and len(testing_panel) > 0, '训练期或最终测试期为空'  # 在窗口无样本时明确停止
assert testing_panel['order_book_id'].nunique() >= 3, '最终测试期不足三家不同公司，无法兑现三项局部解释'  # 按唯一公司而非公司年行数自查对象数量
explanation_positions = testing_panel.reset_index(drop=True).drop_duplicates('order_book_id').index[:3].tolist()  # 固定三家不同公司的首条最终测试观测位置
assert training_panel['modified_opinion'].nunique() == 2 and testing_panel['modified_opinion'].nunique() == 2, '窗口必须同时包含两类标签'  # 避免输出不可定义指标
audit_risk_model = GradientBoostingClassifier(random_state=42)  # 固定模型随机状态
training_weights = compute_sample_weight('balanced', training_panel['modified_opinion'])  # 只从训练期估计类别权重
audit_risk_model.fit(training_panel[feature_names], training_panel['modified_opinion'], sample_weight=training_weights)  # 仅在训练期拟合
testing_probabilities = audit_risk_model.predict_proba(testing_panel[feature_names])[:, 1]  # 生成最终测试期风险概率
testing_predictions = (testing_probabilities >= 0.5).astype(int)  # 使用预先固定的0.5阈值分类
majority_predictions = np.repeat(training_panel['modified_opinion'].mode().iloc[0], len(testing_panel))  # 建立训练期多数类基准

print({'目标年': sorted(audit_financial_panel['target_year'].unique().tolist()), '切分年': int(split_year), '固定阈值': 0.5})  # 输出窗口和确定阈值
print({'样本数量变化': {'点时合并后': len(audit_financial_panel), '训练公司年': len(training_panel), '最终测试公司年': len(testing_panel)}, '公司数': audit_financial_panel['order_book_id'].nunique()})  # 明示小样本证据边界
print({'训练类别': training_panel['modified_opinion'].value_counts().sort_index().to_dict(), '测试类别': testing_panel['modified_opinion'].value_counts().sort_index().to_dict()})  # 输出两窗口类别分布
print('模型/多数类平衡准确率:', balanced_accuracy_score(testing_panel['modified_opinion'], testing_predictions), balanced_accuracy_score(testing_panel['modified_opinion'], majority_predictions))  # 比较排序外的分类表现
print('PR-AUC:', average_precision_score(testing_panel['modified_opinion'], testing_probabilities))  # 报告少数类排序表现
print('混淆矩阵:\n', confusion_matrix(testing_panel['modified_opinion'], testing_predictions))  # 报告误报与漏报
audit_background = training_panel[feature_names].sample(min(200, len(training_panel)), random_state=42)  # 固定训练背景分布
audit_explainer = shap.TreeExplainer(audit_risk_model, audit_background, model_output='raw', feature_perturbation='interventional')  # 在原始模型输出尺度解释条件贡献
audit_shap_values = audit_explainer(testing_panel[feature_names])  # 计算完整最终测试期解释对象
audit_shap_values.feature_names = [audit_feature_display_names[name] for name in feature_names]
print({'SHAP对象': type(audit_shap_values).__name__, '背景样本': len(audit_background), '输出尺度': 'raw', '扰动语义': 'interventional'})  # 报告解释器说明

stability_rows = []  # 收集预先固定种子下的归因摘要
for stability_seed in [7, 42, 2026]:  # 用三组预设种子重复训练而不改测试窗
    stability_model = GradientBoostingClassifier(random_state=stability_seed)  # 只改变模型训练种子
    stability_model.fit(training_panel[feature_names], training_panel['modified_opinion'], sample_weight=training_weights)  # 复用同一训练样本与权重
    stability_background = training_panel[feature_names].sample(min(25, len(training_panel)), random_state=stability_seed)  # 只从训练期重抽背景
    stability_explainer = shap.TreeExplainer(stability_model, stability_background, model_output='raw', feature_perturbation='interventional')  # 保持原始输出尺度
    stability_values = stability_explainer(testing_panel[feature_names]).values  # 对同一六个最终测试观测计算归因
    for feature_position, feature_name in enumerate(feature_names):  # 逐特征记录方向与平均绝对贡献
        stability_rows.append({'随机种子': stability_seed, '财务指标': audit_feature_display_names[feature_name], '平均贡献': stability_values[:, feature_position].mean(), '平均绝对贡献': np.abs(stability_values[:, feature_position]).mean()})
stability_table = pd.DataFrame(stability_rows)  # 形成种子×特征结果表
stability_table['排序'] = stability_table.groupby('随机种子')['平均绝对贡献'].rank(ascending=False, method='dense')
display(stability_table.round(5))  # 如实展示排序或方向是否随种子变化

shap.plots.beeswarm(audit_shap_values, show=False)  # 展示最终测试期特征归因分布
audit_global_figure = plt.gcf()
audit_global_figure.axes[0].set_xlabel('SHAP 贡献值（模型输出尺度）')
if len(audit_global_figure.axes) > 1:
    audit_global_figure.axes[-1].set_ylabel('特征值：低 → 高')
plt.tight_layout()  # 收紧图形边界以适应投影画布
plt.show()  # 输出一幅全局解释图

company_position = explanation_positions[0]  # 固定第一家不同公司的最终测试期解释对象
company_record = testing_panel.iloc[company_position]  # 读取公司、目标年与真实标签
company_raw_output = float(audit_risk_model.decision_function(testing_panel[feature_names].iloc[[company_position]])[0])  # 取得与SHAP相同的原始对数几率尺度
company_shap_total = float(audit_shap_values.base_values[company_position] + audit_shap_values.values[company_position].sum())  # 在raw尺度累加基准与贡献
display(pd.DataFrame({'财务指标': audit_shap_values.feature_names, 'SHAP 贡献值': audit_shap_values.values[company_position]}).round(5))
shap.plots.waterfall(audit_shap_values[company_position], show=False)  # 绘制第一家公司贡献分解
plt.gca().set_xlabel('SHAP 贡献值（模型输出尺度）')
plt.tight_layout()  # 收紧当前瀑布图边界
plt.show()  # 单独输出第一幅局部解释图

company_position = explanation_positions[1]  # 固定第二家不同公司的最终测试期解释对象
company_record = testing_panel.iloc[company_position]  # 读取公司、目标年与真实标签
company_raw_output = float(audit_risk_model.decision_function(testing_panel[feature_names].iloc[[company_position]])[0])  # 取得与SHAP相同的原始对数几率尺度
company_shap_total = float(audit_shap_values.base_values[company_position] + audit_shap_values.values[company_position].sum())  # 在raw尺度累加基准与贡献
display(pd.DataFrame({'财务指标': audit_shap_values.feature_names, 'SHAP 贡献值': audit_shap_values.values[company_position]}).round(5))
shap.plots.waterfall(audit_shap_values[company_position], show=False)  # 绘制第二家公司贡献分解
plt.gca().set_xlabel('SHAP 贡献值（模型输出尺度）')
plt.tight_layout()  # 收紧当前瀑布图边界
plt.show()  # 单独输出第二幅局部解释图

company_position = explanation_positions[2]  # 固定第三家不同公司的最终测试期解释对象
company_record = testing_panel.iloc[company_position]  # 读取公司、目标年与真实标签
company_raw_output = float(audit_risk_model.decision_function(testing_panel[feature_names].iloc[[company_position]])[0])  # 取得与SHAP相同的原始对数几率尺度
company_shap_total = float(audit_shap_values.base_values[company_position] + audit_shap_values.values[company_position].sum())  # 在raw尺度累加基准与贡献
display(pd.DataFrame({'财务指标': audit_shap_values.feature_names, 'SHAP 贡献值': audit_shap_values.values[company_position]}).round(5))
shap.plots.waterfall(audit_shap_values[company_position], show=False)  # 绘制第三家公司贡献分解
plt.gca().set_xlabel('SHAP 贡献值（模型输出尺度）')
plt.tight_layout()  # 收紧当前瀑布图边界
plt.show()  # 单独输出第三幅局部解释图