scipy curve_fit处理大X值拟合错误的原因及拟合方法咨询

问题背景 #

为分析时间趋势，使用scipy curve_fit拟合来自time.time()的X值（量级达16亿，如1663847528.7147126）时，线性拟合频繁出现错误结果，即使提供近似初始参数p0也无法解决，X的量级是导致问题的关键因素。以下是复现问题的代码片段：

import scipy.optimize

def fit_func(x, a, b):
    return a + b * x

y = list(range(5))

x = [1e8 + a for a in range(5)]
print(scipy.optimize.curve_fit(fit_func, x, y, p0=[-x[0], 0]))
# 结果正确：
#   (array([-1.e+08,  1.e+00]), array([[ 0., -0.],
#          [-0.,  0.]]))

x = [1e9 + a for a in range(5)]
print(scipy.optimize.curve_fit(fit_func, x, y, p0=[-x[0], 0.0]))
# 结果不正确：
#   OptimizeWarning: Covariance of the parameters could not be estimated
#   warnings.warn('Covariance of the parameters could not be estimated',
#   (array([-4.53788811e+08,  4.53788812e-01]), array([[inf, inf],
#          [inf, inf]]))

# 为b设置近乎完美的p0可消除警告，但curve_fit仍无法正常工作
print(scipy.optimize.curve_fit(fit_func, x, y, p0=[-x[0], 0.99]))
# 结果不正确：
#   (array([-7.60846335e+10,  7.60846334e+01]), array([[-1.97051972e+19,  1.97051970e+10],
#          [ 1.97051970e+10, -1.97051968e+01]]))
   
# 但设置完全正确的p0时可以正常工作
print(scipy.optimize.curve_fit(fit_func, x, y, p0=[-x[0], 1.0]))
#(array([-1.e+09,  1.e+00]), array([[inf, inf],
#       [inf, inf]]))

测试环境：Python 3.9.6、SciPy 1.7.1、Windows 10。

异常原因 #

核心问题是数值稳定性与参数尺度不匹配：

1. 浮点数精度损失：当X量级达到1e9时，拟合参数a（截距）的量级与b*x的量级相当（如b=1时，b*x为1e9，a约为-1e9）。大数相减的运算会消耗浮点数的有效精度，导致损失函数的计算出现误差，优化算法无法准确迭代。
2. 优化算法对参数尺度敏感：curve_fit默认的Levenberg-Marquardt算法依赖参数的梯度更新，当参数间量级差异过大（如a是1e9级，b是1级），算法的步长调整会出现混乱，无法收敛到正确值。
3. 初始参数的微小偏差被放大：即使初始参数接近正确值，浮点数精度限制会让迭代过程中的微小误差被放大，导致参数更新偏离正确方向，比如设置p0=[-x[0], 0.99]时，计算a + b*x的误差会被X的大尺度放大，最终得到完全错误的参数。

解决方案 #

1. 对X值进行中心化处理 #

这是最直接有效的方法，将X值减去基准值（如第一个样本的X值或X的均值），把X的量级压缩到极小范围：

import scipy.optimize
import numpy as np

def fit_func(x, a, b):
    return a + b * x

y = np.arange(5)
x = np.array([1e9 + a for a in range(5)])
# 中心化X：减去第一个样本值，将X转为[0,1,2,3,4]
x_centered = x - x[0]
# 拟合中心化后的X
params, cov = scipy.optimize.curve_fit(fit_func, x_centered, y, p0=[0, 1.0])
# 转换回原始X对应的参数
a_original = params[0] - params[1] * x[0]
b_original = params[1]
print(f"原始参数：a={a_original:.0f}, b={b_original:.1f}")
# 输出：原始参数：a=-1000000000, b=1.0

2. 使用专用的线性/多项式拟合工具 #

对于线性或多项式拟合，无需使用通用的curve_fit，专用工具内部会优化数值稳定性，效率更高：

• 线性拟合：使用numpy.polyfit或scipy.stats.linregress

import numpy as np
from scipy.stats import linregress

x = np.array([1e9 + a for a in range(5)])
y = np.arange(5)

# numpy.polyfit 一阶拟合（线性）
coeffs = np.polyfit(x, y, 1)
print(f"polyfit结果：b={coeffs[0]}, a={coeffs[1]}")
# 输出：polyfit结果：b=1.0, a=-1000000000.0

# scipy.stats.linregress
result = linregress(x, y)
print(f"linregress结果：b={result.slope}, a={result.intercept}")
# 输出：linregress结果：b=1.0, a=-1000000000.0

• 二阶多项式拟合：使用numpy.polyfit指定阶数为2

# 二阶多项式拟合
coeffs_2 = np.polyfit(x, y, 2)
print(f"二阶多项式系数：x²项={coeffs_2[0]}, x项={coeffs_2[1]}, 常数项={coeffs_2[2]}")
# 输出：二阶多项式系数：x²项=0.0, x项=1.0, 常数项=-1000000000.0

3. 手动缩放参数（若必须用curve_fit） #

如果一定要使用curve_fit，可以对参数进行缩放，让所有参数的量级保持一致：

def fit_func_scaled(x, a_scaled, b):
    # 将a缩放为1e9级的倍数，让a_scaled和b量级相当
    return a_scaled * 1e9 + b * x

x = np.array([1e9 + a for a in range(5)])
y = np.arange(5)
params, cov = scipy.optimize.curve_fit(fit_func_scaled, x, y, p0=[-1.0, 1.0])
a_original = params[0] * 1e9
b_original = params[1]
print(f"缩放拟合结果：a={a_original:.0f}, b={b_original:.1f}")
# 输出：缩放拟合结果：a=-1000000000, b=1.0

总结 #

• 大X值导致拟合异常的核心是浮点数精度损失和参数尺度不匹配。
• 优先选择X中心化或使用专用拟合工具，这两种方法既高效又能保证数值稳定性。
• 若必须使用curve_fit，则需要对X或参数进行缩放处理，确保参数量级一致。

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