44

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.

6
  • This is ... not easy. What have you tried?
    – zellio
    Sep 8, 2011 at 5:59
  • 1
    -1 as way too general. you don't understand how to program, or how to do the algorithm in python??
    – steabert
    Sep 8, 2011 at 6:19
  • Well being a new learner I've thrown myself into the deep end so to speak. I was thinking of using 'for' or 'if' statements in a algorithm. So between numerous ranges of x.
    – Helpless
    Sep 8, 2011 at 6:26
  • IMO, being a newbie is not a good excuse for not even trying—in fact it's exactly the opposite (i.e. an excellent reason to do so).
    – martineau
    Jun 3, 2018 at 22:16

7 Answers 7

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

1
  • 11
    The question asks how to implement the function, not what library provides it. Aug 30, 2019 at 13:47
34

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:

  1. Find the indices of the values in x_values which define an interval containing x. For instance, for x=3 with your example lists, the containing interval would be [x1,x2]=[2.5,3.4], and the indices would be i1=1, i2=2
  2. Calculate the slope on this interval by (y_values[i2]-y_values[i1])/(x_values[i2]-x_values[i1]) (ie dy/dx).
  3. The value at x is now the value at x1 plus the slope multiplied by the distance from x1.

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?

0
20

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. :)

2
  • 6
    I'd find using __call__ instead of __getitem__ to be preferrable in general, its usually an interpolation function.
    – Dave
    May 8, 2015 at 18:30
  • 1
    Does not work with Python3, as indexing a map is not supported anymore
    – Chris_128
    Jan 27, 2020 at 17:03
11
def interpolate(x1: float, x2: float, y1: float, y2: float, x: float):
    """Perform linear interpolation for x between (x1,y1) and (x2,y2) """

    return ((y2 - y1) * x + x2 * y1 - x1 * y2) / (x2 - x1)
1
  • 2
    What a shame. This is the lerp function everyone who searches for wants, in pure Python, even has a docstring—a perfect answer. And I had to scroll to see it. Nov 11, 2021 at 17:48
9

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
4

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 and y_list) your program's behavior may alert you to an issue for values greatly outside x_list. (Linear behavior goes bananas when given non-linear inputs!)

  • Returning the extents of the y_list for Interpolate[x] outside of x_list also means you know the range of your output value. If you extrapolate based on x much, much less than x_list[0] or x much, much greater than x_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])
    
1
  • 1
    I'd find using __call__ instead of __getitem__ to be preferrable in general, its usually an interpolation function.
    – Dave
    May 8, 2015 at 18:31
0

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])
1
  • 1
    I get an error.... TypeError 'Interpolate' object is not callable ?? What is the solution? Dec 10, 2019 at 22:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.