Jak sprawić, by scipy.interpolate dawało ekstrapolowany wynik poza zakresem wejściowym?

Próbuję przenieść program, który używa ręcznie obracanego interpolatora (opracowanego przez kolegę matematyka), aby używał interpolatorów dostarczonych przez scipy. Chciałbym użyć lub owinąć interpolator scipy, aby zachowywał się jak najbliżej starego interpolatora.

Kluczową różnicą między tymi dwiema funkcjami jest to, że w naszym oryginalnym interpolatorze - jeśli wartość wejściowa jest powyżej lub poniżej zakresu wejściowego, nasz oryginalny interpolator ekstrapoluje wynik. Jeśli spróbujesz tego z interpolatorem scipy, podniesie się ValueError. Rozważmy ten program jako przykład:

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)

Czy istnieje rozsądny sposób, aby to zrobić, aby zamiast się rozbijać, ostatnia linia po prostu wykona liniową ekstrapolację, kontynuując gradienty zdefiniowane przez pierwszy i ostatni dwa punkty do nieskończoności.

Zauważ, że w prawdziwym oprogramowaniu tak naprawdę nie używam funkcji exp - to jest tutaj tylko dla ilustracji!

1. Ciągła ekstrapolacja

Możesz użyć interpfunkcji z scipy, ekstrapoluje lewą i prawą wartość jako stałą poza zakresem:

>>> 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. Ekstrapolacja liniowa (lub inna niestandardowa)

Możesz napisać opakowanie wokół funkcji interpolującej, która zajmuje się ekstrapolacją liniową. Na przykład:

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

extrap1dprzyjmuje funkcję interpolacji i zwraca funkcję, która może również ekstrapolować. Możesz go używać w ten sposób:

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

print f_x([9,10])

Wynik:

[ 0.04978707  0.03009069]

Możesz spojrzeć na InterpolatedUnivariateSpline

Oto przykład używający go:

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

Od wersji 0.17.0 SciPy jest nowa opcja dla scipy.interpolate.interp1d, która umożliwia ekstrapolację. Po prostu ustaw w wezwaniu fill_value = 'extrapolate'. Modyfikacja kodu w ten sposób daje:

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)

a wynik to:

0.0497870683679
0.010394302658

A co z scipy.interpolate.splrep (z stopniem 1 i bez wygładzania):

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

Wydaje się, że robi to, co chcesz, ponieważ 34 = 25 + (25-16).

Oto alternatywna metoda, która używa tylko pakietu numpy. Wykorzystuje funkcje tablicowe numpy, więc może być szybsze podczas interpolacji / ekstrapolacji dużych tablic:

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)

Edycja: sugerowana przez Marka Mikofskiego modyfikacja funkcji „ekstrapsu”:

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

Szybsze może być użycie indeksowania logicznego w przypadku dużych zbiorów danych , ponieważ algorytm sprawdza, czy każdy punkt znajduje się poza przedziałem, podczas gdy indeksowanie logiczne umożliwia łatwiejsze i szybsze porównanie.

Na przykład:

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

W moim przypadku przy zestawie danych 300000 punktów oznacza to przyspieszenie z 25,8 do 0,094 sekundy, czyli ponad 250 razy szybciej .

Zrobiłem to, dodając punkt do moich początkowych tablic. W ten sposób unikam definiowania samodzielnie utworzonych funkcji, a ekstrapolacja liniowa (w poniższym przykładzie: prawostronna ekstrapolacja) wygląda dobrze.

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)  

Obawiam się, że o ile wiem, w Scipy nie jest to łatwe. Możesz, o ile jestem dość pewien, że wiesz, wyłączyć błędy granic i wypełnić wszystkie wartości funkcji poza zakresem stałą, ale to naprawdę nie pomaga. Zobacz to pytanie na liście mailingowej, aby uzyskać więcej pomysłów. Może mógłbyś użyć jakiejś funkcji odcinkowej, ale to wydaje się być dużym problemem.

Poniższy kod przedstawia prosty moduł ekstrapolacji. k jest wartością, do której należy ekstrapolować zbiór danych y na podstawie zbioru danych x. numpyJest wymagany moduł.

 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)

Standardowy interpolacja + liniowa ekstrapolacja:

    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)

Nie mam wystarczającej reputacji, aby komentować, ale na wypadek, gdyby ktoś szukał opakowania ekstrapolacji dla liniowej interpolacji 2d z scipy, dostosowałem odpowiedź, która została tu podana dla interpolacji 1d.

def extrap2d(interpolator):
xs = interpolator.x
ys = interpolator.y
zs = interpolator.z
zs = np.reshape(zs, (-1, len(xs)))
def pointwise(x, y):
    if x < xs[0] or y < ys[0]:
        x1_index = np.argmin(np.abs(xs - x))
        x2_index = x1_index + 1
        y1_index = np.argmin(np.abs(ys - y))
        y2_index = y1_index + 1
        x1 = xs[x1_index]
        x2 = xs[x2_index]
        y1 = ys[y1_index]
        y2 = ys[y2_index]
        z11 = zs[x1_index, y1_index]
        z12 = zs[x1_index, y2_index]
        z21 = zs[x2_index, y1_index]
        z22 = zs[x2_index, y2_index]

        return (z11 * (x2 - x) * (y2 - y) +
        z21 * (x - x1) * (y2 - y) +
        z12 * (x2 - x) * (y - y1) +
        z22 * (x - x1) * (y - y1)
       ) / ((x2 - x1) * (y2 - y1) + 0.0)


    elif x > xs[-1] or y > ys[-1]:
        x1_index = np.argmin(np.abs(xs - x))
        x2_index = x1_index - 1
        y1_index = np.argmin(np.abs(ys - y))
        y2_index = y1_index - 1
        x1 = xs[x1_index]
        x2 = xs[x2_index]
        y1 = ys[y1_index]
        y2 = ys[y2_index]
        z11 = zs[x1_index, y1_index]
        z12 = zs[x1_index, y2_index]
        z21 = zs[x2_index, y1_index]
        z22 = zs[x2_index, y2_index]#

        return (z11 * (x2 - x) * (y2 - y) +
        z21 * (x - x1) * (y2 - y) +
        z12 * (x2 - x) * (y - y1) +
        z22 * (x - x1) * (y - y1)
        ) / ((x2 - x1) * (y2 - y1) + 0.0)
    else:
        return interpolator(x, y)
def ufunclike(xs, ys):
    if  isinstance(xs, int) or isinstance(ys, int) or isinstance(xs, np.int32) or isinstance(ys, np.int32):
        res_array = pointwise(xs, ys)
    else:
        res_array = np.zeros((len(xs), len(ys)))
        for x_c in range(len(xs)):
            res_array[x_c, :] = np.array([pointwise(xs[x_c], ys[y_c]) for y_c in range(len(ys))]).T

    return res_array
return ufunclike

Nie komentowałem zbyt wiele i mam świadomość, że kod nie jest super czysty. Jeśli ktoś zauważy błędy, daj mi znać. W moim obecnym przypadku działa bez problemu :)