如何实现线性插值？

Question

假设我有以下数据：

x = [1, 2.5, 3.4, 5.8, 6]
y = [2, 4, 5.8, 4.3, 4]

我希望设计一个函数，使用Python进行线性插值，从1到2.5，从2.5到3.4等等。

我已经尝试查阅这个Python教程，但是仍然无法理解。

- Helpless

---

这并不容易。你尝试过什么？ - zellio

你完全没有理解如何编程，或者如何在Python中实现算法吗？ - steabert

作为一个新手，我已经让自己置身于深水区了。我在考虑在算法中使用“for”或“if”语句，在多个x的范围之间进行操作。 - Helpless

http://en.wikipedia.org/wiki/Linear_least_squares_(mathematics) - warvariuc

在我看来，作为一个新手不去尝试并不是一个好的借口，事实上恰恰相反（即是一个很好的理由去尝试）。 - martineau

如果你只想在两个点之间进行插值，可以使用numpy：https://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html否则，下面的Dave的答案是最好的。 - wordsforthewise

Answer 1

import scipy.interpolate
y_interp = scipy.interpolate.interp1d(x, y)
print y_interp(5.0)

scipy.interpolate.interp1d通过进行线性插值，并可以自定义处理错误条件。

- Dave

---

问题询问如何实现该函数，而不是提供该函数的库。 - Miguel Bartelsman

这样的问题的答案是 - 你不需要重新发明轮子。有些东西你不需要重新发明，因为它们已经被证明是正确的。如果你想要深入了解它们，你应该去源头，比如数学书，并使用最好的工具来尝试它们。你不会在画笔店里寻求绘画课程。 - iantonuk

Answer 2

根据我理解你的问题，你想编写一个函数 y = interpolate(x_values, y_values, x)，它会在一些 x 值上给出相应的 y 值? 那么基本思路如下：

1. 找到包含x的值所定义区间的 x_values 索引。例如，在你的样例列表中，对于 x=3，包含它的区间是 [x1,x2]=[2.5,3.4]，索引将是 i1=1，i2=2。
2. 通过如下式子计算该区间的斜率：(y_values[i2]-y_values[i1])/(x_values[i2]-x_values[i1])（即 dy/dx）。
3. 现在，在 x1 值处的值加上与 x1 的距离乘以斜率即可得到 x 处的值。

此外，你还需要决定如果 x 在 x_values 的区间之外会发生什么情况，这可能是一个错误，或者你可以进行反向插值，假设斜率与第一个/最后一个区间相同。

这有帮助吗？还是你需要更具体的建议？

- carlpett

Answer 3

我想出了一个相当优雅的解决方案（在我看来），所以我忍不住要发帖：

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

我将映射为float，以便在某个区间中如果x1、x2、y1和y2都是整数，则整数除法（Python <= 2.7）不会触发并破坏计算。

在__getitem__中，我利用了self.x_list按升序排序的事实，通过使用bisect_left来（非常）快速地找到比x小的最大元素的索引在self.x_list中。

像这样使用该类：

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]

为了简单起见，我在这里没有处理边界条件。当前情况下，对于 x < 1 的情况，i[x] 会像从 (2.5, 4) 到 (1, 2) 的线段已经延伸至负无穷一样工作，而对于 x == 1 或 x > 6 的情况，i[x] 会引发一个 IndexError。更好的做法是在所有情况下都引发 IndexError，但这留给读者作为练习。 :)

- Lauritz V. Thaulow

---

一般情况下，我更喜欢使用 __call__ 而不是 __getitem__，因为它通常是一个插值函数。 - Dave

无法与Python3一起使用，因为不再支持对映射进行索引。 - Chris_128

发明一个轮子。需要的是 C++ std::lerp 的类似物。 - iantonuk

Answer 4

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)

- jlady

---

真遗憾。这就是每个搜索的人都想要的lerp函数，它是纯Python编写的，甚至还有一个文档字符串 - 完美的答案。而我却不得不滚动才能看到它。 - Ahmed Fasih

Answer 5

在Lauritz的回答基础上，以下是带有以下更改的版本：

• 更新为python3（map对我造成了问题，并且是不必要的）
• 修正边界值下的行为
• 当x越界时引发异常
• 使用__call__而不是__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])

使用示例：

>>> interp = Interpolate([1, 2.5, 3.4, 5.8, 6], [2, 4, 5.8, 4.3, 4])
>>> interp(4)
5.425

- Karl Bartel

Answer 6

不要对两端进行外推，你可以返回 y_list 的范围。大多数情况下，你的应用程序表现良好，Interpolate[x] 将在 x_list 中。在两端外推时产生的（可能是）线性影响可能会让你误以为数据表现良好。

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

- robertberrington

---

一般而言，我发现使用__call__比使用__getitem__更为可取，因为它通常是一个插值函数。 - Dave

Answer 7

你的解决方案在Python 2.7中无效。在检查x元素的顺序时出现了错误。我不得不修改代码才能使其正常工作：

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

- Leo

---

我遇到了一个错误… TypeError: 'Interpolate' 对象不可调用？有什么解决办法吗？ - Rudy Van Drie