高斯过程（一）_高斯过程例题-CSDN博客

本文链接：https://blog.csdn.net/gaofeipaopaotang/article/details/124225012

高斯分布

$\begin{bmatrix} X_1 \\ X_2 \\ \vdots \\X_n \end {bmatrix} \sim \mathcal{N}(\mu, \Sigma)$

称 $X$ 服从高斯分布， $\mu$ 表示n维高斯分布的期望。

$\Sigma = Cov(X_i,X_j) = E[(X_i-\mu_i)(X_j-\mu_j)]$ ，表示n维高斯分布的协方差矩阵。协方差矩阵是一个对称矩阵，矩阵的对角线元素决定 $X$ 对应维度的方差，非对角线元素建模对一个两个维度之间的相关性。

高斯分布具有在条件化和边缘化下闭合的良好代数特性。在条件化和边缘化下关闭意味着这些操作的结果分布也是高斯分布：

$P_{X,Y} = \begin{bmatrix} X \\ Y \end {bmatrix} \sim \mathcal{N}(\mu, \Sigma) = \mathcal{N}(\begin{bmatrix} \mu_X \\\mu_Y \end{bmatrix}, \begin{bmatrix} \Sigma_{XX},\Sigma_{XY} \\ \Sigma_{YX},\Sigma_{YY} \end{bmatrix})$

边缘分布：

$X\sim\mathcal{N}(\mu_X,\Sigma_{XX})$

$Y\sim \mathcal{N}(\mu_Y,\Sigma_{YY})$

条件分布：

$X|Y\sim \mathcal{N}(\mu_X + \Sigma_{XY}\Sigma_{YY}^{-1}(Y-\mu_Y), \Sigma_{XX} - \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX})$

$Y|X\sim \mathcal{N}(\mu_Y+\Sigma_{YX}\Sigma_{XX}^{-1}(X-\mu_X), \Sigma_{YY}-\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY})$

高斯过程

随机过程，比如高斯随机过程，本质上是一组随机变量。或则说，GP 可以看作是一个随机过程，对应于无限维高斯分布。高斯随机过程定义了函数的概率分布，每一个采样得到的样本，对应着一组函数值的实例。

下面通过一个例子说明多元高斯分布于函数之间的关系：

import math
import numpy as np
import matplotlib.pyplot as plt

def kernel(i, j):
    if i == j:
        return 0.5
    else:
        return 0.0

dim = 2
mean = [0 for i in range(dim)]
index = [i + 1 for i in range(dim)]
# 构造协方差矩阵
cov = []
for i in range(dim):
    p = []
    for j in range(dim):
        sigma = kernel(i, j)
        p.append(sigma)
    cov.append(p)
fig, ax = plt.subplots(figsize=(32, 8))
for sample in range(3):
    values = np.random.multivariate_normal(mean, cov, 1).T
    ax.plot(index, values, '-')

请添加图片描述

例子中定义了一个均值为0、方差为0.5、不同维度协方差为0的2元高斯分布
从分布中采样得到3个样本，将样本点的维度作为x轴，维度对应的值作为y轴，因此3个样本点对应3条直线

下面将样本高斯分布的维度扩展到100维(dim = 100)

...
dim = 100
...

![请添加图片描述](

可以看出函数不是很平滑，这是由于我们假设不同维度之间相互独立导致的，因此我们对维度之间的相关性加入一些假设，比如相近的位置之间相关性加到，距离较远的维度之间相关性较小。通过修改协方差举证函数的生成方式引入此假设：

import math
import numpy as np
import matplotlib.pyplot as plt
    
def kernel_rbf(i, j):
    # 新的协方差矩阵生成函数：维度之间不在相互独立，并且距离相近的维度之间相关性较大，
    # 距离较远的维度之间相关性较小
    sigma = 0.5
    l = 4.0
    return math.pow(sigma, 2) * math.exp(0 - math.pow(i - j, 2.0) / (2.0 * math.pow(l, 2.0)))


dim = 100
mean = [0 for i in range(dim)]
index = [i + 1 for i in range(dim)]
cov = []
for i in range(dim):
    p = []
    for j in range(dim):
        sigma = kernel_rbf(i, j)
        p.append(sigma)
    cov.append(p)
fig, ax = plt.subplots(figsize=(32, 8))
for sample in range(3):
    values = np.random.multivariate_normal(mean, cov, 1).T
    ax.plot(index, values, '-')

请添加图片描述

函数图形变得更加的平滑，到此我们对高斯过程和函数 $f$ 之间的对应关系，有了一个直观的认知。

最后我们可将多元高斯分布的维度从正整数扩展到实数 $R$ ，对应于无限维高斯分布，每个样本对应着一个函数，因此高斯过程可以视为函数空间的分布。

...
index = np.linspace(0, 1, dim)
...

请添加图片描述

高斯分布由期望和协方差矩阵唯一确定，不同的均值和协方差矩阵的选择，对应这对函数空间的不同假设。上面例子中采用的协方差生产函数kernel_rbf被称为高斯过程的核，常用的和包括：

$k:R^n\times R^n \rightarrow R, \Sigma = Cov(X,X') = k(t,t')$

$\sigma^2exp(-\frac{||t-t'||^2}{2l^2})$

$\sigma^2exp(-\frac{2sin^2(\pi|t-t'|/p)}{l^2})$

$\sigma_b^2 + \sigma^2(t-c)(t'-c)$

先验分布

在没有观测数据的情况下，采用的高斯分布称为函数空间的先验分布。比如上文中采用的两个分布

均值为0，协方差矩阵为对角矩阵：对应的先验假设是均值为0，不同维度之间相互独立
均值为0，核函数为RBF函数：对应的先验假设是均值为0，距离相近的维度之间相关性较大，距离较远的维度之间相关性较小

后验分布

有了部分观侧数据之后，结合观侧数据和先验假设，由贝叶斯估计，计算得到贝叶斯后验估计 $P_{X|Y}$ ：

$P_{X,Y} = \begin{bmatrix} X \\ Y \end {bmatrix} \sim \mathcal{N}(0, \Sigma) = \mathcal{N}(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \Sigma_{XX}，\Sigma_{XY} \\ \Sigma_{YX}，\Sigma_{YY} + \psi^2 I \end{bmatrix})$

$X|Y\sim \mathcal{N}(\mu_X + \Sigma_{XY}\Sigma_{YY}^{-1}(Y-\mu_Y), \Sigma_{XX} - \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX})$

$=\mathcal{N}(\Sigma_{XY}\Sigma_{YY}^{-1}Y, \Sigma_{XX} - \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX})$

下面继续通过上文的例子，说明如果在先验假设中，结合观侧数据，对函数分布进行改进：

def create_cov(x1, x2):
    cov = []
    for i in range(len(x1)):
        p = []
        for j in range(len(x2)):
            sigma = kernel_rbf(x1[i], x2[j])
            p.append(sigma)
        cov.append(p)
    return np.matrix(cov)
    
test_x = np.linspace(0, 1, dim)

# 加入两个观侧样本
train_x = [0.3, 0.7]
train_y = [0.2, 0.1]

cov_xx = create_cov(test_x, test_x)  # 5 * 5
cov_yy = create_cov(train_x, train_x) # 2 * 2
cov_xy = create_cov(test_x, train_x)  # 5 * 2
cov_yx = create_cov(train_x, test_x)  # 2 * 25

np.matrix(train_y).T # 2 * 1

mean_post = cov_xy * cov_yy.I * np.matrix(train_y).T
cov_post = cov_xx - cov_xy * cov_yy.I * cov_yx

fig, ax = plt.subplots(figsize=(32, 8))
ax.plot(train_x, train_y, '.', markersize=32)

mean_post = np.reshape(mean_post, (1, -1)).tolist()[0]
for r in range(3):
    ys = np.random.multivariate_normal(mean_post, cov_post, 1).T
    ax.plot(test_x, ys, '-')

请添加图片描述

增大采样的次数为100

请添加图片描述