在OpenTURNS中,KrigingAlgorithm
该类可以根据特定输入点的已知输出值估计高斯过程模型的超参数。然后,的getMetamodel
方法KrigingAlgorithm
返回一个对数据进行插值的函数。
首先,我们需要将 Numpy 数组转换coordinates
为observations
OpenTURNSSample
对象:
import openturns as ot
input_train = ot.Sample(coordinates)
output_train = ot.Sample(observations, 1)
数组coordinates
有 shape (6, 2)
,所以它变成了Sample
大小为 6 和维度 2 的 a。数组observations
有 shape (6,)
,这是不明确的:它是Sample
大小为 6 和维度 1 的 a,还是Sample
大小为 1 和维度 6 的 a?Sample
为了澄清这一点,我们在对类构造函数的调用中指定维度 (1) 。
下面,我们定义一个具有恒定趋势函数和平方指数协方差核的高斯过程模型:
inputDimension = 2
basis = ot.ConstantBasisFactory(inputDimension).build()
covariance_kernel = ot.SquaredExponential([1.0]*inputDimension, [1.0])
algo = ot.KrigingAlgorithm(input_train, output_train,
covariance_kernel, basis)
然后我们拟合趋势的值和协方差核的参数(幅度参数和尺度参数),得到一个元模型:
# Fit
algo.run()
result = algo.getResult()
krigingMetamodel = result.getMetaModel()
结果krigingMetamodel
是 a Function
,它以 2DPoint
作为输入并返回 1D Point
。它预测感兴趣的数量。为了说明这一点,让我们构建二维域 [0,1]×[0,1] 并用规则网格对其进行离散化:
# Create the 2D domain
myInterval = ot.Interval([0.0, 0.0], [1.0, 1.0])
# Define the number of interval in each direction of the box
nx = 20
ny = 10
myIndices = [nx - 1, ny - 1]
myMesher = ot.IntervalMesher(myIndices)
myMeshBox = myMesher.build(myInterval)
使用我们krigingMetamodel
来预测该网格上感兴趣的数量所取的值可以通过以下语句完成。我们首先vertices
将网格的 设为Sample
,然后通过对元模型的一次调用来评估预测(这里不需要for
循环):
# Predict
vertices = myMeshBox.getVertices()
predictions = krigingMetamodel(vertices)
为了使用 Matplotlib 查看结果,我们首先必须创建pcolor
函数所需的数据:
# 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
import matplotlib.pyplot as plt
fig = plt.figure()
plt.pcolor(X, Y, predictions_array)
plt.colorbar()
plt.show()
我们看到元模型的预测等于观察到的输入点的观察。
这个元模型是坐标的平滑函数:它的平滑度随着协方差核平滑度的增加而增加,平方指数协方差核恰好是平滑的。