使用高斯过程回归 (GPR) 预测莫纳罗亚数据集上的二氧化碳水平#
本示例基于“机器学习中的高斯过程”[1]的第 5.4.3 节。它说明了使用对数边际似然上的梯度上升进行复杂核工程和超参数优化的示例。数据包括 1958 年至 2001 年间在夏威夷莫纳罗亚天文台收集的月平均大气二氧化碳浓度（以体积百万分率 (ppm) 为单位）。目标是将二氧化碳浓度建模为时间 \(t\) 的函数，并推断 2001 年以后的年份。
参考文献
print(__doc__)
# Authors: Jan Hendrik Metzen <[email protected]>
# Guillaume Lemaitre <[email protected]>
# License: BSD 3 clause
构建数据集#
我们将从收集空气样本的莫纳罗亚天文台获取一个数据集。我们有兴趣估计二氧化碳的浓度并将其外推到未来几年。首先，我们将 OpenML 中可用的原始数据集作为 pandas 数据帧加载。一旦 fetch_openml 添加了对 Polars 的原生支持，这将被 Polars 替换。
from sklearn.datasets import fetch_openml
co2 = fetch_openml(data_id=41187, as_frame=True)
co2.frame.head()
首先，我们处理原始数据帧以创建一个日期列，并将其与二氧化碳列一起选择。
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 月某些日期的二氧化碳浓度。我们可以绘制这些原始信息以更好地理解。
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。因此，总体噪声水平非常低，表明该模型可以很好地解释数据。
脚本总运行时间：（0 分 16.673 秒）
相关示例
