Есть ли простой код для начинающих, где я могу экспериментировать с различными ядрами, используемыми в примере гауссовского процесса в Scikit, научиться знать их функции?

Question

на самом деле я хочу понять, как ядра, используемые в scikit, изучают гауссовский пример, но у меня нет нулевых знаний о том, как работает это ядро ​​и когда его использовать, и я также не получаю никакого примера базового c шаблонного кода, где я могу использовать это ядро один за другим и понять. Частичный код приведен ниже:

X, y = load_mauna_loa_atmospheric_co2()

Ядро с параметрами, указанными в книге GPML

k1 = 66.0**2 * RBF(length_scale=67.0)  # long term smooth rising trend
k2 = 2.4**2 * RBF(length_scale=90.0) \
    * ExpSineSquared(length_scale=1.3, periodicity=1.0)  # seasonal component
# medium term irregularity
k3 = 0.66**2 \
    * RationalQuadratic(length_scale=1.2, alpha=0.78)
k4 = 0.18**2 * RBF(length_scale=0.134) \
    + WhiteKernel(noise_level=0.19**2)  # noise terms
kernel_gpml = k1 + k2 + k3 + k4

gp = GaussianProcessRegressor(kernel=kernel_gpml, alpha=0,
                              optimizer=None, normalize_y=True)
gp.fit(X, y)

print("GPML kernel: %s" % gp.kernel_)
print("Log-marginal-likelihood: %.3f"
      % gp.log_marginal_likelihood(gp.kernel_.theta))

# Kernel with optimized parameters
k1 = 50.0**2 * RBF(length_scale=50.0)  # long term smooth rising trend
k2 = 2.0**2 * RBF(length_scale=100.0) \
    * ExpSineSquared(length_scale=1.0, periodicity=1.0,
                     periodicity_bounds="fixed")  # seasonal component
# medium term irregularities
k3 = 0.5**2 * RationalQuadratic(length_scale=1.0, alpha=1.0)
k4 = 0.1**2 * RBF(length_scale=0.1) \
    + WhiteKernel(noise_level=0.1**2,
                  noise_level_bounds=(1e-3, np.inf))  # noise terms
kernel = k1 + k2 + k3 + k4

gp = GaussianProcessRegressor(kernel=kernel, alpha=0,
                              normalize_y=True)
gp.fit(X, y)

print("\nLearned kernel: %s" % gp.kernel_)
print("Log-marginal-likelihood: %.3f"
      % gp.log_marginal_likelihood(gp.kernel_.theta))

X_ = np.linspace(X.min(), X.max() + 30, 1000)[:, np.newaxis]
y_pred, y_std = gp.predict(X_, return_std=True)

# Illustration
plt.scatter(X, y, c='k')
plt.plot(X_, y_pred)
plt.fill_between(X_[:, 0], y_pred - y_std, y_pred + y_std,
                 alpha=0.5, color='k')
plt.xlim(X_.min(), X_.max())
plt.xlabel("Year")
plt.ylabel(r"CO$_2$ in ppm")
plt.title(r"Atmospheric CO$_2$ concentration at Mauna Loa")
plt.tight_layout()
plt.show()

Answer 1

Все подробности находятся в книге Расмуссена и Уильямса. Пример, который вы показываете, находится в Главе 5 вместе с подробным объяснением всех используемых ядер. Они также показывают много примеров ковариационных функций и соответствующих случайных функций.

Мне неизвестен код для простой визуализации различных ядер, но можно визуализировать популярную квадратную экспоненциальную функцию, которая несколько раз появляется в примере с Мауна-Лоа с разными шкалами длины, как показано ниже:

import numpy as np
import matplotlib.pyplot as plt

def k_se(r,l):
    return np.exp(-r*r/(2*l*l)) 

r = np.arange(0.1,4,0.01)
plt.figure()
for ll in l:
    plt.plot(r,k_se(r,ll),label='length='+str(np.round(ll,1)))
    plt.xlabel('r')
    plt.ylabel('Covariance k(r)')
plt.legend(frameon=False)

Различные ядра для разных масштабов длины выглядят так:

Однако, что более интересно, это рисование случайных функций из гауссовского процесса, который дал ковариационную функцию. Следующий код предназначен не для эффективности или скорости, а для упрощения визуализации этих случайных функций.

def k_se_p(x1, x2, l):
    return np.exp(-((x1-x2)*(x1-x2))/(2*l*l))

def gm(x,l):
    return [[k_se_p(i,j,l) for j in x] for i in x]

x = np.arange(0.1,8,0.01)

Поучительно сначала нарисовать функции из той же шкалы длины:

plt.figure() 
for i in range(5):
    ys = np.random.multivariate_normal(np.zeros(len(x)), gm(x,l[0]))
    if i==0:
        plt.plot(x,ys,color='blue',label='length='+str(np.round(l[0],1)))
    else:
        plt.plot(x,ys,color='blue')
    plt.xlabel('x')
    plt.ylabel('f(x)') 
plt.legend(frameon=False)

Что дает не очень плавные функции:

Увеличенная шкала длины обеспечивает более плавные функции:

plt.figure() 
for i in range(5):
    ys = np.random.multivariate_normal(np.zeros(len(x)), gm(x,l[-1]))
    if i==0:
        plt.plot(x,ys,color='magenta',label='length='+str(np.round(l[-1],1)))
    else:
        plt.plot(x,ys,color='magenta')
    plt.xlabel('x')
    plt.ylabel('f(x)')
plt.legend(frameon=False)

Наконец, мы можем извлечь одну функцию из каждого Длина шкалы и график их вместе:

plt.figure() 
for ll in l:
    ys = np.random.multivariate_normal(np.zeros(len(x)), gm(x,ll))
    plt.plot(x,ys,label='length='+str(np.round(ll,1)))
    plt.xlabel('x')
    plt.ylabel('f(x)') 
plt.legend(frameon=False)