高斯过程回归(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模块实现了高斯过程回归模型，从模型参数、属性和方法等方面介绍该模型，其主要参数包括：

gaussianprocessregressor回归模型的主要属性包括：

gaussianprocessregressor回归模型的常用方法包括：

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