Как реализовать линейную интерполяцию?

Question

Я довольно новичок в программировании и думал, что попробую написать функцию линейной интерполяции.

Скажем, мне даны следующие данные:

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

Я хочу разработать функцию, которая будет интерполировать линейно от 1 до 2,5, от 2,5 до 3,4 и так далее, используя Python.

Я попытался просмотреть 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

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)

Answer 4

Я придумал довольно элегантное решение (ИМХО), поэтому не могу удержаться от публикации:

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 так что целочисленное деление (python <= 2.7) не пойдет и не испортит вещи, если x1, x2, y1 а также y2 все целые числа для некоторого iterval.

В __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 5

Опираясь на ответ Лаурица, вот версия со следующими изменениями

• Обновлен до python3 (карта вызывала у меня проблемы и не нужна)
• Исправлено поведение при краевых значениях
• Вызовите исключение, когда 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

Answer 6

Вместо того, чтобы экстраполировать с концов, вы можете вернуть экстенты 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 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])