使用scipy.optimize.minimize拟合函数和的参数优化异常求助

问题背景 #

我有x、y及y_error三组数据，范围如下：

• x数据：[0,40]
• y数据：[10^-14, 10^-12]
• y_error数据：[10^-15, 10^-14]

拟用两个指数函数的和拟合曲线，模型为：

y = y_e * exp(-b((x/x_e)^(1/n) -1)) + yO * exp(-x / h)
其中 b = 1.9992*n - 0.3271，自由参数为y_e、x_e、n、yO、h。

采用scipy.optimize.minimize结合负对数似然（NLL）方法拟合，代码如下：

对数似然函数 #

def se_log_likelihood(params, x, y, y_error):
    y_e, xe, n, yO, h = params #free parameters which we want to constrain
    b = 1.9992*n - 0.3271
    model = yO * np.exp(-x / h) + y_e * np.exp(-b * ((x/x_e)**(1/n) - 1))
    sigma2 = y_error**2 #variance
    #log_likelihood function 
    log_likelihood = -0.5 * np.sum((I - model)**2 / sigma2 + np.log(2 * np.pi * sigma2))
    return log_likelihood

负对数似然函数 #

nll_se = lambda *args: - se_log_likelihood(*args) #NLL function

优化过程 #

initial = np.array([5*10**-10, 10, 10, 5*10**-12,12]) #initial estimates for our parameters

soln = minimize(nll_se, initial,args=(x_data ,y_data,y_error_data)) 

遇到的问题 #

1. 默认优化方法（L-BFGS-B）下，返回参数与初始值完全一致，优化失败，提示Desired error not necessarily achieved due to precision loss，迭代次数nit=0；
2. 改用Nelder-mead方法后，优化状态显示成功，但拟合结果与预期曲线差距极大。

问题原因分析 #

1. 似然函数存在笔误 #

代码中对数似然函数里的(I - model)是明显错误，应该为(y - model)，这会导致目标函数计算完全偏离真实需求，直接让优化失去意义。

2. 数据与参数尺度不匹配 #

y数据量级在10^-14到10-12，但初始参数y_e设为5e-10，比实际y值大3-5个数量级。小数值的精度损失会让梯度-based优化算法无法有效计算迭代方向，直接导致默认方法停在初始值。

3. 初始值设置严重偏离合理范围 #

• y_e初始值远大于y数据的最大值，导致模型第一项的输出量级和真实数据完全不符；
• n设为10时，(x/x_e)^(1/n)趋近于1，第一项指数项近似为y_e，进一步放大模型和数据的差距；
• 不合理的初始值让无梯度的Nelder-mead方法极易陷入局部最优，得到完全不符合预期的结果。

4. 模型非线性程度高 #

模型包含(x/x_e)^(1/n)这类强非线性项，参数n的变化对模型输出的影响是非线性的，增加了优化难度，无梯度方法很难找到全局最优解。

解决建议 #

1. 修正似然函数笔误 #

将(I - model)替换为(y - model)，确保目标函数计算正确：

log_likelihood = -0.5 * np.sum((y - model)**2 / sigma2 + np.log(2 * np.pi * sigma2))

2. 统一数据与参数尺度 #

对y数据和相关参数做尺度变换，避免小数值导致的精度问题。例如将y、y_error乘以1e12，对应调整y_e、yO的初始值：

# 缩放数据
y_scaled = y_data * 1e12
y_error_scaled = y_error_data * 1e12

# 定义缩放后的似然函数
def se_log_likelihood_scaled(params, x, y_scaled, y_error_scaled):
    y_e_scaled, xe, n, yO_scaled, h = params
    b = 1.9992*n - 0.3271
    # 还原参数到原尺度计算模型后再缩放
    model = (yO_scaled/1e12)*np.exp(-x/h) + (y_e_scaled/1e12)*np.exp(-b*((x/xe)**(1/n)-1))
    model_scaled = model * 1e12
    sigma2 = y_error_scaled**2
    log_likelihood = -0.5 * np.sum((y_scaled - model_scaled)**2 / sigma2 + np.log(2 * np.pi * sigma2))
    return log_likelihood

nll_se_scaled = lambda *args: -se_log_likelihood_scaled(*args)

# 缩放后的初始值（对应原尺度：y_e=0.5e-12，yO=5e-12）
initial_scaled = np.array([0.5, 10, 2, 5, 12])

3. 优化初始值设置 #

根据数据范围和模型特性调整初始值：

• y_e：参考y数据上限设为1e-12左右；
• yO：参考y数据下限设为1e-14左右；
• n：从2-5这类较小值开始尝试，避免指数项过于平缓；
• x_e：设为x数据中间值（如20），让x/x_e范围在0-2之间；
• h：设为x数据范围量级（如20-40），让x/h范围在0-2之间。

4. 选择合适的优化策略 #

• 先用Nelder-mead方法得到较优初始结果，再用梯度-based方法（如L-BFGS-B）做精细优化；
• 给参数添加物理意义的边界约束，避免取值异常：

  bounds = [(1e-15, 1e-11), (1, 50), (0.5, 20), (1e-16, 1e-12), (1, 100)] # (y_e, x_e, n, yO, h)的边界
  soln = minimize(nll_se, initial, args=(x_data, y_data, y_error_data), method='L-BFGS-B', bounds=bounds)
  

5. 分步拟合降低难度 #

先单独拟合yO * exp(-x/h)项，得到yO和h的合理初始值，再加入第一个项进行整体拟合，减少优化变量的耦合度。

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