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

问题描述

我正在尝试移植一个使用手动插值器(由数学家学院开发)的程序来使用 scipy 提供的插值器.我想使用或包装 scipy 插值器，使其行为尽可能接近旧的插值器.

这两个函数之间的一个主要区别在于，在我们的原始内插器中 - 如果输入值高于或低于输入范围，我们的原始内插器将外推结果.如果您使用 scipy 插值器尝试此操作，则会引发 ValueError.以这个程序为例:

将 numpy 导入为 np从 scipy 导入插值x = np.arange(0,10)y = np.exp(-x/3.0)f = interpolate.interp1d(x, y)打印 f(9)print f(11) # 导致 ValueError，因为它大于 max(x)

有没有一种合理的方法可以使它不会崩溃，最后一行将简单地进行线性外推，将第一个和最后两个点定义的梯度延续到无穷大.

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

解决方案

1.常数外推

您可以使用 scipy 的 interp 函数，它将左右值推断为超出范围的常量:

>>>从 scipy 导入 interp、arange、exp>>>x = 范围(0,10)>>>y = exp(-x/3.0)>>>interp([9,10], x, y)数组([0.04978707，0.04978707])

2.线性(或其他自定义)外推

您可以为处理线性外推的插值函数编写一个包装器.例如:

from scipy.interpolate import interp1dfrom scipy import arange, array, expdef extrap1d(插值器):xs = 插值器.xys = interpolator.y逐点定义(x):如果 x xs[-1]:返回 ys[-1]+(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])别的:返回插值器(x)def ufunclike(xs):返回数组(列表(地图(逐点，数组(xs))))返回类似 ufunc 的

extrap1d 接受一个内插函数并返回一个也可以外推的函数.你可以这样使用它:

x = arange(0,10)y = exp(-x/3.0)f_i = interp1d(x, y)f_x = extrap1d(f_i)打印 f_x([9,10])

输出:

[ 0.04978707 0.03009069]

I'm trying to port a program which uses a hand-rolled interpolator (developed by a mathematician colleage) over to use the interpolators provided by scipy. I'd like to use or wrap the scipy interpolator so that it has as close as possible behavior to the old interpolator.

A key difference between the two functions is that in our original interpolator - if the input value is above or below the input range, our original interpolator will extrapolate the result. If you try this with the scipy interpolator it raises a ValueError. Consider this program as an example:

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)

Is there a sensible way to make it so that instead of crashing, the final line will simply do a linear extrapolate, continuing the gradients defined by the first and last two points to infinity.

Note, that in the real software I'm not actually using the exp function - that's here for illustration only!

解决方案

1. Constant extrapolation

You can use interp function from scipy, it extrapolates left and right values as constant beyond the range:

>>> 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. Linear (or other custom) extrapolation

You can write a wrapper around an interpolation function which takes care of linear extrapolation. For example:

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 takes an interpolation function and returns a function which can also extrapolate. And you can use it like this:

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

print f_x([9,10])

Output:

[ 0.04978707  0.03009069]

这篇关于如何使 scipy.interpolate 给出超出输入范围的外推结果?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！