Dependent density regression

在另一个 示例 中，我们展示了如何使用狄利克雷过程执行贝叶斯非参数密度估计。此示例在前一个示例的基础上进行了扩展，说明了依赖密度回归。

正如狄利克雷过程混合模型可以被认为是无限混合模型，它选择活动组件的数量作为推断的一部分一样，依赖密度回归可以被认为是无限的 专家混合，它选择活动专家作为推断的一部分。它们的灵活性和模块化使其成为执行非参数贝叶斯数据分析的强大工具。

import arviz as az
import numpy as np
import pandas as pd
import pymc3 as pm
import seaborn as sns

from IPython.display import HTML
from matplotlib import animation as ani
from matplotlib import pyplot as plt
from theano import tensor as tt

print(f"Running on PyMC3 v{pm.__version__}")

Running on PyMC3 v3.11.2

%config InlineBackend.figure_format = 'retina'
plt.rc("animation", writer="ffmpeg")
blue, *_ = sns.color_palette()
az.style.use("arviz-darkgrid")
SEED = 972915  # from random.org; for reproducibility
np.random.seed(SEED)

我们将使用 Larry Wasserman 的优秀著作《All of Nonparametric Statistics》中的 LIDAR 数据集。我们对数据集进行标准化处理，以提高我们样本的收敛速度。

DATA_URI = "http://www.stat.cmu.edu/~larry/all-of-nonpar/=data/lidar.dat"


def standardize(x):
    return (x - x.mean()) / x.std()


df = pd.read_csv(DATA_URI, sep=r"\s{1,3}", engine="python").assign(
    std_range=lambda df: standardize(df.range), std_logratio=lambda df: standardize(df.logratio)
)

df.head()

	范围	对数比	标准化的范围	标准化的对数比
0	390	-0.050356	-1.717725	0.852467
1	391	-0.060097	-1.707299	0.817981
2	393	-0.041901	-1.686447	0.882398
3	394	-0.050985	-1.676020	0.850240
4	396	-0.059913	-1.655168	0.818631

我们在下面绘制 LIDAR 数据。

fig, ax = plt.subplots(figsize=(8, 6))

ax.scatter(df.std_range, df.std_logratio, color=blue)

ax.set_xticklabels([])
ax.set_xlabel("Standardized range")

ax.set_yticklabels([])
ax.set_ylabel("Standardized log ratio");

此数据集具有两个有趣的属性，使其可用于说明依赖密度回归。

1. 范围和对数比之间的关系是非线性的，但具有局部线性分量。

2. 观测噪声是 异方差 的；也就是说，方差的大小随范围变化。

依赖密度回归背后的直观想法是将问题简化为许多（相关的）密度估计，这些密度估计以预测变量的固定值为条件。以下动画说明了这种直觉。

fig, (scatter_ax, hist_ax) = plt.subplots(ncols=2, figsize=(16, 6))

scatter_ax.scatter(df.std_range, df.std_logratio, color=blue, zorder=2)

scatter_ax.set_xticklabels([])
scatter_ax.set_xlabel("Standardized range")

scatter_ax.set_yticklabels([])
scatter_ax.set_ylabel("Standardized log ratio")

bins = np.linspace(df.std_range.min(), df.std_range.max(), 25)

hist_ax.hist(df.std_logratio, bins=bins, color="k", lw=0, alpha=0.25, label="All data")

hist_ax.set_xticklabels([])
hist_ax.set_xlabel("Standardized log ratio")

hist_ax.set_yticklabels([])
hist_ax.set_ylabel("Frequency")

hist_ax.legend(loc=2)

endpoints = np.linspace(1.05 * df.std_range.min(), 1.05 * df.std_range.max(), 15)

frame_artists = []

for low, high in zip(endpoints[:-1], endpoints[2:]):
    interval = scatter_ax.axvspan(low, high, color="k", alpha=0.5, lw=0, zorder=1)
    *_, bars = hist_ax.hist(
        df[df.std_range.between(low, high)].std_logratio, bins=bins, color="k", lw=0, alpha=0.5
    )

    frame_artists.append((interval,) + tuple(bars))

animation = ani.ArtistAnimation(fig, frame_artists, interval=500, repeat_delay=3000, blit=True)
plt.close()
# prevent the intermediate figure from showing

HTML(animation.to_html5_video())

当我们在左图中用沿 x 轴滑动的窗口切片数据时，窗口中点的 y 值的经验分布在右图中发生变化。这种方法的一个重要方面是，对应于预测变量的相近值的密度估计是相似的。

在前面的示例中，我们看到狄利克雷过程将概率密度估计为具有无限多个分量的混合模型。在正态分量分布的情况下，

\[y \sim \sum_{i = 1}^{\infty} w_i \cdot N(\mu_i, \tau_i^{-1}),\]

其中混合权重 \(w_1, w_2, \ldots\) 由 stick-breaking 过程 生成。

依赖密度回归通过允许混合权重和分量均值随预测变量 \(x\) 的值而变化，从而推广了狄利克雷过程混合模型的这种表示形式。也就是说，

\[y\ |\ x \sim \sum_{i = 1}^{\infty} w_i\ |\ x \cdot N(\mu_i\ |\ x, \tau_i^{-1}).\]

在此示例中，我们将遵循《贝叶斯数据分析》第 23 章，并使用 probit stick-breaking 过程来确定条件混合权重 \(w_i\ |\ x\)。probit stick-breaking 过程首先定义

\[v_i\ |\ x = \Phi(\alpha_i + \beta_i x),\]

其中 \(\Phi\) 是标准正态分布的累积分布函数。然后我们通过将 stick breaking 过程应用于 \(v_i\ |\ x\) 来获得 \(w_i\ |\ x\)。也就是说，

\[w_i\ |\ x = v_i\ |\ x \cdot \prod_{j = 1}^{i - 1} (1 - v_j\ |\ x).\]

对于 LIDAR 数据集，我们使用独立的正态先验 \(\alpha_i \sim N(0, 5^2)\) 和 \(\beta_i \sim N(0, 5^2)\)。我们现在使用 PyMC3 来表达条件混合权重的这个模型。

def norm_cdf(z):
    return 0.5 * (1 + tt.erf(z / np.sqrt(2)))


def stick_breaking(v):
    return v * tt.concatenate(
        [tt.ones_like(v[:, :1]), tt.extra_ops.cumprod(1 - v, axis=1)[:, :-1]], axis=1
    )

N = len(df)
K = 20

std_range = df.std_range.values[:, np.newaxis]
std_logratio = df.std_logratio.values

with pm.Model(coords={"N": np.arange(N), "K": np.arange(K) + 1, "one": [1]}) as model:
    alpha = pm.Normal("alpha", 0.0, 5.0, dims="K")
    beta = pm.Normal("beta", 0.0, 5.0, dims=("one", "K"))
    x = pm.Data("x", std_range)
    v = norm_cdf(alpha + pm.math.dot(x, beta))
    w = pm.Deterministic("w", stick_breaking(v), dims=["N", "K"])

我们将 x 定义为 pm.Data 容器，以便稍后使用 PyMC3 的后验预测功能。

虽然依赖密度回归模型理论上具有无限多个分量，但我们必须将模型截断为有限多个分量（在本例中为 20 个），以便使用 PyMC3 表示它。从模型中抽样后，我们将验证截断是否对我们的结果产生了不当影响。

由于 LIDAR 数据似乎有几个线性分量，我们使用线性模型

\[\begin{split} \begin{align*} \mu_i\ |\ x & \sim \gamma_i + \delta_i x \\ \gamma_i & \sim N(0, 10^2) \\ \delta_i & \sim N(0, 10^2) \end{align*} \end{split}\]

用于条件分量均值。

with model:
    gamma = pm.Normal("gamma", 0.0, 10.0, dims="K")
    delta = pm.Normal("delta", 0.0, 10.0, dims=("one", "K"))
    mu = pm.Deterministic("mu", gamma + pm.math.dot(x, delta))

最后，我们将先验 \(\tau_i \sim \textrm{Gamma}(1, 1)\) 放在分量精度上。

with model:
    tau = pm.Gamma("tau", 1.0, 1.0, dims="K")
    y = pm.Data("y", std_logratio)
    obs = pm.NormalMixture("obs", w, mu, tau=tau, observed=y)

pm.model_to_graphviz(model)

我们现在从依赖密度回归模型中抽样。

SAMPLES = 20000
BURN = 10000

with model:
    step = pm.Metropolis()
    trace = pm.sample(SAMPLES, tune=BURN, step=step, random_seed=SEED, return_inferencedata=True)

100.00% [30000/30000 02:56<00:00 链 0 采样完成，0 次发散]

100.00% [30000/30000 02:50<00:00 链 1 采样完成，0 次发散]

为了验证截断是否对我们的结果产生了不当影响，我们绘制了每个分量的最大后验期望混合权重。（在此模型中，每个点都有每个分量的混合权重，因此我们绘制每个分量在所有数据点上的最大混合权重，以便判断该分量是否对后验产生任何影响。）

fig, ax = plt.subplots(figsize=(8, 6))

max_mixture_weights = trace.posterior["w"].mean(("chain", "draw")).max("N")
ax.bar(max_mixture_weights.coords.to_index(), max_mixture_weights)

ax.set_xlim(1 - 0.5, K + 0.5)
ax.set_xticks(np.arange(0, K, 2) + 1)
ax.set_xlabel("Mixture component")

ax.set_ylabel("Largest posterior expected\nmixture weight");

由于只有三个混合分量对任何数据点具有可观的后验期望权重，我们可以相当肯定截断没有对我们的结果产生不当影响。（如果大多数分量都具有可观的后验期望权重，则截断可能会影响结果，我们将增加分量的数量并再次采样。）

从视觉上看，LIDAR 数据具有三个线性分量是合理的，因此这些后验期望权重似乎很好地识别了数据的结构。我们现在从后验预测分布中抽样，以更好地了解模型的性能。

PP_SAMPLES = 5000

lidar_pp_x = np.linspace(std_range.min() - 0.05, std_range.max() + 0.05, 100)

with model:
    pm.set_data({"x": lidar_pp_x[:, np.newaxis]})
    pp_trace = pm.sample_posterior_predictive(trace, PP_SAMPLES, random_seed=SEED)

/opt/conda/lib/python3.7/site-packages/pymc3/sampling.py:1690: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  "samples parameter is smaller than nchains times ndraws, some draws "

100.00% [5000/5000 03:23<00:00]

下面我们绘制后验期望值和 95% 后验可信区间。

fig, ax = plt.subplots(figsize=(8, 6))

ax.scatter(df.std_range, df.std_logratio, color=blue, zorder=10, label=None)

low, high = np.percentile(pp_trace["obs"], [2.5, 97.5], axis=0)
ax.fill_between(
    lidar_pp_x, low, high, color="k", alpha=0.35, zorder=5, label="95% posterior credible interval"
)

ax.plot(lidar_pp_x, pp_trace["obs"].mean(axis=0), c="k", zorder=6, label="Posterior expected value")

ax.set_xticklabels([])
ax.set_xlabel("Standardized range")

ax.set_yticklabels([])
ax.set_ylabel("Standardized log ratio")

ax.legend(loc=1)
ax.set_title("LIDAR Data");

该模型很好地拟合了数据的线性分量，并且还适应了其异方差性。这种灵活性，以及模块化指定条件混合权重和条件分量密度的能力，使得依赖密度回归成为一种非常有用的非参数贝叶斯模型。

要了解有关依赖密度回归和相关模型的更多信息，请查阅《贝叶斯数据分析》、《贝叶斯非参数数据分析》或《贝叶斯非参数方法》。

此示例首次出现在这里。

作者：Austin Rochford

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

Last updated: Sat Apr 24 2021

Python implementation: CPython
Python version       : 3.7.6
IPython version      : 7.13.0

numpy     : 1.17.5
matplotlib: 3.2.1
pymc3     : 3.11.2
arviz     : 0.11.2
theano    : 1.1.2
pandas    : 1.1.5
seaborn   : 0.10.0

Watermark: 2.2.0