Engee 文档
Notebook

内插法和外推法

内插法是一种技术,可以用样本中现有点之间的点计算出的缺失值来补充数据集。内插法可用于填补数据空白、平滑异常值、预测过程行为等。要估算实验定义域之外的值,需要使用外推法

与各种回归多项式逼近不同的是,内插法的原始数据最好是在规则的网格上。也就是说,必须通过数据类型Range 来指定函数参数。对定义在网格x 上的输出数据y 进行插值的函数$y = F(x)$ ,我们称之为插值函数或插值器。

网格上的多变量插值

只要有足够多的数据点,就可以在任何中间点(包括规则网格或随机抽样点)上进行预测。

In [ ]:
Pkg.add(["Interpolations"])
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)。使用ConstantInterpolationLinearInterpolationCubicSplineInterpolation 命令可以实现完全相同的操作。根据文档,为了方便起见,这些命令被保留在库中。

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,PeriodicReflect 的边界条件。插值算法可将网格周长上的点解释为边界点 (OnGrid()) 或解释为前一行点与假下一行点的中间点 (OnCell()) 。

使用片断线性样条线和常数可以对不规则网格上的数据进行插值。

有关其他参数和插值模式的信息,请参见 文档

外推法

当我们需要预测可用实验数据范围之外的数据时,我们需要在插值参数中添加另一个假设--我们认为函数在可用数据之外的表现如何,这里我们有Throw给出错误)、FlatLinePeriodicReflect 选项,或者我们可以传递一个常数作为参数,其值将在插值过程中被替换。

让我们以一维信号为例演示不同的外推法。

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

让我们通过参数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

当然,同样的方法也可以对高维数据进行插值和外推。例如,让我们在三维空间中放置 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 生态系统中还有几十种软件包,可以在其他问题设置下对数据进行插值处理,它们都列在 本页 上。