如何使scipy.interpolate给出超出input范围的外推结果？
我试图移植一个使用手动内插器（由math家同事开发）的程序来使用scipy提供的内插器。 我想使用或包装scipy插入器，以便它尽可能接近旧插值器的行为。
这两个函数之间的主要区别在于，在我们的原始插补器中 – 如果input值高于或低于input范围，我们的原始插补器将推断结果。 如果你用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.不断的外推
你可以使用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]
你可以看一下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()
从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
什么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）。
这是另一种只使用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)
编辑：标记Mikofskibuild议修改的“extrap”function：
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
使用大数据集的 布尔索引可能会更快，因为algorithm检查每个点是否在区间之外，而布尔索引允许更容易和更快的比较。
例如：
# 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 < fx[0] yo[low] = fy[0] + (xo[low]-fx[0])*(fy[1]-fy[0])/(fx[1]-fx[0]) # Values higher than the maximum "x" are extrapolated at same time high = xo > fx[-1] yo[high] = fy[-1] + (xo[high]-fx[-1])*(fy[-1]-fy[-2])/(fx[-1]-fx[-2]) # Values inside the interpolation range are interpolated directly inside = np.logical_and(xo >= fx[0], xo <= fx[-1]) yo[inside] = f(xo[inside])
在我的情况下，一个300000点的数据集，这意味着从25.8到0.094秒的速度，这是更快的250倍以上 。
我通过给我的初始数组添加一个点来做到这一点。 这样我就避免了自定义函数的定义，并且线性外推（在下面的例子中：右推断）看起来没问题。
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)
恐怕在我看来，Scipy并不容易做到这一点。 你可以，因为我敢肯定，你知道，closures边界错误，并用常量填充范围之外的所有函数值，但这并没有什么帮助。 在邮件列表中查看这个问题 ，获得更多的想法。 也许你可以使用某种分段function，但这似乎是一个主要的痛苦。
标准插值+线性外推：
def interpola(v, x, y): if v <= x[0]: return y[0]+(y[1]-y[0])/(x[1]-x[0])*(vx[0]) elif v >= x[-1]: return y[-2]+(y[-1]-y[-2])/(x[-1]-x[-2])*(vx[-2]) else: f = interp1d(x, y, kind='cubic') return f(v)
下面的代码给你简单的外插模块。 k是数据集y必须根据数据集x外推到的值。 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)
