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

Question

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

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

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

更新：用于此问题并根据Grr's answer更新的Jupyter笔记本是here。

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)，您将无法对它们周围的站点做出任何预测。您需要更改GaussianProcessRegressor 的alpha 参数来解决此问题。这是您可能需要进行网格搜索的内容，因为默认值为 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. 查看高斯过程指南的 Radial Basis Function RBF Kernel 部分，了解有关应用各向异性内核而不是我们上面应用的各向同性内核的一些见解。

评论问题更新

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

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

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

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

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

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 稍微好一点。有趣的是，铜模型的解释方差明显低于铅模型（这实际上也对应于其他两个模型的各个组件的行为）。这一切都非常有趣，并且会引导我们走上许多不同的开发路线，最终得到我们的最终模型。

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