具有高斯过程的多输出空间统计

Question

我最近一直在研究高斯过程。概率多输出的观点在我的领域很有前途。特别是空间统计。但是我遇到了三个问题：

1. 多输出
2. 过拟合和
3. 各向异性。

让我用meuse数据集（来自 R 包sp）运行一个简单的案例研究。

UPDATE：用于这个问题的Jupyter笔记本电脑，并根据更新的格儿的回答，是这里。

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline

meuse = pd.read_csv(filepath_or_buffer='https://gist.githubusercontent.com/essicolo/91a2666f7c5972a91bca763daecdc5ff/raw/056bda04114d55b793469b2ab0097ec01a6d66c6/meuse.csv', sep=',')

例如，我们将重点关注铜和铅。

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(121, aspect=1)
ax1.set_title('Lead')
ax1.scatter(x=meuse.x, y=meuse.y, s=meuse.lead, alpha=0.5, color='grey')

ax2 = fig.add_subplot(122, aspect=1)
ax2.set_title('Copper')
ax2.scatter(x=meuse.x, y=meuse.y, s=meuse.copper, alpha=0.5, color='orange')

事实上，铜和铅的浓度是相关的。

plt.plot(meuse['lead'], meuse['copper'], '.')
plt.xlabel('Lead')
plt.ylabel('Copper')

因此，这是一个多输出问题。

from sklearn.gaussian_process.kernels import RBF
from sklearn.gaussian_process import GaussianProcessRegressor as GPR
reg = GPR(kernel=RBF())
reg.fit(X=meuse[['x', 'y']], y=meuse[['lead', 'copper']])
predicted = reg.predict(meuse[['x', 'y']])

第一个问题：当 y 有多个维度时，内核是否是为相关的多输出构建的？如果没有，我该如何指定内核？

我继续分析以显示第二个问题，过度拟合：

fig = plt.figure(figsize=(12,4))
ax1 = fig.add_subplot(121)
ax1.set_title('Lead')
ax1.set_xlabel('Measured')
ax1.set_ylabel('Predicted')
ax1.plot(meuse.lead, predicted[:,0], '.')

ax2 = fig.add_subplot(122)
ax2.set_title('Copper')
ax2.set_xlabel('Measured')
ax2.set_ylabel('Predicted')
ax2.plot(meuse.copper, predicted[:,1], '.')

我创建了一个 x 和 y 坐标网格，该网格上的所有浓度都被预测为零。

最后，在土壤 3D 中特别出现的最后一个问题是：如何在此类模型中指定各向异性？

Answer 1

首先，您需要拆分数据。训练一个模型，然后在相同的训练数据上进行预测看起来就像您观察到的那样过度拟合，但您没有在任何保留数据上测试您的模型，因此您不知道它在野外的表现如何。尝试sklearn.model_selection.train_test_split像这样拆分您的数据：

X_train, X_test, y_train, y_test = train_test_split(meuse[['x', 'y']], meuse[['lead', 'copper']])

然后你可以训练你的模型。但是，您在那里也有问题。当你按照你的方式训练模型时，你最终会得到一个带有length_scale=1e-05. 基本上，您的模型中没有噪音。使用此设置进行的预测将非常集中在您的输入点 ( X_train)周围，您将无法对其周围的站点进行任何预测。您需要更改 的alpha参数GaussianProcessRegressor来解决此问题。由于默认值为 1e-10，因此您可能需要对其进行网格搜索。例如我使用alpha=0.1.

reg = GPR(RBF(), alpha=0.1)
reg.fit(X_train, y_train)
predicted = reg.predict(X_test)

fig = plt.figure(figsize=(12,4))
ax1 = fig.add_subplot(121)
ax1.set_title('Lead')
ax1.set_xlabel('Measured')
ax1.set_ylabel('Predicted')
ax1.plot(y_test.lead, predicted[:,0], '.')

ax2 = fig.add_subplot(122)
ax2.set_title('Copper')
ax2.set_xlabel('Measured')
ax2.set_ylabel('Predicted')
ax2.plot(y_test.copper, predicted[:,1], '.')

这导致了这个图表：

正如您所看到的，这里没有过拟合问题，实际上这可能是欠拟合的。就像我说的，你需要在这个模型上做一些 GridSearchCV 来根据你的数据提出最佳设置。

所以回答你的问题：

1. 该模型可以很好地处理多输出。

2. 可以通过正确拆分数据或在不同的保留集上进行测试来解决过度拟合问题。

3. 查看高斯过程指南的径向基函数 RBF 内核部分，了解如何应用各向异性内核而不是我们上面应用的各向同性内核。

更新评论中的问题

当您写“模型按原样很好地处理多输出”时，您是说模型“按原样”是为相关目标构建的，还是模型将它们作为独立模型的集合处理得很好？

好问题。根据我对 GaussianProcessRegressor 的了解，我不相信它能够在内部存储多个模型。所以这是一个单一的模型。话虽如此，关于您的问题的有趣之处在于“为相关目标构建”的陈述。在这种情况下，我们的两个目标似乎相当相关（Pearson 相关系数 = 0.818，p=1.25e-38）所以我在这里看到两个问题：

1. 对于相关数据，如果我们为两个目标以及单个目标建立模型，结果将如何比较？

2. 对于非相关数据，上述是否成立？

不幸的是，我们无法在不创建新的“假”数据集的情况下测试第二个问题，这有点超出了我们在这里所做的范围。然而，我们可以很容易地回答第一个问题。使用我们相同的训练/测试分割，我们可以训练两个具有相同超参数的新模型，用于分别预测铅和铜。然后我们可以MultiOutputRegressor使用这两个类来训练 a 。最后将它们与原始模型进行比较。像这样：

reg = GPR(RBF(), alpha=1)
reg.fit(X_train, y_train)
preds = reg.predict(X_test)
reg_lead = GPR(RBF(), alpha=1)
reg_lead.fit(X_train, y_train.lead)
lead_preds = reg_lead.predict(X_test)
reg_cop = GPR(RBF(), alpha=1)
reg_cop.fit(X_train, y_train.copper)
cop_preds = reg_cop.predict(X_test)
multi_reg = MultiOutputRegressor(GPR(RBF(), alpha=1))
multi_reg.fit(X_train, y_train)
multi_preds = multi_reg.predict(X_test)

现在我们有几个模型可以比较。让我们绘制预测图，看看我们得到了什么。

有趣的是，铅预测中没有明显差异，但铜预测中存在一些差异。而这些只存在于原始 GPR 模型和我们的其他模型之间。继续对误差进行更定量的测量，我们可以看到，对于解释方差，原始模型的性能比我们的 MultiOutputRegressor 略好。有趣的是，铜模型的解释方差明显低于铅模型（这实际上也对应于其他两个模型的各个组件的行为）。这一切都非常有趣，并且会引导我们沿着许多不同的开发路线进入最终模型。

我认为这里的重要收获是所有模型迭代似乎都在同一个范围内，并且在这种情况下没有明确的赢家。在这种情况下，您将需要进行一些重要的网格搜索，也许实现各向异性内核和任何其他领域特定的知识都会有所帮助，但正如我们的示例一样，它与有用的模型相去甚远。