Примечание
Перейти в конец, чтобы скачать полный код примера или запустить этот пример в браузере с помощью JupyterLite или Binder.
Прогнозирование уровня CO2 на наборе данных Мауна-Лоа с помощью регрессии гауссовского процесса (GPR)
Этот пример основан на разделе 5.4.3 книги «Gaussian Processes for Machine Learning» [1]. В нём показан пример разработки сложного ядра и оптимизации гиперпараметров с помощью градиентного подъёма по логарифму маргинального правдоподобия. Данные представляют собой среднемесячные концентрации CO2 в атмосфере (в частях на миллион по объёму (ppm)), измеренные в обсерватории Мауна-Лоа на Гавайях в период с 1958 по 2001 год. Цель — смоделировать концентрацию CO2 как функцию времени \(t\) и выполнить экстраполяцию на годы после 2001-го.
Литература
# Authors: The scikit-learn developers # SPDX-License-Identifier: BSD-3-Clause
Создание набора данных
Мы получим набор данных из измерений проб воздуха, собранных в обсерватории Мауна-Лоа. Нас интересует оценка концентрации CO2 и её экстраполяция на последующие годы. Сначала загрузим исходный набор данных из OpenML в виде фрейма данных pandas. Позже мы заменим его на Polars, когда fetch_openml добавит встроенную поддержку этого формата.
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))
Как видно, данные о концентрации CO2 собраны в некоторые дни с марта 1958 года по декабрь 2001 года. Построим график исходных данных, чтобы лучше понять их.
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"
)
В этом примере мы будем предсказывать концентрацию CO2 как функцию даты. Также нас интересует экстраполяция на последующие годы после 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 с большим масштабом длины обеспечивает гладкость этого компонента. При этом возрастающий тренд не задаётся жёстко, чтобы оставить нашей модели некоторую свободу. Конкретные значения масштаба длины и амплитуды являются свободными гиперпараметрами.
from sklearn.gaussian_process.kernels import RBF long_term_trend_kernel = 50.0**2 * RBF(length_scale=50.0)
Сезонная изменчивость описывается периодическим ядром на основе экспоненты квадрата синуса с фиксированным периодом в 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,628 секунды)
Связанные примеры
© 2007–2026 The scikit-learn developers
Licensed under the 3-clause BSD License.
https://scikit-learn.org/1.9/auto_examples/gaussian_process/plot_gpr_co2.html
