Как заставить 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's greater than max(x) 

Есть ли разумный способ сделать так, чтобы вместо сбоя конечная строка просто выполняла линейную экстраполяцию, продолжая градиенты, определяемые первой и последней двумя точками до бесконечности.

Обратите внимание, что в реальном программном обеспечении я фактически не использую функцию exp – это только для иллюстрации!

10 Solutions collect form web for “Как заставить scipy.interpolate дать экстраполированный результат за пределами диапазона ввода?”

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] 

Вы можете взглянуть на 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() 

Начиная с версии 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 

Что относительно 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).

Вот альтернативный метод, который использует только пакет numpy. Он использует функции массива numpy, поэтому может быть быстрее при интерполяции / экстраполяции больших массивов:

 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) 

Изменить: Предположительная модификация Mark Mikofski в функции «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 

Возможно, быстрее использовать булевое индексирование с большими наборами данных , поскольку алгоритм проверяет, находится ли каждая точка вне интервала, тогда как булевское индексирование позволяет более простое и быстрое сравнение.

Например:

 # 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 < fx[0] yo[low] = fy[0] + (xo[low]-fx[0])*(fy[1]-fy[0])/(fx[1]-fx[0]) # Values higher than the maximum "x" are extrapolated at same time high = xo > fx[-1] yo[high] = fy[-1] + (xo[high]-fx[-1])*(fy[-1]-fy[-2])/(fx[-1]-fx[-2]) # Values inside the interpolation range are interpolated directly inside = np.logical_and(xo >= fx[0], xo <= fx[-1]) yo[inside] = f(xo[inside]) 

В моем случае с набором данных 300000 пунктов это означает, что скорость увеличивается с 25,8 до 0,094 секунды, это более чем в 250 раз быстрее .

Я сделал это, добавив точку в свои начальные массивы. Таким образом я избегаю определения самодельных функций, и линейная экстраполяция (в примере ниже: правая экстраполяция) выглядит нормально.

 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) 

Я боюсь, что в Scipy мне нелегко это сделать. Вы можете, поскольку я уверен, что вы знаете, отключите ошибки границ и заполните все значения функций за пределами диапазона константой, но это действительно не помогает. См. Этот вопрос в списке рассылки для некоторых идей. Возможно, вы могли бы использовать какую-то кусочную функцию, но это похоже на серьезную боль.

Стандартная интерполяция + линейная экстраполяция:

  def interpola(v, x, y): if v <= x[0]: return y[0]+(y[1]-y[0])/(x[1]-x[0])*(vx[0]) elif v >= x[-1]: return y[-2]+(y[-1]-y[-2])/(x[-1]-x[-2])*(vx[-2]) else: f = interp1d(x, y, kind='cubic') return f(v) 

Приведенный ниже код дает вам простой модуль экстраполяции. 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) 
Interesting Posts

Что такое эквивалент Python статических переменных внутри функции?

Проверка подлинности Django с пользовательской моделью с идентификатором электронной почты как уникальным ключом

Python: попытка использования формы POST с использованием запросов

Как подключиться / подключить / сгруппировать несколько объектов в Mayavi2

Создание свойств экземпляра класса из словаря?

Python объединяет два списка со всеми возможными перестановками

python: перезапуск цикла

Самый короткий Sudoku Solver в Python – Как это работает?

фильтровать теги django-taggit в Queryset Django

Как проверить путь к классам, где Jython ScriptEngine ищет модуль python?

Проблема с отступом Python?

Как отправлять / помещать данные json в ListSerializer

С socketserver python как передать переменную в конструктор класса обработчика

Передача данных между страницами в функции redirect () в Google App Engine

Поиск с конца файла, бросающего неподдерживаемое исключение

Python - лучший язык программирования в мире.