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

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

时间：2020-08-05 01:03:36

GPR tutorial

1. 高斯过程回归原理1.1 高斯过程1.2 高斯过程回归2. python实现高斯过程回归2.1 参数详解2.2 核函数cookbook2.2 代码模版附录-数学基础知识A1 高斯分布的基本性质A2 贝叶斯框架A3 后验预测分布参考资料

1. 高斯过程回归原理

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

1.1 高斯过程

定义：一个高斯过程是一组随机变量的集合，这组随机变量的每个有限子集构成的联合概率分布都服从多元高斯分布，即：

f∼GP(μ,k)(1−1)f \sim GP(\mu,k) \qquad(1-1) f∼GP(μ,k)(1−1)

其中μ(x)\mu(x)μ(x)和k(x,x′)k(x,x')k(x,x′)分别为随机变量x的均值函数和协方差函数。因此可看出，一个高斯过程实际是由随机变量的均值和协方差函数所决定。

1.2 高斯过程回归

在传统回归模型中，定义Y=f(x)Y=f(x)Y=f(x)，而在高斯过程回归中，设f(x)服从高斯分布。通常假设均值为0，则：

Y=f(x)∼N(0,Σ)(1−2)Y=f(x)\sim N(0,\Sigma) \qquad(1-2)Y=f(x)∼N(0,Σ)(1−2)

其中Σ\SigmaΣ是协方差函数。使用核技巧，令Σij=K(xi,xj)\Sigma_{ij}=K(\mathbf{x}_i,\mathbf{x}_j)Σij=K(xi,xj)，则可由求核函数K来计算协方差函数。如核函数K使用RBF核，则Σij=τe−∥xi−xj∥2σ2\Sigma_{ij}=\tau e^\frac{-\|\mathbf{x}_i-\mathbf{x}_j\|^2}{\sigma^2}Σij=τeσ2−∥xi−xj∥2。如使用多项式核，则Σij=τ(1+xi⊤xj)d\Sigma_{ij}=\tau (1+\mathbf{x}_i^\top \mathbf{x}_j)^dΣij=τ(1+xi⊤xj)d。

因此，协方差函数Σ\SigmaΣ可分解为(K,K∗K∗⊤,K∗∗)\begin{pmatrix} K, K_* \\K_*^\top , K_{**} \end{pmatrix}(K,K∗K∗⊤,K∗∗)。其中K是训练核矩阵，K∗K_*K∗是训练-测试核矩阵，K∗TK_*^TK∗T是测试-训练核矩阵。

隐函数f的条件概率分布可表示为：

f∗∣(Y1=y1,...,Yn=yn,x1,...,xn,xt)∼N(K∗⊤K−1y,K∗∗−K∗⊤K−1K∗)(1−3)f_*|(Y_1=y_1,...,Y_n=y_n,\mathbf{x}_1,...,\mathbf{x}_n,\mathbf{x}_t)\sim \mathcal{N}(K_*^\top K^{-1}y,K_{**}-K_*^\top K^{-1}K_* ) \qquad(1-3)f∗∣(Y1=y1,...,Yn=yn,x1,...,xn,xt)∼N(K∗⊤K−1y,K∗∗−K∗⊤K−1K∗)(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(xi,xj)=exp⁡(−d(xi,xj)22l2)(2−1)k(x_i,x_j)=\exp(-\frac{d(x_i,x_j)^2}{2l^2}) \qquad(2-1)k(xi,xj)=exp(−2l2d(xi,xj)2)(2−1)

式中lll是核函数的长度尺寸，ddd是距离度量函数，这里采用欧式距离度量。

在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(x1,x2)=constant，∀x1,x2(2−2)k(x_1,x_2)=constant， \forall x_1,x_2 \qquad(2-2) k(x1,x2)=constant，∀x1,x2(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 npimport matplotlib.pyplot as pltimport seaborn as snsfrom sklearn.gaussian_process import GaussianProcessRegressorfrom sklearn.gaussian_process.kernels import RBF, ConstantKernel as C, WhiteKernelimport 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 predictionsy_pred_std: std predictionsr2: r2 score of gpr"""gp = GaussianProcessRegressor(kernel=kernel)gp.fit(X_train, y_train) # Instantiated Gaussian regression modelprint("the learned kernel parameters:\t {}".format(gp.kernel_)) # the learned kernel parametersy_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 sizesplt.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 distributionupper, lower = y_pred + 1.96 * y_pred_std, y_pred - 1.96 * y_pred_stdplt.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）条件概率

条件概率分布=联合概率分布边缘概率分布(A−1)条件概率分布=\frac{联合概率分布}{边缘概率分布} \qquad(A-1)条件概率分布=边缘概率分布联合概率分布(A−1) i.e. p(y∣x)=p(x,y)p(x)p(y|x)=\frac{p(x,y)}{p(x)}p(y∣x)=p(x)p(x,y).

（2）乘积法则(product rule)

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

(3) 加和法则(sum rule)

任何边缘概率分布，都可以通过对联合概率分布的积分获得，i.e.

p(y)=∫p(x,y)dxp(y)=\int p(x,y)dxp(y)=∫p(x,y)dx, p(x)=∫p(x,y)dyp(x)=\int p(x,y)dyp(x)=∫p(x,y)dy.

利用加和法则还可以估计条件概率分布，如已知一组随机变量xxx,yyy,zzz构成的联合概率分布p(x,y,z)p(x,y,z)p(x,y,z)，观测到xxx，感兴趣预测yyy,变量zzz是未知的且与待解决的问题无关，如何只用联合概率分布p(x,y,z)p(x,y,z)p(x,y,z)去估计p(y∣x)p(y|x)p(y∣x)呢？

答案就是用加和法则对联合概率分布求积分：

p(y∣x)=p(x,y)p(x)p(y|x)=\frac{p(x,y)}{p(x)}p(y∣x)=p(x)p(x,y)=∫p(x,y,z)dz∫p(x,y,z)dzdy\frac{\int p(x,y,z)dz}{\int p(x,y,z)dzdy}∫p(x,y,z)dzdy∫p(x,y,z)dz

(4)贝叶斯定理(Bayes theorem)

由乘积法则和加和法则可知，任何条件概率分布可由联合概率分布获得。进一步的，利用乘积法则和加和法则，可从条件概率的定义中推导出贝叶斯定理：

p(y∣x)=p(x,y)p(x)=p(x∣y)p(y)p(x)=p(x∣y)p(y)∫p(x∣y)p(y)dy(A−2)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} \qquad(A-2)p(y∣x)=p(x)p(x,y)=p(x)p(x∣y)p(y)=∫p(x∣y)p(y)dyp(x∣y)p(y)(A−2).

从式中可看出，贝叶斯定理定义了当新的信息出现时，不确定性的转化规则，即：

后验概率=似然概率×先验概率置信度(A−3)后验概率=\frac{似然概率\times 先验概率}{置信度} \qquad(A-3)后验概率=置信度似然概率×先验概率(A−3)

(5)贝叶斯框架

统计推断方法有频率框架和贝叶斯框架两大流派，贝叶斯框架基于贝叶斯定理去估计后验概率。应用该框架进行推断的经典描述为，给定一组独立同分布的数据X=(x1,x2,...,xn)X=(x_1,x_2,...,x_n)X=(x1,x2,...,xn)，欲用概率分布p(x∣θ)p(x|\theta)p(x∣θ)估计θ\thetaθ，贝叶斯框架提供了使用先验概率p(θ)p(\theta)p(θ)编码参数θ\thetaθ的不确定度的方法：

p(θ∣X)=∏i=1np(xi∣θ)p(θ)∫∏i=1np(xi∣θ)p(θ)dθ(A−4)p(\theta|X)=\frac{\prod_{i=1}^np(x_i|\theta)p(\theta)}{\int \prod_{i=1}^np(x_i|\theta)p(\theta)d\theta} \qquad(A-4)p(θ∣X)=∫∏i=1np(xi∣θ)p(θ)dθ∏i=1np(xi∣θ)p(θ)(A−4)