高斯过程GaussianProcess
高斯过程的理论知识
非参数方法的基本思想
高斯过程的基本概念
高斯过程的Python实现
使用Numpy手动实现
定义核函数
import numpy as np
def kernel(X1, X2, l=1.0, sigma_f=1.0):
'''
各向同性平方指数内核.
计算点X1与点X2的协方差矩阵.
Args:
X1: ndArray， m个点 (m x d).
X2: ndArray， n个点 (n x d).
返回:
协方差矩阵 (m x n).
'''
sqdist = np.sum(X1**2, 1).reshape(-1, 1) np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
return sigma_f**2 * np.exp(-0.5 / l**2 * sqdist)
定义先验
# 编写绘图函数
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
def plot_gp(mu, cov, X, X_train=None, Y_train=None, samples=[]):
X = X.ravel()
mu = mu.ravel()
uncertainty = 1.96 * np.sqrt(np.diag(cov))
plt.fill_between(X, mu uncertainty, mu - uncertainty, alpha=0.1)
plt.plot(X, mu, label='均值')
for i, sample in enumerate(samples):
plt.plot(X, sample, lw=1, ls='--', label=f'样本 {i 1}')
if X_train is not None:
plt.plot(X_train, Y_train, 'rx')
plt.legend(bbox_to_anchor=(1.04,0.5), loc="center left")
def plot_gp_2D(gx, gy, mu, X_train, Y_train, title, i):
ax = plt.gcf().add_subplot(1, 2, i, projection='3d')
ax.plot_surface(gx, gy, mu.reshape(gx.shape), cmap=cm.coolwarm, linewidth=0, alpha=0.2, antialiased=False)
ax.scatter(X_train[:,0], X_train[:,1], Y_train, c=Y_train, cmap=cm.coolwarm)
ax.set_title(title)
%matplotlib inline
# 有限个输入数据点
X = np.arange(-5, 5, 0.2).reshape(-1, 1)
# 先验的均值与方差
mu = np.zeros(X.shape)
cov = kernel(X, X)
# 从先验分布（多元高斯分布）中抽取样本点
samples = np.random.multivariate_normal(mu.ravel(), cov, 3)
# 画出GP的均值, 置信区间
plot_gp(mu, cov, X, samples=samples)
基于无噪声训练数据进行预测
为了计算充分统计量，即后验预测分布的均值和协方差矩阵，我们用下面代码实现公式（4）和（5）
代码语言：javascript复制# 倒入计算逆矩阵的函数inv()
from numpy.linalg import inv
def posterior_predictive(X_s, X_train, Y_train, l=1.0, sigma_f=1.0, sigma_y=1e-8):
'''
计算后验分布的充分统计量
给定 m 个训练数据 X_train 与 Y_train
给定 n 个新输入数据 inputs X_s.
Args:
X_s: 新输入数据 (n x d).
X_train: 训练输入数据 (m x d).
Y_train: 训练输出数据 (m x 1).
l: 核函数的长度参数.
sigma_f: 核函数的纵向波动参数.
sigma_y: 噪音参数.
返回:
后验均值向量 (n x d) 与协方差矩阵 (n x n).
'''
K = kernel(X_train, X_train, l, sigma_f) sigma_y**2 * np.eye(len(X_train))
K_s = kernel(X_train, X_s, l, sigma_f)
K_ss = kernel(X_s, X_s, l, sigma_f) 1e-8 * np.eye(len(X_s))
K_inv = inv(K)
# 公式 (4)
mu_s = K_s.T.dot(K_inv).dot(Y_train)
# 公式 (5)
cov_s = K_ss - K_s.T.dot(K_inv).dot(K_s)
return mu_s, cov_s
并将它们应用于无噪声训练数据X_train和Y_train。以下示例从后验预测中提取3个样本，并将它们与均值，置信区间和训练数据一起绘制。在无噪声模型中，训练点的方差为 0
代码语言：javascript复制注： 从后验预测分布提取的所有随机函数都经过训练点。
# 无噪音的5个输入数据
X_train = np.array([-4, -3, -2, -1, 1]).reshape(-1, 1)
# y=sin(x)
Y_train = np.sin(X_train)
# 计算后验预测分布的均值向量与协方差矩阵
mu_s, cov_s = posterior_predictive(X, X_train, Y_train)
# 从后验预测分布中抽取3个样本
samples = np.random.multivariate_normal(mu_s.ravel(), cov_s, 3)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train, samples=samples)
基于带有噪音的训练数据进行预测
如果模型中包含一些噪声，则仅对训练点进行近似，并且训练点的方差不为零。
代码语言：javascript复制# 定义噪音参数
noise = 0.4
# 带有噪音的训练数据
X_train = np.arange(-3, 4, 1).reshape(-1, 1)
Y_train = np.sin(X_train) noise * np.random.randn(*X_train.shape)
# 可以对比地看一下带噪音的训练数据与不带噪音的训练数据的区别
plt.figure()
plt.plot(X_train, np.sin(X_train) noise * np.random.randn(*X_train.shape), lw=1, ls='-.', label='有噪音')
plt.plot(X_train, np.sin(X_train) 0.0 * np.random.randn(*X_train.shape), lw=2.5, ls='-', label='无噪音')
plt.plot(X_train, np.sin(X_train) 0.0 * np.random.randn(*X_train.shape), 'rx')
plt.legend(bbox_to_anchor=(1.04,0.5), loc="center left")
plt.show()
# 计算后验预测分布的均值向量以及方差矩阵
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, sigma_y=noise)
# 从后验预测分布中抽取3个样本点
samples = np.random.multivariate_normal(mu_s.ravel(), cov_s, 3)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train, samples=samples)
核函数参数和噪声参数的影响
params = [
(0.3, 1.0, 0.2),
(3.0, 1.0, 0.2),
(1.0, 0.3, 0.2),
(1.0, 3.0, 0.2),
(1.0, 1.0, 0.05),
(1.0, 1.0, 1.5),
]
plt.figure(figsize=(12, 5))
for i, (l, sigma_f, sigma_y) in enumerate(params):
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, l=l,
sigma_f=sigma_f,
sigma_y=sigma_y)
plt.subplot(3, 2, i 1)
plt.subplots_adjust(top=2)
plt.title(f'l = {l}, sigma_f = {sigma_f}, sigma_y = {sigma_y}')
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)
这些参数的最优值可以通过最大化由[1] [3]给出的边际对数似然来得到：
from numpy.linalg import cholesky, det, lstsq
from scipy.optimize import minimize
def nll_fn(X_train, Y_train, noise, naive=True):
'''
基于给定数据X_train和Y_train以及噪声水平
返回一个可以计算负边际对数似然的函数
Args:
X_train: 训练输入 (m x d).
Y_train: 训练输出 (m x 1).
noise: Y_train的噪声水平.
naive: 如果 True 那么使用公式(7)来实现
如果 False 那么使用数值方法来实现
返回:
最小的目标对象
'''
def nll_naive(theta):
# 使用公式(7)来实现
# 与下面的nll_stable的实现相比在数值上不稳定
K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1])
noise**2 * np.eye(len(X_train))
return 0.5 * np.log(det(K))
0.5 * Y_train.T.dot(inv(K).dot(Y_train))
0.5 * len(X_train) * np.log(2*np.pi)
def nll_stable(theta):
# 数值上更稳定，相比于nll_naive
K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1])
noise**2 * np.eye(len(X_train))
L = cholesky(K)
return np.sum(np.log(np.diagonal(L)))
0.5 * Y_train.T.dot(lstsq(L.T, lstsq(L, Y_train)[0])[0])
0.5 * len(X_train) * np.log(2*np.pi)
if naive:
return nll_naive
else:
return nll_stable
# 求解可满足最小化目标函数的参数 l 及 sigma_f.
# 实际上，我们应该使用不同的初始化多次运行最小化，以避免局部最小化，
# 但是为了简单起见，此处将其跳过
res = minimize(nll_fn(X_train, Y_train, noise), [1, 1],
bounds=((1e-5, None), (1e-5, None)),
method='L-BFGS-B')
# 将优化结果存储在全局变量中，以便我们以后可以将其与其他实现的结果进行比较
l_opt, sigma_f_opt = res.x
# 使用优化的核函数参数计算后验预测分布的参数，并绘制结果图
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, l=l_opt, sigma_f=sigma_f_opt, sigma_y=noise)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)
print(l_opt, sigma_f_opt)
0.9872536793237083 0.8613778055591963
更高维的高斯过程
#噪音参数
noise_2D = 0.1
rx, ry = np.arange(-5, 5, 0.3), np.arange(-5, 5, 0.3)
gx, gy = np.meshgrid(rx, rx)
X_2D = np.c_[gx.ravel(), gy.ravel()]
X_2D_train = np.random.uniform(-4, 4, (100, 2))
Y_2D_train = np.sin(0.5 * np.linalg.norm(X_2D_train, axis=1))
noise_2D * np.random.randn(len(X_2D_train))
plt.figure(figsize=(14,7))
mu_s, _ = posterior_predictive(X_2D, X_2D_train, Y_2D_train, sigma_y=noise_2D)
plot_gp_2D(gx, gy, mu_s, X_2D_train, Y_2D_train,
f'参数优化前: l={1.00} sigma_f={1.00}', 1)
res = minimize(nll_fn(X_2D_train, Y_2D_train, noise_2D), [1, 1],
bounds=((1e-5, None), (1e-5, None)),
method='L-BFGS-B')
mu_s, _ = posterior_predictive(X_2D, X_2D_train, Y_2D_train, *res.x, sigma_y=noise_2D)
plot_gp_2D(gx, gy, mu_s, X_2D_train, Y_2D_train,
f'参数优化后: l={res.x[0]:.2f} sigma_f={res.x[1]:.2f}', 2)
使用Scikit-learn实现高斯过程
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF
rbf = ConstantKernel(1.0) * RBF(length_scale=1.0)
gpr = GaussianProcessRegressor(kernel=rbf, alpha=noise**2)
# 1D 训练样本
gpr.fit(X_train, Y_train)
# 计算后验预测分布的均值向量与协方差矩阵
mu_s, cov_s = gpr.predict(X, return_cov=True)
# 获得最优核函数参数
l = gpr.kernel_.k2.get_params()['length_scale']
sigma_f = np.sqrt(gpr.kernel_.k1.get_params()['constant_value'])
# 与前面手写的结果比对
assert(np.isclose(l_opt, l))
assert(np.isclose(sigma_f_opt, sigma_f))
# 绘制结果
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)
小结
从前面我们可以看出，与常见的机器学习模型不同，用高斯过程做预测的方法是直接生成一个后验预测分布（依然是高斯分布）。
这也决定了我们可以不仅仅得到一个“光秃秃”的预测值，还可以得到关于这个预测值的不确定性信息，可以利用这些信息绘制error-bar等等！
从统计学的角度上来看，利用高斯过程模型做预测具有很高的价值。
