高斯过程回归(Kriging)与径向基函数插值

Question

我已经为一个平面图上的传感器之间的温度数据实现了两种类型的插值。由于我对我使用的包的底层过程和数学不是很精通，我发现很难理解为什么他们通过pcolormesh的输出如此不同。

我使用了scipy.interpolate.Rbf和sklearn.gaussian_process。

这些是图片。

Gaussian process regression

Radial Basis Function interpolation

RBF示例看起来与web上的实现完全一样，但是GPR one显示的是这些长线而不是圆形。在Scikit learn的GPR实施中，什么参数将调节这些形状？为什么在探地雷达的温度变化甚至很小的情况下，它们的形状和颜色强度会如此不同？

平面图上的9个传感器(点)均匀分布。

RBF的代码。

# Set X and Y Coordinates for each sensor (pixels)
    days_data['xCoordinate'] = days_data.nodeid.apply(lambda id: createXCoord(id))
    days_data['yCoordinate'] = days_data.nodeid.apply(lambda id: createYCoord(id))

    # Define location of "sensors" on the axes
    x = days_data.xCoordinate.to_list()
    y = days_data.yCoordinate.to_list()
    z = days_data.avgtemperature.to_list() #temperature

    # Use Gaussian function
    rbf_adj = Rbf(x, y, z, function = 'gaussian')

    # Set extent to which colour mesh stretches over
    # the underlying image
    x_fine = np.linspace(0, 1000, 81) #81 - num of samples
    y_fine = np.linspace(0, 700, 81)

    x_grid, y_grid = np.meshgrid(x_fine, y_fine)

    z_grid = rbf_adj(x_grid.ravel(), y_grid.ravel()).reshape(x_grid.shape)

    # Remove the colorbar created by the previous plot, if any
    # To avoid a new colorbar being plotted alongside the previous one each time a different date is selected
    try:
        cb = p.colorbar
        cb.remove()
    except:
        pass

    # plot the pcolor on the Axes. Use alpha to set the transparency
    p=ax.pcolor(x_fine, y_fine, z_grid, alpha=0.3)
    ax.invert_yaxis() #invert Y axis for X and Y to have same starting point

    # Add a colorbar for the pcolor field
    fig.colorbar(p,ax=ax)

探地雷达代码

 # Define location of "sensors" on the axes
    x = days_data.xCoordinate.to_list()
    y = days_data.yCoordinate.to_list()
    z = days_data.avgtemperature.to_list() #temperature

    X = np.array([[a, b] for a, b in zip(x, y)])

    # Set extent to which colour mesh stretches over
    # the underlying image
    x_fine = np.linspace(0, 1000, 81) #81 - num of samples
    y_fine = np.linspace(0, 700, 82)
    X_fine = np.array([[a_fine, b_fine] for a_fine, b_fine in zip(x_fine, y_fine)])

    x_grid, y_grid = np.meshgrid(x_fine, y_fine)

    # Instantiate a Gaussian Process model
    kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))
    gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)

    # Fit to data using Maximum Likelihood Estimation of the parameters
    gp.fit(X, z)
    z_grid, sigma = gp.predict(X_fine, return_std=True)

    # Remove the colorbar created by the previous plot, if any
    # To avoid a new colorbar being plotted alongside the previous one each time a different date is selected
    try:
        cb = p.colorbar
        cb.remove()
    except:
        pass

    # plot the pcolor on the Axes. Use alpha to set the transparency
    p = ax.pcolor(x_grid, y_grid, np.meshgrid(z_grid, y_fine)[0], alpha=0.3)
    ax.invert_yaxis() #invert Y axis for X and Y to have same starting point

    # Add a colorbar for the pcolor field
    fig.colorbar(p,ax=ax)

Answer 1

我猜测高斯过程的尺度参数在x和y方向上有很大的不同，x尺度参数较小，y尺度参数相对较大。这样，x距离较小的两个点的相关性较低，而y距离较小的两个点的相关性较高:这在温度分布中创建了垂直的“带”。

我可以用下面的数据，用OpenTURNS来模拟这一过程。你没有提供温度，所以我只能从图表中猜测出来。

import numpy as np
import openturns as ot
coordinates = ot.Sample([[100.0,100.0],[500.0,100.0],[900.0,100.0], \
                         [100.0,350.0],[500.0,350.0],[900.0,350.0], \
                         [100.0,600.0],[500.0,600.0],[900.0,600.0]])
observations = ot.Sample([25.0,25.0,10.0,20.0,25.0,20.0,15.0,25.0,25.0],1)

# Extract coordinates
x = np.array(coordinates[:,0])
y = np.array(coordinates[:,1])

# Plot the data with a scatter plot and a color map
import matplotlib.pyplot as plt
fig = plt.figure()
plt.scatter(x, y, c=observations, cmap='viridis')
plt.colorbar()
plt.show()

这会产生以下结果：

​

​

可以使用以下脚本来拟合其中的克里金元模型。我使用了平方指数协方差模型。

def fitKriging(coordinates, observations, covarianceModel, basis):
    '''
    Fit the parameters of a kriging metamodel. 
    '''
    algo = ot.KrigingAlgorithm(coordinates, observations, covarianceModel, basis)
    algo.run()
    krigingResult = algo.getResult()
    krigingMetamodel = krigingResult.getMetaModel()
    return krigingResult, krigingMetamodel

inputDimension = 2
basis = ot.ConstantBasisFactory(inputDimension).build()
covarianceModel = ot.SquaredExponential([1.]*inputDimension, [1.0])
krigingResult, krigingMetamodel = fitKriging(coordinates, observations, covarianceModel, basis)

为了绘制来自这个克里金元模型的预测，我使用了以下基于pcolor函数的脚本。

def plotKrigingPredictions(krigingMetamodel):
    '''
    Plot the predictions of a kriging metamodel. 
    '''
    # Create the mesh of the box [0., 1000.] * [0., 700.]
    myInterval = ot.Interval([0., 0.], [1000., 700.])
    # Define the number of interval in each direction of the box
    nx = 20
    ny = 20
    myIndices = [nx-1, ny-1]
    myMesher = ot.IntervalMesher(myIndices)
    myMeshBox = myMesher.build(myInterval)

    # Predict
    vertices = myMeshBox.getVertices()
    predictions = krigingMetamodel(vertices)

    # Format for plot
    X = np.array(vertices[:,0]).reshape((ny,nx))
    Y = np.array(vertices[:,1]).reshape((ny,nx))
    predictions_array = np.array(predictions).reshape((ny,nx))

    # Plot
    plt.figure()
    plt.pcolor(X, Y, predictions_array)
    plt.colorbar()
    plt.show()
    return

plotKrigingPredictions(krigingMetamodel)

这会产生以下结果：

​

​

你可以在你的预测中看到相同的波段。

查看协方差模型解释了原因：

>>> covarianceModel = krigingResult.getCovarianceModel()
>>> print(covarianceModel)
SquaredExponential(scale=[0.167256,1.5929], amplitude=[6.75753])

X标度为0.1672，比y标度1.593小得多。这是因为我们使用了各向异性协方差模型，其中x尺度θ可以与y尺度θ不同。

为了解决这个问题，我们可以使用各向同性协方差模型，其中x尺度保持为y尺度。但是，更简单的方法是将x和y比例设置为给定值，并仅估计振幅参数sigma。

下面的脚本将比例设置为先前估计的平均值，并仅估计sigma。

scales = covarianceModel.getScale()
meanScale = (scales[0]+scales[1])/2.0
covarianceModel.setScale([meanScale]*2)
covarianceModel.setActiveParameter([2]) # Enable sigma (amplitude) only
# Learn amplitude only
krigingResult, krigingMetamodel = fitKriging(coordinates, observations, covarianceModel, basis)
covarianceModel = krigingResult.getCovarianceModel()
print("Covariance model=",covarianceModel)

这将打印：

Covariance model= SquaredExponential(scale=[0.880079,0.880079], amplitude=[6.1472])

最后，下面的脚本：

plotKrigingPredictions(krigingMetamodel)

图：

​

​

由于θ参数在两个方向上是相同的，所以温度现在是球形的，正如我猜你所期望的那样。