如何使scipy.interpolate给出超出输入范围的推断结果?


问题内容

我正在尝试移植一个使用手推插值器(由数学家colleage开发)的程序,以使用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函数-此处仅用于说明!


问题答案:

1.恒定外推

您可以使用interpscipy中的函数,它将左值和右值推断为超出范围的常数:

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