Wie kann man scipy.interpolate zu einem extrapolierten Ergebnis außerhalb des Eingabebereichs führen?

Ich versuche, ein Programm zu portieren, das einen handgerollten Interpolator (entwickelt von einem Mathematikerkollegen) verwendet, um die von scipy bereitgestellten Interpolatoren zu verwenden. Ich möchte den Scipy-Interpolator verwenden oder umschließen, damit er dem alten Interpolator so nahe wie möglich kommt.

Ein wesentlicher Unterschied zwischen den beiden Funktionen besteht darin, dass in unserem ursprünglichen Interpolator - wenn der Eingabewert über oder unter dem Eingabebereich liegt, unser ursprünglicher Interpolator das Ergebnis extrapoliert. Wenn Sie dies mit dem Scipy-Interpolator versuchen, wird a ausgelöst ValueError. Betrachten Sie dieses Programm als Beispiel:

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)

Gibt es eine vernünftige Möglichkeit, dies so zu gestalten, dass die letzte Linie anstelle eines Absturzes einfach eine lineare Extrapolation durchführt und die durch die ersten und letzten beiden Punkte definierten Verläufe bis ins Unendliche fortsetzt.

Beachten Sie, dass ich in der realen Software die exp-Funktion nicht verwende - das dient nur zur Veranschaulichung!

1. Ständige Extrapolation

Sie können die interpFunktion von scipy verwenden. Sie extrapoliert linke und rechte Werte als Konstante außerhalb des Bereichs:

>>> 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. Lineare (oder andere benutzerdefinierte) Extrapolation

Sie können einen Wrapper um eine Interpolationsfunktion schreiben, die sich um die lineare Extrapolation kümmert. Zum Beispiel:

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

extrap1dNimmt eine Interpolationsfunktion und gibt eine Funktion zurück, die auch extrapoliert werden kann. Und Sie können es so verwenden:

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

print f_x([9,10])

Ausgabe:

[ 0.04978707  0.03009069]

Sie können sich InterpolatedUnivariateSpline ansehen

Hier ein Beispiel dafür:

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

Ab SciPy Version 0.17.0 gibt es eine neue Option für scipy.interpolate.interp1d , die eine Extrapolation ermöglicht. Setzen Sie im Aufruf einfach fill_value = 'extrapolate'. Wenn Sie Ihren Code auf diese Weise ändern, erhalten Sie:

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)

und die Ausgabe ist:

0.0497870683679
0.010394302658

Was ist mit scipy.interpolate.splrep (mit Grad 1 und ohne Glättung):

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

Es scheint zu tun, was Sie wollen, da 34 = 25 + (25 - 16).

Hier ist eine alternative Methode, die nur das numpy-Paket verwendet. Es nutzt die Array-Funktionen von numpy und kann daher beim Interpolieren / Extrapolieren großer Arrays schneller sein:

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: Mark Mikofskis vorgeschlagene Modifikation der "extrap" -Funktion:

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

Die Verwendung der booleschen Indizierung mit großen Datenmengen kann schneller sein , da der Algorithmus prüft, ob sich jeder Punkt außerhalb des Intervalls befindet, während die boolesche Indizierung einen einfacheren und schnelleren Vergleich ermöglicht.

Zum Beispiel:

# 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 meinem Fall bedeutet dies bei einem Datensatz von 300000 Punkten eine Geschwindigkeit von 25,8 auf 0,094 Sekunden, was mehr als 250-mal schneller ist .

Ich habe es getan, indem ich meinen anfänglichen Arrays einen Punkt hinzugefügt habe. Auf diese Weise vermeide ich die Definition selbst erstellter Funktionen, und die lineare Extrapolation (im folgenden Beispiel: rechte Extrapolation) sieht in Ordnung aus.

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)  

Ich befürchte, dass es meines Wissens in Scipy nicht einfach ist, dies zu tun. Sie können, wie ich ziemlich sicher bin, die Grenzfehler deaktivieren und alle Funktionswerte außerhalb des Bereichs mit einer Konstanten füllen, aber das hilft nicht wirklich. Weitere Ideen finden Sie in dieser Frage auf der Mailingliste. Vielleicht könnten Sie eine stückweise Funktion verwenden, aber das scheint ein großer Schmerz zu sein.

Der folgende Code gibt Ihnen das einfache Extrapolationsmodul. k ist der Wert, auf den der Datensatz y basierend auf dem Datensatz x extrapoliert werden muss . Das numpyModul wird benötigt.

 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)

Standardinterpolation + lineare Extrapolation:

    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)