How to make scipy.interpolate give an extrapolated result beyond the input range?
Solution 1
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]
Solution 2
As of SciPy version 0.17.0, there is a new option for scipy.interpolate.interp1d that allows extrapolation. Simply set fill_value='extrapolate' in the call. Modifying your code in this way gives:
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)
and the output is:
0.0497870683679
0.010394302658
Solution 3
You can take a look at InterpolatedUnivariateSpline
Here an example using it:
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()
Solution 4
What about scipy.interpolate.splrep (with degree 1 and no smoothing):
>> 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
It seems to do what you want, since 34 = 25 + (25 - 16).
Solution 5
Here's an alternative method that uses only the numpy package. It takes advantage of numpy's array functions, so may be faster when interpolating/extrapolating large arrays:
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)
Edit: Mark Mikofski's suggested modification of the "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
Comments
-
Gordon Pace almost 2 years
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!
-
Gordon Pace about 14 yearsThat's the conclusion I came to, at least with scipy 0.7, however this tutorial written 21 months ago suggests that the interp1d function has a high and low attribute which can be set to "linear", the tutorial is not clear which version of scipy this applies to: projects.scipy.org/scipy/browser/branches/Interpolate1D/docs/…
-
Justin Peel about 14 yearsIt looks like this is part of a branch that hasn't been assimilated into the main version yet so there might still be some problems with it. The current code for this is at projects.scipy.org/scipy/browser/branches/interpolate/… though you might want to scroll to the bottom of the page and click to download it as plain text. I think that this looks promising though I haven't tried it yet myself.
-
Mark Mikofski almost 12 years+1 for an actual example, but you could also use boolean indexing and here
y[x < xp[0]] = fp[0] + (x[x < xp[0]] - xp[0]) / (xp[1] - xp[0]) * (fp[1] - fp[0])
andy[x > xp[-1]] = fp[-1] + (x[x > xp[-1]] - xp[-1]) / (xp[-2] - xp[-1]) * (fp[-2] - fp[-1])
instead ofnp.where
, since theFalse
option,y
doesn't change. -
Stretch over 10 yearsThis is nice, but it does not work if x0 is a float, if y[0] is np.nan, or if y[-1] is np.nan.
-
imbr almost 10 yearsthis is the best answer. That's what i did.
I used k=1 (order)
, so it becomes a linear interpolation, andI used bbox=[xmin-w, xmax+w] where w is my tolerance
-
Félix Adriyel Gagnon-Grenier over 9 yearsHey Federico! If you wonder why you got downvoted, please be aware that when answering questions you need to actually explain how it solves the problem. This answer, as it is, is only a code dump and should have at least a few sentences explaining why and/or how it is usefull. Thanks!
-
a.sam almost 7 yearsIs the extrapolation kind similar to interpolation kind? For example, can we have linear interpolation with nearest point extrapolation?
-
vlmercado almost 7 yearsIf kind='cubic', fill_value='extrapolate' doesn't work.
-
Moot almost 7 years@a.sam: I'm not sure what you mean... Presumably, if you use kind='linear' with fill_value='interpolation' then you get a linear interpolation, and if you use it with fill_value='extrapolation' then you get a linear extrapolation, no?
-
Moot almost 7 years@vlmercado: can you explain in what way it doesn't work? I tried running the above example with the addition of kind='cubic', and it works fine for me.
-
vlmercado almost 7 years@Moot, using scipy 0.18.1, I get the following: ValueError: Extrapolation does not work with kind=spline
-
Moot almost 7 years@vlmercado, indeed, I have SciPy 0.19 and I get an error as well when using kind='spline', so agreed. However, your previous comment was about kind='cubic' not working, but that does seem to work. Also, the documentation for SciPy 0.19 does not mention kind='spline' as a valid option, and also specifies that the options ‘zero’, ‘slinear’, ‘quadratic’ and ‘cubic’ all refer to splines of the respective order, so maybe that takes care of a wanted spline extrapolation?
-
vlmercado almost 7 years@Moot, I used kind='cubic', but the ValueError refers to kind=spline. Here's the complete code: import numpy as np from scipy import interpolate x = np.array([ 21, 63, 126, 252, 504, 756, 1260, 1764, 2520, 5040, 7560]) y = np.array([0.77, 0.93, 1.07, 1.16, 1.28, 1.44, 1.76, 2.02, 2.21, 2.61, 2.88]) tck = interpolate.interp1d(x, y, kind='cubic', fill_value='extrapolate') ValueError: Extrapolation does not work with kind=spline
-
Moot almost 7 years@vlmercado, I tried your code and it does not give any errors for me. I suggest you try upgrading Scipy to 0.19.0 and try again.
-
vlmercado almost 7 years@Moot, I should've qualified my original comment from the start. I shared the error because your initial answer referred to "As of SciPy version 0.17.0..." I hadn't tested on version 0.19.0, but knew it didn't work on version 0.18.0. I don't need to extrapolation feature at this time. Thanks for your effort, however.
-
mostlyWright about 5 yearsIn Python 3.6, I had to add
list
to the return:return array(list(map(pointwise, array(xs))))
to resolve the iterator. -
Wilmer E. Henao over 4 yearsDoes anybody know the method of extrapolation?. You can't just say it extrapolates to something. There are an infinity of options to extrapolate, throw a random number and that is one way to extrapolate. It would be nice to know.
-
Wilmer E. Henao over 4 yearsThis solution is more flexible than the fill_value="extrapolate" option. I was able to adapt the inner function 'pointwise' to my needs, I second the comment above and insert the list when needed. Having said that, sometimes you may just want to have a generator.
-
Moot over 4 years@WilmerE.Henao Presumably the choice for the kind argument sets both the kind of interpolation and the type of extrapolation (the default is 'linear').
-
Yosko over 3 yearsNote that the first solution based on
scipy.interp
is no longer recommanded as it is deprecated and will disappear in SciPy 2.0.0. They recommand usingnumpy.interp
instead but as stated in the question, it won't work here