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


82

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!


2
scipy.interpolate.UnivariateSplinewydaje się ekstrapolować bez problemów.
heltonbiker

Odpowiedzi:


38

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]

1
W Pythonie 3.6 musiałem dodać listdo zwrotu: return array(list(map(pointwise, array(xs)))) aby rozwiązać iterator.
głównie Wright

To rozwiązanie jest bardziej elastyczne niż opcja fill_value = "extrapolate". Udało mi się dostosować wewnętrzną funkcję „punktowo” do moich potrzeb, popieram powyższy komentarz i wstawiam listę w razie potrzeby. Powiedziawszy to, czasami możesz po prostu chcieć mieć generator.
Wilmer E. Henao

1
Zauważ, że pierwsze rozwiązanie oparte na scipy.interpnie jest już zalecane, ponieważ jest przestarzałe i zniknie w SciPy 2.0.0. Zalecają używanie numpy.interpzamiast tego, ale jak stwierdzono w pytaniu, to nie zadziała tutaj
Yosko

86

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

2
to najlepsza odpowiedź. To jest to co zrobiłem. I used k=1 (order), więc staje się interpolacją liniową iI used bbox=[xmin-w, xmax+w] where w is my tolerance
eusoubrasileiro

78

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

Czy rodzaj ekstrapolacji jest podobny do rodzaju interpolacji? Na przykład, czy możemy mieć interpolację liniową z ekstrapolacją najbliższego punktu?
a.sam

Jeśli kind = 'cubic', fill_value = 'extrapolate' nie działa.
vlmercado

@ a.sam: Nie jestem pewien, co masz na myśli ... Przypuszczalnie, jeśli użyjesz kind = 'linear' z fill_value = 'interpolation', otrzymasz liniową interpolację, a jeśli użyjesz jej z fill_value = 'extrapolation' wtedy otrzymujesz liniową ekstrapolację, nie?
Moot

@vlmercado: czy możesz wyjaśnić, w jaki sposób to nie działa? Próbowałem uruchomić powyższy przykład z dodatkiem kind = 'cubic' i dla mnie działa dobrze.
Moot

@Moot, używając scipy 0.18.1, otrzymuję następujące informacje: ValueError: Ekstrapolacja nie działa z kind = spline
vlmercado

8

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


7

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

2
+1 za rzeczywisty przykład, ale możesz również użyć indeksowania logicznego i tutaj y[x < xp[0]] = fp[0] + (x[x < xp[0]] - xp[0]) / (xp[1] - xp[0]) * (fp[1] - fp[0]) i y[x > xp[-1]] = fp[-1] + (x[x > xp[-1]] - xp[-1]) / (xp[-2] - xp[-1]) * (fp[-2] - fp[-1])zamiast np.where, ponieważ Falseopcja ynie zmienia się.
Mark Mikofski

6

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 .


To fajne, ale nie działa, jeśli x0 jest zmiennoprzecinkową, jeśli y [0] to np.nan, lub jeśli y [-1] to np.nan.
Rozciągnij

2

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)  

1

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.


Do takiego wniosku doszedłem, przynajmniej przy scipy 0.7, jednak ten samouczek napisany 21 miesięcy temu sugeruje, że funkcja interp1d ma atrybut high i low, który można ustawić na "linear", samouczek nie jest jasny, która wersja scipy to dotyczy: projects.scipy.org/scipy/browser/branches/Interpolate1D/docs/…
Salim Fadhley

Wygląda na to, że jest to część gałęzi, która nie została jeszcze zasymilowana do głównej wersji, więc mogą nadal występować z nią problemy. Bieżący kod do tego znajduje się na projects.scipy.org/scipy/browser/branches/interpolate/ ... chociaż możesz chcieć przewinąć do dołu strony i kliknąć, aby pobrać go jako zwykły tekst. Myślę, że wygląda to obiecująco, chociaż sam jeszcze tego nie próbowałem.
Justin Peel

1

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)

0

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)

1
Hej Federico! Jeśli zastanawiasz się, dlaczego zostałeś odrzucony, pamiętaj, że odpowiadając na pytania, musisz faktycznie wyjaśnić, w jaki sposób rozwiązuje to problem. Ta odpowiedź w obecnej postaci jest tylko zrzutem kodu i powinna zawierać co najmniej kilka zdań wyjaśniających, dlaczego i / lub jak jest przydatna. Dzięki!
Félix Gagnon-Grenier

0

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

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.