scipy-如何积分线性插值函数？

Question

我有一个函数，它是一组相对较大的数据的插值。我使用线性插值 interp1d，因此有许多非平滑尖锐点，例如 这个。scipy 的 quad 函数会因为这些尖锐点而发出警告。我想知道如何在没有警告的情况下进行积分呢？

谢谢！

---

感谢所有回答。在这里，我总结一下解决方案，以防其他人遇到同样的问题：

1. 像 @Stelios 所做的那样，使用 points 避免警告并获得更准确的结果。
2. 实际上，点数通常大于 quad 的默认限制（limit=50），所以我选择 quad(f_interp, a, b, limit=2*p.shape[0], points=p) 来避免所有这些警告。
3. 如果 a 和 b 不是数据集 x 的起始或结束点，则可以通过 p = x[where(x>=a and x<=b)] 选择点 p。

- Shu Liao

---

欢迎来到 SO。我建议您添加完整的代码，这样其他人可以更轻松地帮助您。 - petezurich

你的原始数据点是否均匀分布？如果是这种情况，那么积分值不会比数据集的缩放总和多很多。 - JohanL

在积分之前，也许可以使用Savitzky-Golay滤波器平滑函数？ - Anon

Answer 1

quad接受一个可选参数，名为points。根据文档：

points : (sequence of floats,ints), 可选

积分区间内可能出现局部积分困难的断点序列（例如奇点、不连续点）。该序列不需要排序。

在您的情况下，“困难”的points恰好是数据点的x坐标。这是一个例子：

import numpy as np 
from scipy.integrate import quad
np.random.seed(123)

# generate random data set 
x = np.arange(0,10)  
y = np.random.rand(10)

# construct a linear interpolation function of the data set 
f_interp = lambda xx: np.interp(xx, x, y)

这里是数据点和 f_interp 的图示： 现在调用 quad：

quad(f_interp,0,9)

返回一系列警告信息以及

(4.89770017785734, 1.3762838395159349e-05)

如果您提供了points参数，即

quad(f_interp,0,9, points = x)

IT问题没有任何警告，结果为

(4.8977001778573435, 5.437539505167948e-14)

这也意味着与之前的调用相比，结果的精度更高。

- Stelios

---

那么如果我想从0.5到9进行积分，x中的点将超出被积函数的边界，quad将返回一个错误消息。在这种情况下我该怎么办？谢谢！ - Shu Liao

@ShuLiao 如果是这样，你应该提供在区间 [0.5, 9] 内对应于 x 的点。这不是很明显吗？ - Stelios

问题在于我担心像 points=x[np.where(x>0.5)] 这样的操作效率，因为这种集成方式在我的代码中会被频繁使用。所以我很好奇是否有其他方法可以避免这种额外的工作。 - Shu Liao

Answer 2

你可以使用scipy.interpolate.InterpolatedUnivariateSpline代替interp1d。这个插值器有integral(a, b)方法，可以计算定积分。

下面是一个例子：

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


# Create some test data.
x = np.linspace(0, np.pi, 21)
np.random.seed(12345)
y = np.sin(1.5*x) + np.random.laplace(scale=0.35, size=len(x))**3

# Create the interpolator.  Use k=1 for linear interpolation.
finterp = InterpolatedUnivariateSpline(x, y, k=1)

# Create a finer mesh of points on which to compute the integral.
xx = np.linspace(x[0], x[-1], 5*len(x))

# Use the interpolator to compute the integral from 0 to t for each
# t in xx.
qq = [finterp.integral(0, t) for t in xx]

# Plot stuff
p = plt.plot(x, y, '.', label='data')
plt.plot(x, y, '-', color=p[0].get_color(), label='linear interpolation')
plt.plot(xx, qq, label='integral of linear interpolation')
plt.grid()
plt.legend(framealpha=1, shadow=True)
plt.show()

故事情节：

（提示：此内容为图片，无法翻译）

- Warren Weckesser

---

谢谢！这确实解决了函数的积分问题。但是，除了函数本身之外，我还想将该函数乘以其他函数进行积分。例如，如何在不使用“quad”的情况下积分函数lambda x: x*finterp(x)[()]？ - Shu Liao