高斯过程（原理和代码实现）

1.高斯过程原理

• 每个点的观测值都是高斯分布，这里面的观测值就是输出$Y$,观测点的组合也符合高斯分布。高斯过程通常可以用来表示一个函数，更具体来说是表示一个函数的分布。高斯过程是非参数化的，针对小样本学习具有很好的效果。参数化的方法把可学习的函数的范围限制死了，无法学习任意类型的函数。而非参数化的方法就没有这个缺点。高斯过程直观来说，两个$x$离得越近，对应的函数值$y$应该相差越小的原理对核函数的参数进行学习。
• 高斯过程中假设每一个观测值都是符合高斯过程，高斯过程是非参数化的，针对小样本学习具有很好的效果。参数化的方法把可学习的函数的范围限制死了，无法学习任意类型的函数。而非参数化的方法就没有这个缺点。高斯过程直观来说，两个x离得越近，对应的函数值y应该相差越小的原理对核函数的参数进行学习。
• 高斯过程主要的两个数值就是方差和均值，利用已知的样本数据构建协方差均值，将需要预测观测值的数据导入计算这些预测观测值的均值方差，然后就能得到预测值，在这个过程中比较重要的是要选取核函数，核函数用来计算映射到高维空间之后内积的一种简便计算方式，将低维线性不可分的数据投影到高维之后变得线性可分的效果。

①高斯分布

  高斯分布（Gaussian distribution），又名正态分布（Normal distribution），正态曲线呈钟型，两头低，中间高，左右对称因其曲线呈钟形，因此人们又经常称之为钟形曲线。
  若随机变量$X$服从一个数学期望为$\mu$、方差为${\sigma }^2$的正态分布，记为$N(\mu ,{\sigma }^2)$。其概率密度函数为正态分布的期望值$\mu$决定了其位置，其标准差$\sigma$决定了分布的幅度。当$\mu =0,\sigma =1$时的正态分布是标准正态分布。其概率密度函数为：

  其中 $\mu$和 $\sigma$分别表示均值和方差，这个概率密度函数曲线画出来就是我们熟悉的钟形曲线，均值和方差唯一地决定了曲线的形状。
  一元高斯分布可以推广到多元高斯分布，前提假设各变量是相互独立的

  其中${\mu }_1,{\mu }_2,...$和${\sigma }_1,{\sigma }_2,...$分别是第 1 维、第二维… 的均值和方差。对上式向量和矩阵表示上式，令


  有：

  带入(2)：

  其中$\mu ϵR^{n\cdot n}$是均值向量，$KϵR^{n\cdot n}$为协方差矩阵，由于我们假设了各维度直接相互独立，因此 $K$是一个对角矩阵。在各维度变量相关时，上式的形式仍然一致，但此时协方差矩阵$K$不再是对角矩阵，只具备半正定和对称的性质。上式通常也简写为

  同样的高斯过程也可以扩充到无限维，这样便形成了高斯过程。

②高斯过程

假设我们在周一到周四每天的 7:00 测试了 4 次心率，如下图中 4 个点，可能的高斯分布如图所示（高瘦的那条）。这是一个一元高斯分布，只有每天 7: 00 的心率这个维度。

现在考虑不仅在每天的 7: 00 测心率（横轴），在 8:00 时也进行测量（纵轴），这个时候变成两个维度（二元高斯分布），如下图所示

更进一步，如果我们在每天的无数个时间点都进行测量，则变成了下图的情况。注意下图中把测量时间作为横轴，则每个颜色的一条线代表一个（无限个时间点的测量）无限维的采样。当对每次对无限维进行采样得到无限多个点时，其实可以理解为我们采样得到了一个函数。

当从函数的视角去看待采样，理解了每次采样无限维相当于采样一个函数之后，原本的概率密度函数不再是点的分布 ，而变成了函数的分布。这个无限元高斯分布即称为高斯过程。高斯过程正式地定义为：对于所有 $x=[x_1,x_2,...,x_n]$，$f(x)=[f(x_1),f(x_2),...,f(x_n)]$都服从多元高斯分布，则称 [公式] 是一个高斯过程，表示为:

这里 $\mu (x):R^n\to R^n$表示均值函数（Mean function），返回各个维度的均值； $K(x,x):R^n\ast R^n\to R^{n\ast n}$为协方差函数 Covariance Function（也叫核函数 Kernel Function）返回两个向量各个维度之间的协方差矩阵。一个高斯过程为一个均值函数和协方差函数唯一地定义，并且一个高斯过程的有限维度的子集都服从一个多元高斯分布（为了方便理解，可以想象二元高斯分布两个维度各自都服从一个高斯分布）。

2.高斯过程代码实现

①源码实现

from scipy.optimize import minimize


from scipy.optimize import minimize


class GPR:

    def __init__(self, optimize=True):
        self.is_fit = False
        self.train_X, self.train_y = None, None
        self.params = {"l": 0.5, "sigma_f": 0.2}
        self.optimize = optimize

    def fit(self, X, y):
        # store train data
        self.train_X = np.asarray(X)
        self.train_y = np.asarray(y)

         # hyper parameters optimization
        def negative_log_likelihood_loss(params):
            self.params["l"], self.params["sigma_f"] = params[0], params[1]
            Kyy = self.kernel(self.train_X, self.train_X) + 1e-8 * np.eye(len(self.train_X))
            loss = 0.5 * self.train_y.T.dot(np.linalg.inv(Kyy)).dot(self.train_y) + 0.5 * np.linalg.slogdet(Kyy)[1] + 0.5 * len(self.train_X) * np.log(2 * np.pi)
            return loss.ravel()

        if self.optimize:
            res = minimize(negative_log_likelihood_loss, [self.params["l"], self.params["sigma_f"]],
                   bounds=((1e-4, 1e4), (1e-4, 1e4)),
                   method='L-BFGS-B')
            self.params["l"], self.params["sigma_f"] = res.x[0], res.x[1]

        self.is_fit = True
def y(x, noise_sigma=0.0):
    x = np.asarray(x)
    y = np.cos(x) + np.random.normal(0, noise_sigma, size=x.shape)
    return y.tolist()

train_X = np.array([3, 1, 4, 5, 9]).reshape(-1, 1)
train_y = y(train_X, noise_sigma=1e-4)
test_X = np.arange(0, 10, 0.1).reshape(-1, 1)

gpr = GPR()
gpr.fit(train_X, train_y)
mu, cov = gpr.predict(test_X)
test_y = mu.ravel()
uncertainty = 1.96 * np.sqrt(np.diag(cov))
plt.figure()
plt.title("l=%.2f sigma_f=%.2f" % (gpr.params["l"], gpr.params["sigma_f"]))
plt.fill_between(test_X.ravel(), test_y + uncertainty, test_y - uncertainty, alpha=0.1)
plt.plot(test_X, test_y, label="predict")
plt.scatter(train_X, train_y, label="train", c="red", marker="x")
plt.legend()

②scikit-learn库实现

scikit-learn可以实现高斯过程回归和高斯过程分类

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF

# fit GPR
kernel = ConstantKernel(constant_value=0.2, constant_value_bounds=(1e-4, 1e4)) * RBF(length_scale=0.5, length_scale_bounds=(1e-4, 1e4))
gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=2)
gpr.fit(train_X, train_y)
mu, cov = gpr.predict(test_X, return_cov=True)
test_y = mu.ravel()
uncertainty = 1.96 * np.sqrt(np.diag(cov))

# plotting
plt.figure()
plt.title("l=%.1f sigma_f=%.1f" % (gpr.kernel_.k2.length_scale, gpr.kernel_.k1.constant_value))
plt.fill_between(test_X.ravel(), test_y + uncertainty, test_y - uncertainty, alpha=0.1)
plt.plot(test_X, test_y, label="predict")
plt.scatter(train_X, train_y, label="train", c="red", marker="x")
plt.legend()