D. Blue-Red Permutation

本文介绍了一种通过调整数组中可变颜色数字实现特定排列的方法。对于可加1的红色数字,期望其能组成接近n的序列;对于可减1的蓝色数字,则期望其组成接近1的序列。通过对两组数字分别排序并检查,验证了数组是否能形成1至n的唯一序列。

1、D. Blue-Red Permutation

题意

给了你 一个大小为n的数组,没个位置都有一个数字和一个颜色,如果这个地方的颜色为红色, 那么这个位置的数字可以+1,如果为蓝色,那么这个地方的颜色可以-1,现在问你是否可以让这个数组中的数字都是1~n,并且只出现过一次。

思路

首先对于能够增加的数字,我们肯定希望他是能够组成……n-3、n-2、n-1、n的,对于能够减的数字,我们是希望它最后是能够组成1、2、3……,那么我们将红色的数字和蓝色的数字分别放进两个容器,单独对他们进行改变,我们对红色的数字进行从小到大进行排序,对于蓝色的数字我们从大到小进行排序,然后单独的对他们进行判断。

#include <bits/stdc++.h>

#define int long long

using namespace std;

const int N = 2e5 + 10;

int n;
int a[N];
string s;

void solve()
{
    cin >> n;
    for (int i = 1; i <= n ; i++) cin >> a[i];
    cin >> s; s = " " + s;

    vector<int> v[2];

    for (int i = 1; i <= n; i ++)
    {
        if (s[i] == 'B') v[0].push_back(a[i]);
        else v[1].push_back(a[i]);
    }

    sort(v[0].begin(), v[0].end());
    sort(v[1].begin(), v[1].end());

    for (int i = 0; i < v[0].size(); i ++)
        if (v[0][i] < i + 1) { cout << "NO" << endl; return ; }
    for (int i = v[1].size() - 1, t = n; i >= 0; i --, t --)
        if (v[1][i] > t) { cout << "NO" << endl; return ;}

    cout << "YES" << endl;
}

signed main()
{
    std::ios::sync_with_stdio(false);

    int t; cin >> t;
    while (t --) solve();

    return 0;
}

## 题目重述 实验数据位于文件 `regress_data1.xlsx` 中,要求完成以下任务: 1. 绘制数据散点图,其中横轴为“人口”,纵轴为“收益”; 2. 构建线性回归模型 $ y = ax + b $,使用最小二乘法求解参数 $ a $ 和 $ b $,采用两种方法: - 使用 NumPy 内置函数; - 使用 Python 循环手动实现; 3. 构建相同模型,以最小平方误差为损失函数,使用梯度下降算法求解参数: - 批量梯度下降(Batch GD); - 随机梯度下降(Stochastic GD); - 小批量梯度下降(Mini-batch GD); 4. 使用 `sklearn` 建立回归模型,并评估其性能(如 $ R^2 $、MSE 等)。 --- ## 给出答案(带注释的完整 Python 代码) ```python import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error, r2_score import warnings # 忽略字体警告和收敛警告 warnings.filterwarnings("ignore") plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 # ----------------------------- # 1. 加载数据 # ----------------------------- # 读取 Excel 文件中的数据 data = pd.read_excel('regress_data1.xlsx') # 提取“人口”和“收益”列 x = data['人口'].values # 形状: (n,) y = data['收益'].values # 形状: (n,) n = len(x) # ----------------------------- # 2. 绘制散点图 # ----------------------------- plt.figure(figsize=(10, 6)) plt.scatter(x, y, color='blue', alpha=0.7, label='数据点') plt.xlabel('人口') plt.ylabel('收益') plt.title('人口 vs 收益 散点图') plt.grid(True, linestyle='--', alpha=0.5) plt.legend() plt.show() # ----------------------------- # 3. 最小二乘法求解 y = ax + b # ----------------------------- # 方法一:使用 NumPy 的 polyfit 函数 a_numpy, b_numpy = np.polyfit(x, y, 1) # 方法二:使用循环手动实现最小二乘法 x_mean = np.mean(x) y_mean = np.mean(y) numerator = 0.0 denominator = 0.0 for i in range(n): numerator += (x[i] - x_mean) * (y[i] - y_mean) denominator += (x[i] - x_mean) ** 2 a_loop = numerator / denominator b_loop = y_mean - a_loop * x_mean print("【最小二乘法结果】") print(f"NumPy 函数法: a = {a_numpy:.4f}, b = {b_numpy:.4f}") print(f"循环实现法: a = {a_loop:.4f}, b = {b_loop:.4f}") # ----------------------------- # 4. 梯度下降算法 # ----------------------------- def add_intercept(X): """添加截距项(x0=1)""" return np.column_stack([np.ones(len(X)), X]) X = add_intercept(x) # (n, 2): 第一列为1,第二列为x y = y.reshape(-1, 1) # (n, 1) # 参数初始化 np.random.seed(42) theta = np.random.randn(2, 1) * 0.1 learning_rate = 0.0001 epochs = 10000 # 批量梯度下降 (BGD) theta_bgd = theta.copy() for _ in range(epochs): gradient = (1/n) * X.T @ (X @ theta_bgd - y) theta_bgd -= learning_rate * gradient a_bgd, b_bgd = theta_bgd[1, 0], theta_bgd[0, 0] # 随机梯度下降 (SGD) theta_sgd = theta.copy() for epoch in range(epochs): for i in range(n): xi = np.array([[1, x[i]]]) # (1, 2) yi = y[i] # (1, 1) gradient = xi.T @ (xi @ theta_sgd - yi) theta_sgd -= learning_rate * gradient a_sgd, b_sgd = theta_sgd[1, 0], theta_sgd[0, 0] # 小批量梯度下降 (MBGD),batch_size=16 theta_mbgd = theta.copy() batch_size = 16 for epoch in range(epochs): indices = np.random.permutation(n) X_shuffled = X[indices] y_shuffled = y[indices] for start in range(0, n, batch_size): end = start + batch_size X_batch = X_shuffled[start:end] y_batch = y_shuffled[start:end] m = len(X_batch) gradient = (1/m) * X_batch.T @ (X_batch @ theta_mbgd - y_batch) theta_mbgd -= learning_rate * gradient a_mbgd, b_mbgd = theta_mbgd[1, 0], theta_mbgd[0, 0] print("\n【梯度下降结果】") print(f"批量梯度下降: a = {a_bgd:.4f}, b = {b_bgd:.4f}") print(f"随机梯度下降: a = {a_sgd:.4f}, b = {b_sgd:.4f}") print(f"小批量梯度下降: a = {a_mbgd:.4f}, b = {b_mbgd:.4f}") # ----------------------------- # 5. 使用 sklearn 建模与评估 # ----------------------------- model_sk = LinearRegression() model_sk.fit(x.reshape(-1, 1), y.ravel()) y_pred_sk = model_sk.predict(x.reshape(-1, 1)) mse = mean_squared_error(y, y_pred_sk) r2 = r2_score(y, y_pred_sk) print("\n【sklearn 模型结果】") print(f"斜率 a = {model_sk.coef_[0]:.4f}, 截距 b = {model_sk.intercept_:.4f}") print(f"均方误差 MSE = {mse:.4f}") print(f"决定系数 R² = {r2:.4f}") # ----------------------------- # 6. 绘制所有回归线对比图 # ----------------------------- plt.figure(figsize=(12, 7)) plt.scatter(x, y, color='blue', alpha=0.6, label='原始数据') x_plot = np.linspace(x.min(), x.max(), 100) plt.plot(x_plot, a_numpy * x_plot + b_numpy, color='green', label='最小二乘法(NumPy)', linewidth=2) plt.plot(x_plot, a_bgd * x_plot + b_bgd, color='orange', linestyle='-.', label='批量梯度下降', linewidth=2) plt.plot(x_plot, model_sk.coef_[0] * x_plot + model_sk.intercept_, color='red', label='Sklearn 模型', linewidth=2) plt.xlabel('人口') plt.ylabel('收益') plt.title('不同方法拟合的线性回归模型对比') plt.legend() plt.grid(True, alpha=0.5) plt.tight_layout() plt.show() 修缮代码使可视化图片的文字可以显示成功
最新发布
10-27
# 导入所需的库 import pandas as pd # 导入pandas库,用于数据处理和分析,特别是DataFrame操作 import numpy as np # 导入numpy库,用于进行数值计算,特别是数组操作 import matplotlib.pyplot as plt # 导入matplotlib的pyplot模块,用于绘制图表 import matplotlib # 导入matplotlib主库,用于更底层的绘图设置 from matplotlib.patches import Patch # 从matplotlib.patches中导入Patch,用于创建图例中的色块 from sklearn.model_selection import train_test_split, GridSearchCV # 从sklearn导入数据划分和网格搜索交叉验证的工具 from sklearn.preprocessing import StandardScaler # 从sklearn导入数据标准化工具 import xgboost as xgb # 【XGBoost 修改】导入XGBoost库,用于构建梯度提升决策树模型 from sklearn.metrics import mean_absolute_error, r2_score, mean_squared_error # 从sklearn导入评估回归模型性能的指标 import joblib # 导入joblib库,用于模型的保存和加载 from scipy.stats import gaussian_kde # 从scipy.stats导入高斯核密度估计,可用于绘制密度图(此脚本中未直接使用,但可能为备用) from sklearn.inspection import PartialDependenceDisplay, partial_dependence # 从sklearn导入部分依赖图的工具(此脚本未使用,改用手动实现) import time # 导入time库,用于计时(此脚本中未直接使用) import os # 导入os库,用于操作系统相关功能,如创建文件夹 from itertools import combinations # 从itertools导入combinations,用于生成组合 import shap # 导入shap库,用于模型解释,计算SHAP值 import warnings # 导入warnings库,用于控制警告信息的显示 from collections import defaultdict # 从collections导入defaultdict,用于创建带有默认值的字典 from PyALE import ale # 导入PyALE库,用于计算和绘制累积局部效应图 from scipy.interpolate import griddata # 从scipy.interpolate导入griddata,用于插值(此脚本中未直接使用) # --- 全局设置 --- # 忽略特定类型的警告,避免在输出中显示不必要的警告信息 warnings.filterwarnings("ignore", category=FutureWarning, module="sklearn.utils._bunch") warnings.filterwarnings("ignore", category=UserWarning) matplotlib.use('TkAgg') # 设置matplotlib的后端,'TkAgg'是一个图形界面后端,确保在某些环境下可以正常显示绘图窗口 import os import tempfile # 设置临时文件夹到一个英文路径 os.environ['JOBLIB_TEMP_FOLDER'] = 'C:/temp' os.environ['TMP'] = 'C:/temp' os.environ['TEMP'] = 'C:/temp' # 确保这个目录存在 if not os.path.exists('C:/temp'): os.makedirs('C:/temp') # --- (函数定义区) --- # 定义一个函数,用于绘制回归模型的拟合效果图 def plot_regression_fit(y_true, y_pred, r2, rmse, mae, data_label_en, title_en, save_path): """ 绘制真实值与预测值的散点图,并显示模型评估指标。 y_true: 真实值 y_pred: 预测值 r2: R-squared值 rmse: 均方根误差 mae: 平均绝对误差 data_label_en: 数据集标签 (如 'Train Set') title_en: 图表标题 save_path: 图片保存路径 """ plt.style.use('seaborn-v0_8-whitegrid') # 使用预设的绘图风格 plt.rcParams['font.family'] = 'sans-serif' # 设置字体为无衬线字体,以获得更好的显示效果 fig, ax = plt.subplots(figsize=(7, 7)) # 创建一个7x7英寸的画布和子图 # 绘制真实值 vs 预测值的散点图 ax.scatter(y_true, y_pred, alpha=0.6, edgecolors='k', label=f'{data_label_en} (n={len(y_true)})') # 计算并设置坐标轴的范围,确保1:1线能完整显示 lims = [np.min([y_true.min(), y_pred.min()]) - 5, np.max([y_true.max(), y_pred.max()]) + 5] # 绘制1:1参考线 (y=x),代表完美预测 ax.plot(lims, lims, 'k--', alpha=0.75, zorder=0, label='1:1 Line (Perfect Fit)') ax.set_aspect('equal') # 设置x和y轴的比例相等 ax.set_xlim(lims) # 设置x轴范围 ax.set_ylim(lims) # 设置y轴范围 y_true_np = np.array(y_true) # 将真实值转换为numpy数组 y_pred_np = np.array(y_pred) # 将预测值转换为numpy数组 m, b = np.polyfit(y_true_np, y_pred_np, 1) # 对散点进行线性拟合,得到斜率m和截距b ax.plot(y_true_np, m * y_true_np + b, 'r-', label='Linear Fit') # 绘制线性拟合线 ax.set_xlabel('True Values (%)', fontsize=12) # 设置x轴标签 ax.set_ylabel('Predicted Values (%)', fontsize=12) # 设置y轴标签 ax.set_title(title_en, fontsize=14, weight='bold') # 设置图表标题 # 准备要在图上显示的评估指标文本 metrics_text = f'$R^2 = {r2:.3f}$\n$RMSE = {rmse:.3f}$\n$MAE = {mae:.3f}$' # 在图的左上角添加文本框,显示评估指标 ax.text(0.05, 0.95, metrics_text, transform=ax.transAxes, fontsize=12, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) ax.legend(loc='lower right') # 在右下角显示图例 fig.savefig(save_path, dpi=200, bbox_inches='tight') # 保存图表到指定路径 plt.close(fig) # 关闭图表,释放内存 # 定义一个函数,用于绘制组合特征重要性图(条形图+甜甜圈图) def plot_importance_combined(df_importance, title, save_path, bar_color='dodgerblue'): """ 绘制特征重要性条形图,并在右下角嵌入一个显示Top-N特征占比的甜甜圈图。 df_importance: 包含'Feature'和'Importance'两列的DataFrame title: 图表标题 save_path: 图片保存路径 bar_color: 条形图的颜色 """ df_importance_sorted = df_importance.sort_values(by='Importance', ascending=True) # 按重要性升序排序 plt.rc("font", family='Microsoft YaHei') # 设置中文字体为微软雅黑,以正确显示中文 fig, ax = plt.subplots(figsize=(14, 10)) # 创建一个14x10英寸的画布和子图 # 绘制水平条形图 bars = ax.barh(df_importance_sorted['Feature'], df_importance_sorted['Importance'], color=bar_color, alpha=0.8) ax.set_title(title, fontsize=18, pad=20) # 设置标题 ax.set_xlabel('重要性得分', fontsize=14) # 设置x轴标签 ax.set_ylabel('特征', fontsize=14) # 设置y轴标签 ax.tick_params(axis='both', which='major', labelsize=12) # 设置刻度标签的大小 ax.grid(axis='x', linestyle='--', alpha=0.6) # 显示x轴方向的网格线 # 在每个条形图旁边显示具体的重要性数值 for bar in bars: width = bar.get_width() ax.text(width, bar.get_y() + bar.get_height() / 2, f' {width:.4f}', va='center', ha='left', fontsize=10) ax.set_xlim(right=ax.get_xlim()[1] * 1.2) # 调整x轴范围,为数值标签留出空间 N_DONUT_FEATURES = 5 # 设置甜甜圈图中要显示的特征数量 if len(df_importance) < N_DONUT_FEATURES: # 如果总特征数小于5,则取全部特征 N_DONUT_FEATURES = len(df_importance) df_desc = df_importance.sort_values(by='Importance', ascending=False) # 按重要性降序排序 top_n_features = df_desc.head(N_DONUT_FEATURES) # 选取最重要的N个特征 donut_feature_names = top_n_features['Feature'].tolist() # 获取这N个特征的名称 # 如果有特征且总重要性大于0,则绘制甜甜圈图 if not top_n_features.empty and top_n_features['Importance'].sum() > 0: total_donut_importance = top_n_features['Importance'].sum() # 计算Top-N特征的重要性总和 donut_percentages = top_n_features['Importance'] / total_donut_importance * 100 # 计算每个特征在Top-N中的百分比 ax_inset = fig.add_axes([0.45, 0.15, 0.28, 0.28]) # 在主图上创建一个嵌入的子图(甜甜圈图的位置) colors = matplotlib.colormaps['tab10'].colors # 获取一组颜色 # 绘制饼图(通过设置width属性使其变为甜甜圈图) wedges, _ = ax_inset.pie( donut_percentages, colors=colors[:len(top_n_features)], startangle=90, counterclock=False, wedgeprops=dict(width=0.45, edgecolor='w') ) # 计算Top-N特征占总特征重要性的比例 subset_importance_ratio = top_n_features['Importance'].sum() / df_importance['Importance'].sum() # 在甜甜圈中心添加文本 ax_inset.text(0, 0, f'Top {N_DONUT_FEATURES} 特征\n占总重要性\n{subset_importance_ratio:.2%}', ha='center', va='center', fontsize=9, linespacing=1.4) label_threshold = 3.0 # 设置标签显示的阈值,小于此值的百分比会用引导线引出 y_text_offsets = {'left': 1.4, 'right': 1.4} # 初始化引导线标签的垂直偏移量 # 为每个扇区添加百分比标签 for i, p in enumerate(wedges): percent = donut_percentages.iloc[i] ang = (p.theta2 - p.theta1) / 2. + p.theta1 # 计算扇区中心角度 y = np.sin(np.deg2rad(ang)) # 计算标签的y坐标 x = np.cos(np.deg2rad(ang)) # 计算标签的x坐标 # 如果百分比小于阈值,使用引导线 if percent < label_threshold and percent > 0: side = 'right' if x > 0 else 'left' # 判断标签在左侧还是右侧 y_pos = y_text_offsets[side] # 获取当前侧的y偏移 y_text_offsets[side] += -0.2 if y > 0 else 0.2 # 更新偏移量,避免标签重叠 connectionstyle = f"angle,angleA=0,angleB={ang}" # 设置引导线样式 # 添加带引导线的注释 ax_inset.annotate(f'{percent:.1f}%', xy=(x, y), xytext=(0.2 * np.sign(x), y_pos), fontsize=9, ha='center', arrowprops=dict(arrowstyle="-", connectionstyle=connectionstyle, relpos=(0.5, 0.5))) # 如果百分比大于阈值,直接在扇区内显示 elif percent > 0: ax_inset.text(x * 0.8, y * 0.8, f'{percent:.1f}%', ha='center', va='center', fontsize=9, fontweight='bold', color='white') # 在甜甜圈图右侧添加图例 ax_inset.legend(wedges, donut_feature_names, loc="center left", bbox_to_anchor=(1.2, 0.5), frameon=False, fontsize=11) plt.savefig(save_path, dpi=720, bbox_inches='tight') # 保存高分辨率图表 plt.close(fig) # 关闭图表,释放内存 # 定义一个函数,用于绘制残差图 def plot_residuals_styled(residuals, y_pred, save_path, title): """ 绘制残差与预测值的关系图,并高亮显示异常值。 residuals: 残差 (真实值 - 预测值) y_pred: 预测值 save_path: 图片保存路径 title: 图表标题 """ plt.style.use('seaborn-v0_8-whitegrid') # 使用预设绘图风格 plt.rc("font", family='Microsoft YaHei') # 设置中文字体 fig, ax = plt.subplots(figsize=(10, 8)) # 创建一个10x8英寸的画布 sd_residuals = np.std(residuals) # 计算残差的标准差 is_outlier = np.abs(residuals) > 2 * sd_residuals # 定义异常值:绝对残差大于2倍标准差 num_outliers = np.sum(is_outlier) # 计算异常值的数量 print(f"在数据集 '{title}' 中, 共找到 {num_outliers} 个异常值 (残差 > 2 S.D.)。") # 打印异常值信息 sd_label = f'S.D. (±{sd_residuals:.2f})' # 准备标准差区间的图例标签 ax.axhspan(-sd_residuals, sd_residuals, color='yellow', alpha=0.3, label=sd_label) # 绘制一个表示一个标准差范围的水平区域 # 绘制正常值的散点图 ax.scatter(y_pred[~is_outlier], residuals[~is_outlier], alpha=0.6, c='green', edgecolors='k', linewidth=0.5, s=50, label='正常值') # 绘制异常值的散点图 ax.scatter(y_pred[is_outlier], residuals[is_outlier], alpha=0.8, c='red', edgecolors='k', linewidth=0.5, s=70, label='异常值 (> 2 S.D.)') ax.axhline(0, color='black', linestyle='--', linewidth=1.5) # 绘制残差为0的参考线 ax.set_title(title, fontsize=16, weight='bold') # 设置标题 ax.set_xlabel('预测值', fontsize=14) # 设置x轴标签 ax.set_ylabel('残差 (真实值 - 预测值)', fontsize=14) # 设置y轴标签 # 设置图表边框样式 for spine in ax.spines.values(): spine.set_visible(True) spine.set_color('black') spine.set_linewidth(1) ax.legend(loc='upper right', fontsize=12) # 显示图例 y_max = np.max(np.abs(residuals)) * 1.2 # 计算y轴的范围 ax.set_ylim(-y_max, y_max) # 设置y轴范围 plt.tight_layout() # 调整布局,防止标签重叠 plt.savefig(save_path, dpi=300, bbox_inches='tight') # 保存图表 plt.close() # 关闭图表 # ========================================================================= # ======================= 新增:手动PDP计算与绘图函数 ======================= # ========================================================================= def manual_pdp_1d(model, X_data, feature_name, grid_resolution=50): """ 手动计算一维部分依赖(PDP)和个体条件期望(ICE)数据。 返回: - grid_values: 特征的网格点。 - pdp_values: 对应的PDP值 (ICE的均值)。 - ice_lines: 所有样本的ICE线数据。 """ # 在特征的最小值和最大值之间生成一系列网格点 grid_values = np.linspace(X_data[feature_name].min(), X_data[feature_name].max(), grid_resolution) # 初始化一个数组来存储每个样本的ICE线数据 ice_lines = np.zeros((len(X_data), grid_resolution)) # 遍历数据集中每一个样本 for i, (_, sample) in enumerate(X_data.iterrows()): # 创建一个临时DataFrame,行数为网格点数,内容为当前样本的重复 temp_df = pd.DataFrame([sample] * grid_resolution) # 将要分析的特征列替换为网格值 temp_df[feature_name] = grid_values # 使用模型进行预测,得到这个样本在不同特征值下的预测结果,即ICE线 ice_lines[i, :] = model.predict(temp_df) # PDP是所有ICE线的平均值,在每个网格点上求均值 pdp_values = np.mean(ice_lines, axis=0) # 返回计算结果 return grid_values, pdp_values, ice_lines def manual_pdp_2d(model, X_data, features_tuple, grid_resolution=20): """ 手动计算二维部分依赖(PDP)数据。 返回: - grid_1: 第一个特征的网格点。 - grid_2: 第二个特征的网格点。 - pdp_values: 二维网格上对应的PDP值。 """ feat1_name, feat2_name = features_tuple # 获取两个特征的名称 # 为第一个特征生成网格点 grid_1 = np.linspace(X_data[feat1_name].min(), X_data[feat1_name].max(), grid_resolution) # 为第二个特征生成网格点 grid_2 = np.linspace(X_data[feat2_name].min(), X_data[feat2_name].max(), grid_resolution) # 初始化一个二维数组来存储PDP值 pdp_values = np.zeros((grid_resolution, grid_resolution)) # 遍历二维网格的每一个点 for i, val1 in enumerate(grid_1): for j, val2 in enumerate(grid_2): # 创建一个原始数据的临时副本 X_temp = X_data.copy() # 将第一个特征的所有值都设为当前的网格点值 X_temp[feat1_name] = val1 # 将第二个特征的所有值都设为当前的网格点值 X_temp[feat2_name] = val2 # 对修改后的整个数据集进行预测 preds = model.predict(X_temp) # 计算预测结果的平均值,作为该网格点的PDP值 pdp_values[j, i] = np.mean(preds) # 返回计算结果 return grid_1, grid_2, pdp_values def plot_3d_scatter_three_features(X_data, y_pred, features_tuple, save_path): """ 绘制三个特征的3D散点图,并用预测值对散点进行着色。 """ feat1_name, feat2_name, feat3_name = features_tuple # 获取三个特征的名称 fig = plt.figure(figsize=(12, 9)) # 创建一个12x9英寸的画布 ax = fig.add_subplot(111, projection='3d') # 添加一个3D子图 # 绘制3D散点图,x,y,z轴分别是三个特征的值,颜色c由模型预测值决定 sc = ax.scatter( X_data[feat1_name], X_data[feat2_name], X_data[feat3_name], c=y_pred, cmap='viridis', s=30, alpha=0.7, edgecolor='k', linewidth=0.5 ) ax.set_xlabel(f'{feat1_name} (标准化值)', fontsize=10, labelpad=10) # 设置x轴标签 ax.set_ylabel(f'{feat2_name} (标准化值)', fontsize=10, labelpad=10) # 设置y轴标签 ax.set_zlabel(f'{feat3_name} (标准化值)', fontsize=10, labelpad=10, rotation=180) # 设置z轴标签 ax.set_title(f'三特征3D散点图\n({feat1_name}, {feat2_name}, {feat3_name})', fontsize=14) # 设置标题 # 添加颜色条,并设置标签 cbar = fig.colorbar(sc, shrink=0.5, aspect=20, label='模型预测值', pad=0.1) ax.view_init(elev=20, azim=45) # 设置3D视图的角度 plt.savefig(save_path, dpi=300) # 保存图表 plt.close(fig) # 关闭图表 print(f"成功绘制 3D 散点图 for {features_tuple}") # 打印成功信息 def plot_3d_pdp_fixed_value(model, X_data, features, save_path, fixed_feature=None, fixed_value=None, grid_resolution=50): """ 绘制三个特征的3D PDP图,其中一个特征被固定在特定值。 """ feature_1, feature_2, feature_3 = features # 获取三个特征名称 if fixed_feature is None: # 如果没有指定要固定的特征 fixed_feature = feature_3 # 默认固定第三个特征 # 找出需要变化的两个特征 varying_features = [f for f in features if f != fixed_feature] varying_feat_1, varying_feat_2 = varying_features[0], varying_features[1] if fixed_value is None: # 如果没有指定固定的值 fixed_value = X_data[fixed_feature].mean() # 默认使用该特征的平均值 # 为两个变化的特征生成网格点 feat1_vals = np.linspace(X_data[varying_feat_1].min(), X_data[varying_feat_1].max(), grid_resolution) feat2_vals = np.linspace(X_data[varying_feat_2].min(), X_data[varying_feat_2].max(), grid_resolution) XX, YY = np.meshgrid(feat1_vals, feat2_vals) # 创建二维网格 # 使用数据集中所有特征的平均值作为背景行,以减少其他特征的影响 background_row = X_data.mean().to_dict() # 创建一个包含网格点组合的DataFrame grid_data = pd.DataFrame(np.c_[XX.ravel(), YY.ravel()], columns=varying_features) # 创建一个用于预测的网格DataFrame,以背景行作为基础 X_grid = pd.DataFrame([background_row] * len(grid_data)) # 将变化的特征列替换为网格值 X_grid[varying_feat_1] = grid_data[varying_feat_1].values X_grid[varying_feat_2] = grid_data[varying_feat_2].values # 将固定的特征列设置为指定的值 X_grid[fixed_feature] = fixed_value X_grid = X_grid[X_data.columns] # 确保列顺序与训练时一致 preds = model.predict(X_grid).reshape(XX.shape) # 进行预测,并重塑为网格形状 plt.rc("font", family='Microsoft YaHei') # 设置中文字体 fig = plt.figure(figsize=(12, 9)) # 创建画布 ax = fig.add_subplot(111, projection='3d') # 创建3D子图 # 绘制3D曲面图 surface = ax.plot_surface(XX, YY, preds, cmap='viridis', alpha=0.9, edgecolor='k', linewidth=0.2) fig.colorbar(surface, ax=ax, shrink=0.6, aspect=20, label='模型预测值') # 添加颜色条 ax.set_xlabel(f'{varying_feat_1} (标准化值)', fontsize=10, labelpad=10) # 设置x轴标签 ax.set_ylabel(f'{varying_feat_2} (标准化值)', fontsize=10, labelpad=10) # 设置y轴标签 ax.set_zlabel('模型预测值', fontsize=10, labelpad=10, rotation=90) # 设置z轴标签 # 设置标题 title_text = f'3D依赖图: {varying_feat_1} vs {varying_feat_2}\n固定 {fixed_feature} = {fixed_value:.3f}' ax.set_title(title_text, fontsize=14) ax.view_init(elev=25, azim=-120) # 设置视角 plt.savefig(save_path, dpi=300) # 保存图片 plt.close(fig) # 关闭图表 print(f"成功绘制 固定值3D PDP for {features},固定 {fixed_feature}") # 打印成功信息 print('-------------------------------------准备数据---------------------------------------') # 从指定的Excel文件中读取数据 df = pd.read_excel(r'D:\巢湖流域.xlsx') y = df.iloc[:, 0] # 提取第一列作为目标变量y x = df.iloc[:, 1:] # 提取从第二列开始的所有列作为特征变量x feature_names_from_df = x.columns.tolist() # 获取特征名称列表 print('-------------------------------------划分数据集---------------------------------------') # 将数据集划分为训练集和测试集,测试集占30%,设置随机种子以保证结果可复现 x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=12) print('-------------------------------------数据标准化---------------------------------------') scaler = StandardScaler() # 实例化一个StandardScaler对象 # 对训练集进行拟合和转换,并将结果转换回DataFrame X_train_scaled_df = pd.DataFrame(scaler.fit_transform(x_train), columns=x_train.columns, index=x_train.index) # 使用在训练集上学习到的参数对测试集进行转换 X_test_scaled_df = pd.DataFrame(scaler.transform(x_test), columns=x_test.columns, index=x_test.index) print('-------------------------------------定义XGBoost模型超参数范围---------------------------------------') # 定义要进行网格搜索的XGBoost超参数范围 xgb_param_grid = { 'n_estimators': range(50, 151, 50), # 树的数量:100, 200, 300 'max_depth': range(3, 7, 1), # 树的最大深度:3, 4, 5, 6 'learning_rate': [0.01, 0.1, 0.2] # 学习率 } print('-------------------------------------搜索最佳超参数---------------------------------------') # 实例化GridSearchCV对象,用于自动寻找最佳超参数组合 gd = GridSearchCV(estimator=xgb.XGBRegressor(objective='reg:squarederror', seed=0), # 使用XGBoost回归器 param_grid=xgb_param_grid, # 指定超参数网格 cv=5, # 使用5折交叉验证 n_jobs=-1, # 使用所有可用的CPU核心进行并行计算 verbose=0) # 不输出详细的搜索过程信息 gd.fit(X_train_scaled_df, y_train) # 在标准化的训练集上执行网格搜索 print('-------------------------------------输出最佳模型---------------------------------------') print("在交叉验证中验证的最好结果:", gd.best_score_) # 打印交叉验证中的最佳得分(R2) print("最好的参数模型:", gd.best_estimator_) # 打印具有最佳参数的模型对象 print("(dict)最佳参数:", gd.best_params_) # 打印最佳参数组合的字典 print('-------------------------------------保存最佳模型---------------------------------------') model_save_dir = r'D:\新建文件夹' # 定义模型保存的目录 os.makedirs(model_save_dir, exist_ok=True) # 创建目录,如果目录已存在则不报错 model_path = os.path.join(model_save_dir, 'XGBoost_model_final.joblib') # 定义模型的完整保存路径 joblib.dump(gd.best_estimator_, model_path) # 将找到的最佳模型保存到文件 print(f"模型已保存至: {model_path}") # 打印保存成功信息 loaded_model = joblib.load(model_path) # 从文件加载模型 print('-------------------------------应用模型--------------------------------------') y_test_pred = loaded_model.predict(X_test_scaled_df) # 使用加载的模型对测试集进行预测 y_train_pred = loaded_model.predict(X_train_scaled_df) # 使用加载的模型对训练集进行预测 #保存数据 X_train_scaled_df.to_excel("xtr.xlsx") X_test_scaled_df.to_excel("xte.xlsx") pd.DataFrame([y_train,y_train_pred]).T.to_excel('ytr.xlsx') pd.DataFrame([y_train,y_test_pred]).T.to_excel('yte.xlsx') print('-------------------------------------训练模型性能---------------------------------------') train_mse = mean_squared_error(y_train, y_train_pred) # 计算训练集的均方误差(MSE) train_rmse = np.sqrt(train_mse) # 计算训练集的均方根误差(RMSE) train_mae = mean_absolute_error(y_train, y_train_pred) # 计算训练集的平均绝对误差(MAE) train_r2 = r2_score(y_train, y_train_pred) # 计算训练集的决定系数(R2) print(f'MSE: {train_mse:.4f}, RMSE: {train_rmse:.4f}, MAE: {train_mae:.4f}, R2: {train_r2:.4f}') print('-------------------------------------验证模型性能---------------------------------------') test_mse = mean_squared_error(y_test, y_test_pred) # 计算测试集的均方误差(MSE) test_rmse = np.sqrt(test_mse) # 计算测试集的均方根误差(RMSE) test_mae = mean_absolute_error(y_test, y_test_pred) # 计算测试集的平均绝对误差(MAE) test_r2 = r2_score(y_test, y_test_pred) # 计算测试集的决定系数(R2) print(f'MSE: {test_mse:.4f}, RMSE: {test_rmse:.4f}, MAE: {test_mae:.4f}, R2: {test_r2:.4f}') print('----------------------------------------结果绘图-----------------------------------------') results_plot_save_dir = r'D:\新建文件夹' # 定义结果图保存的目录 train_path = os.path.join(results_plot_save_dir, 'XGBoost_训练集精度_final.png') # 训练集拟合图的保存路径 test_path = os.path.join(results_plot_save_dir, 'XGBoost_验证集精度_final.png') # 验证集拟合图的保存路径 # 调用函数绘制训练集的拟合图 plot_regression_fit(y_train, y_train_pred, train_r2, train_rmse, train_mae, 'Train Set', 'XGBoost Model Performance (Train Set)', train_path) # 调用函数绘制测试集的拟合图 plot_regression_fit(y_test, y_test_pred, test_r2, test_rmse, test_mae, 'Test Set', 'XGBoost Model Performance (Test Set)', test_path) plt.rcdefaults() # 恢复matplotlib的默认设置 print( '----------------------------------------计算并绘制XGBoost原生特征重要性图-----------------------------------------') importances = loaded_model.feature_importances_ # 获取XGBoost模型内置的特征重要性(基于分裂增益) # 创建一个包含特征名称和重要性分数的DataFrame gbdt_importance_df = pd.DataFrame({'Feature': feature_names_from_df, 'Importance': importances}) save_path_gbdt = os.path.join(results_plot_save_dir, 'XGBoost_特征重要性组合图_final.png') # 定义保存路径 # 调用函数绘制组合特征重要性图 plot_importance_combined(gbdt_importance_df, 'XGBoost模型计算的特征重要性', save_path_gbdt, bar_color='dodgerblue') print( "----------------------------------------计算并绘制Permutation Importance图-----------------------------------------") scores = defaultdict(list) # 创建一个默认值为列表的字典,用于存储每个特征的置换重要性分数 # 遍历每一个特征 for feat_name in feature_names_from_df: X_t = X_test_scaled_df.copy() # 复制一份测试集数据 # 随机打乱当前特征列的顺序 X_t[feat_name] = np.random.permutation(X_t[feat_name].values) # 计算打乱后模型的R2分数 shuff_acc = r2_score(y_test, loaded_model.predict(X_t)) # 计算重要性:(原始R2 - 打乱后R2) / 原始R2,如果原始R2接近0则直接用差值 scores[feat_name].append((test_r2 - shuff_acc) / test_r2 if test_r2 > 1e-6 else test_r2 - shuff_acc) # 对特征按重要性得分从高到低排序 sorted_scores = sorted([(np.mean(score_list), feat) for feat, score_list in scores.items()], reverse=True) perm_feature_names = [feat for _, feat in sorted_scores] # 获取排序后的特征名称 perm_feature_scores = [score for score, _ in sorted_scores] # 获取排序后的重要性分数 # 创建一个包含置换重要性结果的DataFrame perm_importance_df = pd.DataFrame({'Feature': perm_feature_names, 'Importance': perm_feature_scores}) save_path_perm = os.path.join(results_plot_save_dir, 'XGBoost_特征重要性_Permutation_final.png') # 定义保存路径 # 调用函数绘制组合特征重要性图(使用置换重要性数据) plot_importance_combined(perm_importance_df, '特征重要性 (Permutation Importance for XGBoost)', save_path_perm, bar_color='lightcoral') print('----------------------------------------绘制残差分析图-----------------------------------------') train_residuals = y_train - y_train_pred # 计算训练集的残差 test_residuals = y_test - y_test_pred # 计算测试集的残差 train_res_path = os.path.join(results_plot_save_dir, 'XGBoost_训练集残差分析图_final.png') # 训练集残差图保存路径 test_res_path = os.path.join(results_plot_save_dir, 'XGBoost_验证集残差分析图_final.png') # 测试集残差图保存路径 # 调用函数绘制训练集残差图 plot_residuals_styled(train_residuals, y_train_pred, train_res_path, 'XGBoost 训练集残差分析') # 调用函数绘制测试集残差图 plot_residuals_styled(test_residuals, y_test_pred, test_res_path, 'XGBoost 验证集残差分析') # ================================================================================= # ============ 【全新】使用手动计算方法绘制 PDP 和 ICE 相关图 ============ # ================================================================================= print('------------------------开始 PDP 和 ICE 相关绘图 (手动实现)------------------------') # 定义PDP/ICE图的保存目录 pdp_ice_save_dir = os.path.join(results_plot_save_dir, 'XGBoost_PDP_ICE_Plots_final') os.makedirs(pdp_ice_save_dir, exist_ok=True) # 创建目录 # 定义双变量PDP图的保存目录 pdp_2way_save_dir = os.path.join(pdp_ice_save_dir, '2Way_PDP_All_Combinations') os.makedirs(pdp_2way_save_dir, exist_ok=True) # 创建目录 # 定义3D PDP图的保存目录 pdp_3d_save_dir = os.path.join(pdp_ice_save_dir, '3D_PDP_All_Combinations') os.makedirs(pdp_3d_save_dir, exist_ok=True) # 创建目录 n_top_features_for_pdp = 6 # 设置用于PDP分析的最重要特征的数量 if n_top_features_for_pdp > len(feature_names_from_df): # 如果特征总数不足,则取全部特征 n_top_features_for_pdp = len(feature_names_from_df) # 根据XGBoost原生重要性排序,选取最重要的N个特征 top_features_pdp_names = gbdt_importance_df['Feature'].tolist()[:n_top_features_for_pdp] plt.style.use('seaborn-v0_8-whitegrid') # 设置绘图风格 plt.rc("font", family='Microsoft YaHei') # 设置中文字体 # --- 1. 绘制单变量 PDP (含置信区间) 和 ICE 组合图 --- print("\n开始绘制单变量 PDP (含95%置信区间) 和 ICE 组合图...") # 遍历最重要的N个特征 for feature_name in top_features_pdp_names: print(f"正在计算特征 '{feature_name}' 的PDP/ICE数据...") try: # 使用手动编写的函数计算1D PDP和ICE数据 grid_vals, pdp_vals, ice_lines_vals = manual_pdp_1d(loaded_model, X_train_scaled_df, feature_name) # 在每个网格点上计算所有ICE线的标准差,用于构建置信区间 pdp_std = np.std(ice_lines_vals, axis=0) # 开始绘图 fig, ax = plt.subplots(figsize=(10, 8)) # 绘制所有样本的ICE线 (半透明蓝色细线) for ice_line in ice_lines_vals: ax.plot(grid_vals, ice_line, color='tab:blue', alpha=0.05, linewidth=0.5) # 绘制PDP线 (红色虚线),代表平均效应 ax.plot(grid_vals, pdp_vals, color='red', linestyle='--', linewidth=3, label='平均效应 (PDP)') # 绘制95%置信区间 (平均值 ± 1.96 * 标准差) ax.fill_between( grid_vals, pdp_vals - 1.96 * pdp_std, pdp_vals + 1.96 * pdp_std, color='skyblue', alpha=0.4, label='95% 置信区间' ) ax.set_title(f'PDP/ICE 组合图\n特征: {feature_name}', fontsize=16) # 设置标题 ax.set_xlabel(f'{feature_name} (标准化值)', fontsize=12) # 设置x轴标签 ax.set_ylabel('对预测值的依赖性', fontsize=12) # 设置y轴标签 ax.legend() # 显示图例 # 保存图表 plt.savefig(os.path.join(pdp_ice_save_dir, f'XGBoost_Manual_PDP_ICE_{feature_name}.png'), dpi=300, bbox_inches='tight') plt.close(fig) # 关闭图表 print(f"成功绘制特征 '{feature_name}' 的PDP/ICE图。") except Exception as e: print(f"绘制手动 PDP/ICE for {feature_name} 出错: {e}") # 打印错误信息 # --- 2. 绘制双变量 (2D 和 3D) PDP 图 --- print("\n开始绘制双变量 PDP (2D 热力图 和 3D 曲面图)...") if len(top_features_pdp_names) >= 2: # 确保至少有两个特征可以进行组合 # 遍历最重要的N个特征中的所有两两组合 for feat1, feat2 in combinations(top_features_pdp_names, 2): print(f"正在计算特征对 '{feat1}' vs '{feat2}' 的2D PDP数据...") try: # 使用手动编写的函数计算2D PDP数据 grid_x, grid_y, pdp_z = manual_pdp_2d(loaded_model, X_train_scaled_df, (feat1, feat2)) # 创建用于绘图的网格坐标 XX, YY = np.meshgrid(grid_x, grid_y) # 注意:pdp_z的维度可能需要转置以匹配meshgrid的坐标系 ZZ = pdp_z.T # 绘制 2D 热力图 fig_2d, ax_2d = plt.subplots(figsize=(8, 7)) # 使用contourf填充等值线图 c = ax_2d.contourf(XX, YY, ZZ, cmap='viridis', levels=20) fig_2d.colorbar(c, ax=ax_2d, label='部分依赖值') # 添加颜色条 ax_2d.set_title(f'2D PDP: {feat1} vs {feat2}', fontsize=16) # 设置标题 ax_2d.set_xlabel(f'{feat1} (标准化值)', fontsize=12) # 设置x轴标签 ax_2d.set_ylabel(f'{feat2} (标准化值)', fontsize=12) # 设置y轴标签 plt.savefig(os.path.join(pdp_2way_save_dir, f'XGBoost_Manual_PDP_2D_{feat1}_{feat2}.png'), dpi=300) # 保存 plt.close(fig_2d) # 关闭图表 # 绘制 3D 曲面图 fig_3d = plt.figure(figsize=(12, 9)) ax_3d = fig_3d.add_subplot(111, projection='3d') # 创建3D子图 # 绘制3D曲面 surf = ax_3d.plot_surface(XX, YY, ZZ, cmap='viridis', edgecolor='none', antialiased=True) fig_3d.colorbar(surf, shrink=0.5, aspect=20, label='部分依赖值', pad=0.1) # 添加颜色条 ax_3d.set_xlabel(f'{feat1} (标准化值)', fontsize=10, labelpad=10) # x轴标签 ax_3d.set_ylabel(f'{feat2} (标准化值)', fontsize=10, labelpad=10) # y轴标签 ax_3d.set_zlabel('对预测值的依赖性 (PDP)', fontsize=10, labelpad=10, rotation=180) # z轴标签 ax_3d.set_title(f'三维部分依赖图 (3D PDP)\n{feat1} vs {feat2}', fontsize=14) # 标题 ax_3d.view_init(elev=20, azim=45) # 设置视角 plt.savefig(os.path.join(pdp_3d_save_dir, f'XGBoost_Manual_PDP_3D_{feat1}_{feat2}.png'), dpi=300) # 保存 plt.close(fig_3d) # 关闭图表 print(f"成功绘制特征对 '{feat1}' vs '{feat2}' 的2D和3D PDP图。") except Exception as e: print(f"绘制手动 2D/3D PDP for {feat1} & {feat2} 出错: {e}") # 打印错误信息 # --- 绘制三特征3D散点图 --- print("\n开始绘制三特征 (3D) 散点图...") pdp_3d_scatter_save_dir = os.path.join(pdp_ice_save_dir, '3D_Scatter_Three_Features') # 定义保存目录 os.makedirs(pdp_3d_scatter_save_dir, exist_ok=True) # 创建目录 if len(top_features_pdp_names) >= 3: # 确保至少有3个特征 # 最多选择前4个重要特征进行组合,避免组合数过多 n_features_for_3d_scatter = min(len(top_features_pdp_names), 4) # 遍历所有三个特征的组合 for features_tuple in combinations(top_features_pdp_names[:n_features_for_3d_scatter], 3): try: # 定义保存路径 save_path = os.path.join(pdp_3d_scatter_save_dir, f'XGBoost_3D_Scatter_{features_tuple[0]}_{features_tuple[1]}_{features_tuple[2]}.png') # 调用函数绘制3D散点图 plot_3d_scatter_three_features(X_test_scaled_df, y_test_pred, features_tuple, save_path) except Exception as e: print(f"绘制 3D 散点图 for {features_tuple} 出错: {e}") # 打印错误信息 # --- 调用:绘制固定特征值的3D PDP图 --- print("\n开始绘制固定特征值的3D PDP图...") pdp_3d_fixed_save_dir = os.path.join(pdp_ice_save_dir, '3D_PDP_Fixed_Value') # 定义保存目录 os.makedirs(pdp_3d_fixed_save_dir, exist_ok=True) # 创建目录 if len(top_features_pdp_names) >= 3: # 确保至少有3个特征 # 最多选择前4个重要特征进行组合 n_features_for_3d_fixed = min(len(top_features_pdp_names), 4) # 遍历所有三个特征的组合 for features_tuple in combinations(top_features_pdp_names[:n_features_for_3d_fixed], 3): # 对每个组合,轮流固定其中的一个特征 for feature_to_fix in features_tuple: try: features_list = list(features_tuple) # 元组转列表 # 获取另外两个变化的特征 varying_feats = [f for f in features_list if f != feature_to_fix] # 定义保存路径 save_path = os.path.join(pdp_3d_fixed_save_dir, f'XGBoost_3DPDP_{varying_feats[0]}_{varying_feats[1]}_Fix_{feature_to_fix}.png') # 将固定的值设为该特征的中位数 fixed_val = X_train_scaled_df[feature_to_fix].median() # 调用函数绘制固定特征值的3D PDP图 plot_3d_pdp_fixed_value( loaded_model, X_train_scaled_df, features_list, save_path, fixed_feature=feature_to_fix, fixed_value=fixed_val ) except Exception as e: print(f"绘制固定值3D PDP for {features_list} (固定 {feature_to_fix}) 出错: {e}") # 打印错误信息 print('------------------------开始 SHAP 分析------------------------') shap_save_dir = os.path.join(results_plot_save_dir, 'XGBoost_SHAP_Plots_final') # 定义SHAP图的保存目录 os.makedirs(shap_save_dir, exist_ok=True) # 创建目录 explainer = shap.TreeExplainer(loaded_model) # 为树模型创建一个SHAP解释器 shap_values = explainer(X_test_scaled_df) # 计算测试集所有样本的SHAP值 print("\n绘制 SHAP Summary Plot (条形图)...") # 计算每个特征的平均绝对SHAP值,作为其重要性 shap_importance_vals = np.abs(shap_values.values).mean(axis=0) # 创建包含SHAP重要性的DataFrame shap_importance_df = pd.DataFrame({'Feature': X_test_scaled_df.columns, 'Importance': shap_importance_vals}) save_path_shap = os.path.join(shap_save_dir, 'XGBoost_SHAP_特征重要性组合图_final.png') # 定义保存路径 # 调用组合重要性绘图函数,绘制SHAP重要性条形图 plot_importance_combined(shap_importance_df, 'SHAP 特征重要性 (平均绝对SHAP值)', save_path_shap, bar_color='#007bff') print("绘制 SHAP Summary Plot (散点分布图)...") shap.summary_plot(shap_values, X_test_scaled_df, show=False) # 生成SHAP摘要图(散点形式) plt.title('SHAP 特征影响概览 (散点分布)', fontsize=16) # 添加标题 plt.tight_layout() # 调整布局 plt.savefig(os.path.join(shap_save_dir, 'XGBoost_SHAP_summary_scatter.png'), dpi=300, bbox_inches='tight') # 保存 plt.close() # 关闭图表 print("绘制 SHAP Dependence Plots...") shap_dependence_save_dir = os.path.join(shap_save_dir, 'Dependence_Plots') # 定义SHAP依赖图的保存目录 os.makedirs(shap_dependence_save_dir, exist_ok=True) # 创建目录 # 为最重要的N个特征绘制SHAP依赖图 for feature_name in top_features_pdp_names: # 绘制单个特征的依赖图,图中颜色表示交互效应最强的另一个特征 shap.dependence_plot(feature_name, shap_values.values, X_test_scaled_df, interaction_index="auto", show=False) plt.gcf().suptitle(f'SHAP 依赖图: {feature_name}', fontsize=16) # 添加总标题 plt.tight_layout() # 调整布局 plt.savefig(os.path.join(shap_dependence_save_dir, f'XGBoost_SHAP_dependence_{feature_name}.png'), dpi=300, bbox_inches='tight') # 保存图表 plt.close() # 关闭图表 print("绘制 SHAP Waterfall Plot (针对测试集第一个样本)...") plt.figure() # 创建一个新的画布 # 绘制瀑布图,展示单个预测(这里是测试集第一个样本)的SHAP值构成 shap.plots.waterfall(shap_values[0], max_display=15, show=False) plt.title('SHAP Waterfall Plot (测试集样本 0)', fontsize=16) # 添加标题 plt.tight_layout() # 调整布局 plt.savefig(os.path.join(shap_save_dir, 'XGBoost_SHAP_waterfall_sample_0.png'), dpi=300, bbox_inches='tight') # 保存 plt.close() # 关闭图表 print('----------------------------------------SHAP 分析完成-----------------------------------------') print('----------------------------------------开始 ALE 分析-----------------------------------------') ale_save_dir = os.path.join(results_plot_save_dir, 'XGBoost_ALE_Plots_final') # 定义ALE图的保存目录 os.makedirs(ale_save_dir, exist_ok=True) # 创建目录 print(f"ALE 相关图将保存到: {ale_save_dir}") top_features_ale_names = top_features_pdp_names # 使用与PDP相同的最重要特征列表 print(f"\n开始为最重要的 {len(top_features_ale_names)} 个特征绘制一维 ALE 图...") # 生成一组颜色用于绘制不同的ALE图 colors = plt.cm.viridis(np.linspace(0, 0.85, len(top_features_ale_names))) # 遍历最重要的特征 for i, feature_name in enumerate(top_features_ale_names): try: # 使用PyALE库计算并绘制一维ALE图 ale_eff = ale(X=X_train_scaled_df, model=loaded_model, feature=[feature_name], feature_type='continuous', grid_size=50, include_CI=True, C=0.95) fig, ax = plt.gcf(), plt.gca() # 获取当前的图和坐标轴 current_color = colors[i] # 为当前特征选择一个颜色 if ax.lines: # 如果图中有线(ALE主线) ax.lines[0].set_color(current_color) # 设置线的颜色 ax.lines[0].set_linewidth(2.5) # 设置线的宽度 if ax.collections: # 如果图中有集合(置信区间) ax.collections[0].set_facecolor(current_color) # 设置填充颜色 ax.collections[0].set_alpha(0.2) # 设置透明度 ax.set_title(f'累积局部效应 (ALE) - 特征: {feature_name}', fontsize=16) # 设置标题 ax.set_xlabel(f'{feature_name} (标准化值)', fontsize=12) # 设置x轴标签 ax.set_ylabel('ALE (对预测值的影响)', fontsize=12) # 设置y轴标签 plt.tight_layout() # 调整布局 plt.savefig(os.path.join(ale_save_dir, f'XGBoost_ALE_1D_{feature_name}.png'), dpi=300, bbox_inches='tight') # 保存 plt.close(fig) # 关闭图表 except Exception as e: print(f"绘制 1D ALE for {feature_name} 出错: {e}") # 打印错误信息 if plt.get_fignums(): # 如果有未关闭的图表 plt.close('all') # 全部关闭 print(f"\n开始为最重要的特征对绘制二维 ALE 图...") if len(top_features_ale_names) >= 2: # 确保至少有两个特征 # 遍历所有两两特征组合 for feat1_name, feat2_name in combinations(top_features_ale_names, 2): try: # 计算二维ALE效应,但不立即绘图 (plot=False) ale_eff_2d = ale(X=X_train_scaled_df, model=loaded_model, feature=[feat1_name, feat2_name], feature_type='continuous', grid_size=30, plot=False) fig, ax = plt.subplots(figsize=(8, 7)) # 创建画布 # 使用pcolormesh绘制二维ALE热力图 im = ax.pcolormesh(ale_eff_2d.index, ale_eff_2d.columns, ale_eff_2d.values.T, cmap='viridis', shading='auto') fig.colorbar(im, ax=ax, label='ALE (对预测值的影响)') # 添加颜色条 ax.set_title(f'二维 ALE: {feat1_name} vs {feat2_name}', fontsize=16) # 设置标题 ax.set_xlabel(f'{feat1_name} (标准化值)', fontsize=12) # 设置x轴标签 ax.set_ylabel(f'{feat2_name} (标准化值)', fontsize=12) # 设置y轴标签 plt.tight_layout() # 调整布局 plt.savefig(os.path.join(ale_save_dir, f'XGBoost_ALE_2D_{feat1_name}_vs_{feat2_name}.png'), dpi=300, bbox_inches='tight') # 保存图表 plt.close(fig) # 关闭图表 except Exception as e: print(f"绘制 2D ALE for {feat1_name} & {feat2_name} 出错: {e}") # 打印错误信息 if plt.get_fignums(): # 如果有未关闭的图表 plt.close('all') # 全部关闭 print('----------------------------------------ALE 分析完成-----------------------------------------') print('----------------------------------------脚本执行完毕-----------------------------------------') 以上是我对巢湖流域西半湖湖心的国控监测站的2022/4/1 0:00:00 - 2022/7/26 16:00:00时间段的数据做的建模,根据我的代码、数据、结果图、写一篇中文论文,要求可以直接投稿期刊的结果
09-06
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值