GPR tutorial

1. 高斯过程回归原理
- 1.1 高斯过程
- 1.2 高斯过程回归
2. python实现高斯过程回归
- 2.1 参数详解
- 2.2 核函数cookbook
- 2.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模块实现了高斯过程回归模型，从模型参数、属性和方法等方面介绍该模型，其主要参数包括：

参数名	参数含义	备注
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模型预测的R2R^2R2系数，R2R^2R2系数计算方法为1−(((ytrue−ypred)∗∗2).sum()(ytrue−ytrue.mean())∗∗2).sum()1-\frac{(((y_{true} - y_{pred}) 2).sum()}{(y_{true} - y_{true}.mean()) 2).sum()}1−(ytrue−ytrue.mean())∗∗2).sum()(((ytrue−ypred)∗∗2).sum()

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

A3 后验预测分布

考虑一个回归问题：
y=f(x)+ε(A−5)y=f(x)+\varepsilon \qquad(A-5)y=f(x)+ε(A−5)
为从数据中求解概率，给定特定参数w，预测模型可利用后验预测分布（posterior predictive distribution）来求解。
对于数据集D和特征X，根据贝叶斯定理以及概率分布的加和、乘积法则，后验预测分布可表示为：
P(Y∣D,X)=∫wP(Y,w∣D,X)dwP(Y\mid D,X) = \int_{\mathbf{w}}P(Y,\mathbf{w} \mid D,X) d\mathbf{w}P(Y∣D,X)=∫wP(Y,w∣D,X)dw
=∫wP(Y∣w,D,X)P(w∣D)dw(A−6)= \int_{\mathbf{w}} P(Y \mid \mathbf{w}, D,X) P(\mathbf{w} \mid D) d\mathbf{w} \qquad(A-6)=∫wP(Y∣w,D,X)P(w∣D)dw(A−6)
然而上述公式的解析解(closed form)难以计算，但对于具有高斯似然和先验的特殊情况，我们可以求其高斯过程的均值和协方差，即：
P(y∗∣D,x)∼N(μy∗∣D,Σy∗∣D)(A−7)P(y_*\mid D,\mathbf{x}) \sim \mathcal{N}(\mu_{y_*\mid D}, \Sigma_{y_*\mid D}) \qquad(A-7)P(y∗∣D,x)∼N(μy∗∣D,Σy∗∣D)(A−7)
其中均值函数μy∗∣D\mu_{y_*\mid D}μy∗∣D和协方差函数Σy∗∣D\Sigma_{y_*\mid D}Σy∗∣D, 均可由核函数K计算得到：
μy∗∣D=K∗T(K+σ2I)−1y(A−8)\mu_{y_*\mid D} = K_*^T (K+\sigma^2 I)^{-1} y \qquad(A-8)μy∗∣D=K∗T(K+σ2I)−1y(A−8)
Σy∗∣D=K∗∗−K∗T(K+σ2I)−1K∗(A−9)\Sigma_{y_*\mid D} = K_{**} - K_*^T (K+\sigma^2 I)^{-1}K_* \qquad(A-9)Σy∗∣D=K∗∗−K∗T(K+σ2I)−1K∗(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

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

高斯过程回归(Gaussian Process Regression) 粗理解
转自:高斯过程回归 (Gaussian Process Regression) 学习路线 - 知乎一个写高斯过程不错的帖子: 贝叶斯优化 (上) 高斯过程 - 知乎贝叶斯优化 (下) 高斯过程 - ...
视频教程-深度学习原理详解及Python代码实现-深度学习
深度学习原理详解及Python代码实现大学教授,美国归国博士.博士生导师:人工智能公司专家顾问:长期从事人工智能.物联网.大数据研究:已发表学术论文100多篇,授权发明专利10多项白勇 ¥88.0 ...
GBDT（回归树）原理详解与python代码实现
GBDT算法 1.算法原理 2.对数据的要求 3.算法的优缺点 4.算法需要注意的点 5.python代码实现(待更......) 导入相关包读取数据并预处理训练及评估 1.算法原理步骤: 1. ...
随机森林原理详解及python代码实现
随机森林(RF)算法 1.算法原理 2.对数据的要求(无需规范化) 3.算法的优缺点 4.算法需要注意的点 5.python代码实现(待更......) 导入相关包读取数据并预处理(必须处理缺失值) ...
决策树原理详解及python代码实现
决策树算法(信贷中常用来寻找规则) 1.算法原理 1.1 ID3(多叉树分类) 1.2 C4.5(多叉树分类) 1.3 Cart(二叉树分类+回归) 2.ID3.C4.5与Cart比较 3.算法优缺点 ...
KNN算法原理详解及python代码实现
KNN算法算法原理对数据的要求算法的优缺点算法需要注意的点算法实现(python) 算法原理计算待测样本与train_data的距离d并保存数组中对d进行排序,取d最近的k个样本统计样 ...
Dijkstra 路径规划算法原理详解及 Python 代码实现
荷兰数学家 E.W.Dijkstra 于 1959 年提出了 Dijkstra 算法,它是一种适用于非负权值网络的单源最短路径算法,同时也是目前求解最短路径问题的理论上最完备.应用最广的经典算法 ...
google authenticator python_Google Authenticator TOTP原理详解(以Python为例)
http://xsboke.blog.51cto.com 如果有疑问,请点击此处,然后发表评论交流,作者会及时回复(也可以直接在当前文章评论). -------谢谢您的参考,如有疑问,欢迎交流一. ...
DS18B20温度传感器原理详解及例程代码、漏极开路
[常用传感器]DS18B20温度传感器原理详解及例程代码_Z小旋的博客-CSDN博客_ds18b20温度传感器传感器引脚及原理图 DS18B20传感器的引脚及封装图如下: DS18B20一共有三个引 ...
计算机组织与结构poc,CPU漏洞原理详解以及POC代码分享
原标题:CPU漏洞原理详解以及POC代码分享首先,这个漏洞已经公布近一周时间了,看到各大媒体.公众号到处在宣传,本打算不再发布类似信息,但是发现很多媒体的报道达到了一个目的--几乎所有的CPU都有漏 ...

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