如何在 openturns 中获得更好的克里金结果图？

Question

我执行了球面克里金法，但我似乎无法获得好的输出图。 坐标（x 和 y）的范围从大约 51 纬度到大约 6.5 经度 我的观察范围从 -70 到 +10 这是我的代码：

import openturns as ot
import pandas as pd
# your input / output data can be easily formatted as samples for openturns


df = pd.read_csv("kreuzkerpenutm.csv")


inputdata = ot.Sample(df[['x','y']].values)
outputdata = ot.Sample(df[['z']].values)


dimension = 2  # dimension of your input (x,y)
basis = ot.ConstantBasisFactory(dimension).build()
covarianceModel = ot.SphericalModel(dimension)
    
algo = ot.KrigingAlgorithm(inputdata, outputdata, covarianceModel, basis)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()

lower = [-10.0] * 2 # lower bound of the 2D window
upper = [50.0] * 2 # upper bound of the 2D window
graph = metamodel.draw(lower, upper)
graph.setBoundingBox(ot.Interval(lower, upper))
graph.add(ot.Cloud(inputdata)) # overlay a scatter plot of the observation points
graph.setTitle("Kriging metamodel")

# A View object allows us to interact with the underlying matplotlib figure 
from openturns.viewer import View
view = View(graph, legend_kw={'bbox_to_anchor':(1,1), 'loc':"upper left"})
view.getFigure().tight_layout()

这是我的输出：

kriging metamodel graph

我不知道为什么我的图表不会显示我的输入以及我的克里金结果。

感谢您的想法和帮助

Answer 1

如果输入数据未按 [-1,1]^d 进行缩放，则克里金元模型可能无法使用最大似然优化来识别缩放参数。为了对此提供帮助，我们可以：

• 为协方差模型的尺度参数提供更好的起点（这是下面的技巧“A”），
• 设置优化算法的界限，以便搜索参数的区间与手头的数据相对应（这是下面的技巧“B”）。

这就是以下脚本的作用，使用模拟数据而不是 csv 数据文件。在脚本中，我使用 g 函数创建数据，该函数被缩放，以便它产生 [-10, 70] 范围内的结果，就像您的问题一样。请仔细查看设置协方差模型初始值的setScale() 方法：这是优化算法的起点。然后看setOptimizationBounds()方法，它设置了优化算法的界限。

import openturns as ot


dimension = 2  # dimension of your input (x,y)
distribution = ot.ComposedDistribution([ot.Uniform(-10.0, 50.0)] * dimension)
inputdata = distribution.getSample(100)

g = ot.SymbolicFunction(["x", "y"], ["30 + 3.0 * sin(x / 10.0) * (y / 10.0) ^ 2"])

outputdata = g(inputdata)

basis = ot.ConstantBasisFactory(dimension).build()
covarianceModel = ot.SphericalModel(dimension)
covarianceModel.setScale(inputdata.getMax())  # Trick A
algo = ot.KrigingAlgorithm(inputdata, outputdata, covarianceModel, basis)
# Trick B, v2
x_range = inputdata.getMax() - inputdata.getMin()
scale_max_factor = 2.0  # Must be > 1, tune this to match your problem
scale_min_factor = 0.1  # Must be < 1, tune this to match your problem
maximum_scale_bounds = scale_max_factor * x_range
minimum_scale_bounds = scale_min_factor * x_range
scaleOptimizationBounds = ot.Interval(minimum_scale_bounds, maximum_scale_bounds)
algo.setOptimizationBounds(scaleOptimizationBounds)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()
metamodel.setInputDescription(["x", "y"])
metamodel.setOutputDescription(["z"])

lower = [-10.0] * 2  # lower bound of the 2D window
upper = [50.0] * 2  # upper bound of the 2D window
graph = metamodel.draw(lower, upper)
graph.setBoundingBox(ot.Interval(lower, upper))
graph.add(ot.Cloud(inputdata))  # overlay a scatter plot of the observation points
graph.setTitle("Kriging metamodel")

# A View object allows us to interact with the underlying matplotlib figure
from openturns.viewer import View

view = View(graph, legend_kw={"bbox_to_anchor": (1, 1), "loc": "upper left"})
view.getFigure().tight_layout()

前面的脚本生成下图。

还有其他方法可以实现技巧 B。这是 J.Pelamatti 提供的一种方法：

# Trick B, v3
for d in range(X_train.getDimension()):
    dist = scipy.spatial.distance.pdist(X_train[:,d])
    scale_max_factor = 2.0  # Must be > 1, tune this to match your problem
    scale_min_factor = 0.1  # Must be < 1, tune this to match your problem
    maximum_scale_bounds = scale_max_factor * np.max(dist)
    minimum_scale_bounds = scale_min_factor * np.min(dist)

这个话题在this particular thread in OT's forum讨论。

Answer 2

抱歉回复晚了。

您使用的是哪个版本的openturns？ 可能您对（输入）数据进行了嵌入式转换，这使得数据范围大致介于 (-3, 3) 之间（标准缩放）。在这种情况下，克里金结果应该包含转换。 随着最近的openturns implementations，此功能已被删除。

希望这会有所帮助。 干杯