PyMC 多元分布Multivariate Distributions完整指南16 个分布族的定义、参数化与建模实践【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymcPyMC 的pymc/distributions/multivariate.py模块承载了贝叶斯建模中最常用的多元随机变量从基础的多元正态MvNormal、多元 Student-t到 Dirichlet / Multinomial 族、相关性矩阵先验LKJ、矩阵值与 Kronecker 结构化正态再到空间统计中的 CAR / ICAR、非参数模型中的 StickBreakingWeights以及带零和约束的 ZeroSumNormal。本指南以官方 API 文档 docs/source/api/distributions/multivariate.rst 为主线逐一展开每个分布族的数学定义、参数语义、默认变换、源码实现要点与可复制的建模示例帮助读者掌握在 PyMC 中何时选哪个多元分布、如何参数化、如何写进模型的完整方法论。多元分布族总览在 PyMC 的 API 文档体系中多元分布是一个独立的分册其索引位于 docs/source/api/distributions/multivariate.rst并通过 docs/source/api/distributions.rst 的 toctree 与 continuous、discrete、mixture、timeseries 等分册并列。该索引使用 Sphinxautosummary指令配合distribution.rst模板自动生成 16 个分布各自的页面全部通过pymc命名空间导出CAR, Dirichlet, DirichletMultinomial, ICAR, KroneckerNormal, LKJCholeskyCov, LKJCorr, MatrixNormal, Multinomial, MvNormal, MvStudentT, OrderedMultinomial, StickBreakingWeights, Wishart, WishartBartlett, ZeroSumNormal在源码 pymc/distributions/multivariate.py 中__all__列表与上述 16 个名称完全一致说明文档索引与实现保持同步。按建模用途可把这些分布分为六组分组分布典型场景连续多元MvNormal、MvStudentT多维观测、纵向数据、厚尾多元数据成分与计数Dirichlet、Multinomial、DirichletMultinomial、OrderedMultinomial比例/成分数据、多分类计数、有序分类聚合计数协方差/相关矩阵先验LKJCholeskyCov、LKJCorr、Wishart、WishartBartlett为协方差矩阵设定先验常用于多层模型与多元回归矩阵与结构化MatrixNormal、KroneckerNormal矩阵观测、可分离协方差、时空/网格数据空间统计CAR、ICAR区域邻接数据、空间自相关建模贝叶斯空间流行病学特殊约束StickBreakingWeights、ZeroSumNormal非参数混合模型权重、可识别性约束如类别效应哑变量对应的完整测试覆盖位于 tests/distributions/test_multivariate.py其中包含了与 SciPy 对数密度的对照测试类TestMatchesScipy以及为每个分布单独编写的随机采样/对数密度参数化测试类如TestMvNormalCov、TestLKJCorr、TestWishart、TestICAR等是验证各分布实现的直接依据。MvNormal三种协方差参数化与精度矩阵重写定义与参数MvNormal 是多元正态分布概率密度为$$f(x \mid \mu, T) \frac{|T|^{1/2}}{(2\pi)^{k/2}} \exp\left{ -\frac{1}{2} (x-\mu)^{\top} T (x-\mu) \right}$$其支撑集为 $x \in \mathbb{R}^k$均值为 $\mu$方差为 $T^{-1}$。参数语义见 pymc/distributions/multivariate.pymu均值向量tensor_like of floatcov、tau、chol三选一互斥——cov为协方差矩阵、tau为精度矩阵协方差矩阵的逆、chol为协方差矩阵的 Cholesky 分解lower布尔值默认True表示chol为下三角 Cholesky 因子。dist内部调用quaddist_matrix见 pymc/distributions/multivariate.py统一三种参数化若传tau则先求cov inv(tau)若传chol则计算cov chol chol.mT并给chol打上lower_triangular标签以触发 PyTensor 重写。若同时或均未指定tau/cov/chol会抛出ValueError(Incompatible parameterization. Specify exactly one of tau, cov, or chol.)。三种写法的选择从数值稳定性与采样效率考虑PyMC 官方示例见 pymc/distributions/multivariate.py建议给定完整协方差矩阵小维度、直接已知协方差时cov np.array([[1.0, 0.5], [0.5, 2]]) mu np.zeros(2) vals pm.MvNormal(vals, mumu, covcov, shape(5, 2))优先使用 Cholesky 因子尤其是协方差本身是模型参数时配合LKJCholeskyCov见下文mu np.zeros(3) true_cov np.array([[1.0, 0.5, 0.1], [0.5, 2.0, 0.2], [0.1, 0.2, 1.0]]) data np.random.multivariate_normal(mu, true_cov, 10) sd_dist pm.Exponential.dist(1.0, shape3) chol, corr, stds pm.LKJCholeskyCov( chol_cov, n3, eta2, sd_distsd_dist, compute_corrTrue ) vals pm.MvNormal(vals, mumu, cholchol, observeddata)不可观测隐变量时采用非中心化参数化把相关性从均值结构中分离缓解采样退化sd_dist pm.Exponential.dist(1.0, shape3) chol, _, _ pm.LKJCholeskyCov(chol_cov, n3, eta2, sd_distsd_dist, compute_corrTrue) vals_raw pm.Normal(vals_raw, mu0, sigma1, shape(5, 3)) vals pm.Deterministic(vals, pt.dot(chol, vals_raw.T).T)源码级实现精度矩阵的特殊化重写对数密度的计算集中在quaddist_chol见 pymc/distributions/multivariate.py对协方差做 Cholesky 分解后用三角求解计算二次型 $(x-\mu)^{\top}\Sigma^{-1}(x-\mu)$ 与 $\log|\Sigma|$并返回正定性检查标志posdeflogp再通过check_parameters(..., msgposdef covariance)对非正定协方差给出诊断信息。值得注意的一个工程细节当用户用tau精度矩阵定义 MvNormal 时采样与pm.logp之外specialization_ir_rewrites_db中的重写规则mv_normal_to_precision_mv_normal见 pymc/distributions/multivariate.py会把MvNormal(mu, inv(tau))改写为专门的PrecisionMvNormalRV其对数密度直接用 $\delta^{\top}\tau\delta$ 计算而无需先求逆从而获得更高效的 logp。MvStudentT厚尾多元分布MvStudentT 是多元 Student-t 分布适用于存在离群值、尾部比正态更厚的多元数据。其密度为$$f(\mathbf{x}\mid \nu,\mu,\Sigma) \frac{\Gamma[(\nup)/2]}{\Gamma(\nu/2)\nu^{p/2}\pi^{p/2}|\Sigma|^{1/2}}\left[1\frac{1}{\nu}(\mathbf{x}-\mu)^{\top}\Sigma^{-1}(\mathbf{x}-\mu)\right]^{-(\nup)/2}$$参数见 pymc/distributions/multivariate.pynu自由度必须为正标量均值为 $\mu$当 $\nu1$方差为 $\frac{\nu}{\nu-2}\Sigma$当 $\nu2$否则未定义mu均值向量scale尺度矩阵新代码推荐Sigma为旧名称二者只能指定一个同时指定会抛ValueErrortau为精度矩阵chol为尺度矩阵的 Cholesky 因子lower默认True兼容性说明旧参数名cov仍被接受但会触发FutureWarning提示改用scale。实现上MvStudentTRV.rv_op见 pymc/distributions/multivariate.py采用正态/卡方构造法先采样 $Z \sim \mathcal{N}(0, \Sigma)$再采样 $U \sim \chi^2_\nu$输出 $\mu Z / \sqrt{U/\nu}$logp则用gammaln与log1p稳定计算对数密度并检查posdef与nu 0。成分与计数分布Dirichlet、Multinomial 与变体Dirichlet单纯形上的连续分布Dirichlet 分布定义在 K 维单纯形上密度为$$f(\mathbf{x}\mid\mathbf{a}) \frac{\Gamma(\sum_{i1}^k a_i)}{\prod_{i1}^k \Gamma(a_i)}\prod_{i1}^k x_i^{a_i-1}$$支撑为 $x_i \in (0,1)$ 且 $\sum x_i 1$均值为 $a_i / \sum a_i$。唯一参数a为浓度参数须 $a0$类别数由最后一轴长度决定见 pymc/distributions/multivariate.py。Dirichlet 继承SimplexContinuous基类其默认变换通过simplex_cont_transform注册为transforms.simplex见 pymc/distributions/multivariate.pyMCMC 在无约束空间上采样。logp对超出 $[0,1]$ 的值返回 $-\infty$。Multinomial多类别计数Multinomial 把二项分布推广到 K 类结果在 n 次独立试验中每次试验恰好落入 K 个类别之一$x[i]$ 表示第 i 类被观测到的次数$$f(x \mid n, p) \frac{n!}{\prod_{i1}^k x_i!}\prod_{i1}^k p_i^{x_i}$$参数见 pymc/distributions/multivariate.pyn每次复制的总试验次数$n0$整数p各类别概率$0 \le p \le 1$沿最后一轴求和应为 1。一个实用的自动处理当p是常量且最后一轴之和不为 1 时dist会发出UserWarning提示pparameters sum to ... instead of 1.0. They will be automatically rescaled.并自动归一化避免因浮点舍入或手误导致不可逆的报错。logp中同时校验 $0 \le p \le 1$、$\sum p 1$、$n \ge 0$并对非法取值返回 $-\infty$。support_point使用众数近似n*p取整后修正最大项但当期望最大计数很小时可能产生负值此时Assert会提示用户手动提供初始值——这是源码中明确标注的已知局限。DirichletMultinomial边际化成分的复合计数DirichletMultinomial 是Dirichlet 混合的 Multinomial即先 $p \sim \text{Dirichlet}(a)$ 再 $x \sim \text{Multinomial}(n, p)$但密度已对 $p$ 完成边际化$$f(x \mid n, a) \frac{\Gamma(n1)\Gamma(\sum a_k)}{\Gamma(n\sum a_k)}\prod_{k1}^K\frac{\Gamma(x_ka_k)}{\Gamma(x_k1)\Gamma(a_k)}$$参数为n试验次数与aDirichlet 浓度参数。实现上它是一个SymbolicRandomVariable见 pymc/distributions/multivariate.py采样时先抽p ~ Dirichlet(a)再抽Multinomial(n, p)但logp直接使用上面的边际化闭式解。与单独写 DirichletMultinomial 相比边际化降低了采样空间的维度常被用于过度离散overdispersion的计数建模。OrderedMultinomial有序分类的聚合计数回归当响应变量是有序类别如调查问卷的同意程度且观测是按试验聚合的计数向量时OrderedMultinomial是pm.OrderedLogistic的聚合版替代品后者只接受逐条disaggregated数据。它用 sigmoid 把线性预测子 $\eta$ 和 K-1 个切点 $c$ 映射到各类别概率$$f(k \mid \eta, c) \begin{cases} 1-\text{logit}^{-1}(\eta-c_1), k0 \ \text{logit}^{-1}(\eta-c_{k-1})-\text{logit}^{-1}(\eta-c_k), 0kK \ \text{logit}^{-1}(\eta-c_{K-1}), kK \end{cases}$$参数见 pymc/distributions/multivariate.pyeta预测子标量或向量cutpoints长度为 K-1 的切点数组把 $\eta$ 划分为 K 个区间文档明确提示不要把首尾切点显式设为 $\pm\infty$并建议用pm.distributions.transforms.ordered约束切点有序n多项试验总次数compute_p默认True是否把各类别推断概率存入 trace名称形如{name}_probs若关心内存可关闭。官方示例模拟 7 类的选举数据回归true_cum_p np.array([0.1, 0.15, 0.25, 0.50, 0.65, 0.90, 1.0]) true_p np.hstack([true_cum_p[0], true_cum_p[1:] - true_cum_p[:-1]]) fake_elections np.random.multinomial(n1_000, pvalstrue_p, size60) with pm.Model() as model: cutpoints pm.Normal( cutpoints, munp.arange(6) - 2.5, sigma1.5, initvalnp.arange(6) - 2.5, transformpm.distributions.transforms.ordered, ) pm.OrderedMultinomial( results, eta0.0, cutpointscutpoints, nfake_elections.sum(1), observedfake_elections, ) trace pm.sample()实现上_OrderedMultinomial.dist见 pymc/distributions/multivariate.py通过sigmoid(cutpoints - eta)构造累计概率再差分得到各类别概率从而复用Multinomial的完整实现。协方差与相关矩阵先验LKJ 家族与 WishartLKJCholeskyCovCholesky 分解的协方差先验LKJCholeskyCov是对Cholesky 分解后的协方差矩阵定义的先验相关矩阵服从 LKJ 分布见 [1] Lewandowski, Kurowicka and Joe (2009)标准差服从用户任意指定的正数分布。它是协方差先验的首选实现参数如下见 pymc/distributions/multivariate.pyname模型中的变量名etaLKJ 形状参数$eta0$。eta1表示相关矩阵上的均匀分布eta越大越偏向相关性很少接近单位矩阵的矩阵n协方差矩阵的维度$n1$sd_dist标准差分布须用.dist()API 创建的正数标量或向量分布应满足shape[-1]n标量分布会被自动 resize。注意sd_dist会被克隆clone与传入实例相互独立compute_corr默认True返回三元组(chol, corr, stds)为False时只返回 packed Choleskystore_in_trace默认True把corr与stds以{name}_corr、{name}_stds的名字存入后验 trace仅当compute_corrTrue时生效。返回值约定compute_corrTrue时返回chol解包后的 Cholesky 协方差分解、corr相关矩阵、stds标准差否则返回 packed Cholesky 向量。Packed 存储格式Cholesky 因子是下三角矩阵PyMC 将其按行编号压缩为一维数组[[0 - - -] [1 2 - -] [3 4 5 -] [6 7 8 9]]解包可通过pm.math.expand_packed_triangular(packed_cov, lowerTrue)完成helper_deterministics见 pymc/distributions/multivariate.py内部正是这样实现(chol, corr, stds)的拆分的由cov chol chol.T求标准差stds sqrt(diag(cov))再归一化得corr inv_stds[None,:] * cov * inv_stds[:,None]。实现细节在无约束空间中Cholesky 因子除对角线外原样存储对角线元素使用 log 变换保证为正见 pymc/distributions/multivariate.py。logp需要处理第二次变换——把 Cholesky 因子拆成行长度标准差与单位行相关 Cholesky 因子其 Jacobian 行列式按行块对角计算最终形式为 $\det(J_{\phi^{-1}}(U)) \left[\prod_{i2}^N u_{ii}^{i-1} L_{ii}\right]^{-1}$。另外_LKJCholeskyCovRV_logp目前只对常量n与常量eta实现见 pymc/distributions/multivariate.py非常量会抛NotImplementedError。典型用法与MvNormal配合的完整流程见 pymc/distributions/multivariate.pywith pm.Model() as model: # 注意使用 .dist() 访问标准差分布而不是创建新随机变量 sd_dist pm.Exponential.dist(1.0, size10) chol, corr, sigmas pm.LKJCholeskyCov(chol_cov, eta4, n10, sd_distsd_dist) # 只取 packed Cholesky 的写法 # packed_chol pm.LKJCholeskyCov(chol_cov, eta4, n10, sd_distsd_dist, compute_corrFalse) # chol pm.math.expand_packed_triangular(10, packed_chol, lowerTrue) vals pm.MvNormal(vals, munp.zeros(10), cholchol, shape10) # 或变换不相关的正态 vals_raw pm.Normal(vals_raw, mu0, sigma1, shape10) vals2 pt.dot(chol, vals_raw) # 或直接计算协方差矩阵 cov pt.dot(chol, chol.T)LKJCorr相关矩阵先验固定标准差当标准差固定、只有相关矩阵需要建模时LKJCorr更合适。它直接对相关矩阵定义先验eta1对应相关矩阵上的均匀分布$eta \to \infty$ 时逼近单位矩阵。参数n维度$n1$与eta$eta0$见 pymc/distributions/multivariate.py。其默认变换为CholeskyCorrTransform见 pymc/distributions/multivariate.py。with pm.Model() as model: sds 3 * np.ones(10) # 固定的标准差 corr pm.LKJCorr(corr, eta4, n10) vals sds * pm.MvNormal(vals, munp.zeros(10), covcorr, shape10) # 或变换不相关正态 vals_raw pm.Normal(vals_raw, shape10) chol pt.linalg.cholesky(corr) vals2 sds * pt.dot(chol, vals_raw)LKJCorrRV._random_corr_matrix见 pymc/distributions/multivariate.py基于 vine/扩展洋葱法的 Beta-正态构造用pytensor.scan逐步填充矩阵列logp由归一化常数_lkj_normalizing_constant见 pymc/distributions/multivariate.py与对角元素的对数和构成并校验eta 0与行范数等于 1即传入值确为相关矩阵。Wishart对称正定矩阵分布Wishart 分布定义在对称正定SPD矩阵上是一组 i.i.d. 多元正态向量外积和的分布若 $x_i \sim \mathcal{N}(0, V)$$i1,\dots,\nu$则 $X\sum_i x_i x_i^\top \sim \mathrm{Wishart}_p(\nu, V)$。参数见 pymc/distributions/multivariate.pynu自由度须满足 $\nu p-1$V与scale_chol二选一V为 $(p,p)$ SPD 尺度矩阵scale_chol为其下三角 Cholesky 因子$V L L^{\top}$。提供后者可避免在 logp 中重复做 Cholesky 分解两者同时或均不提供会抛ValueError。默认去约束变换为CholeskyCovTransform见 pymc/distributions/multivariate.py即用自由实向量参数化 $X L L^{\top}$、对角线元素取 $\log L_{kk}$因此logp中的cholesky(L L.T)可被重写为L运行时无需再做分解见 pymc/distributions/multivariate.py。采样端采用 Bartlett 分解对角元取自 $\sqrt{\chi^2_{\nu-k}}$严格下三角取标准正态再以scale_chol A组合见 pymc/distributions/multivariate.py。WishartBartlett已废弃的兼容 shimWishartBartlett是历史遗留的Bartlett 分解 Wishart 先验辅助函数见 pymc/distributions/multivariate.py。过去它是 PyMC 中唯一可用于 MCMC 的 Wishart旧版 Wishart 没有去约束变换而新版Wishart已原生支持 Cholesky 参数化与默认变换因此该函数被标记为Deprecated并将在未来版本移除。其参数S尺度矩阵或其 Cholesky 因子、is_cholesky、return_cholesky分别映射到Wishart的V/scale_chol与pt.linalg.cholesky包装传入initval会直接抛NotImplementedError。文档明确建议新代码一律使用pm.Wishart需要is_choleskyTrue时传scale_cholS需要返回 Cholesky 时把pm.Wishart包进pt.linalg.cholesky作为Deterministic。矩阵值与结构化分布MatrixNormal 与 KroneckerNormalMatrixNormal矩阵值正态MatrixNormal 是定义在 $m \times n$ 矩阵上的正态分布其协方差结构可分离为行间协方差与列间协方差的 Kronecker 积密度为$$f(x \mid \mu, U, V) \frac{1}{(2\pi^{mn}|U|^n|V|^m)^{1/2}}\exp\left{-\frac{1}{2}\mathrm{Tr}\left[V^{-1}(x-\mu)^{\top}U^{-1}(x-\mu)\right]\right}$$参数见 pymc/distributions/multivariate.pymu均值数组须与随机变量 X 可广播使mu X的形状为 (M, N)rowcov/rowchol$(M,M)$ 行间协方差矩阵或其 Cholesky 因子二选一rowcov定义了列内方差colcov/colchol$(N,N)$ 列间协方差矩阵或其 Cholesky 因子二选一若rowcov为单位阵则退化为MvNormal中的cov语义。官方示例pymc/distributions/multivariate.pywith pm.Model() as model: colcov np.array([[1.0, 0.5], [0.5, 2]]) rowcov np.array([[1, 0, 0], [0, 4, 0], [0, 0, 16]]) m, n rowcov.shape[0], colcov.shape[0] mu np.zeros((m, n)) vals pm.MatrixNormal(vals, mumu, colcovcolcov, rowcovrowcov)即第 i 行的方差被 $4^i$ 缩放。关键等价关系MatrixNormal等价于MvNormal(mu, np.kron(rowcov, colcov))但计算更快——它利用 Kronecker 积的性质避免构造和求逆 $mn \times mn$ 的大矩阵。一个行被未知常数幂缩放的可学习模型示例with pm.Model() as model: sd_dist pm.HalfCauchy.dist(beta2.5, shape3) colchol, _, _ pm.LKJCholeskyCov(colchol, n3, eta2, sd_distsd_dist) scale pm.LogNormal(scale, munp.log(true_scale), sigma0.5) rowcov pt.diag([scale ** (2 * i) for i in range(m)]) vals pm.MatrixNormal(vals, mumu, colcholcolchol, rowcovrowcov, observeddata)MatrixNormalRV.rng_fn用mu rowchol Z colchol.T采样Z 为标准正态矩阵logp则分两次三角求解计算 $\mathrm{Tr}[V^{-1}\Delta^{\top}U^{-1}\Delta]$避免显式求逆。KroneckerNormalKronecker 结构协方差的多元正态当协方差矩阵本身可写成多个小矩阵的 Kronecker 积 $K \bigotimes K_i$可再加对角噪声 $\sigma^2 I_N$时KroneckerNormal提供比直接使用 $N \times N$ 协方差矩阵更高效的计算。参数见 pymc/distributions/multivariate.pymu均值向量同MvNormalcovs协方差矩阵列表 $[K_1, K_2, \dots]$按给定顺序做 Kronecker 积 $\bigotimes K_i$cholsCholesky 因子列表 $[L_1, L_2, \dots]$满足 $K_i L_i L_i^{\top}$evds特征分解对列表 $[(v_1, Q_1), (v_2, Q_2), \dots]$满足 $K_i Q_i,\mathrm{diag}(v_i),Q_i^{\top}$可由pt.linalg.eigh(K_i)得到sigma高斯白噪声标准差标量默认 0对应协方差 $K \bigotimes K_i \sigma^2 I_N$。covs、chols、evds三者互斥同时传多个抛ValueError。三种方式在内部都会被统一为特征分解形式。示例pymc/distributions/multivariate.pyK1 np.array([[1.0, 0.5], [0.5, 2]]) K2 np.array([[1.0, 0.4, 0.2], [0.4, 2, 0.3], [0.2, 0.3, 1]]) covs [K1, K2] N, mu 6, np.zeros(6) with pm.Model() as model: vals pm.KroneckerNormal(vals, mumu, covscovs, shapeN) chols [np.linalg.cholesky(Ki) for Ki in covs] evds [np.linalg.eigh(Ki) for Ki in covs] with pm.Model() as model: vals2 pm.KroneckerNormal(vals2, mumu, cholschols, shapeN) vals3 pm.KroneckerNormal(vals3, mumu, evdsevds, shapeN) sigma 0.1 with pm.Model() as noise_model: vals pm.KroneckerNormal(vals, mumu, covscovs, sigmasigma, shapeN)实现要点效率增益来自分别分解 $K_1$、$K_2$ 而不是分解大矩阵 $K$即使加了 $\sigma^2 I_N$ 破坏了整体 Kronecker 结构logp仍可借助子矩阵特征分解、用kron_diag组合特征值、kron_dot加速线性求解见 pymc/distributions/multivariate.py。该方法出自 Saatchi (2011) 的结构化高斯过程推断论文适合时空数据与多维网格。空间统计CAR 与 ICARCAR条件自回归CARConditional Autoregression是协方差具有邻接结构的多元正态特例用于区域数据建模密度为$$f(x \mid W, \alpha, \tau) \frac{|T|^{1/2}}{(2\pi)^{k/2}}\exp\left{-\frac{1}{2}(x-\mu)^{\top}T^{-1}(x-\mu)\right}$$其中 $T (\tau D(I-\alpha W))^{-1}$$D \mathrm{diag}(\sum_i W_{ij})$。参数见 pymc/distributions/multivariate.pymu实值均值向量W$(M,M)$ 对称 0/1 邻接矩阵标记区域间是否相邻PyTensor 会尽量转为稀疏矩阵通过as_sparse_or_tensor_variable无法稀疏化时退回稠密变量alpha自回归参数取值 $-1 \alpha 1$越接近 0 相关性越弱越接近 1 自相关越强大多数场景建议把支持限制在 $(0,1)$tau正精度变量控制底层正态变量的尺度。CARRV.rng_fn实现了 Rue (2001) 的高斯马尔可夫随机场快速采样算法构造精度矩阵 $Q \tau(D - \alpha W)$ 后做 Reverse Cuthill-McKee 排序与带状 Cholesky 求解见 pymc/distributions/multivariate.py。logp通过特征值分解 $\mathrm{eigvalsh}(D^{-1/2}WD^{-1/2})$ 计算 $\log\det$并校验 $-1\alpha1$、$\tau0$、$W$ 对称。ICAR内蕴条件自回归先验ICAR 是 CAR 在 $\alpha1$ 时的退化情形主要用于对相邻区域间的空间相关性建模。其对数密度为见 pymc/distributions/multivariate.py$$f(\phi \mid W, \sigma) -\frac{1}{2\sigma^2}\sum_{i\sim j}(\phi_i - \phi_j)^2 - \frac{1}{2}\left(\frac{\sum_i \phi_i}{0.001N}\right)^2 - \ln\sqrt{2\pi} - \ln(0.001N)$$第一项是空间协方差项每个 $\phi_i$ 依据与所有邻居的平方距离被惩罚$i\sim j$ 表示对 $\phi_i$ 所有邻居求和后三项是以 0 为均值、标准差为 $N \times 0.001$ 的正态对数密度施加零和约束——把 $\phi$ 的和惩罚到接近 0。参数W整数 ndarray对称 0/1 邻接矩阵。dist在构造时做四项硬校验ndim2、方阵、对称、仅含 0/1见 pymc/distributions/multivariate.py不符合直接抛ValueErrorsigma标量默认 1$\phi$ 向量的标准差。给sigma设先验得到中心化参数化更推荐的做法是保持默认值用非中心化参数化抽phi ~ ICAR(W)再乘以sigmazero_sum_stdev标量默认 0.001控制零和约束的强度——$\phi$ 的和以 0 为均值、该值为标准差的正态分布。一个关键限制ICARRV.rng_fn直接抛NotImplementedError(Cannot sample from ICAR prior)——ICAR 先验本身无法独立采样只能通过logp参与推断。中心化与非中心化两种写法的官方对照pymc/distributions/multivariate.pyW np.array( [ [0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0], ], ) # 中心化参数化 with pm.Model(): sigma pm.Exponential(sigma, 1) phi pm.ICAR(phi, WW, sigmasigma) mu phi # 非中心化参数化推荐 with pm.Model(): sigma pm.Exponential(sigma, 1) phi pm.ICAR(phi, WW) mu sigma * philogp实现把邻接矩阵转为边列表取tril(W)的下三角逐对计算 $( \phi_i - \phi_j)^2$与零和惩罚项相加。该先验对应 Besag York MolliéBYM空间模型常用于空间流行病学。特殊约束分布StickBreakingWeights 与 ZeroSumNormalStickBreakingWeights截断的棍子破碎权重该分布生成截断的棍子破碎stick-breaking权重$x_k v_k \prod_{\ellk}(1-v_\ell)$$k1,\dots,K$且 $x_{K1} \prod_{\ell1}^{K}(1-v_\ell) 1-\sum_{\ell1}^{K}x_\ell$其中 $v_k \sim \mathrm{Beta}(1,\alpha)$ i.i.d.。密度为$$f(\mathbf{x}\mid\alpha, K) B(1,\alpha)^{-K}x_{K1}^{\alpha}\prod_{k1}^{K1}\left{\sum_{jk}^{K1}x_j\right}^{-1}$$参数见 pymc/distributions/multivariate.pyalpha浓度参数$\alpha0$控制权重分布的均匀程度K要折断的棍子数权重向量长度为 $K1$最后一个权重为 1 减去前 K 个权重之和。均值满足 $\mathbb{E}[x_k] \frac{1}{1\alpha}\left(\frac{\alpha}{1\alpha}\right)^{k-1}$$k1,\dots,K$与 $\mathbb{E}[x_{K1}] \left(\frac{\alpha}{1\alpha}\right)^{K}$。该分布是贝叶斯非参数模型如 Dirichlet 过程混合模型的截断近似中混合权重的核心构件参考 Ishwaran James (2001)。与 Dirichlet 相同它继承SimplexContinuous默认变换为 simplex。ZeroSumNormal沿轴零和约束的正态ZeroSumNormal是一个或多个轴被约束为求和为零的正态分布默认约束最后一个轴数学上等价于$$\mathrm{ZSN}(\sigma) \mathcal{N}\left(0,\ \sigma^2\left(I_K - \tfrac{1}{K}J_K\right)\right), \quad J_{ij}1,\ K\text{约束轴长度}$$参数见 pymc/distributions/multivariate.pysigma尺度参数$\sigma0$实际是底层未约束正态的标准差默认 1且在零和轴上长度不能超过 1否则抛ValueError以保证约束成立n_zerosum_axes整数默认 1指定从最右轴起连续多少个轴被施加零和约束必须 0若想要 0 个约束轴直接用pm.Normaldims/shape维度名或形状与其他分布一致二者与observed至少提供一个否则抛ValueError。典型应用是类别效应的可识别性约束如区域哑变量求和为零。示例pymc/distributions/multivariate.pyCOORDS { regions: [a, b, c], answers: [yes, no, whatever, dont understand question], } with pm.Model(coordsCOORDS) as m: # 零和轴为 answers v pm.ZeroSumNormal(v, dims(regions, answers)) with pm.Model(coordsCOORDS) as m: # 零和轴为 answers 和 regions v pm.ZeroSumNormal(v, dims(regions, answers), n_zerosum_axes2) with pm.Model(coordsCOORDS) as m: # 零和轴为最后两轴 v pm.ZeroSumNormal(v, shape(3, 4, 5), n_zerosum_axes2)实现上ZeroSumNormalRV.rv_op见 pymc/distributions/multivariate.py先按shape抽标准正态再沿指定的n_zerosum_axes逐个减去轴均值logp计算自由度修正后的对数密度自由度等于各零和轴长度减 1 之积并校验各轴均值是否接近 0atol1e-9。默认变换为ZeroSumTransform见 pymc/distributions/multivariate.py。实践建议如何为模型挑选多元分布综合上述源码与文档可提炼出如下选型准则协方差矩阵是模型参数优先LKJCholeskyCovsd_dist指定标准差分布MvNormal(chol...)的非中心化组合标准差固定时才用LKJCorr。协方差已知/可直接给矩阵MvNormal(cov...)追求 logp 效率可用tau但注意 PyTensor 会通过特殊化重写优化精度矩阵路径。数据有厚尾或离群值MvStudentT并注意Sigma已改名scale。成分/比例数据连续单纯形用Dirichlet聚合计数用Multinomial计数且过度离散用DirichletMultinomial有序类别的聚合计数回归用OrderedMultinomial。矩阵观测或可分离协方差MatrixNormal等价于 Kronecker 协方差的MvNormal但更快协方差本身是 Kronecker 积可加对角噪声用KroneckerNormal。区域空间数据CAR或ICAR后者无法直接采样须配sigma做非中心化。权重/可识别性约束非参数混合权重用StickBreakingWeights需要求和为零的效应向量用ZeroSumNormal。Wishart 家族新代码一律使用pm.WishartWishartBartlett已废弃对相关矩阵先验LKJCholeskyCov通常比 Wishart 更易用、数值更稳。所有 16 个分布的实现均位于 pymc/distributions/multivariate.py其行为由 tests/distributions/test_multivariate.py 中的 SciPy 对照与参数化随机测试保障分布索引页面 docs/source/api/distributions/multivariate.rst 与源码__all__保持一致是查阅每个分布详细数学定义、参数说明与完整示例的首选入口。【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考