Dirichlet 过程混合模型用于密度估计#
Dirichlet 过程#
Dirichlet 过程 是分布空间上灵活的概率分布。最广义地说，集合 \( \Omega \) 上的概率分布 \( P \) 是将测度 1 分配给整个空间的[测度](https://en.wikipedia.org/wiki/Measure_(mathematics%29)（\( P(\Omega) = 1 \)）。狄利克雷过程 \( P \sim \textrm{DP}(\alpha, P_0) \) 是一种度量，它具有以下性质：对于 不相交 划分 \( S_1, \ldots, S_n \) 的每个有限 \( \Omega \)，
这里 \( P_0 \) 是空间 \( \Omega \) 上的基本概率测度。精度参数 \( \alpha > 0 \) 控制从 Dirichlet 过程采样的样本与基本测度 \( P_0 \) 的接近程度。当 \( \alpha \to \infty \) 时，来自 Dirichlet 过程的样本接近基本测度 \( P_0 \)。
Dirichlet 过程具有几个使其非常适合 MCMC 模拟的属性。
给定来自 Dirichlet 过程 \( P \sim \textrm{DP}(\alpha, P_0) \) 的 i.i.d. 观测 \( \omega_1, \ldots, \omega_n \) 的后验也是 Dirichlet 过程，其中
\[P\ |\ \omega_1, \ldots, \omega_n \sim \textrm{DP}\left(\alpha + n, \frac{\alpha}{\alpha + n} P_0 + \frac{1}{\alpha + n} \sum_{i = 1}^n \delta_{\omega_i}\right),\]
其中 \( \delta \) 是 狄拉克 delta 测度
新观测的后验预测分布是基本测度和观测的折衷，
\[\omega\ |\ \omega_1, \ldots, \omega_n \sim \frac{\alpha}{\alpha + n} P_0 + \frac{1}{\alpha + n} \sum_{i = 1}^n \delta_{\omega_i}.\]
我们看到先验精度 \( \alpha \) 可以自然地解释为先验样本大小。这种后验预测分布的形式也适用于 Gibbs 抽样。
来自 Dirichlet 过程的样本 \( P \sim \textrm{DP}(\alpha, P_0) \) 以概率 1 是离散的。也就是说，在 \( \Omega \) 中存在元素 \( \omega_1, \omega_2, \ldots \) 和权重 \( \mu_1, \mu_2, \ldots \)，其中 \( \sum_{i = 1}^{\infty} \mu_i = 1 \) 使得
\[P = \sum_{i = 1}^\infty \mu_i \delta_{\omega_i}.\]折棍过程 给出了上述权重 \( \mu_i \) 和样本 \( \omega_i \) 的显式构造，可以直接从中采样。如果 \( \beta_1, \beta_2, \ldots \sim \textrm{Beta}(1, \alpha) \)，则 \( \mu_i = \beta_i \prod_{j = 1}^{i - 1} (1 - \beta_j) \)。这种表示与折棍之间的关系可以说明如下
从长度为 1 的棍子开始。
将棍子分成两部分，第一部分的比例为 \( \mu_1 = \beta_1 \)，第二部分的比例为 \( 1 - \mu_1 \)。
进一步将第二部分分成两部分，第一部分的比例为 \( \beta_2 \)，第二部分的比例为 \( 1 - \beta_2 \)。这根棍子的第一部分的长度为 \( \beta_2 (1 - \beta_1) \)；第二部分的长度为 \( (1 - \beta_1) (1 - \beta_2) \)。
以这种方式继续无限期地折断前一次折断的第二部分。如果 \( \omega_1, \omega_2, \ldots \sim P_0 \)，则
\[P = \sum_{i = 1}^\infty \mu_i \delta_{\omega_i} \sim \textrm{DP}(\alpha, P_0).\]
[建议进一步阅读]：(http://mlg.eng.cam.ac.uk/tutorials/07/ywt.pdf) 和 简要介绍 Dirichlet 过程的其他变体及其应用。
我们可以使用上面的折棍过程轻松地从 Python 中的 Dirichlet 过程进行采样。对于此示例，\( \alpha = 2 \)，基本分布为 \( N(0, 1) \)。
import os
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
import scipy as sp
import seaborn as sns
import xarray as xr
%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
N = 20
K = 30
alpha = 2.0
P0 = sp.stats.norm
我们绘制并绘制来自折棍过程的样本。
beta = sp.stats.beta.rvs(1, alpha, size=(N, K))
w = np.empty_like(beta)
w[:, 0] = beta[:, 0]
w[:, 1:] = beta[:, 1:] * (1 - beta[:, :-1]).cumprod(axis=1)
omega = P0.rvs(size=(N, K))
x_plot = xr.DataArray(np.linspace(-3, 3, 200), dims=["plot"])
sample_cdfs = (w[..., np.newaxis] * np.less.outer(omega, x_plot.values)).sum(axis=1)
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(x_plot, sample_cdfs[0], c="gray", alpha=0.75, label="DP sample CDFs")
ax.plot(x_plot, sample_cdfs[1:].T, c="gray", alpha=0.75)
ax.plot(x_plot, P0.cdf(x_plot), c="k", label="Base CDF")
ax.set_title(rf"$\alpha = {alpha}$")
ax.legend(loc=2);
如上所述，当 \( \alpha \to \infty \) 时，来自 Dirichlet 过程的样本收敛到基本分布。
fig, (l_ax, r_ax) = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=(16, 6))
K = 50
alpha = 10.0
beta = sp.stats.beta.rvs(1, alpha, size=(N, K))
w = np.empty_like(beta)
w[:, 0] = beta[:, 0]
w[:, 1:] = beta[:, 1:] * (1 - beta[:, :-1]).cumprod(axis=1)
omega = P0.rvs(size=(N, K))
sample_cdfs = (w[..., np.newaxis] * np.less.outer(omega, x_plot.values)).sum(axis=1)
l_ax.plot(x_plot, sample_cdfs[0], c="gray", alpha=0.75, label="DP sample CDFs")
l_ax.plot(x_plot, sample_cdfs[1:].T, c="gray", alpha=0.75)
l_ax.plot(x_plot, P0.cdf(x_plot), c="k", label="Base CDF")
l_ax.set_title(rf"$\alpha = {alpha}$")
l_ax.legend(loc=2)
K = 200
alpha = 50.0
beta = sp.stats.beta.rvs(1, alpha, size=(N, K))
w = np.empty_like(beta)
w[:, 0] = beta[:, 0]
w[:, 1:] = beta[:, 1:] * (1 - beta[:, :-1]).cumprod(axis=1)
omega = P0.rvs(size=(N, K))
sample_cdfs = (w[..., np.newaxis] * np.less.outer(omega, x_plot.values)).sum(axis=1)
r_ax.plot(x_plot, sample_cdfs[0], c="gray", alpha=0.75, label="DP sample CDFs")
r_ax.plot(x_plot, sample_cdfs[1:].T, c="gray", alpha=0.75)
r_ax.plot(x_plot, P0.cdf(x_plot), c="k", label="Base CDF")
r_ax.set_title(rf"$\alpha = {alpha}$")
r_ax.legend(loc=2);
Dirichlet 过程混合模型#
对于密度估计任务，来自 Dirichlet 过程的样本的（几乎必然）离散性是一个显着的缺点。这个问题可以通过使用 Dirichlet 过程混合模型进行密度估计来解决。Dirichlet 过程混合模型使用来自参数族 \( \mathcal{F} = \{f_{\theta}\ |\ \theta \in \Theta\} \) 的分量密度，并将混合权重表示为 Dirichlet 过程。如果 \( P_0 \) 是参数空间 \( \Theta \) 上的概率测度，则 Dirichlet 过程混合模型是分层模型
为了说明这个模型，我们模拟了来自 Dirichlet 过程混合模型的抽取，其中 \( \alpha = 2 \)，\( \theta \sim N(0, 1) \)，\( x\ |\ \theta \sim N(\theta, (0.3)^2) \)。
N = 5
K = 30
alpha = 2
P0 = sp.stats.norm
f = lambda x, theta: sp.stats.norm.pdf(x, theta, 0.3)
beta = sp.stats.beta.rvs(1, alpha, size=(N, K))
w = np.empty_like(beta)
w[:, 0] = beta[:, 0]
w[:, 1:] = beta[:, 1:] * (1 - beta[:, :-1]).cumprod(axis=1)
theta = P0.rvs(size=(N, K))
dpm_pdf_components = f(x_plot, theta[..., np.newaxis])
dpm_pdfs = (w[..., np.newaxis] * dpm_pdf_components).sum(axis=1)
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(x_plot, dpm_pdfs.T, c="gray")
ax.set_yticklabels([]);
我们现在专注于单个混合模型，并将其分解为各个（加权）混合分量。
fig, ax = plt.subplots(figsize=(8, 6))
ix = 1
ax.plot(x_plot, dpm_pdfs[ix], c="k", label="Density")
ax.plot(
x_plot,
(w[..., np.newaxis] * dpm_pdf_components)[ix, 0],
"--",
c="k",
label="Mixture components (weighted)",
)
ax.plot(x_plot, (w[..., np.newaxis] * dpm_pdf_components)[ix].T, "--", c="k")
ax.set_yticklabels([])
ax.legend(loc=1);
从这些随机过程中采样很有趣，但是当我们将它们拟合到数据时，这些想法才变得真正有用。Dirichlet 过程的样本的离散性和折棍表示非常适合后验分布的马尔可夫链蒙特卡洛模拟。我们将使用 PyMC 执行此采样。
我们的第一个示例使用 Dirichlet 过程混合模型来估计 老忠实泉 间歇泉在 黄石国家公园 中爆发之间等待时间的密度。
try:
old_faithful_df = pd.read_csv(os.path.join("..", "data", "old_faithful.csv"))
except FileNotFoundError:
old_faithful_df = pd.read_csv(pm.get_data("old_faithful.csv"))
为了方便指定先验，我们标准化了爆发之间的等待时间。
old_faithful_df["std_waiting"] = (
old_faithful_df.waiting - old_faithful_df.waiting.mean()
) / old_faithful_df.waiting.std()
old_faithful_df.head()
| 爆发 | 等待 | 标准化等待时间 | |
|---|---|---|---|
| 0 | 3.600 | 79 | 0.596025 |
| 1 | 1.800 | 54 | -1.242890 |
| 2 | 3.333 | 74 | 0.228242 |
| 3 | 2.283 | 62 | -0.654437 |
| 4 | 4.533 | 85 | 1.037364 |
fig, ax = plt.subplots(figsize=(8, 6))
n_bins = 20
ax.hist(old_faithful_df.std_waiting, bins=n_bins, color="C0", lw=0, alpha=0.5)
ax.set_xlabel("Standardized waiting time between eruptions")
ax.set_ylabel("Number of eruptions");
细心的读者会注意到，我们并没有按照其定义无限期地继续折棍过程，而是将此过程截断为有限次数的折断。显然，在使用 Dirichlet 过程进行计算时，仅需要在内存中存储有限数量的点质量和权重。这种限制并不是非常繁重，因为对于有限数量的观测值，看起来混合模型分量的数量（这些分量为混合模型贡献了不可忽略的质量）的增长速度很可能比样本数量慢。这种直觉可以形式化地表明，为混合模型贡献不可忽略的质量的分量的（预期）数量接近 \( \alpha \log N \)，其中 \( N \) 是样本大小。
有各种巧妙的 Gibbs 抽样 技术用于 Dirichlet 过程，这些技术允许存储的分量数量根据需要增长。随机记忆化 [] 是另一种强大的技术，用于在模拟 Dirichlet 过程的同时仅在内存中存储有限数量的分量。在这个入门示例中，我们采用了不太复杂的方法，即在固定数量 \( K \) 个分量后，简单地截断存储的 Dirichlet 过程分量。 和 为截断提供了理由，表明 \( K > 5 \alpha + 2 \) 最有可能足以捕获几乎所有的混合权重（\( \sum_{i = 1}^{K} w_i > 0.99 \)）。在实践中，我们可以通过检查为混合模型贡献不可忽略的质量的分量数量来验证我们对 Dirichlet 过程的截断近似的适用性。如果在我们的模拟中，所有分量都为混合模型贡献了不可忽略的质量，则我们过早地截断了 Dirichlet 过程。
我们的（截断）Dirichlet 过程混合模型用于标准化等待时间是
请注意，与我们之前的模拟中固定 \( \alpha \) 的值不同，我们在 \( \alpha \) 上指定了先验，以便我们可以从观测值中学习其后验分布。
我们现在使用 PyMC 构建此模型。
N = old_faithful_df.shape[0]
K = 30
def stick_breaking(beta):
portion_remaining = pt.concatenate([[1], pt.extra_ops.cumprod(1 - beta)[:-1]])
return beta * portion_remaining
with pm.Model(coords={"component": np.arange(K), "obs_id": np.arange(N)}) as model:
alpha = pm.Gamma("alpha", 1.0, 1.0)
beta = pm.Beta("beta", 1.0, alpha, dims="component")
w = pm.Deterministic("w", stick_breaking(beta), dims="component")
tau = pm.Gamma("tau", 1.0, 1.0, dims="component")
lambda_ = pm.Gamma("lambda_", 10.0, 1.0, dims="component")
mu = pm.Normal("mu", 0, tau=lambda_ * tau, dims="component")
obs = pm.NormalMixture(
"obs", w, mu, tau=lambda_ * tau, observed=old_faithful_df.std_waiting.values, dims="obs_id"
)
我们使用 NUTS 从模型中采样 1,000 次，并使用 ADVI 初始化。
with model:
trace = pm.sample(
tune=2500,
init="advi",
target_accept=0.975,
random_seed=RANDOM_SEED,
)
Sampling 4 chains for 2_500 tune and 1_000 draw iterations (10_000 + 4_000 draws total) took 127 seconds.
\( \alpha \) 的后验分布高度集中在 0.25 和 1 之间。
az.plot_trace(trace, var_names=["alpha"]);
为了验证截断不会使我们的结果产生偏差，我们绘制了每个分量的后验预期混合权重。
fig, ax = plt.subplots(figsize=(8, 6))
plot_w = np.arange(K) + 1
ax.bar(plot_w - 0.5, trace.posterior["w"].mean(("chain", "draw")), width=1.0, lw=0)
ax.set_xlim(0.5, K)
ax.set_xlabel("Component")
ax.set_ylabel("Posterior expected mixture weight");
我们看到只有三个混合分量具有明显的后验预期权重，因此我们得出结论，将 Dirichlet 过程截断为 30 个分量并没有显着影响我们的估计。
我们现在计算并绘制我们的后验密度估计。
post_pdf_contribs = xr.apply_ufunc(
sp.stats.norm.pdf,
x_plot,
trace.posterior["mu"],
1.0 / np.sqrt(trace.posterior["lambda_"] * trace.posterior["tau"]),
)
post_pdfs = (trace.posterior["w"] * post_pdf_contribs).sum(dim=("component"))
post_pdf_quantiles = post_pdfs.quantile([0.1, 0.9], dim=("chain", "draw"))
fig, ax = plt.subplots(figsize=(8, 6))
n_bins = 20
ax.hist(old_faithful_df.std_waiting.values, bins=n_bins, density=True, color="C0", lw=0, alpha=0.5)
ax.fill_between(
x_plot,
post_pdf_quantiles.sel(quantile=0.1),
post_pdf_quantiles.sel(quantile=0.9),
color="gray",
alpha=0.45,
)
ax.plot(x_plot, post_pdfs.sel(chain=0, draw=0), c="gray", label="Posterior sample densities")
ax.plot(
x_plot,
az.extract(post_pdfs, var_names="x", num_samples=100),
c="gray",
)
ax.plot(x_plot, post_pdfs.mean(dim=("chain", "draw")), c="k", label="Posterior expected density")
ax.set_xlabel("Standardized waiting time between eruptions")
ax.set_yticklabels([])
ax.set_ylabel("Density")
ax.legend(loc=2);
如上所述，我们可以将此密度估计分解为其（加权）混合分量。
fig, ax = plt.subplots(figsize=(8, 6))
n_bins = 20
ax.hist(old_faithful_df.std_waiting.values, bins=n_bins, density=True, color="C0", lw=0, alpha=0.5)
ax.plot(x_plot, post_pdfs.mean(dim=("chain", "draw")), c="k", label="Posterior expected density")
ax.plot(
x_plot,
(trace.posterior["w"] * post_pdf_contribs).mean(dim=("chain", "draw")).sel(component=0),
"--",
c="k",
label="Posterior expected mixture\ncomponents\n(weighted)",
)
ax.plot(
x_plot,
(trace.posterior["w"] * post_pdf_contribs).mean(dim=("chain", "draw")).T,
"--",
c="k",
)
ax.set_xlabel("Standardized waiting time between eruptions")
ax.set_yticklabels([])
ax.set_ylabel("Density")
ax.legend(loc=2);
Dirichlet 过程混合模型在参数分量分布族 \( \{f_{\theta}\ |\ f_{\theta} \in \Theta\} \) 方面非常灵活。我们在下面通过使用泊松分量分布来估计每年太阳黑子的密度来说明这种灵活性。此数据集由 整理，可以下载。
kwargs = dict(sep=";", names=["time", "sunspot.year"], usecols=[0, 1])
try:
sunspot_df = pd.read_csv(os.path.join("..", "data", "sunspot.csv"), **kwargs)
except FileNotFoundError:
sunspot_df = pd.read_csv(pm.get_data("sunspot.csv"), **kwargs)
sunspot_df.head()
| 时间 | 太阳黑子年份 | |
|---|---|---|
| 0 | 1700.5 | 8.3 |
| 1 | 1701.5 | 18.3 |
| 2 | 1702.5 | 26.7 |
| 3 | 1703.5 | 38.3 |
| 4 | 1704.5 | 60.0 |
对于此示例，模型为
K = 50
N = sunspot_df.shape[0]
with pm.Model(coords={"component": np.arange(K), "obs_id": np.arange(N)}) as model:
alpha = pm.Gamma("alpha", 1.0, 1.0)
beta = pm.Beta("beta", 1, alpha, dims="component")
w = pm.Deterministic("w", stick_breaking(beta), dims="component")
# Gamma is conjugate prior to Poisson
lambda_ = pm.Gamma("lambda_", 300.0, 2.0, dims="component")
obs = pm.Mixture(
"obs", w, pm.Poisson.dist(lambda_), observed=sunspot_df["sunspot.year"], dims="obs_id"
)
with model:
trace = pm.sample(
tune=5000,
init="advi",
target_accept=0.95,
random_seed=RANDOM_SEED,
)
Sampling 4 chains for 5_000 tune and 1_000 draw iterations (20_000 + 4_000 draws total) took 508 seconds.
对于太阳黑子模型，\( \alpha \) 的后验分布集中在 0.6 和 1.2 之间，表明我们应该预期比老忠实泉等待时间模型更多的分量为混合模型贡献不可忽略的数量。
az.plot_trace(trace, var_names=["alpha"]);
实际上，我们看到大约有 10 到 15 个混合分量具有明显的后验预期权重。
fig, ax = plt.subplots(figsize=(8, 6))
plot_w = np.arange(K) + 1
ax.bar(plot_w - 0.5, trace.posterior["w"].mean(("chain", "draw")), width=1.0, lw=0)
ax.set_xlim(0.5, K)
ax.set_xlabel("Component")
ax.set_ylabel("Posterior expected mixture weight");
我们现在计算并绘制拟合的密度估计。
x_plot = xr.DataArray(np.arange(250), dims=["plot"])
post_pmf_contribs = xr.apply_ufunc(sp.stats.poisson.pmf, x_plot, trace.posterior["lambda_"])
post_pmfs = (trace.posterior["w"] * post_pmf_contribs).sum(dim=("component"))
post_pmf_quantiles = post_pmfs.quantile([0.025, 0.975], dim=("chain", "draw"))
fig, ax = plt.subplots(figsize=(8, 6))
ax.hist(sunspot_df["sunspot.year"].values, bins=40, density=True, lw=0, alpha=0.75)
ax.fill_between(
x_plot,
post_pmf_quantiles.sel(quantile=0.025),
post_pmf_quantiles.sel(quantile=0.975),
color="gray",
alpha=0.45,
)
ax.plot(x_plot, post_pmfs.sel(chain=0, draw=0), c="gray", label="Posterior sample densities")
ax.plot(
x_plot,
az.extract(post_pmfs, var_names="x", num_samples=100),
c="gray",
)
ax.plot(x_plot, post_pmfs.mean(dim=("chain", "draw")), c="k", label="Posterior expected density")
ax.set_xlabel("Yearly sunspot count")
ax.set_yticklabels([])
ax.legend(loc=1);
同样，我们可以将后验预期密度分解为加权混合密度。
fig, ax = plt.subplots(figsize=(8, 6))
n_bins = 40
ax.hist(sunspot_df["sunspot.year"].values, bins=n_bins, density=True, color="C0", lw=0, alpha=0.75)
ax.plot(x_plot, post_pmfs.mean(dim=("chain", "draw")), c="k", label="Posterior expected density")
ax.plot(
x_plot,
(trace.posterior["w"] * post_pmf_contribs).mean(dim=("chain", "draw")).sel(component=0),
"--",
c="k",
label="Posterior expected mixture\ncomponents\n(weighted)",
)
ax.plot(
x_plot,
(trace.posterior["w"] * post_pmf_contribs).mean(dim=("chain", "draw")).T,
"--",
c="k",
)
ax.set_xlabel("Yearly sunspot count")
ax.set_yticklabels([])
ax.legend(loc=1);
参考文献#
水印#
%load_ext watermark
%watermark -n -u -v -iv -w
Last updated: Mon Feb 06 2023
Python implementation: CPython
Python version : 3.10.9
IPython version : 8.9.0
numpy : 1.24.1
xarray : 2023.1.0
pandas : 1.5.3
matplotlib: 3.6.3
arviz : 0.14.0
pymc : 5.0.1
pytensor : 2.8.11
scipy : 1.10.0
seaborn : 0.12.2
Watermark: 2.3.1