基于I/O设置curve_fit初始猜测值，优化拟合参数匹配已知结果

核心需求 #

编写代码，基于输入X1_data、X2_data和输出Y_list确定预定义函数f_model的参数a1、a2、b1、b2，需满足以下约束：

• 必须使用Levenberg-Marquardt算法，即调用curve_fit(..., method='lm',...)（根据scipy文档，该方法不支持bounds参数）
• 参数需满足不等式：1 < b1 < a1 < a2 < b2
• 函数形式固定为：

  y=((a2+(a1-a2)/(1+((x1/x2/b2)**b1)))*x2)/1000
  

现有实现代码 #

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def f_model(x_data, b1, a1, a2, b2):
    x1, x2 = x_data
    if 1 < b1 < a1 < a2 < b2:
        return 1e-3 * (a2 + (a1 - a2) / (1 + ((x1 / x2 / b2) ** b1))) * x2
    else:
        return np.inf

X1_data = ...
X2_data = ...
Y_list = ...
X_axis = [i + 1 for i in range(X1_data.shape[0])]

for Y_data in Y_list:
    Y_data = np.array(Y_data)
#    if np.min(Y_data) > 1:
#        b1_guess = 1.01
#        a1_guess = np.min(Y_data)
#        a2_guess = np.max(Y_data) * 10
#        b2_guess = np.mean(X1_data / X2_data)
#        guess = [b1_guess, a1_guess, a2_guess, b2_guess]
#    else:
#        guess = [1.01, 100.0, 700.0, 10000.0]
    guess = [3.5, 100.0, 7300.00, 30000.00]

#    bounds = ([1.0, guess[0], guess[1], guess[2]], [np.inf, np.inf, np.inf, np.inf])
    popt, pcov = curve_fit(f_model, (X1_data, X2_data), Y_data, p0=guess, method='lm', maxfev=5000)
    print(f"Guess {guess}")
    print(f"Optimal {popt}")
    optimal_curve = f_model((X1_data, X2_data), *popt)
    plt.scatter(X_axis, Y_data, label='Data')
    plt.plot(X_axis, optimal_curve, label='Fitted Curve')

plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.legend()
plt.show()

当前存在的问题 #

1. 无法根据输入输出合理设置初始猜测值guess，调整初始值会得到多个不同的拟合结果，难以确定最优解
2. 使用已知参数的dataset1验证时，拟合曲线效果符合预期，但无法还原预设的目标参数：b1=3.98; a1=98.4; a2=7337.03; b2=29404.80
3. 需要修改代码，使其既能还原dataset1的目标参数，同时适配dataset2的拟合需求

测试数据集 #

dataset1（2023年12月13日更新Y_list） #

X1_data = np.array([2894648.773, 3157919.698, 3474257.841, 3800193.356, 4156493.087, 4590298.639, 5042658.177])
X2_data = np.array([67.28579923, 69.52170573, 71.68019425, 73.68463411, 75.5179478, 77.1563269, 78.52232434])
Y_list = [[411.1179, 441.5084, 471.0992, 497.1712, 520.8332, 542.983, 560.8622]]
# 目标参数
b1=3.98; a1=98.4; a2=7337.03; b2=29404.80

dataset2 #

X1_data = np.array([240228, 259164, 267785, 300530, 339812])
X2_data = np.array([32.188719, 32.710433, 33.243463, 33.80985, 34.353048])
Y_list = [[23.79969728, 25.34667811, 26.88872834, 30.86784396, 34.04270604],
          [2.70371836, 3.33359208, 3.49283256, 3.9231444, 4.55001784],
          [0.65156594, 0.66139181, 0.50957883, 0.46295474, 0.51407824]]

内容的提问来源于stack exchange，提问作者CLW-II