企业官网建设流程全解析

PyMC 混合分布完全指南：Mixture、NormalMixture 与零膨胀 / Hurdle 分布详解

【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc

本文是 PyMC 概率编程框架中"混合分布家族"的深度技术指南，覆盖文档 docs/source/api/distributions/mixture.rst 所收录的Mixture、NormalMixture、三个零膨胀分布（ZeroInflatedPoisson/ZeroInflatedBinomial/ZeroInflatedNegativeBinomial）与四个 Hurdle 分布（HurdlePoisson/HurdleNegativeBinomial/HurdleGamma/HurdleLogNormal）。读完本文，你将掌握用.dist()API 构造任意分量混合模型、处理子群体异质性、建模"零膨胀"计数数据与"零障碍"（hurdle）数据的完整方法，并理解其底层 logp 计算、参数校验与默认变换机制。

混合分布家族总览

在 PyMC 中，混合分布（mixture distribution）用于描述由多个底层分布按权重叠加而成的随机变量，是建模子群体异质性（subpopulation heterogeneity）的标准工具。其概率质量/密度函数的一般形式为：

$$f(x \mid w, \theta) = \sum_{i=1}^{n} w_i f_i(x \mid \theta_i)$$

其中 $w_i$ 为第 $i$ 个分量的权重（满足 $w_i \ge 0$ 且 $\sum_i w_i = 1$），$f_i$ 为第 $i$ 个分量的分布。该分布的支撑为各分量支撑的并集，均值为各分量均值的加权和 $\sum_i w_i \mu_i$。

PyMC 的混合分布家族全部实现在 pymc/distributions/mixture.py 中，并经由 pymc/distributions/init.py 导出到pm命名空间。API 文档收录的 9 个类按用途可分为三组：

通用Mixture分布

参数说明

pm.Mixture的完整签名如下：

pm.Mixture(name, w, comp_dists, **kwargs)

• w：tensor_like of float，混合权重，要求 $0 \le w \le 1$ 且sum(w) == 1，通常由pm.Dirichlet或pm.Dirichlet的变换版本给出。
• comp_dists：iterable of unnamed distributions或单个批处理（batched）分布。分量必须通过.dist()API 创建（不能直接传入已在模型中注册的随机变量）。若传入单个分布，则其最后一个 size 维度（而非 shape 维度）决定混合分量的个数，例如pm.Poisson.dist(..., size=components)。
• **kwargs：与所有分布一致的通用参数，如shape、observed、dims等。

需要特别注意的是，源码中_BaseMixtureDistribution.dist会对传入的comp_dists做多重校验（见 pymc/distributions/mixture.py 的_BaseMixtureDistribution类）：

1. 所有分量必须是TensorVariable且其owner.op为RandomVariable或SymbolicRandomVariable，否则抛出"Component dist must be a distribution created via the .dist() API"错误；
2. 分量不能是模型中已注册的变量（通过check_dist_not_registered检查）；
3. 所有分量的支撑维度（ndim_supp）必须一致；
4. 若分量超过一个，则要么全部为连续类型、要么全部为离散类型，混用会抛出ValueError。

分量构造的两种等价形式

Mixture的分量既可以用分量列表给出，也可以用单个批处理分布给出，两者在数学上等价，但后者得益于向量化，计算效率更高。以下为源码 docstring 中的官方示例。

例 1：两个 Poisson 分量的混合

import numpy as np import pymc as pm with pm.Model() as model: w = pm.Dirichlet("w", a=np.array([1, 1])) # 2 个混合权重 lam1 = pm.Exponential("lam1", lam=1) lam2 = pm.Exponential("lam2", lam=1) # 只需 logp，因此用 .dist() 创建分量而不向模型添加 RV # 以下两种形式等价，但第二种受益于向量化 components = [ pm.Poisson.dist(mu=lam1), pm.Poisson.dist(mu=lam2), ] # shape=(2,) 表示 2 个混合分量 components = pm.Poisson.dist(mu=pm.math.stack([lam1, lam2]), shape=(2,)) like = pm.Mixture("like", w=w, comp_dists=components, observed=data)

例 2：Normal 与 StudentT 的混合

with pm.Model() as model: w = pm.Dirichlet("w", a=np.array([1, 1])) # 2 个混合权重 mu = pm.Normal("mu", 0, 1) components = [ pm.Normal.dist(mu=mu, sigma=1), pm.StudentT.dist(nu=4, mu=mu, sigma=1), ] like = pm.Mixture("like", w=w, comp_dists=components, observed=data)

例 3：5×3 的 Normal 混合（批量混合）

with pm.Model() as model: # w 是 5 个独立的长度为 3 的权重向量 # 若 shape 为 (3,)，则权重会在 5 个复制维度间共享 w = pm.Dirichlet("w", a=np.ones(3), shape=(5, 3)) # 3 个混合分量各有一个独立的均值 mu = pm.Normal("mu", mu=np.arange(3), sigma=1, shape=3) # 两种形式等价，第二种受益于向量化 components = [ pm.Normal.dist(mu=mu[0], sigma=1, shape=(5,)), pm.Normal.dist(mu=mu[1], sigma=1, shape=(5,)), pm.Normal.dist(mu=mu[2], sigma=1, shape=(5,)), ] components = pm.Normal.dist(mu=mu, sigma=1, shape=(5, 3)) # 混合结果是长度为 5 的数组 # 每个元素可视为 3 个不同均值分量的独立标量混合 like = pm.Mixture("like", w=w, comp_dists=components, observed=data)

例 4：多变量分量的混合（Dirichlet）

with pm.Model() as model: w = pm.Dirichlet("w", a=np.ones(2)) # 2 个混合权重 # 两种形式等价，第二种受益于向量化 components = [ pm.Dirichlet.dist(a=[1, 10, 100], shape=(3,)), pm.Dirichlet.dist(a=[100, 10, 1], shape=(3,)), ] components = pm.Dirichlet.dist(a=[[1, 10, 100], [100, 10, 1]], shape=(2, 3)) # 混合结果是长度为 3 的数组 # 每个元素只来自两个核心 Dirichlet 分量之一 like = pm.Mixture("like", w=w, comp_dists=components, observed=data)

示例 4 说明Mixture不仅支持标量分量，也支持多元分量——源码中分量堆叠轴（mixture axis）会根据分量的支撑维度ndim_supp动态计算（mix_axis = -ndim_supp - 1），因此 MvNormal、Dirichlet 等多变量分布同样可以充当分量。

底层 logp 计算原理

Mixture的对数概率由 pymc/distributions/mixture.py 中的mixture_logprob函数（注册于MixtureRV）实现，其核心是一个数值稳定的加权求和：

mix_logp = pt.logsumexp(pt.log(weights) + components_logp, axis=-1)

其中components_logp是各分量在该取值处的 logp 在混合轴上堆叠的结果。使用logsumexp而非朴素地先算log(w_i) + log(f_i(x))再直接相加，可以避免下溢并保持梯度稳定。随后通过check_parameters校验权重约束（0 <= w <= 1且sum(w) == 1），违反约束时返回-inf并携带错误信息"0 <= weights <= 1, sum(weights) == 1"。

此外，MixtureRV还注册了：

• mixture_logcdf：累积分布函数的对数，同样以logsumexp(log(weights) + components_logcdf, axis=-1)计算，并做相同的权重约束校验；
• mixture_support_point：初始支持点，取各分量支持点的加权和（离散分量会四舍五入）；
• mixture_default_transform：默认变换逻辑。只有当所有分量的默认变换类型一致且属于白名单（CholeskyCovPacked、CircularTransform、IntervalTransform、LogTransform、LogExpM1、LogOddsTransform、Ordered、SimplexTransform、SumTo1）时才自动施加变换；否则发出MixtureTransformWarning警告，提示"未找到安全的默认变换，若合适可指定自定义变换以提高采样效率"。对于IntervalTransform，还会进一步校验各分量区间变换的向后表达式是否等价（例如不允许Interval(0, 1)与Interval(0, 2)混用）。

NormalMixture：高斯混合的便捷封装

pm.NormalMixture是高斯混合模型（GMM）的开箱即用封装，其密度为：

$$f(x \mid w, \mu, \sigma^2) = \sum_{i=1}^{n} w_i N(x \mid \mu_i, \sigma^2_i)$$

支撑为 $x \in \mathbb{R}$，均值为 $\sum_i w_i \mu_i$，方差为 $\sum_i w_i (\sigma^2_i + \mu_i^2) - \left(\sum_i w_i \mu_i\right)^2$。

参数与签名

pm.NormalMixture(name, w, mu, sigma=None, tau=None, **kwargs) pm.NormalMixture.dist(w, mu, sigma=None, tau=None, **kwargs)

注意：sigma与tau只需传入其一，不可同时传入。源码通过get_tau_sigma(tau=tau, sigma=sigma)统一解析，其内部实现为Mixture(name, w, Normal.dist(mu, sigma=sigma), **kwargs)——即NormalMixture本质上是Mixture与单个批处理Normal.dist的组合，最后一个 size 维度对应分量个数。

官方示例：三分量高斯混合

n_components = 3 with pm.Model() as gauss_mix: μ = pm.Normal( "μ", mu=data.mean(), sigma=10, shape=n_components, transform=pm.distributions.transforms.ordered, # 施加排序变换以解决标签交换问题 initval=[1, 2, 3], ) σ = pm.HalfNormal("σ", sigma=10, shape=n_components) weights = pm.Dirichlet("w", np.ones(n_components)) y = pm.NormalMixture("y", w=weights, mu=μ, sigma=σ, observed=data)

示例中为均值向量施加了ordered排序变换并给出initval=[1, 2, 3]，这是处理高斯混合模型**标签交换（label switching）**问题的常用手段，可让后验分布保持可识别。

零膨胀分布：Zero-Inflated 系列

零膨胀分布用于建模计数数据中零的比例高于基础分布（如 Poisson）所预期的情形。零膨胀的含义是：数据的一部分"额外零"来自一个独立于计数过程本身的机制（例如未患病的人群永远不会产生发病计数），因此零的密度为 $(1-\psi)$ 与基础分布在 0 处的概率之和。

ZeroInflatedPoisson

概率质量函数为：

$$f(x \mid \psi, \mu) = \begin{cases} (1-\psi) + \psi e^{-\mu}, & \text{if } x = 0 \ \psi \frac{e^{-\mu}\mu^x}{x!}, & \text{if } x=1,2,3,\ldots \end{cases}$$

支撑为 $x \in \mathbb{N}_0$，均值为 $\psi\mu$，方差为 $\mu + \frac{1-\psi}{\psi}\mu^2$（呈过度离散）。

用法：pm.ZeroInflatedPoisson("y", psi=psi, mu=mu, observed=data)。

ZeroInflatedBinomial

$$f(x \mid \psi, n, p) = \begin{cases} (1-\psi) + \psi (1-p)^{n}, & \text{if } x = 0 \ \psi {n \choose x} p^x (1-p)^{n-x}, & \text{if } x=1,2,\ldots,n \end{cases}$$

支撑为 $x \in \mathbb{N}_0$，均值为 $\psi n p$，方差为 $(1-\psi) n p [1 - p(1 - \psi n)]$。

ZeroInflatedNegativeBinomial

$$f(x \mid \psi, \mu, \alpha) = \begin{cases} (1-\psi) + \psi \left( \frac{\alpha}{\alpha+\mu} \right)^\alpha, & \text{if } x = 0 \ \psi \frac{\Gamma(x+\alpha)}{x! \Gamma(\alpha)} \left( \frac{\alpha}{\mu+\alpha} \right)^\alpha \left( \frac{\mu}{\mu+\alpha} \right)^x, & \text{if } x=1,2,3,\ldots \end{cases}$$

支撑为 $x \in \mathbb{N}_0$，均值为 $\psi\mu$，方差为 $\psi \left(\frac{\mu^2}{\alpha}\right) + \psi \mu + \psi \mu^2 - \psi^2 \mu^2$。

该分布支持两种参数化方式，两者之间由以下关系连接：

$$\mu = \frac{n(1-p)}{p}, \qquad \alpha = n$$

mu/alpha与p/n两组参数任选其一即可：pm.ZeroInflatedNegativeBinomial("y", psi=psi, mu=mu, alpha=alpha, observed=data)或pm.ZeroInflatedNegativeBinomial("y", psi=psi, p=p, n=n, observed=data)。

零膨胀的实现机制

从源码看，三个零膨胀分布都经由私有辅助函数_zero_inflated_mixture构造，其本质是Mixture的特例：权重被硬编码为[1 - psi, psi]，两个分量分别为DiracDelta.dist(0)（确定性地产生 0 的点质量分布）和对应的基础分布：

weights = pt.stack([1 - nonzero_p, nonzero_p], axis=-1) comp_dists = [ DiracDelta.dist(0), # 零部分：额外的零 nonzero_dist, # 非零部分：Poisson / Binomial / NegativeBinomial ]

因此"零膨胀"在模型层面被翻译为一个二分量混合：要么以概率 $1-\psi$ 来自确定为零的DiracDelta分量，要么以概率 $\psi$ 来自基础计数分布。测试文件 tests/distributions/test_mixture.py 的TestZeroInflatedMixture类覆盖了各分布在不同psi、size组合下的随机采样、logp 与 logcdf 正确性验证。

Hurdle 分布：零障碍（两阶段）模型

Hurdle（又称"零障碍"）模型与零膨胀模型的关键区别在于：hurdle 模型中的零不是"膨胀"出来的，而是来自一个完全独立的机制。建模思想是"先跨过一道门槛（hurdle）决定结果是否为零，跨过后再决定非零部分的具体取值"。其一般形式为：

$$f(x \mid \psi, \theta) = \begin{cases} 1 - \psi, & \text{if } x = 0 \ \psi \frac{\text{PDF}(x \mid \theta)}{1 - \text{CDF}(\epsilon \mid \theta)}, & \text{if } x \ge 1 \end{cases}$$

其中 $\epsilon$ 为机器精度（machine precision），非零部分的密度被截断分布（排除零）重新归一化。直观理解：若基础分布本身在 0 处已有概率质量（如离散 Poisson），hurdle 模型会把"自然的零"也一并排除，交由独立的零过程负责，从而让 $\psi$ 精确刻画"非零发生的概率"。

PyMC 提供四个 Hurdle 分布，全部通过_Hurdle基类与_Hurdle._create构造。

HurdlePoisson

$$f(x \mid \psi, \mu) = \begin{cases} 1 - \psi, & \text{if } x = 0 \ \psi \frac{\text{PoissonPDF}(x \mid \mu)}{1 - \text{PoissonCDF}(0 \mid \mu)}, & \text{if } x=1,2,3,\ldots \end{cases}$$

参数：psi（期望来自 Poisson 过程的样本比例，$0 < \psi < 1$）、mu（期望事件数，$\mu \ge 0$）。

HurdleNegativeBinomial

$$f(x \mid \psi, \mu, \alpha) = \begin{cases} 1 - \psi, & \text{if } x = 0 \ \psi \frac{\text{NegativeBinomialPDF}(x \mid \mu, \alpha)}{1 - \text{NegativeBinomialCDF}(0 \mid \mu, \alpha)}, & \text{if } x=1,2,3,\ldots \end{cases}$$

参数：psi、mu（均值，$\mu > 0$）、alpha（形状参数，$\alpha > 0$），同样支持p/n备选参数化。

HurdleGamma

$$f(x \mid \psi, \alpha, \beta) = \begin{cases} 1 - \psi, & \text{if } x = 0 \ \psi \frac{\text{GammaPDF}(x \mid \alpha, \beta)}{1 - \text{GammaCDF}(\epsilon \mid \alpha, \beta)}, & \text{if } x \ge 1 \end{cases}$$

参数：psi、alpha（形状参数，$\alpha > 0$）、beta（速率参数，$\beta > 0$），或使用备选的mu（$\mu > 0$）与sigma（$\sigma > 0$）参数化。

⚠️ 重要限制：HurdleGamma无法用 MCMC 方法正确采样（需要专用步进采样器），因此它只能用作观测变量（observed），或者仅用前向方法（sample_prior_predictive、sample_posterior_predictive）采样。这一限制在源码 docstring 中以.. warning::明确标出。

HurdleLogNormal

$$f(x \mid \psi, \mu, \sigma) = \begin{cases} 1 - \psi, & \text{if } x = 0 \ \psi \frac{\text{LogNormalPDF}(x \mid \mu, \sigma)}{1 - \text{LogNormalCDF}(\epsilon \mid \mu, \sigma)}, & \text{if } x \ge 1 \end{cases}$$

参数：psi、mu（位置参数，默认 0）、sigma（标准差，$\sigma > 0$，默认 1）、tau（尺度参数，$\tau > 0$，默认 1）。sigma与tau二选一即可。

⚠️ 重要限制：与HurdleGamma相同，HurdleLogNormal也不能用 MCMC 正确采样，仅适用于作为观测变量或前向采样。

Hurdle 的实现机制

从源码_Hurdle._create可以看出，Hurdle 分布的构造逻辑与零膨胀分布同构但有一处关键差异——对离散非零分布进行截断以排除零：

1. 若非零分布是离散的（dtype以"int"开头），则用Truncated.dist(nonzero_dist, lower=1, max_n_steps=10_000)将其截断，使零从基础分布中被彻底排除（连续分布在 0 处概率质量为零，无需截断即可直接使用；max_n_steps为截断随机采样时的最大迭代步数）；
2. 权重硬编码为[1 - psi, psi]，分量仍为DiracDelta.dist(0)与截断后的非零分布；
3. 若非零分布类型既非int也非float，抛出ValueError。

Hurdle 的 logp 由marginal_hurdle_logprob（注册于_HurdleRV）实现。值得关注的是，其内部为规避 NaN 梯度做了精细处理：pt.where会同时求值两个分支，而对连续分布（如 Gamma）计算logp(dist, 0)会得到-inf进而污染梯度，因此先用pt.switch(pt.eq(value, 0), 1.0, value)将零值替换为合法取值 1.0 再计算非零分支的 logp，最后才用pt.where按"是否为零"选择 $-\log(1-\psi)$ 或 $\log(\psi)+\text{logp}$。校验条件为 $0 \le \psi \le 1$。

采样、验证与实战建议

测试验证与可靠性

混合分布家族的正确性在 tests/distributions/test_mixture.py（共 1756 行）中有系统性保障，包括：

• TestMixture：单分量/多分量在不同 size、权重（含确定性权重如[1, 0]）下的形状断言、随机采样正确性与 logp 对照验证；
• TestNormalMixture：与scipy统计分布对照的 logp/logcdf 一致性检验；
• TestZeroInflatedMixture：零膨胀三个分布在多种psi取值下的采样与概率验证；
• TestHurdleDistributions：四个 Hurdle 分布从psi=0.05到psi=0.9的参数化采样验证，并断言其owner.op确为_Hurdle类型。

选型与使用建议

1. 自定义任意分量的混合：用pm.Mixture，分量必须经.dist()创建；优先用单个批处理分布（shape=(n_components,)）以利用向量化。
2. 高斯混合模型：直接用pm.NormalMixture，并配合ordered变换缓解标签交换问题。
3. 计数数据的"额外零"：若零来自独立机制（如未暴露人群），用零膨胀系列（ZIP、ZIB、ZINB）；若零与非零由完全不同的过程决定（先决定"是否发生"，再决定"发生多少"），用 Hurdle 系列。
4. 连续 Hurdle 建模：HurdleGamma/HurdleLogNormal只能作为观测变量或用于前向采样，不要试图用 MCMC 推断其作为潜变量时的后验。
5. 权重约束：无论哪种混合，权重都必须满足 $0 \le w \le 1$ 且归一化，违反时 logp 返回-inf；建议始终以pm.Dirichlet或pm.Dirichlet派生变量作为权重来源。

相关资源

• API 文档入口：docs/source/api/distributions/mixture.rst，与连续/离散/多元分布同属 docs/source/api/distributions.rst 的toctree；
• 核心实现：pymc/distributions/mixture.py（_BaseMixtureRV、Mixture、NormalMixture、零膨胀与 Hurdle 系列、mixture_logprob、mixture_logcdf、mixture_default_transform）；
• 公共导出：pymc/distributions/init.py；
• 测试与行为契约：tests/distributions/test_mixture.py。

【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc