具有高斯过程的多输出空间统计
Multi-output spatial statistics with gaussian processes
我最近一直在研究高斯过程。概率多输出的观点在我的领域很有前途。特别是空间统计。但是我遇到了三个问题：
- 多输出
- 过拟合和
- 各向异性。
让我 运行 使用 meuse 数据集（来自 R 包 sp）进行一个简单的案例研究。
UPDATE：本题使用的Jupyter notebook，根据, is 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 土壤中特别出现的最后一个问题：我如何在此类模型中指定 各向异性？
首先，您需要拆分数据。训练一个模型然后根据相同的训练数据进行预测看起来就像您观察到的过度拟合，但是您没有在任何保留数据上测试您的模型，因此您不知道它在野外的表现如何。尝试使用 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 操作，以便根据您的数据得出最佳设置。
所以回答你的问题：
该模型可以很好地处理多输出。
可以通过适当拆分数据或在不同的保留集上进行测试来解决过度拟合问题。
查看高斯过程指南的 Radial Basis Function RBF Kernel 部分，了解应用各向异性核而不是我们上面应用的各向同性核的一些见解。
更新评论中的问题
当你写 "the model handles multi output quite well as is" 时，你是说模型 "as is" 是为相关目标构建的，还是说模型作为独立模型的集合处理得很好?
好问题。根据我对 GaussianProcessRegressor 的了解，我不相信它能够在内部存储多个模型。所以这是一个单一的模型。话虽如此，您的问题的有趣之处在于陈述 "built for correlated targets"。在这种情况下，我们的两个目标似乎确实相当相关（皮尔逊相关系数 = 0.818，p=1.25e-38）所以我真的在这里看到两个问题：
对于相关数据，如果我们为两个目标以及单个目标建立模型，结果将如何比较？
对于不相关的数据，上述是否成立？
不幸的是，我们无法在不创建新的 "fake" 数据集的情况下测试第二个问题，这有点超出了我们在这里所做的范围。然而，我们可以很容易地回答第一个问题。使用我们相同的 train/test 拆分，我们可以训练两个具有相同超参数的新模型来分别预测铅和铜。然后我们可以使用 类 训练一个 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 稍微好一点。有趣的是，铜模型的解释方差明显低于铅模型（这实际上也对应于其他两个模型的各个组件的行为）。这一切都非常有趣，并且会引导我们沿着许多不同的开发路线走向我们的最终模型。
我认为这里的重要收获是所有模型迭代似乎都在同一个球场，并且在这种情况下没有明确的赢家。在这种情况下，您将需要进行一些重要的网格搜索，也许实施各向异性内核和任何其他领域特定知识会有所帮助，但我们的示例与有用的模型相去甚远。
我最近一直在研究高斯过程。概率多输出的观点在我的领域很有前途。特别是空间统计。但是我遇到了三个问题：
- 多输出
- 过拟合和
- 各向异性。
让我 运行 使用 meuse 数据集（来自 R 包 sp）进行一个简单的案例研究。
UPDATE：本题使用的Jupyter notebook，根据
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 土壤中特别出现的最后一个问题：我如何在此类模型中指定 各向异性？
首先，您需要拆分数据。训练一个模型然后根据相同的训练数据进行预测看起来就像您观察到的过度拟合，但是您没有在任何保留数据上测试您的模型，因此您不知道它在野外的表现如何。尝试使用 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 操作，以便根据您的数据得出最佳设置。
所以回答你的问题：
该模型可以很好地处理多输出。
可以通过适当拆分数据或在不同的保留集上进行测试来解决过度拟合问题。
查看高斯过程指南的 Radial Basis Function RBF Kernel 部分，了解应用各向异性核而不是我们上面应用的各向同性核的一些见解。
更新评论中的问题
当你写 "the model handles multi output quite well as is" 时，你是说模型 "as is" 是为相关目标构建的，还是说模型作为独立模型的集合处理得很好?
好问题。根据我对 GaussianProcessRegressor 的了解，我不相信它能够在内部存储多个模型。所以这是一个单一的模型。话虽如此，您的问题的有趣之处在于陈述 "built for correlated targets"。在这种情况下，我们的两个目标似乎确实相当相关（皮尔逊相关系数 = 0.818，p=1.25e-38）所以我真的在这里看到两个问题：
对于相关数据，如果我们为两个目标以及单个目标建立模型，结果将如何比较？
对于不相关的数据，上述是否成立？
不幸的是，我们无法在不创建新的 "fake" 数据集的情况下测试第二个问题，这有点超出了我们在这里所做的范围。然而，我们可以很容易地回答第一个问题。使用我们相同的 train/test 拆分，我们可以训练两个具有相同超参数的新模型来分别预测铅和铜。然后我们可以使用 类 训练一个 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 稍微好一点。有趣的是，铜模型的解释方差明显低于铅模型（这实际上也对应于其他两个模型的各个组件的行为）。这一切都非常有趣，并且会引导我们沿着许多不同的开发路线走向我们的最终模型。
我认为这里的重要收获是所有模型迭代似乎都在同一个球场，并且在这种情况下没有明确的赢家。在这种情况下，您将需要进行一些重要的网格搜索，也许实施各向异性内核和任何其他领域特定知识会有所帮助，但我们的示例与有用的模型相去甚远。