如何用leastsq实现最小二乘多项式拟合？附实验与分析数据

嘿，我来一步步教你怎么用scipy.optimize.leastsq完成多项式的最小二乘拟合，同时对比实验值和分析值，生成优化后的拟合曲线。

1. 准备工作：导入依赖库与整理数据 #

首先得把需要的库导全，你给出的代码里缺了leastsq所在的模块，补上就行：

from scipy.integrate import quad
import pylab as py
import numpy as np
from scipy.optimize import leastsq  # 关键模块，别忘啦！

然后把你的数据整理成numpy数组（方便后续计算），补全y_data和y1_data的完整内容：

# 温度数据
x_data = np.array([329.74403, 329.21733, 328.73927, 328.25111, 327.75969, 
                   327.26852, 326.7746, 326.28142, 325.78471, 325.28635, 324.78976])
# 分析值（替换成你的完整数据）
y_data = np.array([...])  
# 实验值（补全省略部分，这里给个示例）
y1_data = np.array([1.55e-06, 1.82e-06, 1.93e-06, 1.17e-06, 1.93e-06, 
                    1.79e-06, 1.31e-06, 1.75e-06, 1.68e-06, 1.52e-06, 1.45e-06])

2. 定义多项式模型与残差函数 #

leastsq的核心是残差函数——也就是拟合值和真实值（这里是实验值y1_data）的差值。我们先以3阶多项式为例，定义模型和残差：

# 3阶多项式模型：y = a0 + a1*x + a2*x² + a3*x³
def poly_model(params, x):
    a0, a1, a2, a3 = params
    return a0 + a1*x + a2*x**2 + a3*x**3

# 残差函数：拟合值 - 实验值（leastsq会最小化残差的平方和）
def residual_func(params, x, y):
    return poly_model(params, x) - y

如果需要其他阶数（比如2阶、4阶），直接修改poly_model里的参数和多项式项就行。

3. 执行最小二乘拟合 #

给拟合参数一个初始猜测值（比如全1），然后调用leastsq得到优化后的系数：

# 初始化参数：3阶多项式需要4个参数，初始值设为[1,1,1,1]
initial_params = np.array([1, 1, 1, 1])

# 执行拟合：传入残差函数、初始参数、以及额外的x和y数据
optimized_params, _ = leastsq(residual_func, initial_params, args=(x_data, y1_data))

optimized_params就是我们最终得到的多项式系数啦。

4. 生成拟合曲线并可视化对比 #

现在用优化后的参数生成平滑的拟合曲线，把实验值、分析值、拟合曲线画在一起直观对比：

# 生成更密集的x点，让拟合曲线更平滑
x_fit = np.linspace(min(x_data), max(x_data), 100)
# 计算拟合值
y_fit = poly_model(optimized_params, x_fit)

# 绘图
py.figure(figsize=(10, 6))
# 实验值散点
py.scatter(x_data, y1_data, label='实验值', color='#ff4444', marker='o', s=60)
# 分析值散点
py.scatter(x_data, y_data, label='分析值', color='#33b5e5', marker='s', s=60)
# 拟合曲线
py.plot(x_fit, y_fit, label=f'3阶多项式拟合曲线', color='#99cc00', linewidth=2)

# 设置图表样式
py.xlabel('温度', fontsize=12)
py.ylabel('数值', fontsize=12)
py.title('实验值、分析值与拟合曲线对比', fontsize=14)
py.legend(fontsize=10)
py.grid(True, alpha=0.3)
py.show()

5. 进阶技巧：适配任意阶数的多项式 #

如果想快速切换多项式阶数，可以把阶数做成变量，写个通用版的函数：

# 通用n阶多项式模型
def poly_model_general(params, x, n):
    return sum(params[i] * x**i for i in range(n+1))

# 通用残差函数
def residual_func_general(params, x, y, n):
    return poly_model_general(params, x, n) - y

# 比如用2阶多项式拟合
n_order = 2
initial_params = np.ones(n_order + 1)  # 阶数+1个参数
optimized_params, _ = leastsq(residual_func_general, initial_params, args=(x_data, y1_data, n_order))

小提醒 #

• 尽量用numpy数组存储数据，比普通列表计算效率高；
• 初始参数别太离谱，比如如果你的y值是1e-6量级，初始参数可以设成[1e-6, 0, 0, 0]这类更贴合的数值，避免拟合不收敛；
• 如果拟合效果不好，可以尝试调整多项式阶数，或者检查数据是否有异常值。

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