Как заставить scipy.interpolate дать экстраполированный результат за пределами диапазона ввода?
Я пытаюсь перенести программу, которая использует ручной ролл-интерполятор (разработанный математиком-коллегой) для использования интерполяторов, предоставляемых scipy. Я хотел бы использовать или обрезать scipy-интерполятор так, чтобы он был как можно ближе к старому интерполятору.
Ключевое различие между двумя функциями заключается в том, что в нашем исходном интерполяторе - если входное значение выше или ниже диапазона ввода, наш исходный интерполятор экстраполирует результат. Если вы попробуете это с помощью scipy-интерполятора, он поднимает значение ValueError. Рассмотрим эту программу как пример:
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 greater than max(x)
Есть ли разумный способ сделать так, чтобы вместо сбоя конечная строка просто выполнила линейную экстраполяцию, продолжая градиенты, определяемые первой и последней двумя точками до бесконечности.
Обратите внимание, что в реальном программном обеспечении я фактически не использую функцию exp - это только для иллюстрации!
Ответы
Ответ 1
1. Постоянная экстраполяция
Вы можете использовать функцию interp из scipy, она экстраполирует значения слева и справа как постоянные за пределы диапазона:
>>> 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. Линейная (или другая обычная) экстраполяция
Вы можете написать обертку вокруг функции интерполяции, которая заботится о линейной экстраполяции. Например:
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 принимает функцию интерполяции и возвращает функцию, которая также может экстраполироваться. И вы можете использовать его следующим образом:
x = arange(0,10)
y = exp(-x/3.0)
f_i = interp1d(x, y)
f_x = extrap1d(f_i)
print f_x([9,10])
Вывод:
[ 0.04978707 0.03009069]
Ответ 2
Вы можете посмотреть InterpolatedUnivariateSpline
Вот пример использования:
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()
Ответ 3
Начиная с версии SciPy 0.17.0, существует новая опция scipy.interpolate.interp1d, которая позволяет экстраполяцию. Просто установите fill_value = 'extrapolate' в вызове. Модификация кода таким образом дает:
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)
а выход:
0.0497870683679
0.010394302658
Ответ 4
Как насчет scipy.interpolate.splrep(со степенью 1 и без сглаживания):
>> 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
Кажется, что вы делаете то, что хотите, поскольку 34 = 25 + (25 - 16).
Ответ 5
Здесь альтернативный метод, который использует только пакет numpy. Он использует функции numpy array, поэтому может быть быстрее при интерполяции/экстраполяции больших массивов:
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)
Изменить: Марк Микофски предложил модификацию функции "extrap":
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
Быстрее использовать булево индексирование с большими наборами данных, так как алгоритм проверяет, находится ли каждая точка вне интервала, тогда как булевское индексирование позволяет легче и быстрее сравнение.
Например:
# 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])
В моем случае с набором данных 300000 пунктов это означает, что скорость увеличивается с 25,8 до 0,094 секунды, это более чем в 250 раз быстрее.
Ответ 7
Я сделал это, добавив точку в мои начальные массивы. Таким образом я избегаю определения самодельных функций, и линейная экстраполяция (в приведенном ниже примере: правильная экстраполяция) выглядит нормально.
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)
Ответ 8
Я боюсь, что в Scipy нелегко сделать это, насколько мне известно. Вы можете, поскольку я уверен, что вы знаете, отключите ошибки границ и заполните все значения функций за пределами диапазона константой, но это действительно не помогает. См. этот вопрос в списке рассылки для некоторых идей. Возможно, вы могли бы использовать какую-то кусочную функцию, но это похоже на большую боль.
Ответ 9
В приведенном ниже коде вы получите простой модуль экстраполяции. k - значение, на которое набор данных y должен быть экстраполирован на основе набора данных x. Требуется модуль numpy.
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)
Ответ 10
Стандартная интерполяция + линейная экстраполяция:
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)
