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_{i1}^{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任意分量连续/离散的自定义加权混合便捷封装NormalMixture高斯混合模型GMM直接传均值/标准差零膨胀分布ZeroInflatedPoisson、ZeroInflatedBinomial、ZeroInflatedNegativeBinomial计数数据中额外的零来自独立机制Hurdle 分布HurdlePoisson、HurdleNegativeBinomial、HurdleGamma、HurdleLogNormal零与非零部分由完全不同过程产生的两阶段模型通用Mixture分布参数说明pm.Mixture的完整签名如下pm.Mixture(name, w, comp_dists, **kwargs)wtensor_like of float混合权重要求 $0 \le w \le 1$ 且sum(w) 1通常由pm.Dirichlet或pm.Dirichlet的变换版本给出。comp_distsiterable of unnamed distributions或单个批处理batched分布。分量必须通过.dist()API 创建不能直接传入已在模型中注册的随机变量。若传入单个分布则其最后一个 size 维度而非 shape 维度决定混合分量的个数例如pm.Poisson.dist(..., sizecomponents)。**kwargs与所有分布一致的通用参数如shape、observed、dims等。需要特别注意的是源码中_BaseMixtureDistribution.dist会对传入的comp_dists做多重校验见 pymc/distributions/mixture.py 的_BaseMixtureDistribution类所有分量必须是TensorVariable且其owner.op为RandomVariable或SymbolicRandomVariable否则抛出Component dist must be a distribution created via the .dist() API错误分量不能是模型中已注册的变量通过check_dist_not_registered检查所有分量的支撑维度ndim_supp必须一致若分量超过一个则要么全部为连续类型、要么全部为离散类型混用会抛出ValueError。分量构造的两种等价形式Mixture的分量既可以用分量列表给出也可以用单个批处理分布给出两者在数学上等价但后者得益于向量化计算效率更高。以下为源码 docstring 中的官方示例。例 1两个 Poisson 分量的混合import numpy as np import pymc as pm with pm.Model() as model: w pm.Dirichlet(w, anp.array([1, 1])) # 2 个混合权重 lam1 pm.Exponential(lam1, lam1) lam2 pm.Exponential(lam2, lam1) # 只需 logp因此用 .dist() 创建分量而不向模型添加 RV # 以下两种形式等价但第二种受益于向量化 components [ pm.Poisson.dist(mulam1), pm.Poisson.dist(mulam2), ] # shape(2,) 表示 2 个混合分量 components pm.Poisson.dist(mupm.math.stack([lam1, lam2]), shape(2,)) like pm.Mixture(like, ww, comp_distscomponents, observeddata)例 2Normal 与 StudentT 的混合with pm.Model() as model: w pm.Dirichlet(w, anp.array([1, 1])) # 2 个混合权重 mu pm.Normal(mu, 0, 1) components [ pm.Normal.dist(mumu, sigma1), pm.StudentT.dist(nu4, mumu, sigma1), ] like pm.Mixture(like, ww, comp_distscomponents, observeddata)例 35×3 的 Normal 混合批量混合with pm.Model() as model: # w 是 5 个独立的长度为 3 的权重向量 # 若 shape 为 (3,)则权重会在 5 个复制维度间共享 w pm.Dirichlet(w, anp.ones(3), shape(5, 3)) # 3 个混合分量各有一个独立的均值 mu pm.Normal(mu, munp.arange(3), sigma1, shape3) # 两种形式等价第二种受益于向量化 components [ pm.Normal.dist(mumu[0], sigma1, shape(5,)), pm.Normal.dist(mumu[1], sigma1, shape(5,)), pm.Normal.dist(mumu[2], sigma1, shape(5,)), ] components pm.Normal.dist(mumu, sigma1, shape(5, 3)) # 混合结果是长度为 5 的数组 # 每个元素可视为 3 个不同均值分量的独立标量混合 like pm.Mixture(like, ww, comp_distscomponents, observeddata)例 4多变量分量的混合Dirichletwith pm.Model() as model: w pm.Dirichlet(w, anp.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, ww, comp_distscomponents, observeddata)示例 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_{i1}^{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, sigmaNone, tauNone, **kwargs) pm.NormalMixture.dist(w, mu, sigmaNone, tauNone, **kwargs)参数说明w混合权重$0 \le w \le 1$mu各分量的均值sigma各分量的标准差tau各分量的精度precision注意sigma与tau只需传入其一不可同时传入。源码通过get_tau_sigma(tautau, sigmasigma)统一解析其内部实现为Mixture(name, w, Normal.dist(mu, sigmasigma), **kwargs)——即NormalMixture本质上是Mixture与单个批处理Normal.dist的组合最后一个 size 维度对应分量个数。官方示例三分量高斯混合n_components 3 with pm.Model() as gauss_mix: μ pm.Normal( μ, mudata.mean(), sigma10, shapen_components, transformpm.distributions.transforms.ordered, # 施加排序变换以解决标签交换问题 initval[1, 2, 3], ) σ pm.HalfNormal(σ, sigma10, shapen_components) weights pm.Dirichlet(w, np.ones(n_components)) y pm.NormalMixture(y, wweights, muμ, sigmaσ, observeddata)示例中为均值向量施加了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 } x1,2,3,\ldots \end{cases}$$支撑为 $x \in \mathbb{N}_0$均值为 $\psi\mu$方差为 $\mu \frac{1-\psi}{\psi}\mu^2$呈过度离散。参数说明psi期望来自 Poisson 过程的样本比例$0 \psi 1$mu给定时间区间内的期望事件数$\mu \ge 0$用法pm.ZeroInflatedPoisson(y, psipsi, mumu, observeddata)。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 } x1,2,\ldots,n \end{cases}$$支撑为 $x \in \mathbb{N}_0$均值为 $\psi n p$方差为 $(1-\psi) n p [1 - p(1 - \psi n)]$。参数说明psi期望来自 Binomial 过程的样本比例$0 \psi 1$n伯努利试验次数$n \ge 0$p单次试验成功概率$0 p 1$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 } x1,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$$参数说明psi期望来自 NegativeBinomial 过程的样本比例$0 \psi 1$muPoisson 参数均值$\mu 0$alphaGamma 参数形状$\alpha 0$p备选参数化单次试验成功概率$0 p 1$n备选参数化目标成功次数$n 0$mu/alpha与p/n两组参数任选其一即可pm.ZeroInflatedNegativeBinomial(y, psipsi, mumu, alphaalpha, observeddata)或pm.ZeroInflatedNegativeBinomial(y, psipsi, pp, nn, observeddata)。零膨胀的实现机制从源码看三个零膨胀分布都经由私有辅助函数_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 处已有概率质量如离散 Poissonhurdle 模型会把自然的零也一并排除交由独立的零过程负责从而让 $\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 } x1,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 } x1,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 分布的构造逻辑与零膨胀分布同构但有一处关键差异——对离散非零分布进行截断以排除零若非零分布是离散的dtype以int开头则用Truncated.dist(nonzero_dist, lower1, max_n_steps10_000)将其截断使零从基础分布中被彻底排除连续分布在 0 处概率质量为零无需截断即可直接使用max_n_steps为截断随机采样时的最大迭代步数权重硬编码为[1 - psi, psi]分量仍为DiracDelta.dist(0)与截断后的非零分布若非零分布类型既非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 分布从psi0.05到psi0.9的参数化采样验证并断言其owner.op确为_Hurdle类型。选型与使用建议自定义任意分量的混合用pm.Mixture分量必须经.dist()创建优先用单个批处理分布shape(n_components,)以利用向量化。高斯混合模型直接用pm.NormalMixture并配合ordered变换缓解标签交换问题。计数数据的额外零若零来自独立机制如未暴露人群用零膨胀系列ZIP、ZIB、ZINB若零与非零由完全不同的过程决定先决定是否发生再决定发生多少用 Hurdle 系列。连续 Hurdle 建模HurdleGamma/HurdleLogNormal只能作为观测变量或用于前向采样不要试图用 MCMC 推断其作为潜变量时的后验。权重约束无论哪种混合权重都必须满足 $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创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考