如何制作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

你可以看看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 2

从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

Answer 3

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

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

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多倍.