Note
Go to the end to download the full example code. or to run this example in your browser via Binder
使用高斯过程回归（GPR）对莫纳罗亚数据集的CO2水平进行预测#
本示例基于《高斯过程机器学习》第5.4.3节[1]。它展示了复杂核工程和使用梯度上升法对对数边际似然进行超参数优化的示例。数据包括1958年至2001年在夏威夷莫纳罗亚天文台收集的每月平均大气CO2浓度（以体积百万分之一（ppm）为单位）。目标是将CO2浓度建模为时间 \(t\) 的函数，并对2001年以后的年份进行外推。
参考文献
print(__doc__)
# 作者：scikit-learn 开发者
# SPDX-License-Identifier：BSD-3-Clause
构建数据集#
我们将从毛纳罗亚天文台获取一个收集空气样本的数据集。我们感兴趣的是估算二氧化碳的浓度并将其外推到未来几年。首先，我们加载在OpenML上可用的原始数据集，并将其作为pandas数据框。这将在 fetch_openml 添加对Polars的原生支持后被替换。
from sklearn.datasets import fetch_openml
co2 = fetch_openml(data_id=41187, as_frame=True)
co2.frame.head()
首先，我们处理原始数据框以创建一个日期列，并选择它和CO2列。
import polars as pl
co2_data = pl.DataFrame(co2.frame[["year", "month", "day", "co2"]]).select(
pl.date("year", "month", "day"), "co2"
)
co2_data.head()
co2_data["date"].min(), co2_data["date"].max()
(datetime.date(1958, 3, 29), datetime.date(2001, 12, 29))
我们可以看到，从1958年3月到2001年12月的某些天，我们得到了CO2浓度数据。我们可以绘制这些原始信息，以便更好地理解。
import matplotlib.pyplot as plt
plt.plot(co2_data["date"], co2_data["co2"])
plt.xlabel("date")
plt.ylabel("CO$_2$ concentration (ppm)")
_ = plt.title("Raw air samples measurements from the Mauna Loa Observatory")
我们将通过取每月平均值来预处理数据集，并删除未收集测量值的月份。这样的处理将对数据产生平滑效果。
co2_data = (
co2_data.sort(by="date")
.group_by_dynamic("date", every="1mo")
.agg(pl.col("co2").mean())
.drop_nulls()
)
plt.plot(co2_data["date"], co2_data["co2"])
plt.xlabel("date")
plt.ylabel("Monthly average of CO$_2$ concentration (ppm)")
_ = plt.title(
"Monthly average of air samples measurements\nfrom the Mauna Loa Observatory"
)
本示例的目的是根据日期预测二氧化碳浓度。我们还对2001年之后的未来年份进行外推感兴趣。
作为第一步，我们将划分数据和目标进行估计。由于数据是日期，我们会将其转换为数值。
X = co2_data.select(
pl.col("date").dt.year() + pl.col("date").dt.month() / 12
).to_numpy()
y = co2_data["co2"].to_numpy()
设计合适的内核#
为了设计用于高斯过程的核函数，我们可以对手头的数据做出一些假设。我们观察到它们具有以下几个特征：长期上升趋势、明显的季节性变化和一些较小的不规则性。我们可以使用不同的合适核函数来捕捉这些特征。
首先，可以使用具有大尺度参数的径向基函数（RBF）核来拟合长期上升趋势。具有大尺度的RBF核会使这一部分变得平滑。没有强制要求趋势性增加，以便为我们的模型提供一定的自由度。具体的尺度和幅度是自由超参数。
季节性变化由具有固定周期为1年的周期指数正弦平方核解释。该周期成分的长度尺度控制其平滑度，是一个自由参数。为了允许偏离精确周期性的衰减，采用了与RBF核的乘积。该RBF成分的长度尺度控制衰减时间，也是另一个自由参数。这种类型的核也被称为局部周期核。
from sklearn.gaussian_process.kernels import ExpSineSquared
seasonal_kernel = (
2.0**2
* RBF(length_scale=100.0)
* ExpSineSquared(length_scale=1.0, periodicity=1.0, periodicity_bounds="fixed")
)
小的不规则性可以通过有理二次核成分来解释，其长度尺度和量化长度尺度扩散性的 alpha 参数需要确定。有理二次核等价于具有多个长度尺度的 RBF 核，并且能够更好地适应不同的不规则性。
from sklearn.gaussian_process.kernels import RationalQuadratic
irregularities_kernel = 0.5**2 * RationalQuadratic(length_scale=1.0, alpha=1.0)
最终，数据集中的噪声可以通过一个包含RBF核贡献的核来解释，该核可以解释局部天气现象等相关噪声成分，以及一个用于白噪声的白核贡献。相对幅度和RBF的长度尺度是进一步的自由参数。
from sklearn.gaussian_process.kernels import WhiteKernel
noise_kernel = 0.1**2 * RBF(length_scale=0.1) + WhiteKernel(
noise_level=0.1**2, noise_level_bounds=(1e-5, 1e5)
)
因此，我们的最终核是所有先前核的总和。
co2_kernel = (
long_term_trend_kernel + seasonal_kernel + irregularities_kernel + noise_kernel
)
co2_kernel
50**2 * RBF(length_scale=50) + 2**2 * RBF(length_scale=100) * ExpSineSquared(length_scale=1, periodicity=1) + 0.5**2 * RationalQuadratic(alpha=1, length_scale=1) + 0.1**2 * RBF(length_scale=0.1) + WhiteKernel(noise_level=0.01)
模型拟合与外推#
现在，我们准备使用高斯过程回归器并拟合可用数据。为了遵循文献中的示例，我们将从目标值中减去均值。我们本可以使用 normalize_y=True 。然而，这样做还会缩放目标值（将 y 除以其标准差）。因此，不同核的超参数将具有不同的含义，因为它们将不再以百万分率（ppm）表示。
from sklearn.gaussian_process import GaussianProcessRegressor
y_mean = y.mean()
gaussian_process = GaussianProcessRegressor(kernel=co2_kernel, normalize_y=False)
gaussian_process.fit(X, y - y_mean)
现在，我们将使用高斯过程进行预测：
训练数据用于检查拟合优度；
未来数据用于查看模型的外推结果。
因此，我们从1958年到当前月份创建合成数据。此外，我们还需要添加在训练期间计算的减去的均值。
import datetime
import numpy as np
today = datetime.datetime.now()
current_month = today.year + today.month / 12
X_test = np.linspace(start=1958, stop=current_month, num=1_000).reshape(-1, 1)
mean_y_pred, std_y_pred = gaussian_process.predict(X_test, return_std=True)
mean_y_pred += y_mean
plt.plot(X, y, color="black", linestyle="dashed", label="Measurements")
plt.plot(X_test, mean_y_pred, color="tab:blue", alpha=0.4, label="Gaussian process")
plt.fill_between(
X_test.ravel(),
mean_y_pred - std_y_pred,
mean_y_pred + std_y_pred,
color="tab:blue",
alpha=0.2,
)
plt.legend()
plt.xlabel("Year")
plt.ylabel("Monthly average of CO$_2$ concentration (ppm)")
_ = plt.title(
"Monthly average of air samples measurements\nfrom the Mauna Loa Observatory"
)
我们的拟合模型能够很好地拟合先前的数据，并且有信心外推到未来的年份。
核超参数的解释#
现在，我们可以看看核函数的超参数。
gaussian_process.kernel_
44.8**2 * RBF(length_scale=51.6) + 2.64**2 * RBF(length_scale=91.5) * ExpSineSquared(length_scale=1.48, periodicity=1) + 0.536**2 * RationalQuadratic(alpha=2.89, length_scale=0.968) + 0.188**2 * RBF(length_scale=0.122) + WhiteKernel(noise_level=0.0367)
因此，大部分目标信号（减去均值后）可以通过一个约45 ppm的长期上升趋势和约52年的长度尺度来解释。周期性成分的振幅约为2.6 ppm，衰减时间约为90年，长度尺度约为1.5。长衰减时间表明我们有一个非常接近季节性周期的成分。相关噪声的振幅约为0.2 ppm，长度尺度约为0.12年，白噪声贡献约为0.04 ppm。因此，总体噪声水平非常小，表明数据可以很好地用模型解释。
Total running time of the script: (0 minutes 2.635 seconds)
Related examples
