首页 > 解决方案 > 是否有适合初学者的简单代码,我可以在其中试验 scikit 中高斯过程示例中使用的差异内核来学习了解它们的功能?

问题描述

实际上,我想了解 scikit learn gaussian example 中使用的内核,但我对这些内核的行为方式以及何时使用哪个内核知之甚少,而且我也没有获得任何示例基本模板代码,我可以在其中一一使用这些内核并理解。部分代码如下:

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()

标签: scikit-learnkernelgaussian

解决方案


所有细节都在拉斯穆森和威廉姆斯的书中。您展示的示例在第 5 章中,并详细说明了所有使用的内核。它们还展示了协方差函数和相应随机函数的许多示例。

我不知道有一种代码可以简单地可视化各种内核,但可以可视化流行的平方指数函数,该函数在 Mauna Loa 示例中以不同的长度尺度出现多次,如下所示:

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)

在此处输入图像描述


推荐阅读