企业官网建设流程全解析

PyMC 概率计算 API 完全指南：logp、logcdf、icdf 与条件对数概率深入解析

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

本文是 PyMC 概率编程框架中「Probability（概率计算）」API 的深度技术指南，核心围绕 docs/source/api/logprob.rst 文档所定义的五个核心函数展开：pm.logp、pm.logcdf、pm.icdf（位于pymc顶层命名空间）以及pymc.logprob.conditional_logp、pymc.logprob.transformed_conditional_logp（位于pymc.logprob子模块）。读完本文，你将掌握如何手工构造随机变量的对数概率/对数分布函数/逆分布函数计算图，理解条件对数概率如何实现分层模型的因子化联合概率，并能借助仓库源码理解其背后的 IR（中间表示）图重写与singledispatch派发机制，为自定义分布、诊断与模型调试提供精确的计算工具。

一、概述：PyMC 概率计算 API 的定位

在贝叶斯建模流程中，除了通过pm.Model建模并调用pm.sample采样之外，我们经常需要直接对随机变量的概率量进行求值——例如计算某个观测值的对数似然、分位数、尾部概率，或者提取分层模型中每个变量对应的条件对数概率项。PyMC 在顶层pymc命名空间（见 pymc/init.py 通过from pymc.logprob import *导出）和pymc.logprob子模块（见 pymc/logprob/init.py）中提供了这组完整的 API：

所有函数都返回PyTensor 符号计算图（Variable），而非直接返回数值。这一设计使其可以与eval()（调试用）或pm.compile（编译为可重复调用的函数）无缝衔接，并且能够被进一步嵌入更大的计算图中参与自动微分。

从模块 docstring（pymc/logprob/init.py）可以看到，这一整套 API 的核心使命是：「Conversion of PyMC graphs into logp graphs」——将 PyMC 图转换为 logp 图。

二、pm.logp：构造随机变量的对数概率计算图

pm.logp(rv, value, warn_rvs=True, **kwargs)接受两个核心参数（实现见 pymc/logprob/basic.py）：

• rv：一个随机变量（RandomVariable或其可测变换，如pt.exp(rv)）；
• value：与rv同类型的值变量（Variable或tensor_like）。如果传入的不是Variable，会被pt.as_tensor_variable(value, dtype=rv.dtype)自动转换（basic.py）；
• warn_rvs：默认True。当 logp 图中出现其他未替换的随机变量时发出警告，因为此时这些 RV 会被克隆、与原始变量身份不一致，可能导致后续无法替换为观测值（basic.py）；
• 返回值是一个 logp 符号Variable；若无法推导会抛出RuntimeError。

2.1 基本用法：单变量对数概率

这是 docstring 中最直接的示例（basic.py），演示两种求值方式：

import pymc as pm import pytensor.tensor as pt mu = pt.scalar("mu") rv = pm.Normal.dist(mu, 1.0) value = pt.scalar("value") rv_logp = pm.logp(rv, value) # 方式一：用 .eval() 调试 print(rv_logp.eval({value: 0.9, mu: 0.0})) # -1.32393853 # 方式二：编译为函数反复求值 rv_logp_fn = pm.compile([value, mu], rv_logp) print(rv_logp_fn(value=0.9, mu=0.0)) # -1.32393853

这里的pm.Normal.dist(mu, 1.0)创建的是一个自由随机变量（dist不进入任何模型上下文），而pm.logp把它与取值点value绑定，构造出标准正态分布在 0.9 处的对数密度。

2.2 随机变量的变换：自动推导密度

pm.logp的独特之处在于，当rv是某个随机变量的确定性变换（例如对数正态pt.exp(rv)）时，它会自动推导变换后的密度（basic.py）：

mu = pt.scalar("mu") rv = pm.Normal.dist(mu, 1.0) exp_rv = pt.exp(rv) # 对数正态变量 value = pt.scalar("value") exp_rv_logp = pm.logp(exp_rv, value) print(exp_rv_logp.eval({value: 0.9, mu: 0.0})) # -0.81912844 exp_rv_logp_fn = pm.compile([value, mu], exp_rv_logp) print(exp_rv_logp_fn(value=0.9, mu:0.0) if False else exp_rv_logp_fn(value=0.9, mu=0.0)) # -0.81912844

2.3 在CustomDist中定义自定义分布的对数概率

pm.logp还是自定义分布的标准构建块（basic.py）：

import pymc as pm import pytensor.tensor as pt def normal_logp(value, mu, sigma): return pm.logp(pm.Normal.dist(mu, sigma), value) with pm.Model() as model: mu = pm.Normal("mu") sigma = pm.HalfNormal("sigma") pm.CustomDist("x", mu, sigma, logp=normal_logp)

2.4 底层原理：helper 派发与 IR 图重建

从源码看，pm.logp内部先尝试直接调用_logprob_helper(rv, value, **kwargs)（basic.py）；如果抛出NotImplementedError，则通过construct_ir_fgraph({rv: value})将原始图转换为 IR 中间表示图（basic.py），取出其中的(ir_rv, ir_value)配对后再调用_logprob_helper，最后用cleanup_ir清理图并可选地检查残留 RV 发出警告。

_logprob_helper（pymc/logprob/abstract.py）则调用基于singledispatch的_logprob分发器，按rv.owner.op的具体类型路由到对应分布的密度实现，并自动给结果命名（如x_logprob）。该分发器支持用户注册新的分布实现：

@singledispatch def _logprob(op, values, *inputs, **kwargs): raise NotImplementedError(f"Logprob method not implemented for {op}")

这正是整个pymc/logprob/目录下各文件（binary、censoring、arithmetic、mixture、linalg、scan、switch、tensor、order、cumsum、checks、transforms等，见 pymc/logprob/init.py）通过「注册 rewrite」来扩展可测操作覆盖面的基础。

三、pm.logcdf：对数累积分布函数

pm.logcdf(rv, value, warn_rvs=True)（basic.py）返回随机变量在给定取值处的对数 CDF（即 log P(X ≤ value)）。对数化处理在尾部概率极小时比直接计算 CDF 更稳定，也是后续logccdf的基础。

3.1 基本用法

import pymc as pm import pytensor.tensor as pt mu = pt.scalar("mu") rv = pm.Normal.dist(mu, 1.0) value = pt.scalar("value") rv_logcdf = pm.logcdf(rv, value) print(rv_logcdf.eval({value: 0.9, mu: 0.0})) # -0.2034146 rv_logcdf_fn = pm.compile([value, mu], rv_logcdf) print(rv_logcdf_fn(value=0.9, mu=0.0)) # -0.2034146

3.2 变换变量的 logcdf

与logp一致，pm.logcdf也能处理随机变量的确定性变换（basic.py）：

exp_rv = pt.exp(rv) # 对数正态 exp_rv_logcdf = pm.logcdf(exp_rv, value) print(exp_rv_logcdf.eval({value: 0.9, mu: 0.0})) # -0.78078813

3.3 在CustomDist中定义 logcdf

def normal_logcdf(value, mu, sigma): return pm.logcdf(pm.Normal.dist(mu, sigma), value) with pm.Model() as model: mu = pm.Normal("mu") sigma = pm.HalfNormal("sigma") pm.CustomDist("x", mu, sigma, logcdf=normal_logcdf)

3.4 底层的_logcdf分发器

pm.logcdf的实现路径与logp完全平行：先尝试_logcdf_helper(rv, value)，失败则走construct_ir_fgraph→_logcdf_helper→cleanup_ir（basic.py）。_logcdf同样是一个singledispatch分发器（abstract.py），默认抛NotImplementedError，由各分布模块按op类型注册具体实现。

补充说明：仓库中还提供了pm.logccdf（对数互补 CDF，即对数生存函数 log(1 - CDF(x))，basic.py）。它优先使用分布注册的数值稳定实现，否则回退到log1mexp(logcdf)以避免尾部精度损失（abstract.py）。该函数已随__all__导出（pymc/logprob/init.py），但未列入本文所依据的 API 文档页面。

四、pm.icdf：逆累积分布函数（分位数函数）

pm.icdf(rv, value, warn_rvs=True)（basic.py）返回随机变量在给定累积概率（0~1 区间）处的逆 CDF，即分位数。注意与logp/logcdf不同，这里的value表示概率而非取值点，其 dtype 允许与rv不同（转换时使用dtype="floatX"，见 basic.py）。

import pymc as pm import pytensor.tensor as pt mu = pt.scalar("mu") rv = pm.Normal.dist(mu, 1.0) value = pt.scalar("value") # 概率值 rv_icdf = pm.icdf(rv, value) # 标准正态分布的 0.9 分位数 print(rv_icdf.eval({value: 0.9, mu: 0.0})) # 1.28155157 rv_icdf_fn = pm.compile([value, mu], rv_icdf) print(rv_icdf_fn(value=0.9, mu=0.0)) # 1.28155157

同样支持随机变量的变换（basic.py）：

exp_rv = pt.exp(rv) # 对数正态的 0.9 分位数 exp_rv_icdf = pm.icdf(exp_rv, value) print(exp_rv_icdf.eval({value: 0.9, mu: 0.0})) # 3.60222448

其底层同样通过_icdf_helper→_icdf分发器实现（abstract.py）。对于离散分布，icdf也具备实现——tests/logprob/test_basic.py中的test_icdf_discrete测试用例验证了这一点。

五、pymc.logprob.conditional_logp：因子化联合对数概率

这是整个 logprob API 中最核心也最强大的函数。conditional_logp(rv_values, warn_rvs=True, ir_rewriter=None, extra_rewrites=None, **kwargs)接收一个dict，将每个随机变量映射到其值变量，返回dict，将每个值变量映射到其对应的条件对数概率项，且这些项之和等于联合对数概率（basic.py）。

5.1 动机：分层模型的条件概率

docstring 给出了一个经典的分层结构示例（basic.py）：

import pytensor.tensor as pt sigma2_rv = pt.random.invgamma(0.5, 0.5) Y_rv = pt.random.normal(0, pt.sqrt(sigma2_rv))

其对应的分层模型为：

σ² ~ InvGamma(0.5, 0.5) Y ~ N(0, σ²)

• 若只对Y_rv指定值变量y_vv = pt.scalar("y")，则conditional_logp({Y_rv: y_vv})得到的是条件对数概率 log p(Y | σ²)(y | s²)，其中sigma2_rv仍是随机变量；
• 若同时对sigma2_rv指定值变量s2_vv = pt.scalar("s2")，则conditional_logp({Y_rv: y_vv, sigma2_rv: s2_vv})会返回两个变量的条件对数概率项，其和即为联合对数概率：

log p(Y, σ²)(y, s²) = log p(Y|σ²)(y | s²) + log p(σ²)(s²)

5.2 返回结构与参数语义

• rv_values: dict[Variable, Variable]：随机变量 → 值变量（测度空间输入参数）的映射，定义了一个联合概率图；
• warn_rvs=True：当 logp 图中出现未在rv_values中指定值变量的随机变量时发出警告；
• ir_rewriter：用于生成可测变量 IR 的 rewriter；
• extra_rewrites：附加应用的重写（如重参数化、变换等），transformed_conditional_logp正是通过该参数注入值变换的；
• 返回值values_to_logps: dict：每个值变量 → 从对应随机变量推导出的条件对数概率项，顺序与输入一致；
• 若某个值变量的 logprob 项无法推导，会抛出RuntimeError并列出缺失项；若同一个值变量被赋了多个 logprob 项，则抛出ValueError（basic.py）。

5.3 实现机制：IR 图上的拓扑遍历

从实现看（basic.py），conditional_logp的执行流程是：

1. construct_ir_fgraph(rv_values, ir_rewriter=ir_rewriter)构造 IR 函数图；
2. 若传入extra_rewrites，先对其执行rewrite(fgraph)；
3. 预先把所有值变量及其祖先加入replacements映射（避免后续克隆破坏身份）；
4. 按fgraph.toposort()拓扑序遍历图中的每个节点，仅处理MeasurableOp类型的节点，通过get_related_valued_nodes找到关联的「可测变量-值变量」配对；
5. 将可测变量替换为对应的值变量后，调用_logprob(node.op, node_values, *node_inputs, **kwargs)生成该节点的条件对数概率项；
6. 最终用cleanup_ir统一清理，保证输出顺序与输入一致。

这里的MeasurableOp（abstract.py）定义了「输出可被赋予测度/对数概率的操作」，RandomVariable是它最基础的注册类型；而ValuedRV（abstract.py）则代表「可测变量与其取值」的配对 (Y, y)，它的引入既标识了多个相互依赖可测变量之间的条件点，又防止了跨条件点的自动重写。

5.4 在模型与测试中的印证

这一 API 正是pm.Model内部计算对数似然的基础。tests/logprob/test_basic.py中的test_factorized_joint_logprob_basic、test_factorized_joint_logprob_multi_obs、test_factorized_joint_logprob_diff_dims、test_hierarchical_logp、test_hierarchical_obs_logp、test_warn_rvs_conditional_logp等测试用例，覆盖了基本因子化、多观测、不同维度以及分层模型等多种场景，验证了「各项之和等于联合对数概率」的核心性质。

六、pymc.logprob.transformed_conditional_logp：带变换的联合对数概率

transformed_conditional_logp(rvs, *, rvs_to_values, rvs_to_transforms, jacobian=True, **kwargs)是conditional_logp的薄封装（basic.py），专门用于在计算联合对数概率的同时处理值变换及其雅可比修正：

• rvs：需要返回 logprob 项的随机变量序列；
• rvs_to_values：随机变量 → 值变量的完整映射（所有映射都必须提供）；
• rvs_to_transforms：随机变量 → 变换对象（Transform）的映射；
• jacobian=True：是否在变换后的 logprob 中计入雅可比行列式的对数修正；
• 返回值：仅包含rvs对应项的 logprob 列表，顺序与rvs一致，多余的中间值变量会被过滤掉。

6.1 内部原理：TransformValuesRewrite 与雅可比修正

实现中，该函数把非空的变换映射构造成TransformValuesRewrite(values_to_transforms)作为extra_rewrites传入conditional_logp（basic.py）。

TransformValuesRewrite（pymc/logprob/transform_value.py）负责在 IR 图中注入两类特殊算子：

• TransformedValue：一个 no-op，把原始值与其变换后的版本配对；
• TransformedValueRV：标识「值被变换过」的随机变量，其_logprob注册实现（transform_value.py）会先调用底层分布的_logprob，然后按transform.log_jac_det计算对数雅可比行列式并叠加。当多变量变换（如 Simplex、Ordered）作用于单变量分布、导致雅可比维度小于 logp 维度时，还会对 logp 的尾部维度求和以对齐；而维度相反的不支持场景（单变量变换作用于多变量分布）则会抛出NotImplementedError。

这正是pm.Model中「变换空间采样」的核心支撑：模型在变换后的无约束空间上做 HMC/NUTS 采样时，必须通过这种方式把联合对数概率从原始空间折算到变换空间。

6.2 与conditional_logp的差异总结

七、实践要点与常见陷阱

1. 返回的是符号图而非数值：所有函数都返回 PyTensorVariable，必须用.eval({...})或pm.compile([...], expr)求值；直接print只会看到计算图结构。这是 docstring 中反复强调的两种标准求值方式。
2. value的类型与 dtype：logp/logcdf会把非Variable输入按rv.dtype转换（basic.py、basic.py），而icdf按floatX转换（basic.py），因为分位数函数的输入是概率而非观测值。
3. warn_rvs警告的含义：当 logp 图中残留了未绑定值变量的随机变量时（例如分层模型只绑定了部分变量），PyMC 会警告这些 RV 是被克隆的、与原始对象身份不一致。解决方式是先用model.replace_rvs_by_values替换，或改用conditional_logp显式声明所有变量的值变量（basic.py）。
4. conditional_logp要求映射完整：rvs_to_values与rvs_to_transforms映射缺一不可；transformed_conditional_logp会在图中发现未替换的随机变量时直接抛出ValueError，提示可能是混用了不同模型的变量，或CustomDist/Interval 变换函数引用了非局部变量（basic.py）。
5. 自定义分布注册：如果希望某个自定义操作也能参与logp/logcdf/icdf推导，可以在_logprob/_logcdf/_icdf这三个singledispatch分发器上注册新的实现（abstract.py），并在pymc/logprob/__init__.py的「Add rewrites to the DBs」区域（pymc/logprob/init.py）引入对应模块以完成 rewrite 注册。

八、测试验证与进一步探索

仓库的测试套件是理解这些 API 行为边界的最佳参考：

• tests/logprob/test_basic.py：覆盖logp/logcdf/icdf/conditional_logp的核心行为，包括联合对数概率的因子化（test_factorized_joint_logprob_basic）、多观测（test_factorized_joint_logprob_multi_obs）、不同维度（test_factorized_joint_logprob_diff_dims）、分层模型（test_hierarchical_logp、test_hierarchical_obs_logp）、RV 残留警告（test_warn_rvs_conditional_logp）以及离散分布 icdf（test_icdf_discrete）；
• tests/logprob/test_abstract.py：验证_logcdf_helper等底层 helper 的行为；
• tests/logprob/test_composite_logprob.py：组合表达式（如混合、比较运算）的对数概率推导；
• pymc/logprob/ 目录下的binary.py、censoring.py、arithmetic.py、mixture.py、linalg.py、scan.py、switch.py、tensor.py、order.py、cumsum.py、checks.py、transforms.py则展示了各类可测操作（位运算比较、审查、混合、矩阵运算、scan 循环、变换等）是如何通过注册 rewrite 融入统一推导框架的。

通过本文介绍的五个 API，你可以脱离采样流程、直接以编程方式操纵 PyMC 的概率计算内核：从单变量对数密度，到任意变换的分布函数，再到分层模型的完整因子化联合概率——这套计算图既服务于pm.Model内部的推理引擎，也向开发者开放，成为自定义分布、似然诊断与贝叶斯工具开发的可靠地基。

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