Comment faire en sorte que scipy.interpolate donne un résultat extrapolé au-delà de la plage d'entrée ?

Question

Je tente de porter un programme qui utilise un interpolateur fait maison (développé par un collègue mathématicien) pour utiliser les interpolateurs fournis par scipy. Je voudrais utiliser ou envelopper l'interpolateur scipy pour qu'il ait un comportement aussi proche que possible de l'ancien interpolateur.

Une différence clé entre les deux fonctions est que dans notre interpolateur original - si la valeur d'entrée est au-dessus ou en dessous de la plage d'entrée, notre interpolateur original extrapolera le résultat. Si vous essayez cela avec l'interpolateur scipy, une ValueError est lancée. Considérez ce programme comme exemple:

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) # Cause une ValueError, car c'est plus grand que max(x)

Y a-t-il un moyen sensé de faire en sorte que, au lieu de planter, la dernière ligne fasse simplement une extrapolation linéaire, prolongeant les gradients définis par les deux premiers et les deux derniers points jusqu'à l'infini.

Notez que dans le logiciel réel, je n'utilise pas réellement la fonction exp - elle est là pour illustration seulement!

Answer 1

À partir de la version 0.17.0 de SciPy, il existe une nouvelle option pour scipy.interpolate.interp1d qui permet l'extrapolation. Il suffit de définir fill_value='extrapolate' dans l'appel. En modifiant votre code de cette manière :

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)

le résultat est le suivant :

0.0497870683679
0.010394302658

Answer 2

Vous pouvez jeter un œil à InterpolatedUnivariateSpline

Voici un exemple d'utilisation:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline

# valeurs données
xi = np.array([0.2, 0.5, 0.7, 0.9])
yi = np.array([0.3, -0.1, 0.2, 0.1])
# positions pour interp/extrapoler
x = np.linspace(0, 1, 50)
# ordre du spline : 1 linéaire, 2 quadratique, 3 cubique ...
order = 1
# faire l'interp/extrapolation
s = InterpolatedUnivariateSpline(xi, yi, k=order)
y = s(x)

# exemple montrant l'interpolation pour linéaire, quadratique et cubique
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 3

1. Extrapolation constante

Vous pouvez utiliser la fonction interp de scipy, elle extrapole les valeurs de gauche et de droite de manière constante au-delà de la plage :

>>> 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. Extrapolation linéaire (ou personnalisée)

Vous pouvez écrire une fonction qui enveloppe une fonction d'interpolation et qui gère 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(list(map(pointwise, array(xs))))

    return ufunclike

extrap1d prend une fonction d'interpolation et renvoie une fonction qui peut également 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])

Sortie :

[ 0.04978707  0.03009069]

Answer 4

Que dire de scipy.interpolate.splrep (avec un degré de 1 et aucun lissage) :

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

Il semble faire ce que vous voulez, puisque 34 = 25 + (25 - 16).

Answer 5

Voici une méthode alternative qui utilise uniquement le package numpy. Elle tire parti des fonctions de tableau de numpy, il peut donc être plus rapide lors de l'interpolation/extrapolation de grands tableaux :

import numpy as np

def extrap(x, xp, yp):
    """fonction np.interp avec extrapolation linéaire"""
    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)

Modification suggérée par Mark Mikofski de la fonction "extrap" :

def extrap(x, xp, yp):
    """fonction np.interp avec extrapolation linéaire"""
    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