How to make scipy.interpolate give an extrapolated result beyond the input range?

Question

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!

Answer 1

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

Answer 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