Способность регрессии на основе гауссовских процессов (GPR) оценивать уровень шума в данных

Генерация данных

Мы будем работать в условиях, где X содержит один признак. Создадим функцию, которая генерирует целевую переменную для предсказания. Добавим возможность зашумить сгенерированную целевую переменную.

import numpy as np


def target_generator(X, add_noise=False):
    target = 0.5 + np.sin(3 * X)
    if add_noise:
        rng = np.random.RandomState(1)
        target += rng.normal(0, 0.3, size=target.shape)
    return target.squeeze()

Посмотрим на генератор целевой переменной, не добавляя шум, чтобы увидеть сигнал, который мы хотим предсказать.

X = np.linspace(0, 5, num=80).reshape(-1, 1)
y = target_generator(X, add_noise=False)

import matplotlib.pyplot as plt

plt.plot(X, y, label="Expected signal")
plt.legend()
plt.xlabel("X")
_ = plt.ylabel("y")

Целевая переменная получается преобразованием входных данных X с помощью синусоидальной функции. Теперь сгенерируем несколько зашумлённых обучающих выборок. Чтобы проиллюстрировать уровень шума, построим истинный сигнал вместе с зашумлёнными обучающими выборками.

rng = np.random.RandomState(0)
X_train = rng.uniform(0, 5, size=20).reshape(-1, 1)
y_train = target_generator(X_train, add_noise=True)

plt.plot(X, y, label="Expected signal")
plt.scatter(
    x=X_train[:, 0],
    y=y_train,
    color="black",
    alpha=0.4,
    label="Observations",
)
plt.legend()
plt.xlabel("X")
_ = plt.ylabel("y")

Оптимизация гиперпараметров ядра в GPR

Теперь создадим GaussianProcessRegressor с аддитивным ядром, состоящим из ядер RBF и WhiteKernel. Ядро WhiteKernel позволит оценить уровень шума в данных, а ядро RBF будет использоваться для моделирования нелинейной зависимости между данными и целевой переменной.

Однако мы покажем, что пространство гиперпараметров содержит несколько локальных минимумов. Это подчеркнёт важность начальных значений гиперпараметров.

Создадим модель с ядром, которому задан высокий уровень шума и большой масштаб длины. Такая модель будет объяснять все вариации в данных шумом.

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

kernel = 1.0 * RBF(length_scale=1e1, length_scale_bounds=(1e-2, 1e3)) + WhiteKernel(
    noise_level=1, noise_level_bounds=(1e-10, 1e1)
)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0)
gpr.fit(X_train, y_train)
y_mean, y_std = gpr.predict(X, return_std=True)

/home/circleci/project/sklearn/gaussian_process/kernels.py:455: ConvergenceWarning: The optimal value found for dimension 0 of parameter k1__k2__length_scale is close to the specified upper bound 1000.0. Increasing the bound and calling fit again may find a better value.
  warnings.warn(

plt.plot(X, y, label="Expected signal")
plt.scatter(x=X_train[:, 0], y=y_train, color="black", alpha=0.4, label="Observations")
plt.errorbar(X, y_mean, y_std, label="Posterior mean ± std")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
_ = plt.title(
    (
        f"Initial: {kernel}\nOptimum: {gpr.kernel_}\nLog-Marginal-Likelihood: "
        f"{gpr.log_marginal_likelihood(gpr.kernel_.theta)}"
    ),
    fontsize=8,
)

Мы видим, что найденное оптимальное ядро по-прежнему имеет высокий уровень шума и ещё больший масштаб длины. Масштаб длины достигает максимальной границы, заданной для этого параметра, в результате чего появляется предупреждение.

Что ещё важнее, модель не даёт полезных предсказаний: среднее предсказание выглядит постоянным и не следует ожидаемому сигналу без шума.

Теперь инициализируем ядро RBF большим начальным значением length_scale, а ядро WhiteKernel — меньшим начальным уровнем шума, оставив границы параметров без изменений.

kernel = 1.0 * RBF(length_scale=1e-1, length_scale_bounds=(1e-2, 1e3)) + WhiteKernel(
    noise_level=1e-2, noise_level_bounds=(1e-10, 1e1)
)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0)
gpr.fit(X_train, y_train)
y_mean, y_std = gpr.predict(X, return_std=True)

plt.plot(X, y, label="Expected signal")
plt.scatter(x=X_train[:, 0], y=y_train, color="black", alpha=0.4, label="Observations")
plt.errorbar(X, y_mean, y_std, label="Posterior mean ± std")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
_ = plt.title(
    (
        f"Initial: {kernel}\nOptimum: {gpr.kernel_}\nLog-Marginal-Likelihood: "
        f"{gpr.log_marginal_likelihood(gpr.kernel_.theta)}"
    ),
    fontsize=8,
)

Сначала отметим, что предсказания новой модели точнее, чем предсказания предыдущей: новая модель способна оценить функциональную зависимость без шума.

Рассмотрев гиперпараметры ядра, мы видим, что найденная наилучшая комбинация имеет меньший уровень шума и меньший масштаб длины, чем у первой модели.

Чтобы получить представление о локальных минимумах, можно изучить отрицательное логарифмическое маргинальное правдоподобие (LML) GaussianProcessRegressor при разных гиперпараметрах.

from matplotlib.colors import LogNorm

length_scale = np.logspace(-2, 4, num=80)
noise_level = np.logspace(-2, 1, num=80)
length_scale_grid, noise_level_grid = np.meshgrid(length_scale, noise_level)

log_marginal_likelihood = [
    gpr.log_marginal_likelihood(theta=np.log([0.36, scale, noise]))
    for scale, noise in zip(length_scale_grid.ravel(), noise_level_grid.ravel())
]
log_marginal_likelihood = np.reshape(log_marginal_likelihood, noise_level_grid.shape)

vmin, vmax = (-log_marginal_likelihood).min(), 50
level = np.around(np.logspace(np.log10(vmin), np.log10(vmax), num=20), decimals=1)
plt.contour(
    length_scale_grid,
    noise_level_grid,
    -log_marginal_likelihood,
    levels=level,
    norm=LogNorm(vmin=vmin, vmax=vmax),
)
plt.colorbar()
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Length-scale")
plt.ylabel("Noise-level")
plt.title("Negative log-marginal-likelihood")
plt.show()

Мы видим два локальных минимума, соответствующих ранее найденным комбинациям гиперпараметров. В зависимости от начальных значений гиперпараметров оптимизация на основе градиента может как сойтись к наилучшей модели, так и не сойтись. Поэтому важно повторять оптимизацию несколько раз с разными начальными значениями. Это можно сделать, задав параметр n_restarts_optimizer класса GaussianProcessRegressor.

Попробуем снова обучить модель с неудачными начальными значениями, но на этот раз выполним 10 случайных перезапусков.

kernel = 1.0 * RBF(length_scale=1e1, length_scale_bounds=(1e-2, 1e3)) + WhiteKernel(
    noise_level=1, noise_level_bounds=(1e-10, 1e1)
)
gpr = GaussianProcessRegressor(
    kernel=kernel, alpha=0.0, n_restarts_optimizer=10, random_state=0
)
gpr.fit(X_train, y_train)
y_mean, y_std = gpr.predict(X, return_std=True)

plt.plot(X, y, label="Expected signal")
plt.scatter(x=X_train[:, 0], y=y_train, color="black", alpha=0.4, label="Observations")
plt.errorbar(X, y_mean, y_std, label="Posterior mean ± std")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
_ = plt.title(
    (
        f"Initial: {kernel}\nOptimum: {gpr.kernel_}\nLog-Marginal-Likelihood: "
        f"{gpr.log_marginal_likelihood(gpr.kernel_.theta)}"
    ),
    fontsize=8,
)

Как мы и надеялись, случайные перезапуски позволяют оптимизации найти наилучший набор гиперпараметров, несмотря на неудачные начальные значения.

Общее время выполнения скрипта: (0 минут 4.420 секунды)

Download Jupyter notebook: plot_gpr_noisy.ipynb

Download Python source code: plot_gpr_noisy.py

Download zipped: plot_gpr_noisy.zip

Связанные примеры

Иллюстрация априорного и апостериорного гауссовских процессов для разных ядер

Регрессия на основе гауссовских процессов: простой вводный пример

Вероятностные предсказания при классификации на основе гауссовских процессов (GPC)

Прогнозирование уровня CO2 на наборе данных Мауна-Лоа с помощью регрессии на основе гауссовских процессов (GPR)