期刊图片复现|Python实现非线性极值约束分析及阈值响应特征分析

# =========================================================================================# ====================================== 1. 环境设置 =======================================# =========================================================================================import osimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom matplotlib.colors import LinearSegmentedColormap

# =========================================================================================# ======================================2.颜色库=======================================# =========================================================================================COLOR_SCHEMES = {    1: ['#0000CD', '#FF0000', '#2E8B57', '#9400D3'],}

# =========================================================================================# ======================================3.曲线拟合函数=======================================# =========================================================================================def fit_best_curve(x, y, x_smooth):    models = {}  # 用于存储模型    n = len(x)  # 数据点数量    # 线性回归    p1 = np.polyfit(x, y, 1)    models['Linear'] = (np.poly1d(p1)(x), np.poly1d(p1)(x_smooth), 1)    # 二次多项式    p2 = np.polyfit(x, y, 2)    models['Poly2'] = (np.poly1d(p2)(x), np.poly1d(p2)(x_smooth), 2)    # 三次多项式    p3 = np.polyfit(x, y, 3)    models['Poly3'] = (np.poly1d(p3)(x), np.poly1d(p3)(x_smooth), 3)    has_non_positive = np.min(x) <= 0  # 检查是否存在小于等于0的值    best_r2 = -float('inf')  # 初始化最高R2    best_y_smooth = None  # 初始化最佳平滑y值    best_p_val = 1.0  # 初始化最佳P值    y_mean = np.mean(y)  # y的均值    ss_tot = np.sum((y - y_mean) ** 2) + 1e-8  # 计算总平方和，加极小值防分母为零        # 如果当前模型R2更好        if r2 > best_r2:            best_r2 = r2  # 更新最佳R2            best_y_smooth = y_pred_smooth  # 更新对应的预测曲线Y值            best_p_val = p_val  # 更新对应P值    return best_y_smooth, best_r2, best_p_val  # 返回平滑预测值、R2和P值

# =========================================================================================# ======================================4.寻找曲线拐点函数=======================================# =========================================================================================def find_inflection_threshold(x, y):    y_double_prime = np.gradient(y_prime) / (dx + 1e-8)  # 计算二阶导数    curvature = np.abs(y_double_prime) / (1 + y_prime ** 2) ** 1.5  # 计算曲率    margin = int(len(x) * 0.1)  # 忽略前后10%的数据以避免边界效应    valid_curvature = curvature[margin:-margin]  # 提取中间部分的曲率    max_idx = margin + np.argmax(valid_curvature)  # 寻找最大曲率位置对应原始数据的索引    return max_idx  # 返回拐点索引

# =========================================================================================# ======================================6.散点图绘制函数=======================================# =========================================================================================def draw_scatter(ax, df, feature_col, xlabel_text, analysis_res, scheme_id):    scatter_pt_c, scatter_line_c, _, _ = COLOR_SCHEMES[scheme_id]  # 获取散点及拟合线颜色    # 绘制散点    r2 = analysis_res['r2']  # 获取分析结果中的R2    p_val = analysis_res['p_value']  # P值    # 设置显著性标记    if p_val < 0.001:        p_text = r'p < 0.001'    elif p_val < 0.01:        p_text = r'p < 0.01'    elif p_val < 0.05:        p_text = r'p < 0.05'    else:        p_text = f'p={p_val:.3f}'    # 添加文本标注    ax.text(0.65,  # x            0.88,  # y            f'$R^2={r2:.4f}$',  # 文本            fontsize=22,  # 字体大小            transform=ax.transAxes)  # 坐标系    ax.text(0.65,  # x            0.78,  # y            f'${p_text}$',  # 文本            fontsize=22,  # 大小            transform=ax.transAxes)  # 坐标系    apply_axis_styles(ax, df, feature_col, xlabel_text)  # 设置轴和刻度

# =========================================================================================# ======================================7.热图绘制函数=======================================# =========================================================================================def draw_heatmap(fig, ax, df, feature_col, xlabel_text, scheme_id):    scatter_pt_c, scatter_line_c, _, _ = COLOR_SCHEMES[scheme_id]  # 提取配色    # 创建渐变色    custom_cmap = LinearSegmentedColormap.from_list(f'cmap_{scheme_id}',                                                    ['#101030', scatter_pt_c, '#f5f5f5', scatter_line_c])    if len(x_data) > 10000:        idx = np.random.choice(len(x_data), 10000, replace=False)  # 无放回抽样        x_sample = x_data[idx]  # x轴抽样数据        y_sample = y_data[idx]  # y轴抽样数据    else:        x_sample = x_data  # x        y_sample = y_data  # y    # 绘制颜色条    cbar = fig.colorbar(h, cax=cax, format='%.4f')    # 颜色条顶部标题    cbar.ax.set_title('Density', fontsize=16, pad=15)    cbar.ax.tick_params(labelsize=14)  # 颜色条刻度标注    apply_axis_styles(ax, df, feature_col, xlabel_text)  # 设置轴和刻度

# =========================================================================================# ======================================8.RCS绘制函数======================================# =========================================================================================def draw_rcs(ax, df, feature_col, xlabel_text, analysis_res, scheme_id):    _, _, rcs_fill_c, rcs_pt_c = COLOR_SCHEMES[scheme_id]  # 提取配色    x_smooth = analysis_res['x_smooth']  # 平滑x    y_bottom = 0.0  # 曲线下面积的底部参考线    im.set_clip_path(patch)  # 把之前创建的曲线路径作为掩膜实现曲线下渐变填充    # 绘制RCS主平滑线    ax.plot(x_smooth,  # x            y_rcs,  # y            color='#2C2C2C',  # 颜色            linestyle='--',  # 虚线            linewidth=1.5)  # 粗细    x_min_data, x_max_data = df[feature_col].min(), df[feature_col].max()  # 真实数据x极值    x_pad = (x_max_data - x_min_data) * 0.05  # 两侧留白    y_pad = (1.0 - 0.0) * 0.05  # y两侧六百    # 文本    ax.text(threshold_x,  # x            threshold_y - 0.1,  # y            f'({threshold_x:.2f},'  # 文本x            f' {threshold_y:.4f})',  # 文本y            fontsize=16,  # 大小            ha='center',  # 水平            va='top')  # 垂直    apply_axis_styles(ax, df, feature_col, xlabel_text)  # 设置轴和刻度

# =========================================================================================# ======================================9.执行部分=======================================# =========================================================================================if __name__ == '__main__':    all_analysis_results = {}  # 存储各特征分析结果    # 遍历特征    for col, xlabel in feature_configs:        num_bins = 100  # 分箱数        x_data = df_real[col]  # 取该特征数据        bins = np.linspace(x_data.min(), x_data.max(), num_bins + 1)  # 生成等间距的分箱区间        bin_centers = []  # 存储分箱中点的列表        boundary_y = []  # 存储边界Y值(98分位数)的列表        # 遍历分箱        for i in range(num_bins):            mask = (x_data >= bins[i]) & (x_data < bins[i + 1])  # 生成当前分箱的数据掩码            y_in_bin = df_real['LERI'][mask]  # 提取落入当前分箱的目标值            if len(y_in_bin) > 10:                bin_centers.append(0.5 * (bins[i] + bins[i + 1]))  # 计算当前分箱中心坐标                boundary_y.append(np.percentile(y_in_bin, 98))  # 计算箱内数据的第98百分位数作为边界    # 是否批量绘图    plot_all = True    schemes_to_run = list(COLOR_SCHEMES.keys()) if plot_all else [1]  # 获取配色    # 遍历配色进行绘图    for scheme_id in schemes_to_run:        print(f'正在绘制并保存方案：[{scheme_id}]')        fig_combo.tight_layout(w_pad=0.0, h_pad=2)  # 组图布局调整        # 保存组图        fig_combo.savefig(os.path.join(BASE_DIR, f'Combined_AllFeatures_scheme_{scheme_id}.png'), dpi=300,                          bbox_inches='tight')        plt.close(fig_combo)  # 关闭

df_real = pd.read_excel(r'data.xlsx') 

2.设置绘图结果的保存路径，执行部分：

BASE_DIR = r'分析'  # 保存路径DIR_SCATTER = os.path.join(BASE_DIR, 'Scatter')  # 散点图保存路径DIR_HEATMAP = os.path.join(BASE_DIR, 'Heatmap')  # 热力图保存路径DIR_RCS = os.path.join(BASE_DIR, 'RCS')  # RCS图保存路径

feature_configs = [    ('Elevation', 'Elevation/m'),    ('Temperature', 'Temperature/°C'),    ('Precipitation', 'Precipitation/mm'),    ('WindSpeed', 'Wind speed(m/s)')]

4.设置是否要进行批量绘图，执行部分：

plot_all = True