Линейная интерполяция в Python

Вывод практических алгоритмов

Линейная интерполяция в Python - одна строка кода

Использование `numpy` для упрощения техники, которая должна быть в арсенале каждого

Введение

Линейная интерполяция создает непрерывную функцию из дискретных данных. Это фундаментальный строительный блок для алгоритма градиентного спуска, который используется при обучении практически всем методам машинного обучения. Глубокое обучение, при всей его сложности, не могло бы работать без него. Здесь мы рассмотрим теоретические основы построения алгоритма линейной интерполяции, а затем реализуем практический алгоритм и применим его к конкретному примеру.

Вывод

Если вы предпочитаете пропустить вывод, прокрутите вниз до раздела с заголовком The One-liner для кульминации. В противном случае продолжайте читать! Это не займет так много времени ...

В предыдущей статье мы обсуждали основы линейной регрессии с использованием некоторой линейной алгебры. В частности, у нас был такой результат для определения коэффициентов регрессии:

И у нас также было это уравнение, чтобы делать прогнозы с этими коэффициентами:

Если мы подставим первое уравнение во второе, мы получим прямую оценку предсказанных значений:

И все эти X термины образуют нечто, называемое матрицей проекции, также называемой матрицей шляпы:

Матрица проекции преобразует (или проецирует) наблюдаемые значения Y в прогнозируемые значения Ŷ. Для этого он должен содержать всю информацию, необходимую для прогнозирования любого произвольного значения Y.

Таким образом, кажется, что если мы знаем H, мы знаем, как предсказать любое значение! Чтобы проиллюстрировать это, давайте создадим игрушечный пример и посмотрим, как выглядит матрица шляп. Мы будем использовать две точки данных, которые лежат в координатах (0, 1) и (5, 2).

Разработка

Продемонстрируем на Python:

import numpy as np
np.set_printoptions(suppress=True)  # no scientific formatting
# np.r_ and np.c_ are short-hands for concatenation of arrays
# np.r_ performs along the first-axis (rows)
# np.c_ performs along the second-axis (columns)
# our x-values
x = np.c_[1., 1.,], [0., 5.]]
# our y-values
y = np.r_[1., 2.]
print(x)
array([[1., 0.],
       [1., 5.]])

Пока все довольно просто. Мы можем найти коэффициенты регрессии, используя первое уравнение:

# @ is the numpy operator for matrix multiplication
c = np.linalg.inv(x.T @ x) @ x.T @ Y
print(c)
array([1., 0.2])

И мы видим наклон 1/5 и точку пересечения 1, как и ожидалось. Матрица проекции следующая:

# @ is the numpy operator for matrix multiplication
h = x @ np.linalg.inv(x.T @ x) @ x.T
print(h)
array([[1., 0.],
       [0., 1.]])

Напомним, что когда мы умножаем H на Y, мы получаем прогнозы Ŷ . Поскольку у нас есть единицы на диагонали и нули в других местах, в этом случае H на самом деле является единичной матрицей, поэтому мы просто возвращаем прогнозы (1, 2) для фактического y -значения (1, 2). Это не тот случай, когда наши точки не образуют прямую линию, как почти в любой задаче линейной регрессии. Однако для интерполяции мы всегда можем сделать это предположение, поскольку есть только две точки, ограничивающие желаемую точку, в которой мы хотим выполнить интерполяцию!

Интерполяция

А как насчет фактической интерполяции? Если мы вспомним второе уравнение, мы можем предоставить любые значения для X, чтобы создать прогноз для значений. Давайте перепишем нашу матрицу проекции, используя этот факт, где X₀ представляет наши фактические значения, а X₁ представляет значение, при котором мы хотим интерполировать. Мы будем использовать W, чтобы представить это как матрицу весов, по причинам, которые мы вскоре увидим:

Что произойдет, если мы попробуем разные значения X₀ между значениями X₁, которые равны 0 и 5?

for i in np.linspace(0, 5, 6):
    w = np.c_[1., i] @ np.linalg.inv(x.T @ x) @ x.T
    y_hat = w @ y
    print(f'x0 = {i:.1f}, w = {w}, y_hat = {y_hat}')
x0 = 0.0, w = [[1.0 0.0]], y_hat = [1.0]
x0 = 1.0, w = [[0.8 0.2]], y_hat = [1.2]
x0 = 2.0, w = [[0.6 0.4]], y_hat = [1.4]
x0 = 3.0, w = [[0.4 0.6]], y_hat = [1.6]
x0 = 4.0, w = [[0.2 0.8]], y_hat = [1.8]
x0 = 5.0, w = [[0.0 1.0]], y_hat = [2.0]

Происходит кое-что интересное! В каждом случае сумма весовой матрицы равна единице. Когда мы умножаем эту весовую матрицу на значения Y, мы получаем средневзвешенное значение двух значений. Первый вес пропорционален расстоянию интерполированного X₀ между двумя значениями X₁, а второй вес - это просто дополнение к первому весу. Другими словами:

Однострочный

Это дает нам линейную интерполяцию в одну строку:

new_y = np.c_[1., new_x] @ np.linalg.inv(x.T @ x) @ x.T @ y

Конечно, это немного уловка. Мы должны точно знать два значения в исходном массиве значений x, между которыми находится наше новое интерполированное значение x. Нам нужна функция для определения индексов этих двух значений. К счастью, numpy содержит именно такую функцию: np.searchsorted. Давайте воспользуемся им, чтобы превратить линейную алгебру в векторизованную функцию. Вместо того, чтобы вычислять матрицу весов для каждого интерполированного значения, мы просто сохраним соответствующую информацию: сам вес и индексы значений X (которые будут также по индексам значений Y при вычислении интерполированных значений).

from typing import NamedTuple
class Coefficient(NamedTuple):
    lo: np.ndarray  # int  
    hi: np.ndarray  # int
    w: np.ndarray  # float
def interp_coef(x0: np.ndarray, x: np.ndarray) -> Coefficient:
    # find the indices into the original array
    hi = np.minimum(len(x0) - 1, np.searchsorted(x0, x, 'right'))
    lo = np.maximum(0, hi - 1)
    
    # calculate the distance within the range
    d_left = x - x0[lo]
    d_right = x0[hi] - x
    d_total = d_left + d_right
    # weights are the proportional distance
    w = d_right / d_total
    # correction if we're outside the range of the array
    w[np.isinf(w)] = 0.0
    # return the information contained by the projection matrices
    return Coefficient(lo, hi, w)
def interp_with(y0: np.ndarray, coef: Coefficient) -> np.ndarray:
    return coef.w * y0[coef.lo] + (1 - coef.w) * y0[coef.hi]

Этот алгоритм удобен тем, что нам нужно вычислить веса только один раз. Затем с помощью весов мы можем интерполировать любой вектор значений y, который совпадает с нашим вектором значений x.

Пример приложения

В инженерных задачах моделируемые нами физические процессы часто являются экспоненциальными или степенными функциями. Когда мы собираем данные для этих процессов, мы можем делать это через равные промежутки времени, разнесенные по логам. Многие методы интерпретации предназначены для работы с четными интервалами времени, разнесенными на логарифмы. Однако с некоторыми численными методами удобнее работать даже с линейно разнесенными значениями. Поэтому возможность трансформировать туда-сюда удобна. Линейная интерполяция - один из инструментов для выполнения этого преобразования.

Фактически, мы не ограничиваемся линейной интерполяцией. Мы можем сначала линеаризовать наши значения, а затем выполнить интерполяцию. Давайте посмотрим на пример расхода q и давления p в радиальной системе. . Эти данные были созданы синтетически с равномерно разнесенными по логарифмической шкале выборками, и когда мы наносим их на логарифмический график, мы видим отчетливое степенное поведение для большинства данных. В более поздние времена скорость потока становится экспоненциальной, но это не должно сильно влиять на нашу интерполяцию.

Давайте выполним интерполяцию, используя массив равномерно распределенных значений. Мы собираемся принять степенную функцию в качестве основы для нашей интерполяции, применив сначала логарифмическое преобразование к нашим значениям. Мы также обратим интерполяцию, чтобы убедиться, что она выполнена правильно:

# construct out linear time series
tDi = np.linspace(1e-3, 1e4, 100_001)
# interpolate from log to linear
coef = interp_coef(np.log(time), np.log(time_interp))
pi = np.exp(interp_with(np.log(pressure), coef))
qi = np.exp(interp_with(np.log(rate), coef))
# reverse from linear to log
coef_i = interp_coef(np.log(time_interp), np.log(time))
pr = np.exp(interp_with(np.log(pi), coef_i))
qr = np.exp(interp_with(np.log(qi), coef_i))

Резюме

Линейные методы - важный инструмент для манипулирования данными. Преобразование переменных, например, с использованием np.log, как показано в примере, позволяет нам применять эти методы в ситуациях, когда мы не ожидаем, что данные будут линейными.

Хотя мы разработали и реализовали наш собственный алгоритм, стандартный алгоритм действительно существует: scipy.interpolated.interp1d. Однако для использования в некоторых пакетах, например при использовании pymc3, может потребоваться реализовать что-то подобное для себя.

Существуют также более общие методы интерполяции. Один из них - B-сплайны, которые мы обсудим в следующей статье этой серии. B-сплайны также могут выполнять линейную интерполяцию, а также квадратичную, кубическую и т. Д. B-сплайны работают, расширяя концепции, которые мы здесь разработали. А именно, весовая вектор расширяется до матрицы весов. Вместо того, чтобы писать функцию для поиска индексов для каждого интерполированного значения, мы можем просто вычислить интерполированные значения y с матричным умножением коэффициентов и новых значений x.

Линейная интерполяция в Python - одна строка кода

Вывод практических алгоритмов