注意
跳转到末尾 下载完整的示例代码。或通过 JupyterLite 或 Binder 在浏览器中运行此示例
使用高斯过程回归 (GPR) 预测蒙纳罗亚数据集上的二氧化碳水平#
此示例基于《机器学习中的高斯过程》[1] 的第 5.4.3 节。它展示了一个复杂的核工程和超参数优化示例，使用对数-边际-似然上的梯度上升。数据包含 1958 年至 2001 年间在夏威夷蒙纳罗亚天文台收集的每月平均大气二氧化碳浓度（体积百万分率 (ppm)）。目标是建立二氧化碳浓度随时间 \(t\) 变化的模型，并对 2001 年之后的年份进行外推。
参考文献
print(__doc__)
# Authors: The scikit-learn developers
# 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 月某些日期的二氧化碳浓度。我们可以绘制这些原始信息以更好地理解。
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 分 4.369 秒)
相关示例
