Scipy Lagrange 除以零答案

【问题标题】：Scipy Lagrange division by zeroScipy Lagrange 除以零
【发布时间】：2021-09-30 20:49:48
【问题描述】：

我正在尝试使用 Numpy 的拉格朗日插值法对一组有序对进行插值；我以前做过这件事，没有发生任何事故。然而，这一次，我不断得到“除以零误差”，并且插值多项式出现了无限系数。我知道由于拉格朗日方法的内部工作原理，数据点不能重复，并且它们不会重复。

这是我的代码和有问题的有序对，采用 numpy 矢量格式。代码：

x = out["x"].round(decimals=3)
x = np.array(x)
y = out["y"].round(decimals=3)
y = np.array(y)
print(x)
print(y)
pol = lagrange(x,y)
print(pol)

有序对：

[273.324 285.579 309.292 279.573 297.427 290.681 276.621 293.586 283.463
 284.674 273.904 288.064 280.125 294.269 288.51  285.898 273.419 273.023
 281.754 281.546 283.21  303.399 297.392 293.359 306.404 356.285 302.487
 280.586 299.487 302.487]
[ 0.     5.414  6.202  0.     9.331 11.52   0.    10.495  5.439  4.709
  0.     4.916  0.    10.508  6.736  5.25   0.     0.     6.53   4.305
  5.124  6.753 10.175 10.545  5.98   9.147 11.137  0.     8.764  9.57 ]

非常感谢。

【问题讨论】：

您是否注意到他们警告不要使用超过 20 分？您可以想象创建 19 次多项式的问题。
作为旁注，四舍五入会破坏数据，并且可能是您最终得到重复数据的方式。您应该始终保持完整分辨率，并且仅在您打印数据时进行四舍五入。毕竟，只有人类有处理多个小数位的问题。
感谢您的意见。我不知道不要使用它超过 20 点的警告。它确实工作不到 20 点。

标签： python numpy

【解决方案1】：

为什么拉格朗日插值对您不起作用。

您的数组x 中有两次值302.487。 IE。你确实重复了。

为什么拉格朗日插值不是您想要的。

正如 Tim Roberts 所指出的，拉格朗日插值实际上并不是针对 20 点进行的。问题是高度多项式倾向于过拟合。查看wikipedia article of overfitting.中的以下示例

图 2. 将噪声（大致线性）数据拟合到线性函数和多项式函数。尽管多项式函数是完美拟合的，但可以预期线性函数具有更好的泛化能力：如果将这两个函数用于推断拟合数据之外，线性函数应该可以做出更好的预测。

替代回归

至少有两个有效的替代方案。其中之一是维基百科文章中推荐的内容。如果您知道您的数据来自哪种类型的函数，请使用回归以将该类型的函数拟合到数据中。在上面的例子中，这是一个线性函数。如果你想这样做，请查看scipy's curve fit。

替代样条插值

另一种选择是样条插值。再次来自Spline Interpolation的维基百科文章

样条插值不是一次将单个高次多项式拟合到所有值，而是将低次多项式拟合到值的小子集，例如，在每对十个点之间拟合九个三次多项式，而不是为所有这些拟合一个单一的 10 次多项式。样条插值通常优于多项式插值，因为即使在样条中使用低次多项式时，插值误差也可以很小。样条插值也避免了龙格现象的问题，即使用高次多项式插值时点之间会发生振荡。

我只想指出两个小的技术细节。第一点是您需要订购积分，所以我为您做了。并且两个 scipy 的 UnivariateSpline 有一个平滑参数 s 你需要选择。如果你把它选小，它会像你习惯使用拉格朗日插值一样坚持数据，但如果你把它变大，它会变得更平滑，并希望能更好地概括。下面我选择了 2 个不同的值供您查看，但您可能应该自己尝试一下。我包括了一个很小的，所以你会看到它可以做你习惯于拉格朗日插值的事情，但不推荐它。此外，您可能应该使用更多数据，对其进行预处理等。但这不是问题所在。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

idx = np.argsort(x)
x = x[idx]
y = y[idx]

for s in [10,60]:
    t = np.linspace(np.min(x), np.max(x), 10**4)
    f = UnivariateSpline(x,y, s=s)

    plt.scatter(x,y)
    plt.plot(t,f(t))
    plt.title(f'{s=}')
    plt.show()

【讨论】：