Comment faire de scipy.interpoler donner une extrapolation de résultats au-delà de l'entrée de gamme?

1. Constante extrapolation

Vous pouvez utiliser interp fonction de scipy, il extrapole à gauche et à droite des valeurs comme constante au-delà de la portée:

>>> 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. Linéaire (ou d'autres personnalisé) extrapolation

Vous pouvez écrire un wrapper autour d'une fonction d'interpolation qui prend soin de l'extrapolation linéaire. Par exemple:

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 prend une fonction d'interpolation et retourne une fonction qui peut aussi extrapoler. Et vous pouvez l'utiliser comme ceci:

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

print f_x([9,10])

De sortie:

[ 0.04978707  0.03009069]

Que de SciPy version 0.17.0, il y a une nouvelle option pour scipy.interpoler.interp1d qui permet l'extrapolation. Réglez simplement fill_value='extrapoler" dans l'appel. La modification de votre code de cette manière donne:

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)

et la sortie est:

0.0497870683679
0.010394302658

Voici une autre méthode qui utilise uniquement le numpy paquet. Il faut profiter de numpy du tableau de fonctions, de manière peut-être plus rapide lorsque l'interpolation/extrapolation de grands tableaux:

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: Marque Mikofski de la modification proposée de la "extrap" fonction:

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

Je l'ai fait par l'ajout d'un point pour mes premiers tableaux. De cette façon, j'ai éviter de définir des self-made fonctions, et l'extrapolation linéaire (dans l'exemple ci-dessous: droit de l'extrapolation) semble 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)  

Standard interpoler + linéaire extrapoler:

    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)