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

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

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 by sastanin for How to make scipy.interpolate give an extrapolated result beyond the input range?

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

Answer by Justin Peel for How to make scipy.interpolate give an extrapolated result beyond the input range?

I'm afraid that there is no easy to do this in Scipy to my knowledge. You can, as I'm fairly sure that you are aware, turn off the bounds errors and fill all function values beyond the range with a constant, but that doesn't really help. See this question on the mailing list for some more ideas. Maybe you could use some kind of piecewise function, but that seems like a major pain.

Answer by Giovanni for How to make scipy.interpolate give an extrapolated result beyond the input range?

I did it by adding a point to my initial arrays. In this way I avoid defining self-made functions, and the linear extrapolation (in the example below: right extrapolation) looks ok.

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)    

Answer by ryggyr for How to make scipy.interpolate give an extrapolated result beyond the input range?

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(xxp[-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[-1]]= yp[-1] + (x[x>xp[-1]]-xp[-1])*(yp[-1]-yp[-2])/(xp[-1]-xp[-2])      return y  

Answer by Joma for How to make scipy.interpolate give an extrapolated result beyond the input range?

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 by subnivean for How to make scipy.interpolate give an extrapolated result beyond the input range?

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

Answer by I?igo Hernez Corres for How to make scipy.interpolate give an extrapolated result beyond the input range?

It may be faster to use boolean indexing with large datasets, since the algorithm checks if every point is in outside the interval, whereas boolean indexing allows an easier and faster comparison.

For example:

# 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])  

In my case, with a data set of 300000 points, this means an speed up from 25.8 to 0.094 seconds, this is more than 250 times faster.

Answer by Kashyap Narasimha for How to make scipy.interpolate give an extrapolated result beyond the input range?

The below code gives u the simple extrapolation module. 'k' is the value to which the data set 'y' has to be extrapolated based on the data set 'x'. 'numpy' module is needed

 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)  

Answer by Federico Gonzlez Tamez for How to make scipy.interpolate give an extrapolated result beyond the input range?

Standard interpolate + linear extrapolate:

    def interpola(v, x, y):          if v <= x[0]:              return y[0]+(y[1]-y[0])/(x[1]-x[0])*(v-x[0])          elif v >= x[-1]:              return y[-2]+(y[-1]-y[-2])/(x[-1]-x[-2])*(v-x[-2])          else:              f = interp1d(x, y, kind='cubic')               return f(v)  

Answer by Moot for How to make scipy.interpolate give an extrapolated result beyond the input range?

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  

(Note: I would have just added this as a comment instead of an answer due to the short scope, but I don't have the 50 required points)

Fatal error: Call to a member function getElementsByTagName() on a non-object in D:\XAMPP INSTALLASTION\xampp\htdocs\endunpratama9i\www-stackoverflow-info-proses.php on line 72