How to implement linear interpolation?
-
27-10-2019 - |
문제
Say I am given data as follows:
x = [1, 2.5, 3.4, 5.8, 6]
y = [2, 4, 5.8, 4.3, 4]
I want to design a function that will interpolate linearly between 1
and 2.5
, 2.5
to 3.4
, and so on using Python.
I have tried looking through this Python tutorial, but I am still unable to get my head around it.
해결책
As I understand your question, you want to write some function y = interpolate(x_values, y_values, x)
, which will give you the y
value at some x
? The basic idea then follows these steps:
- Find the indices of the values in
x_values
which define an interval containingx
. For instance, forx=3
with your example lists, the containing interval would be[x1,x2]=[2.5,3.4]
, and the indices would bei1=1
,i2=2
- Calculate the slope on this interval by
(y_values[i2]-y_values[i1])/(x_values[i2]-x_values[i1])
(iedy/dx
). - The value at
x
is now the value atx1
plus the slope multiplied by the distance fromx1
.
You will additionally need to decide what happens if x
is outside the interval of x_values
, either it's an error, or you could interpolate "backwards", assuming the slope is the same as the first/last interval.
Did this help, or did you need more specific advice?
다른 팁
import scipy.interpolate
y_interp = scipy.interpolate.interp1d(x, y)
print y_interp(5.0)
scipy.interpolate.interp1d
does linear interpolation by and can be customized to handle error conditions.
I thought up a rather elegant solution (IMHO), so I can't resist posting it:
from bisect import bisect_left
class Interpolate(object):
def __init__(self, x_list, y_list):
if any(y - x <= 0 for x, y in zip(x_list, x_list[1:])):
raise ValueError("x_list must be in strictly ascending order!")
x_list = self.x_list = map(float, x_list)
y_list = self.y_list = map(float, y_list)
intervals = zip(x_list, x_list[1:], y_list, y_list[1:])
self.slopes = [(y2 - y1)/(x2 - x1) for x1, x2, y1, y2 in intervals]
def __getitem__(self, x):
i = bisect_left(self.x_list, x) - 1
return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])
I map to float
so that integer division (python <= 2.7) won't kick in and ruin things if x1
, x2
, y1
and y2
are all integers for some iterval.
In __getitem__
I'm taking advantage of the fact that self.x_list is sorted in ascending order by using bisect_left
to (very) quickly find the index of the largest element smaller than x
in self.x_list
.
Use the class like this:
i = Interpolate([1, 2.5, 3.4, 5.8, 6], [2, 4, 5.8, 4.3, 4])
# Get the interpolated value at x = 4:
y = i[4]
I've not dealt with the border conditions at all here, for simplicity. As it is, i[x]
for x < 1
will work as if the line from (2.5, 4) to (1, 2) had been extended to minus infinity, while i[x]
for x == 1
or x > 6
will raise an IndexError
. Better would be to raise an IndexError in all cases, but this is left as an exercise for the reader. :)
Instead of extrapolating off the ends, you could return the extents of the y_list
. Most of the time your application is well behaved, and the Interpolate[x]
will be in the x_list
. The (presumably) linear affects of extrapolating off the ends may mislead you to believe that your data is well behaved.
Returning a non-linear result (bounded by the contents of
x_list
andy_list
) your program's behavior may alert you to an issue for values greatly outsidex_list
. (Linear behavior goes bananas when given non-linear inputs!)Returning the extents of the
y_list
forInterpolate[x]
outside ofx_list
also means you know the range of your output value. If you extrapolate based onx
much, much less thanx_list[0]
orx
much, much greater thanx_list[-1]
, your return result could be outside of the range of values you expected.def __getitem__(self, x): if x <= self.x_list[0]: return self.y_list[0] elif x >= self.x_list[-1]: return self.y_list[-1] else: i = bisect_left(self.x_list, x) - 1 return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])
Your solution did not work in Python 2.7. There was an error while checking for the order of the x elements. I had to change to code to this to get it to work:
from bisect import bisect_left
class Interpolate(object):
def __init__(self, x_list, y_list):
if any([y - x <= 0 for x, y in zip(x_list, x_list[1:])]):
raise ValueError("x_list must be in strictly ascending order!")
x_list = self.x_list = map(float, x_list)
y_list = self.y_list = map(float, y_list)
intervals = zip(x_list, x_list[1:], y_list, y_list[1:])
self.slopes = [(y2 - y1)/(x2 - x1) for x1, x2, y1, y2 in intervals]
def __getitem__(self, x):
i = bisect_left(self.x_list, x) - 1
return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])
Building on Lauritz` answer, here's a version with the following changes
- Updated to python3 (the map was causing problems for me and is unnecessary)
- Fixed behavior at edge values
- Raise exception when x is out of bounds
- Use
__call__
instead of__getitem__
from bisect import bisect_right
class Interpolate:
def __init__(self, x_list, y_list):
if any(y - x <= 0 for x, y in zip(x_list, x_list[1:])):
raise ValueError("x_list must be in strictly ascending order!")
self.x_list = x_list
self.y_list = y_list
intervals = zip(x_list, x_list[1:], y_list, y_list[1:])
self.slopes = [(y2 - y1) / (x2 - x1) for x1, x2, y1, y2 in intervals]
def __call__(self, x):
if not (self.x_list[0] <= x <= self.x_list[-1]):
raise ValueError("x out of bounds!")
if x == self.x_list[-1]:
return self.y_list[-1]
i = bisect_right(self.x_list, x) - 1
return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])
Example usage:
>>> interp = Interpolate([1, 2.5, 3.4, 5.8, 6], [2, 4, 5.8, 4.3, 4])
>>> interp(4)
5.425