scipy curve_fit返回alpha_下限值问题的解决方法问询

使用scipy.optimize.curve_fit拟合自定义的correction_function时，参数alpha_始终返回设定的下限值，调整边界、拟合方法、max_nfev等参数后仍无改善。相关代码如下：

import numpy as np
from scipy import optimize

def correction_function(z, A, z_c, alpha_):
    return (A * np.power(z, 2) * np.exp(-np.power((z/z_c), alpha_)))

x_data = [0.0025000000000000001, 0.0074999999999999997, 0.012500000000000001, 0.017499999999999998, 0.022499999999999999, 0.0275, 0.032500000000000001, 0.037500000000000006, 0.042500000000000003, 0.047500000000000001, 0.052500000000000005, 0.057500000000000002, 0.0625, 0.067500000000000004, 0.072500000000000009, 0.077499999999999999, 0.082500000000000004, 0.087500000000000008, 0.092499999999999999, 0.097500000000000003, 0.10250000000000001, 0.1075, 0.1125, 0.11750000000000001, 0.1225, 0.1275, 0.13250000000000001, 0.13750000000000001, 0.14250000000000002, 0.14749999999999999, 0.1525, 0.1575, 0.16250000000000001, 0.16750000000000001, 0.17250000000000001, 0.17750000000000002, 0.1825, 0.1875, 0.1925, 0.19750000000000001, 0.20250000000000001, 0.20750000000000002, 0.21249999999999999, 0.2175, 0.2225, 0.22750000000000001, 0.23250000000000001, 0.23750000000000002, 0.24249999999999999, 0.2475, 0.2525, 0.25750000000000001, 0.26250000000000001, 0.26750000000000002, 0.27250000000000002, 0.27750000000000002, 0.28250000000000003, 0.28750000000000003, 0.29249999999999998, 0.29749999999999999, 0.30249999999999999, 0.3075, 0.3125, 0.3175, 0.32250000000000001, 0.32750000000000001, 0.33250000000000002, 0.33750000000000002, 0.34250000000000003, 0.34750000000000003, 0.35250000000000004, 0.35749999999999998, 0.36249999999999999, 0.36749999999999999, 0.3725, 0.3775, 0.38250000000000001, 0.38750000000000001, 0.39250000000000002, 0.39750000000000002]

y_data = [0.0, 2714.0, 6116.0, 7371.0, 9134.0, 9134.0, 8445.0, 8257.0, 8446.0, 9257.0, 11076.0, 12533.0, 12761.0, 11939.0, 12002.0, 12343.0, 13220.0, 14164.0, 14584.0, 15273.0, 15179.0, 14264.0, 14603.0, 15420.0, 15361.0, 15480.0, 16446.0, 17508.0, 19001.0, 19043.0, 18968.0, 18596.0, 18686.0, 18480.0, 19462.0, 19672.0, 20841.0, 22097.0, 23246.0, 24085.0, 24126.0, 23526.0, 22624.0, 22491.0, 21883.0, 21290.0, 21063.0, 19879.0, 18824.0, 17181.0, 15809.0, 15745.0, 17262.0, 19663.0, 22068.0, 22657.0, 22121.0, 20647.0, 20100.0, 19418.0, 19477.0, 19453.0, 20282.0, 20706.0, 21027.0, 20365.0, 19868.0, 18350.0, 17041.0, 15142.0, 13702.0, 12106.0, 10833.0, 10004.0, 9283.0, 8700.0, 8607.0, 8500.0, 8861.0, 9132.0]

# 注：原代码存在变量名错误，xdata/ydata应为x_data/y_data
params, params_covariance = optimize.curve_fit(correction_function, x_data, y_data, p0 = [2.0e6, 0.2, 1.5], bounds=((0.1e6, 0.15, 1.4), (10.0e6, 0.25, 1.6)))

• 修正基础变量名错误：原代码中curve_fit调用时误用了未定义的xdata和ydata，需替换为实际定义的x_data和y_data，否则会直接触发运行错误。

• 核心问题：模型与数据不匹配：观察y_data的分布，其呈现双峰形态，但当前的correction_function是单峰函数（先上升后下降的单峰曲线），单峰模型无法拟合双峰数据，因此拟合算法会将alpha_推到边界，试图让模型尽可能贴近数据，这是根本原因。

• 调整模型结构适配双峰数据：将模型修改为两个单峰函数的叠加，以适配双峰分布：

  def correction_function(z, A1, z_c1, alpha_1, A2, z_c2, alpha_2):
      peak1 = A1 * np.power(z, 2) * np.exp(-np.power((z/z_c1), alpha_1))
      peak2 = A2 * np.power(z, 2) * np.exp(-np.power((z/z_c2), alpha_2))
      return peak1 + peak2
  
  # 对应调整初始猜测和边界
  p0 = [2.0e6, 0.2, 1.5, 1.5e6, 0.3, 1.5]
  bounds = ((0.1e6, 0.15, 0.5, 0.1e6, 0.25, 0.5), (10.0e6, 0.25, 5, 10.0e6, 0.35, 5))
  params, params_covariance = optimize.curve_fit(correction_function, x_data, y_data, p0=p0, bounds=bounds)
  

• 放宽参数边界并优化初始值（单峰模型尝试）：如果必须使用单峰模型，先大幅放宽alpha_的边界范围，同时调整A和z_c的初始值，让模型先贴近数据的主要峰值：

  params, params_covariance = optimize.curve_fit(correction_function, x_data, y_data, 
                                                 p0=[3.0e6, 0.25, 2.0], 
                                                 bounds=((0.1e6, 0.1, 0.5), (10.0e6, 0.4, 5)),
                                                 max_nfev=10000)
  

• 添加拟合权重：如果数据两端的低y值区域噪声较大，可以通过sigma参数降低其权重，避免拟合被干扰：

  # 给低y值数据加更大的sigma（即更低权重）
  sigma = np.where(y_data < 5000, 10, 1)
  params, params_covariance = optimize.curve_fit(correction_function, x_data, y_data, 
                                                 p0=[2.0e6, 0.2, 1.5], 
                                                 bounds=((0.1e6, 0.15, 1.4), (10.0e6, 0.25, 1.6)),
                                                 sigma=sigma)
  

内容的提问来源于stack exchange，提问作者Vitor Silva