GP-Circular

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
import theano.tensor as tt

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 42
np.random.seed(RANDOM_SEED)
az.style.use("arviz-darkgrid")

圆形域对于高斯过程是一个挑战。

• 假定周期性模式，但使用原语很难捕获

• 对于圆形域 \([0, \pi)\) 如何建模 \(\pi-\varepsilon\) 和 \(\varepsilon\) 之间的相关性，真实距离是 \(2\varepsilon\)，但如果您只是将其视为非圆形 \((\pi-\varepsilon) - \varepsilon\)，则计算方式不同

• 对于正确计算的距离，我们需要验证我们获得的核是否是正定的

需要另一种方法。

在以下论文中，Weinland 函数用于解决问题，并确保圆形域（不仅仅是圆形域）上的正定核。

\[ W_c(t) = \left(1 + \tau \frac{t}{c}\right)\left(1-\frac{t}{c}\right)_+^\tau \]

其中 \(c\) 是 \(t\) 的最大值，而 \(\tau\ge 4\) 是一些正数

圆上测地线距离（弧长）的核看起来像

\[ k_g(x, y) = W_\pi(\text{dist}_{\mathit{geo}}(x, y)) \]

简而言之，您可以认为

• \(t\) 是时间，它从 \(0\) 运行到 \(24\)，然后返回到 \(0\)

• \(c\) 是任何时间戳之间的最大距离，这里将是 \(12\)

• \(\tau\) 与相关强度成正比。 让我们看看有多少！

在 python 中，Weinland 函数的实现如下

def weinland(t, c, tau=4):
    return (1 + tau * t / c) * np.clip(1 - t / c, 0, np.inf) ** tau

我们还需要圆形域上距离的实现

def angular_distance(x, y, c):
    # https://stackoverflow.com/questions/1878907/the-smallest-difference-between-2-angles
    return (x - y + c) % (c * 2) - c

C = np.pi
x = np.linspace(0, C)

让我们可视化 Weinland 函数是什么，以及它如何影响核

plt.figure(figsize=(16, 9))
for tau in range(4, 10):
    plt.plot(x, weinland(x, C, tau), label=f"tau={tau}")
plt.legend()
plt.ylabel("K(x, y)")
plt.xlabel("dist");

正如我们所见，\(\tau\) 越高，样本的相关性越小

此外，让我们验证我们的圆形距离函数是否按预期工作

plt.plot(
    np.linspace(0, 10 * np.pi, 1000),
    abs(angular_distance(np.linspace(0, 10 * np.pi, 1000), 1.5, C)),
)
plt.ylabel(r"$\operatorname{dist}_{geo}(1.5, x)$")
plt.xlabel("$x$");

在 pymc3 中，我们将使用 pm.gp.cov.Circular 来建模圆形函数

angles = np.linspace(0, 2 * np.pi)
observed = dict(x=np.random.uniform(0, np.pi * 2, size=5), y=np.random.randn(5) + 4)


def plot_kernel_results(Kernel):
    """
    To check for many kernels we leave it as a parameter
    """
    with pm.Model() as model:
        cov = Kernel()
        gp = pm.gp.Marginal(pm.gp.mean.Constant(4), cov_func=cov)
        lik = gp.marginal_likelihood("x_obs", X=observed["x"][:, None], y=observed["y"], noise=0.2)
        mp = pm.find_MAP()
        # actual functions
        y_sampled = gp.conditional("y", angles[:, None])
        # GP predictions (mu, cov)
        y_pred = gp.predict(angles[:, None], point=mp)
        trace = pm.sample_posterior_predictive([mp], var_names=["y"], samples=100)
    plt.figure(figsize=(9, 9))
    paths = plt.polar(angles, trace["y"].T, color="b", alpha=0.05)
    plt.scatter(observed["x"], observed["y"], color="r", alpha=1, label="observations")
    plt.polar(angles, y_pred[0], color="black")
    plt.fill_between(
        angles,
        y_pred[0] - np.diag(y_pred[1]) ** 0.5,
        y_pred[0] + np.diag(y_pred[1]) ** 0.5,
        color="gray",
        alpha=0.5,
        label=r"$\mu\pm\sigma$",
    )
    plt.fill_between(
        angles,
        y_pred[0] - np.diag(y_pred[1]) ** 0.5 * 3,
        y_pred[0] + np.diag(y_pred[1]) ** 0.5 * 3,
        color="gray",
        alpha=0.25,
        label=r"$\mu\pm3\sigma$",
    )
    plt.legend()

def circular():
    tau = pm.Deterministic("τ", pm.Gamma("_τ", alpha=2, beta=1) + 4)
    cov = pm.gp.cov.Circular(1, period=2 * np.pi, tau=tau)
    return cov

plot_kernel_results(circular)

100.00% [7/7 00:00<00:00 logp = -6.9811, ||grad|| = 0.94385]

100.00% [100/100 00:01<00:00]

另一种解决方案是周期核。

注意:

• 在周期核中，控制点之间相关性的关键参数是 ls

• 在圆形核中，它是 tau，添加 ls 参数没有意义，因为它会抵消

基本上，这些核之间几乎没有区别，只是建模相关性的方式不同。

def periodic():
    ls = pm.Gamma("ℓ", alpha=2, beta=1)
    return pm.gp.cov.Periodic(1, 2 * np.pi, ls=ls)

plot_kernel_results(periodic)

100.00% [13/13 00:00<00:00 logp = -7.1213, ||grad|| = 0.0066278]

100.00% [100/100 00:00<00:00]

从模拟中，我们看到圆形核导致更不确定的后验。

让我们看看指数核如何失败

def rbf():
    ls = pm.Gamma("ℓ", alpha=2, beta=1)
    return pm.gp.cov.Exponential(1, ls=ls)

plot_kernel_results(rbf)

100.00% [6/6 00:00<00:00 logp = -7.4099, ||grad|| = 1.8782]

100.00% [100/100 00:01<00:00]

结果看起来与我们使用圆形核得到的结果相似，但变化点 \(0^\circ\) 未考虑在内

结论

• 一旦您坚信函数应平滑地穿过循环边界，请使用圆形/周期核

• 周期核与圆形核一样好，只是后者允许更多不确定性

• RBF 核不是正确的选择

%load_ext watermark
%watermark -n -u -v -iv -w

pymc3 3.9.3
numpy 1.19.0
arviz 0.9.0
last updated: Sat Oct 10 2020 

CPython 3.8.6
IPython 7.17.0
watermark 2.0.2