wake-up-neo.com

Wie lässt sich scipy.interpolate außerhalb des Eingabebereichs extrapolieren

Ich versuche, ein Programm zu portieren, das einen handgerollten Interpolator (entwickelt von einem Kollegen eines Mathematikers) verwendet, um die von scipy bereitgestellten Interpolatoren zu verwenden. Ich möchte den Scipy-Interpolator so verwenden oder umschließen, dass er dem alten Interpolator möglichst nahe 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 extrapolieren. Wenn Sie dies mit dem Scipy-Interpolator versuchen, wird eine ValueError ausgelöst. 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 sinnvolle Möglichkeit, so zu machen, dass die letzte Linie nicht abstürzt, sondern eine lineare Extrapolation durchführt, wobei die durch die ersten und letzten beiden Punkte definierten Gradienten auf unendlich gesetzt werden. 

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

65
Salim Fadhley

1. Konstante Hochrechnung

Sie können die interp-Funktion 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 für die lineare Extrapolation sorgt. 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(map(pointwise, array(xs)))

    return ufunclike

extrap1d nimmt eine Interpolationsfunktion und gibt eine Funktion zurück, die auch extrapolieren kann. Und du kannst es so benutzen:

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]
32
sastanin

Sie können einen Blick auf InterpolatedUnivariateSpline werfen.

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

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

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

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

8
subnivean

Hier ist eine alternative Methode, die nur das Paket numpy verwendet. Es nutzt die Array-Funktionen von Numpy und kann daher schneller interpolieren/extrapolieren, wenn große Arrays interpoliert werden:

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 Änderung 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
6
ryggyr

Die Verwendung von boolescher Indexierung mit großen Datensätzen ist möglicherweise schneller, da der Algorithmus überprüft, ob jeder Punkt außerhalb des Intervalls liegt, während die boolesche Indexierung 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 Beschleunigung von 25,8 auf 0,094 Sekunden, dies ist mehr als 250-mal schneller .

Ich habe es getan, indem ich meinen ursprünglichen Arrays einen Punkt hinzugefügt habe. Auf diese Weise vermeide ich, selbst erstellte Funktionen zu definieren, und die lineare Extrapolation (im Beispiel unten: 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)  
2
Giovanni

Ich fürchte, dass es in Scipy meines Wissens nicht einfach ist. Sie können, da bin ich mir ziemlich sicher, dass Sie sich bewusst sind, die Begrenzungsfehler ausschalten und alle Funktionswerte außerhalb des Bereichs mit einer Konstanten füllen. Siehe diese Frage auf der Mailingliste für weitere Ideen. Vielleicht könnten Sie eine Art stückweise Funktion verwenden, aber das scheint ein schwerer Schmerz zu sein.

1
Justin Peel

Standardinterpolation + lineares Extrapolat:

    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)

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 Modul numpy ist erforderlich.

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