绘制 GPR 二氧化碳图

简介

本实验是一个循序渐进的教程，介绍如何使用 scikit-learn 包对莫纳罗亚二氧化碳数据进行高斯过程回归（GPR）。数据包括 1958 年至 2001 年期间在夏威夷莫纳罗亚天文台收集的每月平均大气二氧化碳浓度。目标是将二氧化碳浓度建模为时间的函数，并将其外推到 2001 年以后的年份。

虚拟机提示

虚拟机启动完成后，点击左上角切换到“笔记本”标签，以访问 Jupyter Notebook 进行练习。

有时，你可能需要等待几秒钟让 Jupyter Notebook 完成加载。由于 Jupyter Notebook 的限制，操作验证无法自动化。

如果你在学习过程中遇到问题，随时向 Labby 提问。课程结束后提供反馈，我们将立即为你解决问题。

Skills Graph

构建数据集

第一步是通过从收集空气样本的莫纳罗亚天文台获取数据来构建数据集。我们感兴趣的是估计二氧化碳浓度并将其外推到未来几年。我们加载 OpenML 中可用的原始数据集，并通过计算每月平均值和删除未收集测量值的月份来预处理数据集。

from sklearn.datasets import fetch_openml
import pandas as pd

co2 = fetch_openml(data_id=41187, as_frame=True, parser="pandas")
co2_data = co2.frame
co2_data["date"] = pd.to_datetime(co2_data[["year", "month", "day"]])
co2_data = co2_data[["date", "co2"]].set_index("date")
co2_data = co2_data.resample("M").mean().dropna(axis="index", how="any")

X = (co2_data.index.year + co2_data.index.month / 12).to_numpy().reshape(-1, 1)
y = co2_data["co2"].to_numpy()

设计合适的核函数

为了设计与我们的高斯过程一起使用的核函数，我们可以对手头的数据做出一些假设。我们观察到它们有几个特征：我们看到一个长期上升趋势、一个明显的季节性变化以及一些较小的不规则性。我们可以使用不同的合适核函数来捕捉这些特征。

from sklearn.gaussian_process.kernels import RBF, ExpSineSquared, RationalQuadratic, WhiteKernel

long_term_trend_kernel = 50.0**2 * RBF(length_scale=50.0)
seasonal_kernel = (
    2.0**2
    * RBF(length_scale=100.0)
    * ExpSineSquared(length_scale=1.0, periodicity=1.0, periodicity_bounds="fixed")
)
irregularities_kernel = 0.5**2 * RationalQuadratic(length_scale=1.0, alpha=1.0)
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
)

模型拟合与外推

现在，我们准备好使用高斯过程回归器并拟合可用数据。为了遵循文献中的示例，我们将从目标值中减去均值。我们创建从 1958 年到当前月份的合成数据，并使用高斯过程对训练数据进行预测以检查拟合优度，对未来数据进行预测以查看模型的外推效果。

from sklearn.gaussian_process import GaussianProcessRegressor
import datetime
import numpy as np
import matplotlib.pyplot as plt

y_mean = y.mean()
gaussian_process = GaussianProcessRegressor(kernel=co2_kernel, normalize_y=False)
gaussian_process.fit(X, y - y_mean)

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"
)
plt.show()

核函数超参数的解读

现在，我们可以查看一下核函数的超参数。

gaussian_process.kernel_

因此，在减去均值后，大部分目标信号可以由一个约 45 ppm 的长期上升趋势和一个约 52 年的长度尺度来解释。周期性分量的幅度约为 2.6 ppm，衰减时间约为 90 年，长度尺度约为 1.5。较长的衰减时间表明我们有一个非常接近季节性周期的分量。相关噪声的幅度约为 0.2 ppm，长度尺度约为 0.12 年，白噪声贡献约为 0.04 ppm。因此，整体噪声水平非常小，这表明数据可以被模型很好地解释。

总结

本实验展示了如何使用 scikit-learn 包对莫纳罗亚的二氧化碳数据进行高斯过程回归（GPR）。数据包括 1958 年至 2001 年期间在夏威夷莫纳罗亚天文台收集的每月平均大气二氧化碳浓度。目标是将二氧化碳浓度建模为时间的函数，并将其外推到 2001 年之后的年份。我们构建了数据集，设计了合适的核函数，拟合了模型，并进行了外推。最后，我们解读了核函数的超参数。