高斯过程回归(Gaussian process regression)原理详解及python代码实战

1. 高斯过程回归原理

高斯过程回归(Gaussian process regression,GPR)是一个随机过程按时间或空间索引的随机变量集合这些随机变量的每个有限集合都服从多元正态分布即它们的每个有限线性组合都是正态分布。高斯过程的分布是所有这些无限多随机变量的联合概率分布。

1.1 高斯过程

定义一个高斯过程是一组随机变量的集合这组随机变量的每个有限子集构成的联合概率分布都服从多元高斯分布即
$$f\sim GP(\mu ,k)(1-1)$$
其中$\mu (x)$和$k(x,x^{\prime })$分别为随机变量x的均值函数和协方差函数。因此可看出一个高斯过程实际是由随机变量的均值和协方差函数所决定。

1.2 高斯过程回归

在传统回归模型中定义$Y=f(x)$而在高斯过程回归中设f(x)服从高斯分布。通常假设均值为0则
$$Y=f(x)\sim N(0,\Sigma )(1-2)$$
其中$\Sigma$是协方差函数。使用核技巧令${\Sigma }_{ij}=K({\mathbf{x}}_i,{\mathbf{x}}_j)$则可由求核函数K来计算协方差函数。如核函数K使用RBF核则${\Sigma }_{ij}=\tau e^{\frac{-\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }^2}{{\sigma }^2}}$。如使用多项式核则${\Sigma }_{ij}=\tau (1+{\mathbf{x}}_i^{\top }{\mathbf{x}}_j{)}^d$。
因此协方差函数$\Sigma$可分解为$\left(\begin{array}{c}K,K_{\ast }\\ K_{\ast }^{\top },K_{\ast \ast }\end{array}\right)$。其中K是训练核矩阵$K_{\ast }$是训练-测试核矩阵$K_{\ast }^T$是测试-训练核矩阵。
隐函数f的条件概率分布可表示为
$$f_{\ast }|(Y_1=y_1,...,Y_n=y_n,{\mathbf{x}}_1,...,{\mathbf{x}}_n,{\mathbf{x}}_t)\sim \mathcal{N}(K_{\ast }^{\top }{K}^{-1}y,K_{\ast \ast }-K_{\ast }^{\top }{K}^{-1}{K}_{\ast })(1-3)$$

---

2. python实现高斯过程回归

2.1 参数详解

基于机器学习库sklearn实现高斯过程回归。sklearn中GaussianProcessRegressor模块实现了高斯过程回归模型从模型参数、属性和方法等方面介绍该模型其主要参数包括

参数名	参数含义	备注
kernel	核函数形式的高斯过程的协方差函数	常用核函数有RBF、ConstantKernel核函数的常见超参数有核的长度尺寸、长度尺寸的上下限
alpha	在模型拟合过程中加入核矩阵对角线位置的值	1确保计算值形成正定矩阵防止拟合过程中出现潜在的数值问题2也可以解释为训练观测值上附加高斯测量噪声的方差3如果alpha参数传递的是一个数组它必须与用于拟合的数据具有相同的条目数并用作依赖于数据点的噪声级
random_state	随机状态数	决定用于初始化中心的随机数的生成在多次函数调用时指定此参数保证可复现性

GaussianProcessRegressor回归模型的主要属性包括

属性名称	尺寸
X_train_	array-like of shape (n_samples, n_features) or list of object
y_train_	array-like of shape (n_samples,) or (n_samples, n_targets)

GaussianProcessRegressor回归模型的常用方法包括

方法名称	参数	返回值
predict(X, return_std=False, return_cov=False)	X是高斯过程要拟合的样本点	用GPR模型进行预测返回样本点预测概率分布的均值、标准差和预测联合概率分布的协方差
score(X, y, sample_weight=None)	X是测试样本, y是X对应的真值	返回GPR模型预测的$R^2$系数$R^2$系数计算方法为$$1-\frac{(((y_{true}-y_{pred})\ast \ast 2).sum()}{(y_{true}-y_{true}.mean())\ast \ast 2).sum()}$$

2.2 核函数cookbook

核函数在sklearn.gaussian_process.kernels模块中常用的核函数有

• RBF核函数Radial basis function kernel
  RBF核函数又被称为平方指数核其计算方式为
  $$k(x_i,x_j)=\exp (-\frac{d(x_i,x_j{)}^2}{2l^2})(2-1)$$
  式中$l$是核函数的长度尺寸$d$是距离度量函数这里采用欧式距离度量。

在sklearn中RBF函数有两个参数

RBF(length_scale=1.0, length_scale_bounds=(1e-05, 100000.0))
# length_scale:核函数的长度尺寸float or ndarray of shape (n_features,), default=1.0
# length_scale_bounds核函数长度尺寸的上下限若设为'fixed'则无法在超参数调整期间修改核函数长度尺寸。

• 常数核ConstantKernel
  常数核可以作为乘积核(product kernel)的一部分用于缩放其他因子核的大小也可以用作和核的一个部分用于修改高斯过程的平均值。其数学形式表示为
  $$k(x_1,x_2)=constant\forall x_1,x_2(2-2)$$
  同样在skelearn中有ConstantKernel函数两个参数,含义与RBF的参数类似。

kernel = RBF() + ConstantKernel(constant_value=2,constant_value_bounds=(1e-5, 1e5))
# 作用等价于kernel = RBF() + 2

2.2 代码模版

导入库

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C, WhiteKernel
import sklearn.gaussian_process.kernels as k

训练GPR模型

def gpr_regressor(X_train,y_train,X_test,y_test,kernel=C(constant_value=1) * RBF(length_scale=1, length_scale_bounds=(1e-2, 1e2))):
    """
    gpr model for regression
    :param X_train: (n_samples, n_features)
    :param y_train: (n_samples,)
    :param X_test: (n_samples, n_features)
    :param y_test: (n_samples,)
    :param kernel: kernel of gpr
    :return:
        y_pred: mean predictions
        y_pred_std: std predictions
        r2: r2 score of gpr
    """
    gp = GaussianProcessRegressor(kernel=kernel)
    gp.fit(X_train, y_train)  # Instantiated Gaussian regression model
    print("the learned kernel parameters:\t {}".format(gp.kernel_))  # the learned kernel parameters
    y_pred, y_pred_std = gp.predict(X_test, return_std=True)
    r2 = gp.score(X_test, y_test)
    print('r2 coefficient is {:.2f}'.format(r2))
    return y_pred,y_pred_std,r2

预测结果可视化

def plot_errorbar_gpr(y_pred,y_pred_std,r2,y_test):
    """
    plot errorbar for gpr predictions
    :param y_pred:
    :param y_pred_std:
    :param r2:
    :param y_test: one-dimension
    :return:
    """

    plt.errorbar(x=y_test, y=y_pred, yerr=y_pred_std, fmt="o", label="Samples", markersize=5,color='#2698eb')
    #x, y define the data locations, xerr, yerr define the errorbar sizes
    plt.xlabel("ground true")
    plt.ylabel("predicted ")
    plt.title("Gaussian process regression, R2=%.2f" % (r2))
    print("finished!")

def plot_intervel_gpr(y_pred,y_pred_std,r2,X_test):
    """
    plot confidence interval for gpr predictions
    :param y_pred:
    :param y_pred_std:
    :param r2:
     :param X_test: should be one-dimension shape
    :return:
    """
    # 1.96 sigma = 95% confidence interval for a normal distribution
    upper, lower = y_pred + 1.96 * y_pred_std, y_pred - 1.96 * y_pred_std
    plt.plot(X_test, y_pred, label="GPR", ls="-")
    plt.fill_between(X_test, y1=upper, y2=lower, alpha=0.2, label="95% confidence", color="#2698eb")
    plt.legend(ncol=4, fontsize=12)
    plt.title("Gaussian process regression, R2=%.2f" %(r2))
    print("finished!")

预测概率区间图如下

以标准差为尺度的误差线图如下

---

附录-数学基础知识

这里列出了高斯过程回归涉及到的数学基础知识方便大家参考。

A1 高斯分布的基本性质

高斯分布的四大属性标准化Normalization、边缘化(Marginalization)、可加性(Summation)、条件性(Conditioning)具体数学表示如下图所示

A2 贝叶斯框架

贝叶斯框架的基础概念包括条件概率、乘积法则、加和法则、贝叶斯定理等。
1条件概率
$$条件概率分布=\frac{联合概率分布}{边缘概率分布}(A-1)$$ i.e. $p(y|x)=\frac{p(x,y)}{p(x)}$.

2乘积法则(product rule)
任何联合概率分布都可以表示为一维条件概率分布的乘积i.e. $p(x,y,z)=p(x|y,z)p(y|z)p(z)$.

(3) 加和法则(sum rule)
任何边缘概率分布都可以通过对联合概率分布的积分获得i.e.
$p(y)=\int p(x,y)dx$, $p(x)=\int p(x,y)dy$.

利用加和法则还可以估计条件概率分布如已知一组随机变量$x$,$y$,$z$构成的联合概率分布$p(x,y,z)$观测到$x$感兴趣预测$y$,变量$z$是未知的且与待解决的问题无关如何只用联合概率分布$p(x,y,z)$去估计$p(y|x)$呢
答案就是用加和法则对联合概率分布求积分
$p(y|x)=\frac{p(x,y)}{p(x)}$=$\frac{\int p(x,y,z)dz}{\int p(x,y,z)dzdy}$

(4)贝叶斯定理(Bayes theorem)
由乘积法则和加和法则可知任何条件概率分布可由联合概率分布获得。进一步的利用乘积法则和加和法则可从条件概率的定义中推导出贝叶斯定理
$$p(y|x)=\frac{p(x,y)}{p(x)}=\frac{p(x|y)p(y)}{p(x)}=\frac{p(x|y)p(y)}{\int p(x|y)p(y)dy}(A-2)$$.
从式中可看出贝叶斯定理定义了当新的信息出现时不确定性的转化规则即
$$后验概率=\frac{似然概率\times 先验概率}{置信度}(A-3)$$
(5)贝叶斯框架
统计推断方法有频率框架和贝叶斯框架两大流派贝叶斯框架基于贝叶斯定理去估计后验概率。应用该框架进行推断的经典描述为给定一组独立同分布的数据$X=(x_1,x_2,...,x_n)$欲用概率分布$p(x|\theta )$估计$\theta$贝叶斯框架提供了使用先验概率$p(\theta )$编码参数$\theta$的不确定度的方法
$$p(\theta |X)=\frac{{\prod }_{i=1}^{n}p(x_i|\theta )p(\theta )}{\int {\prod }_{i=1}^{n}p(x_i|\theta )p(\theta )d\theta }(A-4)$$

A3 后验预测分布

考虑一个回归问题
$$y=f(x)+\epsilon (A-5)$$
为从数据中求解概率给定特定参数w预测模型可利用后验预测分布posterior predictive distribution来求解。
对于数据集D和特征X根据贝叶斯定理以及概率分布的加和、乘积法则后验预测分布可表示为
$$P(Y\mid D,X)={\int }_{\mathbf{w}}P(Y,\mathbf{w}\mid D,X)d\mathbf{w}$$
$$={\int }_{\mathbf{w}}P(Y\mid \mathbf{w},D,X)P(\mathbf{w}\mid D)d\mathbf{w}(A-6)$$
然而上述公式的解析解(closed form)难以计算但对于具有高斯似然和先验的特殊情况我们可以求其高斯过程的均值和协方差即
$$P(y_{\ast }\mid D,\mathbf{x})\sim \mathcal{N}({\mu }_{y_{\ast }\mid D},{\Sigma }_{y_{\ast }\mid D})(A-7)$$
其中均值函数${\mu }_{y_{\ast }\mid D}$和协方差函数${\Sigma }_{y_{\ast }\mid D}$, 均可由核函数K计算得到
$${\mu }_{y_{\ast }\mid D}=K_{\ast }^T(K+{\sigma }^{2}I{)}^{-1}y(A-8)$$
$${\Sigma }_{y_{\ast }\mid D}=K_{\ast \ast }-K_{\ast }^T(K+{\sigma }^{2}I{)}^{-1}{K}_{\ast }(A-9)$$

参考资料

[1]cornell 高斯过程lecture
[2]https://cosmiccoding.com.au/tutorials/gaussian_processes
[3]sklearn GaussianProcessRegressor文档
[4].高斯过程可视化展示网页
[5].http://gaussianprocess.org/gpml/chapters/RW.pdf
[6].https://towardsdatascience.com/getting-started-with-gaussian-process-regression-modeling-47e7982b534d