Python中基于上百组观测值求解三次函数系数的方法

三次函数$f(x)=ax^3+bx2+cx+d$仅包含4个未知参数，只有输入4组无误差的精确数据时，才存在能完全匹配所有点的精确解析解。你手上的上百组观测数据是带测量噪声的，方程数量远多于未知数数量，属于超定方程组，不存在能让所有等式严格成立的解析解。sympy是符号计算库，仅支持适定方程组的精确求解，遇到这种超定场景自然会返回空列表。

对于多组带噪声的观测数据求多项式参数，标准做法是使用最小二乘法拟合，核心逻辑是找到一组a、b、c、d参数，让所有数据点的模型预测值与实际观测值的误差平方和最小，不需要严格经过每个数据点。
最简便的实现方式是直接使用numpy库内置的多项式拟合函数，不需要手动实现最小二乘逻辑：

import numpy as np

# 原始观测数据
xs = [28.0, 29.0, 12.0, 12.0, 42.0, 35.0, 28.0, 30.0, 32.0, 46.0, 18.0, 28.0, 28.0, 64.0, 
38.0, 18.0, 49.0, 37.0, 25.0, 24.0, 42.0, 50.0, 12.0, 64.0, 23.0, 35.0, 22.0, 16.0, 44.0, 
77.0, 26.0, 44.0, 38.0, 37.0, 45.0, 42.0, 24.0, 42.0, 12.0, 46.0, 12.0, 26.0, 37.0, 15.0, 
67.0, 36.0, 43.0, 36.0, 45.0, 82.0,
44.0, 30.0, 33.0, 51.0, 50.0]

fxs = [59.5833333333333, 59.5833333333333, 10.0, 10.0, 47.0833333333333, 51.2499999999999, 
34.5833333333333, 88.75, 63.7499999999999, 34.5833333333333, 51.2499999999999, 10.0, 
63.7499999999999, 51.0, 59.5833333333333,47.0833333333333, 49.5625, 43.5624999999999, 
63.7499999999999, 10.0, 76.25, 47.0833333333333,10.0, 51.2499999999999,47.0833333333333,10.0, 
35.0, 51.2499999999999, 76.25, 100.0, 51.2499999999999, 59.5833333333333, 63.7499999999999, 
76.25, 100.0, 51.2499999999999, 10.0, 22.5, 10.0, 88.75, 10.0, 59.5833333333333, 
47.0833333333333, 34.5833333333333, 51.2499999999999, 63.7499999999999,63.7499999999999, 10.0, 
76.25, 62.1249999999999, 47.0833333333333, 10.0, 76.25, 47.0833333333333, 88.75]

# 三次多项式拟合，deg=3对应最高次为3次
a, b, c, d = np.polyfit(xs, fxs, deg=3)
print(f"拟合得到的参数：a={a:.6f}, b={b:.6f}, c={c:.6f}, d={d:.6f}")

# 如需预测新x对应的f(x)，可调用poly1d生成可直接调用的模型
model = np.poly1d((a,b,c,d))
# 例：预测x=30时的f(x)，直接调用model(30)即可

• 原始数据里存在多组x相同但f(x)不同的矛盾样本（比如x=12对应多个不同的f(x)值），这种情况本来就不可能存在能同时满足所有点的函数，用最小二乘拟合是唯一合理的处理方式。
• 如果后续需要更高次的多项式拟合，只需要修改deg参数即可，不需要改动其他逻辑。
• 如果需要评估拟合效果，可以计算所有点的预测值和真实值的均方误差，判断模型是否符合预期。

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