Как заставить 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 – это только для иллюстрации!

  • построение автомасштабированных подзаговоров с фиксированными пределами в matplotlib
  • Как вычислить вероятность значения, заданного списком выборок из дистрибутива в Python?
  • Причина, почему numpy rollaxis настолько запутанна?
  • Does Conda заменяет необходимость в virtualenv?
  • Линейная интерполяция Python 4D на прямоугольной сетке
  • Линейная подгонка, включая все ошибки с NumPy / SciPy
  • Ошибка при вызове scikit-learn с использованием сборки AMD64 для Scipy в Windows
  • Как быстро выполнить установку наименьших квадратов по множеству наборов данных?
  • 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) 
    Python - лучший язык программирования в мире.