scipy.interpolate가 입력 범위를 초과하는 외삽 결과를 제공하도록 만드는 방법은 무엇입니까?

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 함수를 사용하지 않습니다. 여기에서는 설명 용으로 만 사용합니다!

---

1. 상수 외삽

interpscipy의 함수를 사용할 수 있으며 왼쪽 및 오른쪽 값을 범위를 초과하는 상수로 외삽합니다.

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

---

It may be faster to use boolean indexing with large datasets, since the algorithm checks if every point is in outside the interval, whereas boolean indexing allows an easier and faster comparison.

For example:

# 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])

In my case, with a data set of 300000 points, this means an speed up from 25.8 to 0.094 seconds, this is more than 250 times faster.

---

I did it by adding a point to my initial arrays. In this way I avoid defining self-made functions, and the linear extrapolation (in the example below: right extrapolation) looks ok.

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)  

---

I'm afraid that there is no easy to do this in Scipy to my knowledge. You can, as I'm fairly sure that you are aware, turn off the bounds errors and fill all function values beyond the range with a constant, but that doesn't really help. See this question on the mailing list for some more ideas. Maybe you could use some kind of piecewise function, but that seems like a major pain.

---

Standard interpolate + linear extrapolate:

    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)

---

The below code gives you the simple extrapolation module. k is the value to which the data set y has to be extrapolated based on the data set x. The numpy module is required.

 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)

참고URL : https://stackoverflow.com/questions/2745329/how-to-make-scipy-interpolate-give-an-extrapolated-result-beyond-the-input-range