期刊图片复现|Python进行约束线分析

方法原理

这份Python代码实现了一套完整且高度自动化的双约束线分析方法，通常应用于生态学、环境科学或地理空间分析领域。核心目的在于探究多个环境因子（如温度、降水、坡度、NDVI等）对某个目标响应变量的非线性限制效应。与探讨“平均趋势”的常规回归分析不同，约束线分析关注的是数据分布的“轮廓边缘”，用以揭示当其他潜在环境条件处于最适宜状态时，单一自变量对因变量的最大承载潜力（上限）和最小限制底线（下限）。为了清晰且极具深度地展示其中复杂的数学与统计原理，以下从数据降噪、模型优选、极值定位以及不确定性评估等多个核心维度进行剖析：

1. 数据分箱与分位数特征提取

从原始数据中提取出真正能够代表上下边界的数据点。根据输入自变量的极小值和极大值，将其定义域均匀切割成60个等距区间。同时对于每个区间设定了15的门槛。只有当某个分箱区间内的数据点数量大于或等于15个时，该区间才被认为具有统计学上的代表性，从而有效避免了数据稀疏区域可能导致的伪边界现象。对于满足条件的有效区间，代码计算其X值的中位数作为该区间的代表横坐标。在纵坐标方向上，采用95%分位数作为上限，5%分位数作为下限。

2. 自适应非线性多项式拟合

获取了上下边界后，采用Levenberg-Marquardt（LM）非线性优化算法进行多项式回归拟合。尝试二次多项式和三次多项式的双重拟合测试。在评估拟合效果时，算法引入了校正决定系数（Adjusted R2）来进行模型评估。只有当三次拟合的校正R2严格大于二次拟合的校正R2，并且该值大于0时，才会采纳三次多项式曲线；否则选用二次曲线，在极端数据下退化为常数模型。

3. 极值点求解与阈值定位

拟合出平滑的边界曲线后，通过多项式对象的deriv()方法获取拟合曲线的一阶导函数，并利用roots属性求解导数为0的根。算法找到的上限极大值点代表了在给定环境因子下，系统理论上能达到的最高产出量；而下限极小值点则代表了系统运行的最脆弱点或环境胁迫的绝对底线。

4. Bootstrap自助法重抽样与置信区间构建

利用Bootstrap技术，从原始数据点中进行有放回的随机抽样，生成大小相同但包含随机重复与缺失项的全新观测集。针对每一次生成的虚拟样本，都会完整复现一遍“等距分箱中位数与分位数提取自适应多项式优选”的整个流程。循环执行1000次。提取这1000次成功模拟所产生的Y值阵列，并截取每个X节点上的2.5%和97.5%分位数，构建出95%置信区间。

5. 基准线交点计算与可视化

设定一条基准线，用于预警因变量何时会跌破安全阈值。将下限拟合多项式整体向下平移，随后寻找这个新多项式等于0的实数根。计算出了为了维持系统底线安全，环境自变量所必须跨过或保持的临界红线坐标。生成多子图画布。最底层是灰色原始散点云以呈现密度，浅层是半透明的上下置信带，中层是提取出的代表性散点与拟合边界线，最顶层则是极值点坐标、阈值垂直虚线以及基准刻度交点。同时，内置60套配色方案，实现了满足顶级SCI学术期刊出版要求的高标准数据可视化。

# =========================================================================================# ====================================== 1. 环境设置 =======================================# =========================================================================================import numpy as npimport matplotlib.pyplot as pltfrom scipy.optimize import curve_fitimport pandas as pdimport stringimport matplotlibfrom scipy import statsfrom sklearn.metrics import mean_squared_errorfrom sklearn.metrics import r2_scoreplt.rcParams['font.family'] = 'serif'plt.rcParams['font.serif'] = ['Times New Roman']plt.rcParams['axes.unicode_minus'] = Falsematplotlib.rcParams['pdf.fonttype'] = 42matplotlib.rcParams['ps.fonttype'] = 42plt.rcParams['mathtext.fontset'] = 'stix'

# =========================================================================================# ======================================2.颜色库=======================================# =========================================================================================COLOR_SCHEMES = {    1: {'UPPER': '#E53935', 'LOWER': '#1E88E5', 'GREY': '#B0BEC5'},}

# =========================================================================================# ======================================3.多项式函数=======================================# =========================================================================================def poly2(x, a, b, c):    return a * x ** 2 + b * x + cdef poly3(x, a, b, c, d):    return a * x ** 3 + b * x ** 2 + c * x + d

# =========================================================================================# ======================================4.将系数格式化函数=======================================# =========================================================================================def format_formula(coeffs):    deg = len(coeffs) - 1  # 根据系数的个数计算多项式的最高次幂    terms = []  # 用于存储公式中的每一项文本    # 遍历系数列表，获取索引和系数数值    for i, c in enumerate(coeffs):        elif power == 3:  # 三次项            terms.append(f"{c_str}x³")  # 添加到列表中    return "y=" + "+".join(terms).replace("+-", "-")  # 返回公式

# =========================================================================================# ======================================5.多项式拟合函数及选择函数=======================================# =========================================================================================def lm_fit_and_select(x_data, y_data):    n = len(x_data)#样本数    #避免点数过少无法计算自由度    if n <= 4:        poly_func = np.poly1d([np.mean(y_data)])        return poly_func, "y = const", 0.0, 0    if adj_r2_3 > adj_r2_2 and adj_r2_3 > 0:        poly_func = np.poly1d(popt3)  # 使用三次拟合的参数生成多项式对象        return poly_func, format_formula(popt3), r2_3, 3  # 返回多项式对象、公式文本、原始R2展示、阶数3    elif popt2 is not None:        poly_func = np.poly1d(popt2)  # 使用二次拟合的参数生成多项式对象        return poly_func, format_formula(popt2), max(r2_2, 0), 2  # 返回多项式对象、公式文本、原始R2展示、阶数2    else:  # 如果失败        poly_func = np.poly1d([np.mean(y_data)])  # 生成一个只包y数据平均值的常数多项式        return poly_func, "y = const", 0.0, 0  # 返回常数多项式、占位文本、R2、阶数0

# =========================================================================================# ===============================6.散点分箱、提取分位数并拟合上下边界的函数==================# =========================================================================================def process_and_fit_quantiles(x, y, bins=60):    bin_edges = np.linspace(np.min(x), np.max(x), bins + 1)  # 根据给定的区间数生成x的等距分割边界    valid_x, q_up_list, q_low_list = [], [], []  # 用来保存有效的x值、95%分位数和5%分位数    # 遍历每一个分箱区间    for i in range(bins):    fit_up, form_up, r2_up, deg_up = lm_fit_and_select(valid_x, q_up_list)  # 多项式拟合挑选    fit_low, form_low, r2_low, deg_low = lm_fit_and_select(valid_x, q_low_list)  # 多项式拟合挑选    # 返回提取的点阵以及上下两条拟合曲线的所有信息，包含新增的阶数    return valid_x, q_up_list, q_low_list, fit_up, form_up, r2_up, deg_up, fit_low, form_low, r2_low, deg_low

# =========================================================================================# ======================================7.寻找多项式在特定区间内极大或极小值点的函数=======================================# =========================================================================================def find_extreme_point(poly_func, x_min, x_max, find_max=True):    deriv = poly_func.deriv()  # 对传入的多项式求一阶导数    roots = deriv.roots  # 解导数为0的点，即驻点    y_vals = poly_func(valid_roots)  # 计算这些内部有效极值根对应的y值    idx = np.argmax(y_vals) if find_max else np.argmin(y_vals)  # 如果要求找极大值，取y最大值的索引；否则取最小值的索引    return valid_roots[idx], y_vals[idx]  # 返回找到的目标极值点的横坐标和纵坐标

# =========================================================================================# ===================== 8.Bootstrap自助法置信区间计算函数 =====================# =========================================================================================def calculate_bootstrap_ci(raw_x, raw_y, x_smooth, is_upper=True, bins=60, n_bootstraps=1000):    n_points = len(raw_x)  #原始数据总数        # 生成一个与原始数据总数相同的随机索引数组，replace=True表示有放回的随机抽样        indices = np.random.choice(n_points,size=n_points,replace=True)        x_boot = raw_x[indices]  #本次重抽样的X数据集        y_boot = raw_y[indices]  #本次重抽样的Y数据集        #根据当前重抽样出的X数据集的极值，生成等间距的分箱边界        bin_edges = np.linspace(np.min(x_boot), np.max(x_boot), bins + 1)        valid_x, q_list = [], []  #用于存储有效分箱的中心X坐标和对应提取的分位数Y值        valid_x = np.array(valid_x)        q_list = np.array(q_list)    #判断是否至少有一次成功的拟合    if valid_sims > 0:        #删除空行        simulated_y = simulated_y[:valid_sims]        #置信区间下界        ci_lower = np.percentile(simulated_y, 2.5, axis=0)        #置信区间上界        ci_upper = np.percentile(simulated_y, 97.5, axis=0)    else:        #没有任何一次重抽样能成功拟合出曲线。生成nan数组        ci_lower = np.full_like(x_smooth, np.nan)        ci_upper = np.full_like(x_smooth, np.nan)    #返回置信区间下界和上界数组    return ci_lower, ci_upper

# =========================================================================================# ======================================9.双约束曲线图绘制函数=======================================# =========================================================================================def create_dual_constraint_plot(df, scheme_id=1):    scheme = COLOR_SCHEMES.get(scheme_id, COLOR_SCHEMES[1])  # 提取配色方案    num_features = len(feature_cols)  # 获取特征总数量    num_cols = 3  # 组合图列    num_rows = (num_features + num_cols - 1) // num_cols  # 行    fig, axes = plt.subplots(num_rows,  # 行                             num_cols,  # 列                             figsize=(16, 5.5 * num_rows))  # 大画布长宽    # 若子图多于1个将二维轴矩阵展平为一维方便遍历，否则套成列表    axes = np.array(axes).flatten() if num_features > 1 else [axes]    # 调整布局    plt.subplots_adjust(left=0.08,  # 左侧留白边缘                        right=0.98,  # 右侧留白边缘                        bottom=0.12,  # 底部留白边缘                        top=0.92,  # 顶部留白边缘                        wspace=0.25,  # 列间距                        hspace=0.35)  # 行间距        #接收函数返回值，获取阶数        val_x, q_up, q_low, fit_up, form_up, r2_up, deg_up, fit_low, form_low, r2_low, deg_low = process_and_fit_quantiles( feature_X, target_Y)        # 绘制背景散点        ax.scatter(feature_X,  # X                   target_Y,  # Y                   color=COLOR_GREY,  # 颜色                   alpha=0.15,  # 透明度                   s=8,  # 大小        # 绘制上限拟合线        ax.plot(x_smooth,  # x                y_smooth_up,  # y                color=COLOR_UPPER,  # 颜色                linewidth=2,  # 线宽                zorder=4)  # 层        # 绘制下限拟合线        ax.plot(x_smooth,  # x                y_smooth_low,  # y                color=COLOR_LOWER,  # 颜色                linewidth=2,  # 线宽                zorder=4)  # 层        #绘制上限置信区间        ax.fill_between(x_smooth,  # 拟合曲线                        ci_up_lower,  #下界                        ci_up_upper,  #上界                        color=COLOR_UPPER,  # 颜色                        alpha=0.15,  # 透明度                        zorder=2)  # 层        #绘制下限置信区间        ax.fill_between(x_smooth,# 拟合曲线                        ci_low_lower,  #下界                        ci_low_upper,  #上界                        color=COLOR_LOWER,  # 颜色                        alpha=0.15,  # 透明度                        zorder=2)  # 层                   color='black',  # 颜色                   linewidth=1.5,  # 线宽                   zorder=2)  # 层        # 顶部标题        ax.set_title(title,  # 文本                     fontsize=14,  # 大小                     fontweight='bold')  # 加粗        # X轴标题        ax.set_xlabel(xlabel, fontsize=12)        # Y轴标题

# =========================================================================================# =====================================10.执行部分======================================# =========================================================================================if __name__ == '__main__':    excel_path = r'data.xlsx'  # 原始数据路径    df = pd.read_excel(excel_path)  # 读取数据    # 用于决定是批量生成所有配色的图表，还是只生成单张特定配色的图表    plot_all = True    if plot_all:        for i in range(1, 61):            create_dual_constraint_plot(df, scheme_id=i)    else:        target_scheme = 22        create_dual_constraint_plot(df, scheme_id=target_scheme)

plot_all = True

2.单次绘图所使用的配色方案，最底下：

target_scheme = 22

excel_path = r'data.xlsx'  # 原始数据路径

4.设置分箱的数，第六部分的函数：

def process_and_fit_quantiles(x, y, bins=60):

5.设置分箱的数，抽样次数，第八部分的函数：

def calculate_bootstrap_ci(raw_x, raw_y, x_smooth, is_upper=True, bins=60, n_bootstraps=1000):

6.修改基准线的值，第九部分的函数：

BASELINE_Y = 557.5  # 基准水平线Y的恒定值

7.定义特征名映射，第九部分的函数：

label_map = {    'TEM': 'Temperature (°C)',    'PRE': 'Precipitation (mm)',    'SLP': 'Slope (°)',    'NDVI': 'Normalized Difference Vegetation Index',    'POP': 'Population (People/km²)',    'GDP': 'Gross Domestic Product (RMB yuan/km²)'}

plt.savefig(fr'plot{scheme_id}.png', dpi=300, bbox_inches='tight')