高斯过程#

GP 基础#

有时，模型中的未知参数或变量不是一个标量值或固定长度的向量，而是一个函数。高斯过程（GP）可以用作先验概率分布，其支持度覆盖连续函数的空间。函数 $f(x)$ 上的 GP 先验通常写作，

\[ \begin{align}\begin{aligned}f(x) \sim \mathcal{GP}(m(x), \, k(x, x')) \,.\\f(x) 服从高斯过程 \mathcal{GP}(m(x), \, k(x, x')) \,.\end{aligned}\end{align} \]

函数值被建模为从多元正态分布中抽取的样本，该分布由均值函数 $m(x)$ 和协方差函数 $k(x, x')$ 参数化。高斯过程由于多元正态分布的边缘化和条件化特性，成为函数先验的方便选择。通常，在推断步骤中评估 $f(x)$ 的边缘分布。然后使用条件分布来预测新点 $x_*$ 处的函数值 $f(x_*)$。

联合分布 $f(x)$ 和 $f(x_*)$ 是多元正态分布，

\[\begin{split}\begin{bmatrix} f(x) \ f(x_*) \ \end{bmatrix} \sim \text{N}\left( \begin{bmatrix} m(x) \ m(x_*) \ \end{bmatrix} \,, \begin{bmatrix} k(x,x') & k(x_*, x) \\ k(x_*, x) & k(x_*, x_*') \ \end{bmatrix} \right) \,.\end{split}\]

从联合分布出发，可以得到 $f(x)$ 的边缘分布，即 $\text{N}(m(x),\, k(x, x'))$。条件分布为

\[ \begin{align}\begin{aligned}f(x_*) \mid f(x) \sim \text{N}\left( k(x_*, x) k(x, x)^{-1} [f(x) - m(x)] + m(x_*) ,\, k(x_*, x_*) - k(x, x_*) k(x, x)^{-1} k(x, x_*) \right) \,.\\f(x_*) \mid f(x) \sim \text{N}\left( k(x_*, x) k(x, x)^{-1} [f(x) - m(x)] + m(x_*) ,\, k(x_*, x_*) - k(x, x_*) k(x, x)^{-1} k(x, x_*) \right) \,.\end{aligned}\end{align} \]

备注

关于GPs的更多信息，请参阅Rasmussen & Williams所著的《Gaussian Processes for Machine Learning <http://www.gaussianprocess.org/gpml/>`_》，或D. Mackay的`这篇介绍 <http://www.inference.org.uk/mackay/gpB.pdf>`_。

PyMC 是一个非常适合使用全贝叶斯高斯过程模型的工作环境。PyMC 中的高斯过程具有清晰的语法且高度可组合，并且包含许多预定义的协方差函数（或核）、均值函数以及几种高斯过程实现。高斯过程被视为可以在更大或层次模型中使用的分布，而不仅仅是作为独立的回归模型。

均值和协方差函数#

那些使用过 GPy 或 GPflow Python 包的人会发现构造均值和协方差函数的语法有些熟悉。当首次实例化时，均值和协方差函数被参数化，但尚未给出它们的输入。协方差函数还必须提供输入矩阵的维度数量，以及一个索引列表，该列表指示它们将操作哪些维度。这种设计的原因是为了使协方差函数可以由其他协方差函数的组合构成。

例如，要构建一个在表示三个预测变量的三列矩阵的第二列和第三列上操作的指数二次协方差函数:

ls = [2, 5] # the lengthscales
cov_func = pm.gp.cov.ExpQuad(input_dim=3, ls=ls, active_dims=[1, 2])

这里，ls 参数，即长度尺度，是二维的，允许第二和第三维度具有不同的长度尺度。我们必须指定 input_dim，即 X 的总列数，以及 active_dims，即协方差函数将作用于哪些列或维度，原因是 cov_func 实际上还没有看到输入数据。active_dims 参数是可选的，默认情况下为输入矩阵的所有列。

PyMC 中的协方差函数紧密遵循核的代数规则，这允许用户将协方差函数组合成新的函数，例如：

两个协方差函数的和也是一个协方差函数:

cov_func = pm.gp.cov.ExpQuad(...) + pm.gp.cov.ExpQuad(...)

两个协方差函数的乘积也是一个协方差函数:

cov_func = pm.gp.cov.ExpQuad(...) * pm.gp.cov.Periodic(...)

协方差函数与标量的乘积（或和）是一个协方差函数:
```
cov_func = eta**2 * pm.gp.cov.Matern32(...)
```

协方差函数定义后，现在可以通过调用 cov_func(x, x)`（或 :code:`mean_func(x)）来评估它。由于 PyMC 是基于 PyTensor 构建的，因此定义和实验非标准的协方差和均值函数相对容易。更多信息请查看关于协方差函数的教程。

GP 实现#

PyMC 包含几种 GP 实现，包括边际和潜变量模型，以及一些快速近似。它们的使用都遵循类似的模式：首先，用均值函数和协方差函数实例化一个 GP。然后，GP 对象可以相加，使得函数的特征可以被仔细建模和分离。最后，在 GP 对象上调用 prior、marginal_likelihood 或 conditional 方法之一，以实际构建表示函数先验的 PyMC 随机变量。

使用 gp.Latent 作为示例，首先指定 GP 的语法是:

gp = pm.gp.Latent(mean_func, cov_func)

第一个参数是均值函数，第二个是协方差函数。我们已经创建了GP对象，但我们还没有明确它将作为哪个函数的先验，输入是什么，或者它将基于哪些参数进行条件化。

备注

The gp.Marginal 类及其类似类没有 prior 方法。相反，它们有一个 marginal_likelihood 方法，该方法的使用方式类似，但有额外的必需参数，例如观测数据、噪声或其他，具体取决于实现。请参阅笔记本以获取示例。conditional 方法的工作方式类似。

调用 prior 方法将创建一个 PyMC 随机变量，该变量表示潜在函数 $f(x) = \mathbf{f}$:

f = gp.prior("f", X)

f 是一个可以在 PyMC 模型中像其他类型的随机变量一样使用的随机变量。第一个参数是我们正在为其放置先验的函数的随机变量的名称。第二个参数是先验所在的函数的输入，X。输入通常是已知并存在于数据中的，但它们也可以是 PyMC 随机变量。如果输入是 PyTensor 张量或 PyMC 随机变量，则需要给出 shape。

通常在这个时候，会对模型进行推断。conditional 方法在任意 $x_*$ 输入点上创建潜在函数的条件分布，即 $f(x_*)$。为了构建条件分布，我们编写如下:

f_star = gp.conditional("f_star", X_star)

加性高斯过程#

PyMC 中的 GP 实现构建得非常易于定义加性 GP 并从单个 GP 组件中采样。我们可以写:

gp1 = pm.gp.Marginal(mean_func1, cov_func1)
gp2 = pm.gp.Marginal(mean_func2, cov_func2)
gp3 = gp1 + gp2

GP 对象必须具有相同的类型，gp.Marginal 不能添加到 gp.Latent。

考虑两个独立的高斯过程分布函数，$f_1(x) \sim \mathcal{GP}\left(m_1(x),\, k_1(x, x')\right)$ 和 $f_2(x) \sim \mathcal{GP}\left( m_2(x),\, k_2(x, x')\right)$。$f_1,\, f_1^*,\, f_2,\, f_2^*,\, f_1 + f_2$ 和 $f_1^* + f_2^*$ 的联合分布为

\[\begin{split}\begin{bmatrix} f_1 \ f_1^* \ f_2 \ f_2^* \ f_1 + f_2 \ f_1^* + f_2^* \end{bmatrix} \sim \text{N}\left( \begin{bmatrix} m_1 \ m_1^* \ m_2 \ m_2^* \ m_1 + m_2 \ m_1^* + m_2^* \end{bmatrix} \,,\begin{bmatrix} K_1 & K_1^* & 0 & 0 & K_1 & K_1^* \\ K_1^{*^T} & K_1^{**} & 0 & 0 & K_1^* & K_1^{**} \\ 0 & 0 & K_2 & K_2^* & K_2 & K_2^{*} \\ 0 & 0 & K_2^{*^T} & K_2^{**} & K_2^{*} & K_2^{**} \\ K_1 & K_1^{*} & K_2 & K_2^{*} & K_1 + K_2 & K_1^{*} + K_2^{*} \\ K_1^{*^T} & K_1^{**} & K_2^{*^T} & K_2^{**} & K_1^{*^T}+K_2^{*^T} & K_1^{**}+K_2^{**} \end{bmatrix} \right) \,.\end{split}\]

使用联合分布来获得 $f_1^*$ 的条件分布，其中 $f_1 + f_2$ 的贡献已被剔除，我们得到

\[$f_1^* \mid f_1 + f_2 \sim \text{N}\left( m_1^* + K_1^{*^T}(K_1 + K_2)^{-1}\left[f_1 + f_2 - m_1 - m_2\right] \,,\, K_1^{**} - K_1^{*^T}(K_1 + K_2)^{-1}K_1^* \right) \,.\]

这些方程展示了如何将GP模型分解为各个组成部分，以了解每个部分对数据的贡献。更多信息，请查看 David Duvenaud的博士论文。

PyMC 中的 GP 对象会自动跟踪这些边际。以下代码草图展示了如何定义 $f_2^*$ 的条件分布。我们在示例中使用了 gp.Marginal，但对于其他实现也是如此。第一个代码块拟合了 GP 先验。为了简洁起见，我们将 $f_1 + f_2$ 表示为 $f$:

with pm.Model() as model:
    gp1 = pm.gp.Marginal(mean_func1, cov_func1)
    gp2 = pm.gp.Marginal(mean_func2, cov_func2)

    # gp represents f1 + f2.
    gp = gp1 + gp2

    f = gp.marginal_likelihood("f", X, y, sigma)

    idata = pm.sample(1000)

要构建 gp1 或 gp2 的条件分布，我们还需要包括额外的参数，X、y 和 sigma:

with model:
    # conditional distributions of f1 and f2
    f1_star = gp1.conditional("f1_star", X_star,
                              given={"X": X, "y": y, "sigma": sigma, "gp": gp})
    f2_star = gp2.conditional("f2_star", X_star,
                              given={"X": X, "y": y, "sigma": sigma, "gp": gp})

    # conditional of f1 + f2, `given` not required
    f_star = gp.conditional("f_star", X_star)

第二个块生成条件分布。注意，对于 $f1$ 和 $f2$ 的条件分布需要额外的参数，而对于 $f$ 则不需要。这是因为当在 gp 上调用 .marginal_likelihood 时，这些参数会被缓存。

备注

在构建条件时，必须提供额外的参数 X、y、sigma 和 gp 作为名为 given 的字典！

由于未调用 gp1 或 gp2 的边缘似然方法，因此需要为它们的条件提供所需的输入。与先验类似，f_star、f1_star 和 f2_star 现在是随机变量，可以像 PyMC 中的任何其他随机变量一样使用。

查看笔记本以获取在 PyMC 中使用 GP 功能的详细演示。