scipy curve_fit 对于大 X 值不正确答案

【问题标题】：scipy curve_fit incorrect for large X valuesscipy curve_fit 对于大 X 值不正确
【发布时间】：2022-09-22 23:55:27
【问题描述】：

为了确定一段时间内的趋势，我使用scipy curve_fit 和来自time.time() 的X 值，例如1663847528.7147126（16 亿）。进行线性插值有时会产生错误的结果，并且提供近似的初始 p0 值无济于事。我发现 X 的大小是这个错误的关键因素，我想知道为什么？

这是一个简单的 sn-p 显示工作和非工作 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]))
# Result is correct:
#   (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]))
# Result is not correct:
#   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]]))

Almost perfect p0 for b removes the warning but still curve_fit doesn\'t work
print(scipy.optimize.curve_fit(fit_func, x, y, p0=[-x[0], 0.99]))
# Result is not correct:
#   (array([-7.60846335e+10,  7.60846334e+01]), array([[-1.97051972e+19,  1.97051970e+10],
#          [ 1.97051970e+10, -1.97051968e+01]]))
   
# ...but perfect p0 works
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]]))

作为一个附带问题，也许有一种更有效的线性拟合方法？不过，有时我想找到二阶多项式拟合。

在 Windows 10 下使用 Python 3.9.6 和 SciPy 1.7.1 进行测试。

拟合过程对比例敏感。规范化可能是您需要的。

标签： python scipy scipy-optimize

【解决方案1】：

如果您只需要计算线性拟合，我相信 curve_fit 不是必需的，我也将使用 SciPy 中的 linregress 函数：

>>> from scipy import stats

>>> y = list(range(5))

>>> x = [1e8 + a for a in range(5)]
>>> stats.linregress(x, y)
LinregressResult(slope=1.0, intercept=-100000000.0, rvalue=1.0, pvalue=1.2004217548761408e-30, stderr=0.0, intercept_stderr=0.0)

>>> x2 = [1e9 + a for a in range(5)]
>>> stats.linregress(x2, y)
LinregressResult(slope=1.0, intercept=-1000000000.0, rvalue=1.0, pvalue=1.2004217548761408e-30, stderr=0.0, intercept_stderr=0.0)

一般来说，如果您需要多项式拟合，我会使用 NumPy polyfit。

【讨论】：

【解决方案2】：

根本原因

你面临两个问题：

装配过程对比例敏感。这意味着在特定变量上选择的单位（例如，µA 而不是 kA）可以人为地阻止算法正确收敛（例如，一个变量比另一个变量大几个数量级并主导回归）；
浮点算术错误。当从1e8 切换到1e9 时，当这种错误占主导地位时，您就会遇到严重的问题。

第二个是非常重要的实现。假设您仅限于 8 位有效数字表示，那么 1 000 000 000 和 1 000 000 001 是相同的数字，因为它们都仅限于本文 1.0000000e9 并且我们无法准确表示 1.0000000_e9 需要多一位数字 (@987654327 @)。这就是您的第二个示例失败的原因。

此外，您正在使用非线性最小二乘算法来解决线性最小二乘问题，但这与您的问题无关。

你有两个解决方案：

在执行计算时提高机器精度；
标准化您的问题。

我会选择第二个，因为它更通用。

正常化

为了缓解这两个问题，一个常见的解决方案是标准化。在您的情况下，一个简单的标准化就足够了：

import numpy as np
import scipy.optimize

y = np.arange(5)
x = 1e9 + y

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

xm = np.mean(x)         # 1000000002.0
xs = np.std(x)          # 1.4142135623730951

result = scipy.optimize.curve_fit(fit_func, (x - xm)/xs, y)

# (array([2.        , 1.41421356]),
# array([[0., 0.],
#        [0., 0.]]))

# Back transformation:
a = result[0][1]/xs                    # 1.0
b = result[0][0] - xm*result[0][1]/xs  # -1000000000.0

或者使用sklearn 接口得到同样的结果：

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.linear_model import LinearRegression

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("regressor", LinearRegression())
])

pipe.fit(x.reshape(-1, 1), y)

pipe.named_steps["scaler"].mean_          # array([1.e+09])
pipe.named_steps["scaler"].scale_         # array([1.41421356])
pipe.named_steps["regressor"].coef_       # array([1.41421356])
pipe.named_steps["regressor"].intercept_  # 2.0

反向变换

实际上，当归一化拟合结果时，然后用归一化变量表示。要获得所需的拟合参数，您只需做一些数学运算即可将回归参数转换回原始变量比例。

简单地写下并解决转换：

 y = x'*a' + b'
x' = (x - m)/s
 y = x*a + b

这为您提供了以下解决方案：

a = a'/s
b = b' - m/s*a'

【讨论】：