大X值的scipy curve_fit不正确

ddrv8njm  于 2022-11-10  发布在  其他
关注(0)|答案(2)|浏览(219)

为了确定一段时间内的趋势,我使用scipy curve_fittime.time()中的X值,例如1663847528.7147126(16亿)。进行线性插值有时会产生错误的结果,提供近似的初始p0值也没有帮助。我发现X的大小是造成这种错误的关键因素,我想知道为什么?
下面是一个简单的代码片段,显示了工作和非工作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进行测试。

b4wnujal

b4wnujal1#

根源

您面临两个问题:

  • 拟合过程是尺度敏感的,这意味着在特定变量上选择的单位(例如μA而不是kA)可能会人为地阻止算法正确收敛(例如,一个变量比另一个变量大几个数量级,并主导回归);
  • 浮点算术错误。当从1e8切换到1e9时,您刚刚达到这种错误占主导地位的幅度。

第二个问题非常重要,我们假设你只能用8位有效数字来表示,那么1 000 000 0001 000 000 001是相同的数字,因为它们都被限制为1.0000000e9,我们不能准确地表示1.0000000_e9,因为1.0000000_e9需要多一位数字(_)。这就是第二个例子失败的原因。
此外,您正在使用非线性最小二乘算法来解决线性最小二乘问题,这也在某种程度上与您的问题有关。
您有三种解决方案:

  • 标准化;
  • 规范和变更方法/算法;
  • 提高机床精度。

我会选择第一个,因为它更通用,第二个是由@blunova提出的,完全有意义,后者可能是一个固有的限制。

标准化

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

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'

精度附录

Numpy的默认浮点精度为float64,大约有15个有效数字:

x.dtype                            # dtype('float64')
np.finfo(np.float64).precision     # 15

但是scipy.curve_fit依赖于scipy.least_square,后者利用平方度量来驱动优化。
在不深入研究细节的情况下,我怀疑这就是问题发生的地方,当处理所有接近1e9的值时,您将达到浮点算术错误占主导地位的阈值。
因此,你所达到的1e9的阈值与变量x上的数字之间的区别无关(float64有足够的精度,使它几乎完全不同),而是与求解时对它的使用有关:

minimize F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)
subject to lb <= x <= ub`

您还可以在其签名中检查容差是否约为80倍宽:

scipy.optimize.least_squares(fun, x0, jac='2-point', bounds=(- inf, inf),
    method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0,
    loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, 
    tr_options={}, jac_sparsity=None, max_nfev=None, verbose=0,
    args=(), kwargs={})

这可能会让您调整算法,在达到收敛之前添加额外的步骤(如果是这样的话),但这不会取代或击败规范化的有用性。

方法比较

scipy.stats.linregress方法的有趣之处在于设计中处理的尺度容差。该方法使用变量归一化和纯线性代数以及数值稳定性技巧(参见TINY变量)来解决LS问题,即使是在有问题的条件下。
这当然与scipy.optimize.curve_fit方法形成对比,scipy.optimize.curve_fit方法是作为优化梯度下降算法实现的NLLS解算器(参见Levenberg–Marquardt algorithm)。
如果你坚持线性最小二乘问题(线性的参数而不是变量,所以二阶多项式是LLS),那么LLS可能是一个更简单的选择,因为它可以为你处理归一化。

kwvwclae

kwvwclae2#

如果您只需要计算线性拟合,我认为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

相关问题