Документация Engee
Notebook

Интерполяция и экстраполяция

Интерполяция – это прием, который позволяет дополнить набор данных недостающими значениями, рассчитанными в точках, находящихся между имеющимися точками в выборке. При помощи интерполяции можно заполнять пропуски в данных, сглаживать выбросы, прогнозировать поведение процесса, и так далее. Чтобы оценить значения за пределами области определения эксперимента используются методы экстраполяции.

В отличие от разных видов регрессии или полиномиальной аппроксимации, для интерполяции желательно, чтобы исходные данные располагались на регулярной сетке. То есть аргументы функции должны быть заданы через тип данных Range. Функцию $y = F(x)$, осуществляющую интерполяцию выходных данных y заданных на сетке x, мы будем называть интерполирующей функцией или интерполянтом.

Многомерная интерполяция на сетке

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

In [ ]:
using Interpolations, Plots, Random
Random.seed!(12);
gr( fmt=:png, size=(500,400) );

Предположим, что наши эмпирические данные образуют матрицу 4 на 4:

In [ ]:
xs = range( 1, 10, length=4 );
ys = range( 1, 10, length=4 );
zs = rand( length(xs), length(ys) );

Создав нужный нам интерполянт (например, через кубический сплайн), мы можем вычислить его значение в любой точке в пределах диапазона входных данных:

In [ ]:
itp = cubic_spline_interpolation( (xs,ys), zs );
itp( 2, 2 )
Out[0]:
0.7444301350128467

Эта модель нам позволит рассчитать выходные значения в любой точке внутри пространства аргументов, ограниченного периметром, заданным сеткой, на которой определены эмпирические данные.

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

In [ ]:
scatter( repeat(xs, inner=length(ys)),
         repeat(ys, outer=length(xs)),
         vec(reshape( zs, (1,:))), markercolor=:red, label="Данные" )

xs_new = range(1, 10, length=30)
ys_new = range(1, 10, length=30)
surface!( xs_new, ys_new, [itp(x,y) for x in xs_new, y in ys_new],
          colorbar = false, alpha = 0.5, c=:viridis )
Out[0]:
No description has been provided for this image

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

Таким же образом можно выполнить кусочно-линейную интерполяцию (linear_interpolation) или интерполяцию по "методу ближайших соседей" (constant_interpolation). Точно такие же действия можно реализовать при помощи команд ConstantInterpolation, LinearInterpolation и CubicSplineInterpolation. Согласно документации, они оставлены в библиотеке ради удобства.

In [ ]:
itp0 = constant_interpolation( (xs,ys), zs ); # Константа
itp1 = linear_interpolation( (xs,ys), zs ); # Линейная интерполяция

plot(
    (scatter( repeat(xs, inner=length(ys)), repeat(ys, outer=length(xs)), vec(reshape( zs, (1,:))), markercolor=:red );
     surface!( xs_new, ys_new, itp0(xs_new, ys_new), colorbar = false, alpha = 0.5, c=:viridis, legend=:none )
     ),
    (scatter( repeat(xs, inner=length(ys)), repeat(ys, outer=length(xs)), vec(reshape( zs, (1,:))), markercolor=:red );
     surface!( xs_new, ys_new, itp1(xs_new, ys_new), colorbar = false, alpha = 0.5, c=:viridis, legend=:none )
     ),
    size=(1000,400)
    )
Out[0]:
No description has been provided for this image

Вид функции интерполяции определяется степенью полинома, задающего сплайн: Constant, Linear, Quadratic, или Cubic. До сих пор мы пользовались удобным, упрощенным синтаксисом этих команд. Но есть и более сложный синтаксис:

In [ ]:
# Вычисляем интерполяцию квадратичным сплайном,
# указав граничное условие – отражение (Reflect)
# и интерпретируя точки на периферии как центроиды (OnCell)
# а также дополнительно масштабируем полученный интерполянт, который при таком синтаксисе поначалу имеет единичный диапазон
sitp4 = scale( interpolate(zs, BSpline(Quadratic(Reflect(OnCell())))), (xs, ys,))

# Линейная интерполяция по одной оси, ступенчатая – по другой оси
# аргумент Gridded означает, что точки на периметре интерпретируются как OnGrid()
itp5 = interpolate( (xs, ys), zs, (Gridded(Linear()),Gridded(Constant())) )

plot(
    (scatter( repeat(xs, inner=length(ys)), repeat(ys, outer=length(xs)), vec(reshape( zs, (1,:))), markercolor=:red );
     surface!( xs_new, ys_new, sitp4(xs_new, ys_new), colorbar = false, alpha = 0.5, c=:viridis, legend=:none )
    ),
    
    (scatter( repeat(xs, inner=length(ys)), repeat(ys, outer=length(xs)), vec(reshape( zs, (1,:))), markercolor=:red );
     surface!( xs_new, ys_new, itp5(xs_new, ys_new), colorbar = false, alpha = 0.5, c=:viridis, legend=:none )
    ),
    size=(1000,400)
)
Out[0]:
No description has been provided for this image

Не все комбинации условий могут дать результат. К примеру, квадратичная и кубическая интерполяция требуют указания граничных условий из перечня Flat, Line (то же, что и Natural), Free, Periodic или Reflect. Точки на периметре сетки могут быть интерпретированы алгоритмом интерполяции как граничные (OnGrid()) или как располагающиеся на полпути между предшествующим рядом и фиктивным следующим рядом точек (OnCell()).

Интерполяцию данных, заданных на нерегурярной сетке, можно выполнять при помощи кусочно-линейного сплайна и при помощи констант.

Информацию о дополнительных параметрах и режимах интерполяции можно почерпнуть в документации.

Экстраполяция

При необходимости прогнозировать данные за пределами диапазона имеющихся экспериментальных данных к параметрам интерполяции нужно добавить еще одну гипотезу – как, по нашему мнению, ведет себя функция за пределами имеющихся данных, и тут у нас есть варианты Throw (выдать ошибку), Flat, Line, Periodic и Reflect, либо можно передать в качестве аргумента константу, значение которой будет подставляться при интерполяции.

Продемонстрируем разные методы экстраполяции на примере одномерного сигнала.

In [ ]:
x = 1:5;
y = [2,7,3,8,6];

Зададим разные граничные условия (boundary conditions) через аргумент extrapolation_bc.

In [ ]:
new_x = range(-2, 8, length=100)

plot(
    plot(x, y, title="Исходные данные"),
    plot(new_x, cubic_spline_interpolation(x, y, extrapolation_bc=Line())(new_x), title="Line()"),
    plot(new_x, cubic_spline_interpolation(x, y, extrapolation_bc=Flat())(new_x), title="Flat()"),
    plot(new_x, cubic_spline_interpolation(x, y, extrapolation_bc=Periodic())(new_x), title="Periodic()"),
    plot(new_x, cubic_spline_interpolation(x, y, extrapolation_bc=Reflect())(new_x), title="Reflect()"),
    plot(new_x, cubic_spline_interpolation(x, y, extrapolation_bc=5)(new_x), title="fill(5)"),
    size=(1000,400), layout=(2,3)
)
Out[0]:
No description has been provided for this image

Разумеется, эти же методы позволяют осуществлять интерполяцию и экстраполяцию данных более высокой размерности. Например, разместим в 3-мерном пространстве 64 точки данных и построим гладкую функцию на их основе.

In [ ]:
a = range( 1, 10, length=4 );
b = range( 1, 10, length=4 );
c = range( 1, 10, length=4 );
d = rand( length(a), length(b), length(c) );

# 3D-интерполяция (легко записать в одну строку)
itp = extrapolate(
          scale(
              interpolate( d, BSpline(
                                  Quadratic(
                                      Reflect(OnCell())
                                  )
                              )
                ), (a,b,c,)
      ), Flat() )

# Задаем новую сетку
new_a = new_b = -2:0.1:13; new_c = 1:10;

vPlots = [ contour( new_a, new_b, itp(new_a, new_b, ic), colorbar = false, alpha = 0.5, c=:viridis, legend=:none, title="c = $ic", titlefont = font(8)) for ic in new_c ] 
plot!( vPlots..., layout=(2,:), size=(1000,400) )
Out[0]:
No description has been provided for this image

Заключение

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

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