¿Cómo hacer scipy.interpolate dar un resultado extrapolado más allá del rango de entrada?

Estoy tratando de trasladar un progtwig que utiliza un interpolador manual (desarrollado por un colega matemático) para usar los interpoladores provistos por scipy. Me gustaría usar o envolver el interpolador scipy para que tenga el mayor comportamiento posible con el interpolador anterior.

Una diferencia clave entre las dos funciones es que en nuestro interpolador original: si el valor de entrada está por encima o por debajo del rango de entrada, nuestro interpolador original extrapolará el resultado. Si lo intentas con el interpolador scipy, genera un ValueError . Considere este progtwig como ejemplo:

 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) 

¿Existe una forma sensata de hacerlo de modo que, en lugar de estrellarse, la línea final simplemente realice una extrapolación lineal, continuando los gradientes definidos por el primero y los últimos dos puntos hasta el infinito?

Tenga en cuenta que en el software real no estoy usando la función exp. ¡Eso está aquí solo como ilustración!

1. Extrapolación constante.

Puede usar la función interp de scipy, extrapola valores de izquierda y derecha como constantes más allá del rango:

 >>> 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. Extrapolación lineal (u otra costumbre).

Puede escribir una envoltura alrededor de una función de interpolación que se encargue de la extrapolación lineal. Por ejemplo:

 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 toma una función de interpolación y devuelve una función que también puede extrapolar. Y puedes usarlo así:

 x = arange(0,10) y = exp(-x/3.0) f_i = interp1d(x, y) f_x = extrap1d(f_i) print f_x([9,10]) 

Salida:

 [ 0.04978707 0.03009069] 

Puedes echar un vistazo a InterpolatedUnivariateSpline

Aquí un ejemplo usándolo:

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

A partir de la versión 0.17.0 de SciPy, existe una nueva opción para scipy.interpolate.interp1d que permite la extrapolación. Simplemente configure fill_value = ‘extrapolar’ en la llamada. Modificar su código de esta manera le da:

 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) 

y la salida es:

 0.0497870683679 0.010394302658 

¿Qué pasa con scipy.interpolate.splrep (con grado 1 y sin suavizado):

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

Parece hacer lo que quieres, ya que 34 = 25 + (25 – 16).

Aquí hay un método alternativo que usa solo el paquete numpy. Aprovecha las funciones de matriz de Numpy, por lo que puede ser más rápido al interpolar / extrapolar matrices grandes:

 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) 

Edición: la modificación sugerida por Mark Mikofski de la función “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[-1]]= yp[-1] + (x[x>xp[-1]]-xp[-1])*(yp[-1]-yp[-2])/(xp[-1]-xp[-2]) return y 

Puede ser más rápido utilizar la indexación booleana con grandes conjuntos de datos , ya que el algoritmo verifica si cada punto está fuera del intervalo, mientras que la indexación booleana permite una comparación más fácil y rápida.

Por ejemplo:

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

En mi caso, con un conjunto de datos de 300000 puntos, esto significa una aceleración de 25.8 a 0.094 segundos, esto es más de 250 veces más rápido .

Lo hice agregando un punto a mis arreglos iniciales. De esta manera, evito definir funciones hechas por nosotros mismos, y la extrapolación lineal (en el siguiente ejemplo: extrapolación correcta) se ve bien.

 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) 

Me temo que no hay nada fácil de hacer en Scipy, que yo sepa. Puede, como estoy bastante seguro de que está al tanto, desactivar los errores de límites y rellenar todos los valores de funciones más allá del rango con una constante, pero eso realmente no ayuda. Vea esta pregunta en la lista de correo para obtener más ideas. Tal vez podrías usar algún tipo de función por partes, pero eso parece ser un gran dolor.

Interpolación estándar + extrapolación lineal:

  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) 

El siguiente código le da el módulo de extrapolación simple. k es el valor al que se debe extrapolar el conjunto de datos y en función del conjunto de datos x. Se numpy módulo 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)