如何让scipy.interpolate在输入范围之外给出外推结果？

Question

我试图将一个使用自定义插值器（由数学家同事开发）的程序移植到使用scipy提供的插值器上。 我想使用或包装scipy插值器，使其行为尽可能接近旧插值器。

两个函数之间的关键区别在于，在我们的原始插值器中 - 如果输入值高于或低于输入范围，我们的原始插值器将对结果进行外推。 如果您尝试使用scipy插值器执行此操作，它会引发 ValueError 。 请考虑以下示例程序：

import numpy as np
from scipy import interpolate

x = np.arange(0,10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y)

print f(9)
print f(11) # Causes ValueError, because it's greater than max(x)

有没有一种明智的方法，使得最后一行不会崩溃，而是通过线性外推来延续由前两个点和最后两个点定义的渐变到无限远处。

请注意，在实际软件中，我并没有使用exp函数 - 这只是为了举例说明！

- Salim Fadhley

---

scipy.interpolate.UnivariateSpline 似乎可以很好地进行外推。 - heltonbiker

Answer 1

从SciPy版本0.17.0开始，scipy.interpolate.interp1d有一个新的选项，允许外推。只需在调用中设置fill_value='extrapolate'。以这种方式修改您的代码如下：

import numpy as np
from scipy import interpolate

x = np.arange(0,10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y, fill_value='extrapolate')

print f(9)
print f(11)

输出结果为：

0.0497870683679
0.010394302658

- Moot

---

外推类型是否类似于内插类型？例如，我们可以使用最近点外推进行线性内插吗？ - a.sam

如果 kind='cubic'，fill_value='extrapolate' 不起作用。 - vlmercado

@a.sam：我不确定你的意思是什么...假设你使用kind='linear'和fill_value='interpolation'，那么你会得到线性插值；如果你使用它并且fill_value='extrapolation'，那么你会得到线性外推，对吧？ - Moot

@vlmercado，确实，我有SciPy 0.19版本，当使用kind='spline'时也会出现错误，所以同意。然而，您之前的评论是关于kind='cubic'无法工作，但似乎它确实可以工作。此外，SciPy 0.19的文档没有将kind='spline'列为有效选项，并且还指定了选项‘zero’、‘slinear’、‘quadratic’和‘cubic’都是相应阶数的样条插值，因此这可能已经满足了所需的样条外推？ - Moot

@Moot，我使用了kind='cubic'，但是ValueError指的是kind=spline。以下是完整代码：import numpy as np from scipy import interpolatex = np.array([ 21, 63, 126, 252, 504, 756, 1260, 1764, 2520, 5040, 7560]) y = np.array([0.77, 0.93, 1.07, 1.16, 1.28, 1.44, 1.76, 2.02, 2.21, 2.61, 2.88]) tck = interpolate.interp1d(x, y, kind='cubic', fill_value='extrapolate')ValueError: Extrapolation does not work with kind=spline - vlmercado

Answer 2

你可以看一下 InterpolatedUnivariateSpline

这里是使用它的一个例子：

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

# given values
xi = np.array([0.2, 0.5, 0.7, 0.9])
yi = np.array([0.3, -0.1, 0.2, 0.1])
# positions to inter/extrapolate
x = np.linspace(0, 1, 50)
# spline order: 1 linear, 2 quadratic, 3 cubic ... 
order = 1
# do inter/extrapolation
s = InterpolatedUnivariateSpline(xi, yi, k=order)
y = s(x)

# example showing the interpolation for linear, quadratic and cubic interpolation
plt.figure()
plt.plot(xi, yi)
for order in range(1, 4):
    s = InterpolatedUnivariateSpline(xi, yi, k=order)
    y = s(x)
    plt.plot(x, y)
plt.show()

- Joma

---

这是最佳答案。这就是我所做的。我使用了k=1（次序），因此它变成了线性插值，并且我使用了bbox=[xmin-w, xmax+w]，其中w是我的公差。 - imbr

Answer 3

1. 常值外推

您可以使用scipy中的interp函数，它会将范围外的左右值外推为常数:

>>> from scipy import interp, arange, exp
>>> x = arange(0,10)
>>> y = exp(-x/3.0)
>>> interp([9,10], x, y)
array([ 0.04978707,  0.04978707])

2. 线性（或其他自定义）外推

你可以编写一个包装器函数，用于处理线性外推的插值函数。例如：

from scipy.interpolate import interp1d
from scipy import arange, array, exp

def extrap1d(interpolator):
    xs = interpolator.x
    ys = interpolator.y

    def pointwise(x):
        if x < xs[0]:
            return ys[0]+(x-xs[0])*(ys[1]-ys[0])/(xs[1]-xs[0])
        elif x > xs[-1]:
            return ys[-1]+(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])
        else:
            return interpolator(x)

    def ufunclike(xs):
        return array(list(map(pointwise, array(xs))))

    return ufunclike

extrap1d接受一个插值函数并返回一个能够进行外推的函数。您可以像这样使用它:

x = arange(0,10)
y = exp(-x/3.0)
f_i = interp1d(x, y)
f_x = extrap1d(f_i)

print f_x([9,10])

输出:

[ 0.04978707  0.03009069]

- sastanin

---

在Python 3.6中，我不得不添加list到返回语句中：return array(list(map(pointwise, array(xs))))以解决迭代器问题。 - mostlyWright

这个解决方案比fill_value="extrapolate"选项更灵活。我能够根据自己的需求调整内部函数“pointwise”，我赞同上面的评论并在需要时插入列表。话虽如此，有时您可能只想要一个生成器。 - Wilmer E. Henao

请注意，基于scipy.interp的第一个解决方案已不再推荐使用，因为它已被弃用并将在SciPy 2.0.0中消失。他们建议改用numpy.interp，但正如问题中所述，这里无法使用。 - Yosko

Answer 4

那么scipy.interpolate.splrep呢（使用1阶和无平滑）：

>> tck = scipy.interpolate.splrep([1, 2, 3, 4, 5], [1, 4, 9, 16, 25], k=1, s=0)
>> scipy.interpolate.splev(6, tck)
34.0

看起来它可以做到你想要的，因为34 = 25 + (25 - 16)。

- subnivean

Answer 5

这里有一种只使用numpy包的替代方法。它利用了numpy的数组函数，因此在插值/外推大型数组时可能更快：

import numpy as np

def extrap(x, xp, yp):
    """np.interp function with linear extrapolation"""
    y = np.interp(x, xp, yp)
    y = np.where(x<xp[0], yp[0]+(x-xp[0])*(yp[0]-yp[1])/(xp[0]-xp[1]), y)
    y = np.where(x>xp[-1], yp[-1]+(x-xp[-1])*(yp[-1]-yp[-2])/(xp[-1]-xp[-2]), y)
    return y

x = np.arange(0,10)
y = np.exp(-x/3.0)
xtest = np.array((8.5,9.5))

print np.exp(-xtest/3.0)
print np.interp(xtest, x, y)
print extrap(xtest, x, y)

编辑：Mark Mikofski建议修改“extrap”函数：

def extrap(x, xp, yp):
    """np.interp function with linear extrapolation"""
    y = np.interp(x, xp, yp)
    y[x < xp[0]] = yp[0] + (x[x<xp[0]]-xp[0]) * (yp[0]-yp[1]) / (xp[0]-xp[1])
    y[x > xp[-1]]= yp[-1] + (x[x>xp[-1]]-xp[-1])*(yp[-1]-yp[-2])/(xp[-1]-xp[-2])
    return y

- ryggyr

---

你可以使用布尔索引和这里y[x < xp[0]] = fp[0] + (x[x < xp[0]] - xp[0]) / (xp[1] - xp[0]) * (fp[1] - fp[0])以及y[x > xp[-1]] = fp[-1] + (x[x > xp[-1]] - xp[-1]) / (xp[-2] - xp[-1]) * (fp[-2] - fp[-1])代替np.where，因为在False选项下，y不会改变。+1 如果提供实际示例。 - Mark Mikofski

Answer 6

当处理大型数据集时，使用布尔索引可能会更快，因为该算法检查每个点是否在区间之外，而布尔索引允许更轻松和更快速的比较。

例如：

# Necessary modules
import numpy as np
from scipy.interpolate import interp1d

# Original data
x = np.arange(0,10)
y = np.exp(-x/3.0)

# Interpolator class
f = interp1d(x, y)

# Output range (quite large)
xo = np.arange(0, 10, 0.001)

# Boolean indexing approach

# Generate an empty output array for "y" values
yo = np.empty_like(xo)

# Values lower than the minimum "x" are extrapolated at the same time
low = xo < f.x[0]
yo[low] =  f.y[0] + (xo[low]-f.x[0])*(f.y[1]-f.y[0])/(f.x[1]-f.x[0])

# Values higher than the maximum "x" are extrapolated at same time
high = xo > f.x[-1]
yo[high] = f.y[-1] + (xo[high]-f.x[-1])*(f.y[-1]-f.y[-2])/(f.x[-1]-f.x[-2])

# Values inside the interpolation range are interpolated directly
inside = np.logical_and(xo >= f.x[0], xo <= f.x[-1])
yo[inside] = f(xo[inside])

在我的情况下，使用300000个数据点，这意味着从25.8秒加速到0.094秒，这是快了250多倍。

- Iñigo Hernáez Corres

---

这很好，但是如果x0是浮点数，y [0]是np.nan或y [-1]是np.nan，则不起作用。 - Stretch

Answer 7

我通过在我的初始数组中添加一个点来实现。这样我就避免了定义自制函数，线性外推（在下面的示例中：向右外推）看起来还不错。

import numpy as np  
from scipy import interp as itp  

xnew = np.linspace(0,1,51)  
x1=xold[-2]  
x2=xold[-1]  
y1=yold[-2]  
y2=yold[-1]  
right_val=y1+(xnew[-1]-x1)*(y2-y1)/(x2-x1)  
x=np.append(xold,xnew[-1])  
y=np.append(yold,right_val)  
f = itp(xnew,x,y)  

- Giovanni

Answer 8

我没有足够的声望来发表评论，但是如果有人正在寻找一个用于scipy线性二维插值的外推包装器，我已经根据这里给出的1d插值答案进行了修改。

def extrap2d(interpolator):
xs = interpolator.x
ys = interpolator.y
zs = interpolator.z
zs = np.reshape(zs, (-1, len(xs)))
def pointwise(x, y):
    if x < xs[0] or y < ys[0]:
        x1_index = np.argmin(np.abs(xs - x))
        x2_index = x1_index + 1
        y1_index = np.argmin(np.abs(ys - y))
        y2_index = y1_index + 1
        x1 = xs[x1_index]
        x2 = xs[x2_index]
        y1 = ys[y1_index]
        y2 = ys[y2_index]
        z11 = zs[x1_index, y1_index]
        z12 = zs[x1_index, y2_index]
        z21 = zs[x2_index, y1_index]
        z22 = zs[x2_index, y2_index]

        return (z11 * (x2 - x) * (y2 - y) +
        z21 * (x - x1) * (y2 - y) +
        z12 * (x2 - x) * (y - y1) +
        z22 * (x - x1) * (y - y1)
       ) / ((x2 - x1) * (y2 - y1) + 0.0)


    elif x > xs[-1] or y > ys[-1]:
        x1_index = np.argmin(np.abs(xs - x))
        x2_index = x1_index - 1
        y1_index = np.argmin(np.abs(ys - y))
        y2_index = y1_index - 1
        x1 = xs[x1_index]
        x2 = xs[x2_index]
        y1 = ys[y1_index]
        y2 = ys[y2_index]
        z11 = zs[x1_index, y1_index]
        z12 = zs[x1_index, y2_index]
        z21 = zs[x2_index, y1_index]
        z22 = zs[x2_index, y2_index]#

        return (z11 * (x2 - x) * (y2 - y) +
        z21 * (x - x1) * (y2 - y) +
        z12 * (x2 - x) * (y - y1) +
        z22 * (x - x1) * (y - y1)
        ) / ((x2 - x1) * (y2 - y1) + 0.0)
    else:
        return interpolator(x, y)
def ufunclike(xs, ys):
    if  isinstance(xs, int) or isinstance(ys, int) or isinstance(xs, np.int32) or isinstance(ys, np.int32):
        res_array = pointwise(xs, ys)
    else:
        res_array = np.zeros((len(xs), len(ys)))
        for x_c in range(len(xs)):
            res_array[x_c, :] = np.array([pointwise(xs[x_c], ys[y_c]) for y_c in range(len(ys))]).T

    return res_array
return ufunclike

我没有写很多注释，也知道代码不是特别干净。如果有人发现任何错误，请告诉我。在我的当前使用情况下，它可以正常运行 :)

- KRS

---

非常有用的代码，谢谢 - 但是我认为存在一个错误。当len(xs)！= len(ys)时，reshape无法正常工作。我认为使用unravel_index可能会解决这个问题。 - 16Aghnar

Answer 9

据我所知，Scipy中没有简单的方法来实现这一点。你可以关闭边界错误并用一个常数填充范围之外的所有函数值，但这并没有什么帮助。请参见邮件列表上的this question以获取更多想法。也许你可以使用某种分段函数，但那似乎很麻烦。

- Justin Peel

---

这是我得出的结论，至少在scipy 0.7中是这样的。然而，21个月前编写的这篇教程表明interp1d函数具有可以设置为“linear”的高和低属性，但该教程不清楚适用于哪个版本的scipy：http://projects.scipy.org/scipy/browser/branches/Interpolate1D/docs/tutorial.rst?rev=4591 - Salim Fadhley

看起来这是一个尚未合并到主版本中的分支的一部分，因此可能仍然存在一些问题。当前的代码位于http://projects.scipy.org/scipy/browser/branches/interpolate/interpolate1d.py，但您可能需要滚动到页面底部并单击以下载纯文本格式。我认为这看起来很有前途，尽管我自己还没有尝试过。 - Justin Peel

Answer 10

下面的代码为您提供了简单的外推模块。 k 是基于数据集 x 对数据集 y 进行外推的值。需要使用numpy 模块。

 def extrapol(k,x,y):
        xm=np.mean(x);
        ym=np.mean(y);
        sumnr=0;
        sumdr=0;
        length=len(x);
        for i in range(0,length):
            sumnr=sumnr+((x[i]-xm)*(y[i]-ym));
            sumdr=sumdr+((x[i]-xm)*(x[i]-xm));

        m=sumnr/sumdr;
        c=ym-(m*xm);
        return((m*k)+c)

- Kashyap Narasimha