为了确定一段时间内的趋势,我使用scipy curve_fit
和time.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进行测试。
2条答案
按热度按时间b4wnujal1#
根源
您面临两个问题:
1e8
切换到1e9
时,您刚刚达到这种错误占主导地位的幅度。第二个问题非常重要,我们假设你只能用8位有效数字来表示,那么
1 000 000 000
和1 000 000 001
是相同的数字,因为它们都被限制为1.0000000e9
,我们不能准确地表示1.0000000_e9
,因为1.0000000_e9
需要多一位数字(_
)。这就是第二个例子失败的原因。此外,您正在使用非线性最小二乘算法来解决线性最小二乘问题,这也在某种程度上与您的问题有关。
您有三种解决方案:
我会选择第一个,因为它更通用,第二个是由
@blunova
提出的,完全有意义,后者可能是一个固有的限制。标准化
为了缓解这两个问题,一个常见的解决方案是标准化。在您的情况下,简单的标准化就足够了:
或者使用
sklearn
接口得到相同的结果:反向转换
实际上,当归一化时,拟合结果是以归一化变量的形式表示的。要获得所需的拟合参数,您只需要做一点数学运算,将回归参数转换回原始变量标度。
只需写下并求解转换:
这将为您提供以下解决方案:
精度附录
Numpy的默认浮点精度为
float64
,大约有15个有效数字:但是
scipy.curve_fit
依赖于scipy.least_square
,后者利用平方度量来驱动优化。在不深入研究细节的情况下,我怀疑这就是问题发生的地方,当处理所有接近
1e9
的值时,您将达到浮点算术错误占主导地位的阈值。因此,你所达到的
1e9
的阈值与变量x
上的数字之间的区别无关(float64
有足够的精度,使它几乎完全不同),而是与求解时对它的使用有关:您还可以在其签名中检查容差是否约为80倍宽:
这可能会让您调整算法,在达到收敛之前添加额外的步骤(如果是这样的话),但这不会取代或击败规范化的有用性。
方法比较
scipy.stats.linregress
方法的有趣之处在于设计中处理的尺度容差。该方法使用变量归一化和纯线性代数以及数值稳定性技巧(参见TINY
变量)来解决LS问题,即使是在有问题的条件下。这当然与
scipy.optimize.curve_fit
方法形成对比,scipy.optimize.curve_fit
方法是作为优化梯度下降算法实现的NLLS解算器(参见Levenberg–Marquardt algorithm)。如果你坚持线性最小二乘问题(线性的参数而不是变量,所以二阶多项式是LLS),那么LLS可能是一个更简单的选择,因为它可以为你处理归一化。
kwvwclae2#
如果您只需要计算线性拟合,我认为
curve_fit
是不必要的,我也会使用SciPy中的linregress
函数:一般来说,如果需要多项式拟合,我会使用NumPy polyfit。