如何使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函数-这里只是为了说明！

Answer 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(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]

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()

Answer 3

从SciPy version 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