Как заставить scipy.interpolate дать экстраполируемый результат вне входного диапазона?
Я пытаюсь портировать программу, которая использует скрученный вручную интерполятор (разработанный математиком colleage) для использования интерполяторов, обеспеченных scipy. Я хотел бы использовать или перенести scipy интерполятор так, чтобы он имел максимально близко поведение к старому интерполятору.
Основное отличие между двумя функциями - то, что в нашем исходном интерполяторе - если входное значение выше или ниже входного диапазона, наш исходный интерполятор будет экстраполировать результат. При попытке этого scipy интерполятором, он повышает a 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's greater than max(x)
Есть ли разумный способ сделать его так, чтобы вместо катастрофического отказа, заключительная строка просто сделала, линейное экстраполирует, продолжая градиенты, определенные первыми и последними двумя точками к бесконечности.
Обратите внимание на то, что в реальном программном обеспечении я на самом деле не использую функцию exp - это здесь для иллюстрации только!
2 ответа
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])
Output:
[ 0.04978707 0.03009069]
Боюсь, насколько мне известно, сделать это в Scipy непросто. Вы можете, как я почти уверен, что вы знаете, отключить ошибки границ и заполнить все значения функций, выходящие за пределы диапазона, константой, но это на самом деле не помогает. См. этот вопрос в списке рассылки для получения дополнительных идей. Возможно, вы могли бы использовать какую-то кусочную функцию, но это кажется большой проблемой.