高斯过程：潜在变量实现#

The gp.Latent 类是高斯过程的直接实现，没有近似。给定均值和协方差函数，我们可以在函数 \(f(x)\) 上放置一个先验，

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

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

The `.prior` 方法#

The 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}\]

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

The `.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

分类

标签

高斯过程：潜在变量实现#

The `.prior` 方法#

The `.conditional` 方法#

示例 1：带有学生-T 分布噪声的回归#

在 PyMC 中编写模型代码#

结果#

使用`.conditional`进行预测#

示例 2：分类#

作者#

水印#

分类

标签

高斯过程：潜在变量实现#

The .prior 方法#

The .conditional 方法#

示例 1：带有学生-T 分布噪声的回归#

在 PyMC 中编写模型代码#

结果#

使用.conditional进行预测#

示例 2：分类#

作者#

水印#

The `.prior` 方法#

The `.conditional` 方法#

使用`.conditional`进行预测#