高斯过程：隐变量实现#

gp.Latent 类是高斯过程的直接实现，无需近似。给定均值和协方差函数，我们可以对函数 \(f(x)\) 设置先验，

\[ f(x) \sim \mathcal{GP}(m(x),\, k(x, x')) \,. \]

它被称为“隐变量”，因为 GP 本身作为隐变量包含在模型中，它不像 gp.Marginal 那样被边缘化处理。与 gp.Latent 不同，您不会在 gp.Marginal 的迹中找到来自 GP 后验的样本。这是 GP 最直接的实现，因为它不假设特定的似然函数或数据或协方差矩阵中的结构。

`.prior` 方法#

prior 方法将多元正态先验分布添加到 PyMC 模型中，分布对象是 GP 函数值向量 \(\mathbf{f}\)，

\[ \mathbf{f} \sim \text{MvNormal}(\mathbf{m}_{x},\, \mathbf{K}_{xx}) \,, \]

其中向量 \(\mathbf{m}_x\) 和矩阵 \(\mathbf{K}_{xx}\) 是在输入 \(x\) 上评估的均值向量和协方差矩阵。默认情况下，PyMC 通过使用协方差矩阵的 Cholesky 因子旋转 f 上的先验来进行重参数化。这通过减少变换后的随机变量 v 的后验协方差来改进采样。重参数化模型为：

\[\begin{split} \begin{aligned} \mathbf{v} \sim \text{N}(0, 1)& \\ \mathbf{L} = \text{Cholesky}(\mathbf{K}_{xx})& \\ \mathbf{f} = \mathbf{m}_{x} + \mathbf{Lv} \\ \end{aligned} \end{split}\]

有关此重参数化的更多信息，请参阅关于从多元分布中抽取值的部分。

`.conditional` 方法#

conditional 方法实现了函数值的预测分布，这些函数值不是原始数据集的一部分。此分布为：

\[ \mathbf{f}_* \mid \mathbf{f} \sim \text{MvNormal} \left( \mathbf{m}_* + \mathbf{K}_{*x}\mathbf{K}_{xx}^{-1} \mathbf{f} ,\, \mathbf{K}_{**} - \mathbf{K}_{*x}\mathbf{K}_{xx}^{-1}\mathbf{K}_{x*} \right) \]

使用我们上面定义的相同 gp 对象，我们可以通过以下方式构造具有此分布的随机变量：

# vector of new X points we want to predict the function at
X_star = np.linspace(0, 2, 100)[:, None]

with latent_gp_model:
    f_star = gp.conditional("f_star", X_star)

示例 2：分类#

首先，我们使用 GP 生成一些遵循伯努利分布的数据，其中 \(p\)，即获得 1 而不是 0 的概率是 \(x\) 的函数。我重置了种子并添加了更多伪造数据点，因为模型可能难以辨别 0.5 附近的变异，而观测值很少。

# reset the random seed for the new example
RANDOM_SEED = 8888
rng = np.random.default_rng(RANDOM_SEED)

# number of data points
n = 300

# x locations
x = np.linspace(0, 10, n)

# true covariance
ell_true = 0.5
eta_true = 1.0
cov_func = eta_true**2 * pm.gp.cov.ExpQuad(1, ell_true)
K = cov_func(x[:, None]).eval()

# zero mean function
mean = np.zeros(n)

# sample from the gp prior
f_true = pm.draw(pm.MvNormal.dist(mu=mean, cov=K), 1, random_seed=rng)

# Sample the GP through the likelihood
y = pm.Bernoulli.dist(p=pm.math.invlogit(f_true)).eval()

fig = plt.figure(figsize=(10, 4))
ax = fig.gca()

ax.plot(x, pm.math.invlogit(f_true).eval(), "dodgerblue", lw=3, label="True rate")
# add some noise to y to make the points in the plot more visible
ax.plot(x, y + np.random.randn(n) * 0.01, "kx", ms=6, label="Observed data")

ax.set_xlabel("X")
ax.set_ylabel("y")
ax.set_xlim([0, 11])
plt.legend(loc=(0.35, 0.65), frameon=True);

../_images/7cbdd65e8062052021b15738c704a2958553a62b2adb54315b394fe300a02943.png

with pm.Model() as model:
    ell = pm.InverseGamma("ell", mu=1.0, sigma=0.5)
    eta = pm.Exponential("eta", lam=1.0)
    cov = eta**2 * pm.gp.cov.ExpQuad(1, ell)

    gp = pm.gp.Latent(cov_func=cov)
    f = gp.prior("f", X=x[:, None])

    # logit link and Bernoulli likelihood
    p = pm.Deterministic("p", pm.math.invlogit(f))
    y_ = pm.Bernoulli("y", p=p, observed=y)

    idata = pm.sample(1000, chains=2, cores=2, nuts_sampler="numpyro")

We recommend running at least 4 chains for robust computation of convergence diagnostics

# check Rhat, values above 1 may indicate convergence issues
n_nonconverged = int(np.sum(az.rhat(idata)[["eta", "ell", "f_rotated_"]].to_array() > 1.03).values)
if n_nonconverged == 0:
    print("No Rhat values above 1.03, \N{check mark}")
else:
    print(f"The MCMC chains for {n_nonconverged} RVs appear not to have converged.")

No Rhat values above 1.03, ✓

ax = az.plot_pair(
    idata,
    var_names=["eta", "ell"],
    kind=["kde", "scatter"],
    scatter_kwargs={"color": "darkslategray", "alpha": 0.4},
    gridsize=25,
    divergences=True,
)

ax.axvline(x=eta_true, color="dodgerblue")
ax.axhline(y=ell_true, color="dodgerblue");

../_images/d442fa2b2bbf8af2362dabefac7ea3039dcbb72068b5b3d30ebc99176f755c44.png

n_pred = 200
X_new = np.linspace(0, 12, n_pred)[:, None]

with model:
    f_pred = gp.conditional("f_pred", X_new, jitter=1e-4)
    p_pred = pm.Deterministic("p_pred", pm.math.invlogit(f_pred))

with model:
    idata.extend(pm.sample_posterior_predictive(idata.posterior, var_names=["f_pred", "p_pred"]))

Sampling: [f_pred]

# plot the results
fig = plt.figure(figsize=(10, 4))
ax = fig.gca()

# plot the samples from the gp posterior with samples and shading
p_pred = az.extract(idata.posterior_predictive, var_names="p_pred").transpose("sample", ...)
plot_gp_dist(ax, p_pred, X_new)

# plot the data (with some jitter) and the true latent function
plt.plot(x, pm.math.invlogit(f_true).eval(), "dodgerblue", lw=3, label="True f")
plt.plot(
    x,
    y + np.random.randn(y.shape[0]) * 0.01,
    "kx",
    ms=6,
    alpha=0.5,
    label="Observed data",
)

# axis labels and title
plt.xlabel("X")
plt.ylabel("True f(x)")
plt.xlim([0, 12])
plt.title("Posterior distribution over $f(x)$ at the observed values")
plt.legend(loc=(0.32, 0.65), frameon=True);

../_images/349b027dc4df637d4a7867a8945a47149339694a42841034355b699f00a6e48a.png

作者#

由 Bill Engels 于 2017 年创建 (pymc#1674)
由 Colin Caroll 于 2019 年重新执行 (pymc#3397)
由 Bill Engels 于 2022 年 9 月更新为 V4 版本 (pymc-examples#237)
由 Chris Fonnesbeck 于 2023 年 7 月更新为 V5 版本 (pymc-examples#549)
由 Alexandre Andorra 于 2024 年 5 月更新

水印#

%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,aeppl,xarray

Last updated: Mon May 27 2024

Python implementation: CPython
Python version       : 3.12.2
IPython version      : 8.22.2

pytensor: 2.20.0
aeppl   : not installed
xarray  : 2024.3.0

matplotlib: 3.8.3
numpy     : 1.26.4
pymc      : 5.15.0+14.gfd11cf012
arviz     : 0.17.1

Watermark: 2.4.3

许可声明#

此示例库中的所有笔记本均根据 MIT 许可证提供，该许可证允许修改和再分发以用于任何用途，前提是保留版权和许可声明。

引用 PyMC 示例#

要引用此笔记本，请使用 Zenodo 为 pymc-examples 存储库提供的 DOI。

重要提示

许多笔记本改编自其他来源：博客、书籍……在这种情况下，您也应该引用原始来源。

另请记住引用您的代码使用的相关库。

这是一个 bibtex 中的引用模板

@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}

渲染后可能如下所示

分类

标签

高斯过程：隐变量实现#

`.prior` 方法#

`.conditional` 方法#

示例 1：具有 Student-T 分布噪声的回归#

在 PyMC 中编写模型代码#

结果#

使用 `.conditional` 进行预测#

示例 2：分类#

作者#

水印#

许可声明#

引用 PyMC 示例#

分类

标签

高斯过程：隐变量实现#

.prior 方法#

.conditional 方法#

示例 1：具有 Student-T 分布噪声的回归#

在 PyMC 中编写模型代码#

结果#

使用 .conditional 进行预测#

示例 2：分类#

作者#

水印#

许可声明#

引用 PyMC 示例#

`.prior` 方法#

`.conditional` 方法#

使用 `.conditional` 进行预测#