如何实现线性插值？

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之间线性插值，依此类推。

我试过通过The Python Tutorial，但我仍然无法理解它。

Answer 1

据我了解你的问题，你想写一些函数y = interpolate(x_values, y_values, x)，它会给你一些y的x值？然后基本思路遵循以下步骤：

1. 找到x_values中值的索引，它定义了包含x的区间。例如，对于带有示例列表的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. x的值现在是x1的值加上斜率乘以距离x1的距离。

如果x超出x_values的间隔，你还需要决定会发生什么，或者它是一个错误，或者你可以“向后”插值，假设斜率与第一个/最后一个间隔相同。

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

Answer 2

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

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

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都是某些iterval的整数，整数除法（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]

为简单起见，我根本没有处理边境条件。事实上，i[x]的x < 1将起作用，好像从（2.5,4）到（1,2）的线已经延伸到负无穷大，而i[x]或x == 1的x > 6将提升IndexError。更好的是在所有情况下都会引发IndexError，但这仍然是读者的练习。 :)

Answer 4

您可以返回y_list的范围，而不是对端点进行外推。大多数时候你的应用程序表现良好，而Interpolate[x]将在x_list中。 （推测）推断末端的线性影响可能会误导您相信您的数据表现良好。

• 返回非线性结果（以x_list和y_list的内容为界），您的程序的行为可能会提醒您在x_list之外的值很大的问题。 （当给出非线性输入时，线性行为会变成香蕉！）
• 在y_list之外返回Interpolate[x]的x_list范围也意味着你知道输出值的范围。如果基于x进行外推很多，远小于x_list[0]或x，远远大于x_list[-1]，那么您的返回结果可能超出了您预期的值范围。

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

Answer 5

您的解决方案在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])