PyMC 概率编程内部机制详解从 Distribution 到 logp、模型上下文与推理的开发者指南【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc本篇指南围绕 PyMC 仓库中的开发者文档 docs/source/contributing/developer_guide.md 展开深入剖析概率编程语言PPL在 PyMC 中的设计与实现分布类如何构建随机变量、logp如何被求值、pm.Model上下文管理器如何充当磁带机、以及这些机制如何支撑 MCMC、变分推断VI与前向采样。读完你将理解with pm.Model():背后发生了什么并能在当前仓库源码中定位对应的实现。PyMC 是一个基于 PyTensor 的 Python 贝叶斯统计建模包。本文面向希望理解其内部设计的开发者读者将掌握随机变量的符号图表示、条件对数概率函数的编译、CompoundStep的采样循环以及 VI 模块中 Operator/Approximation/TestFunction 的抽象关系。Distribution概率分布的类层次结构PyMC 中的概率分布实现为继承自pymc.Continuous或pymc.Discrete的类二者又共同继承自pymc.Distribution由后者定义高层 API。以当前仓库源码 pymc/distributions/distribution.py 为例模块导出了Continuous、Discrete、Distribution、DiracDelta、SymbolicRandomVariable等核心符号__all__ [ Continuous, DiracDelta, Discrete, Distribution, SymbolicRandomVariable, ]对于希望实现新分布如自定义似然或新概率族的开发者仓库提供了专门的实现指南 docs/source/contributing/implementing_distribution.md其中详细说明了rv_op的注册、logp分发函数的编写与测试要求。从源码结构看DistributionMeta元类会在类定义时读取rv_op属性若它是一个RandomVariable则自动确定rv_type并在后续根据rv_type是否为SymbolicRandomVariable创建相应的分发函数dispatch functions。也就是说一个分布类的核心工作就是声明它对应的 PyTensor 随机算子Op以及该算子参数的取值方式。反射机制调用分布构造器时发生了什么PyMC 中模型变量通过调用概率分布类并传入参数来定义其记号与数学记号一一对应z Normal(z, 0, 5)等价于数学表达式$$ z \sim \text{Normal}(0, 5) $$该调用发生在pm.Model()的上下文中模型会拦截部分信息例如已知维度。一次构造调用返回的是一个 PyTensorTensorVariable——它是模型变量及其输入依赖图的符号表示。底层而言变量是通过Distribution.distAPI 创建的该 API 会调用与分布对应的 PyTensorRandomVariableOp。RandomVariableOp 的核心思想是创建可以与概率分布属性相关联的符号变量TensorVariable。例如进入符号计算图的RandomVariableOp 与分布的随机数生成器、概率质量/密度函数相关联。当我们创建z为 $\text{Normal}(0, 5)$ 随机变量后可以获得对应的RandomVariableOp 实例句柄with pm.Model(): z pm.Normal(z, 0, 5) print(type(z.owner.op)) # pytensor.tensor.random.basic.NormalRV isinstance(z.owner.op, pytensor.tensor.random.basic.RandomVariable) # True因为NormalRV可以与 Normal 分布的概率密度函数关联我们就能通过特殊的pm.logp函数对它求值。pm.logp在仓库中的实现位于 pymc/logprob/basic.py签名核心为logp(rv, value, ...)with pm.Model(): z pm.Normal(z, 0, 5) symbolic pm.logp(z, 2.5) numeric symbolic.eval() # array(-2.65337645)也可以手算验证$$ \begin{aligned} pdf_{\mathcal{N}}(\mu, \sigma, x) \frac{1}{\sigma \sqrt{2 \pi}} \exp^{- 0.5 (\frac{x - \mu}{\sigma})^2} \ pdf_{\mathcal{N}}(0, 5, 2.5) 0.070413 \ ln(0.070413) -2.6533 \end{aligned} $$在概率编程语境下这一机制使 PyMC 及其后端 PyTensor 能够创建并求值计算图以计算例如对数先验或对数似然值。关于distAPI当前源码 pymc/distributions/distribution.py 给出了清晰的参数契约shape新 RV 各维度大小的元组、return_next_rng为True时返回(next_rng, rv)元组、size/dtype等关键字会转发给 PyTensor 的 RV Op。注意dist不支持initval与dims参数会分别抛出TypeError与NotImplementedError且shape与size不能同时传入。PyMC 与其他 PPL 的比较在 PyMC 的模型上下文中随机变量本质上就是 PyTensor 张量可以像 NumPy 数组一样参与各种运算。这与 TensorFlow ProbabilityTFP和 Pyro 不同——后两者需要更显式地在随机变量与张量之间转换。考虑如下模型$$ \begin{aligned} z \sim \mathcal{N}(0, 5) \ x \sim \mathcal{N}(z, 1) \ \end{aligned} $$PyMCwith pm.Model() as model: z pm.Normal(z, mu0., sigma5.) # pytensor.tensor.var.TensorVariable x pm.Normal(x, muz, sigma1., observed5.) # pytensor.tensor.var.TensorVariable # z2.5 的对数先验 pm.logp(z, 2.5).eval() # -2.65337645 # 在给定 z2.5 下 x 的条件对数概率 x.logp({z: 2.5}) # -4.0439386 # 整个模型的对数概率 model.logp({z: 2.5}) # -6.6973152注上述x.logp(...)、model.logp(...)的调用风格取自开发者文档中的早期示例。在当前版本中Model.logp返回的是 logp 计算图而非立即求值详见 pymc/model/core.py实际计算需配合compile_logp或initial_point使用。示例中的数值用于展示语义等价性。TensorFlow Probabilityimport tensorflow.compat.v1 as tf from tensorflow_probability import distributions as tfd with tf.Session() as sess: z_dist tfd.Normal(loc0., scale5.) # class tfp.python.distributions.normal.Normal z z_dist.sample() # class tensorflow.python.framework.ops.Tensor x tfd.Normal(locz, scale1.).log_prob(5.) # class tensorflow.python.framework.ops.Tensor model_logp z_dist.log_prob(z) x print(sess.run(x, feed_dict{z: 2.5})) # -4.0439386 print(sess.run(model_logp, feed_dict{z: 2.5})) # -6.6973152Pyroz_dist dist.Normal(loc0., scale5.) # class pyro.distributions.torch.Normal z pyro.sample(z, z_dist) # class torch.Tensor # 重置/指定 z 的值 z.data torch.tensor(2.5) x dist.Normal(locz, scale1.).log_prob(5.) # class torch.Tensor model_logp z_dist.log_prob(z) x x # -4.0439386 model_logp # -6.6973152三者的关键差异PyMC 的随机变量天然就是符号张量无需像 TFP 那样通过Session/feed_dict注入值也无需像 Pyro 那样直接修改.data属性来指定值。logp 函数幕后logp本身很直接——它是每个分布内部的一个 PyTensor 函数。开发者文档给出了如下示意签名文档明确标注代码块已过时仅用于理解思想def logp(self, value): # 获取参数 param1, param2, ... self.params1, self.params2, ... # 计算对数似然函数所有输入都是或可转换为PyTensor 张量 total_log_prob f(param1, param2, ..., value) return total_log_prob在logp方法中参数和值要么是 PyTensor 张量要么可以转换为张量。便利之处在于 logp 的求值本身也表示为张量RV.logpt当我们把不同的logp链接起来例如对所有RVs.logpt求和得到模型总 logp时依赖关系由 PyTensor 在构建与编译图时自动处理。由于编译后的函数依赖图中已有的节点每当需要生成一个以新输入张量为参数的新函数时要么重新生成带正确依赖的图要么通过编辑现有图来替换节点。PyMC 在需要时采用第二种方法即pytensor.clone_replace()。在pm.Model()上下文中分布会自动变成带分布属性的张量即 PyMC 随机变量。要获取自由随机变量free_RV的 logp只需对它自身求值logp()# self 是附着了分布的 pytensor.tensor self.logp_sum_unscaledt distribution.logp_sum(self) self.logp_nojac_unscaledt distribution.logp_nojac(self)对于观测随机变量observed RV则在数据上求 logpself.logp_sum_unscaledt distribution.logp_sum(data) self.logp_nojac_unscaledt distribution.logp_nojac(data)当前仓库中这一机制已演化为更精细的transformed_conditional_logp见 pymc/logprob/basic.py与conditional_logp同文件它们处理 RV 到 value 变量的替换、变换transform与雅可比修正并由 pymc/model/core.py 中的Model.logp统一调用。Model 上下文与随机变量with pm.Model() ...是 PyMC 模型语言的签名式语法。本质上Python 上下文管理器见 PEP 343的行为是with EXPR as VAR: USERCODE粗略等价于VAR EXPR VAR.__enter__() try: USERCODE finally: VAR.__exit__()那么在with pm.Model() as model: ...块内除了初始的model pm.Model()之外还发生了什么随机变量的三分类从上面的讨论可知在 Model 上下文中调用pm.Normal(x, ...)会返回一个随机变量。因此向模型添加随机变量有两种等价方式with pm.Model() as m: x pm.Normal(x, mu0., sigma1.) print(type(x)) # class pytensor.tensor.var.TensorVariable print(m.free_RVs) # [x] print(logpt(x, 5.0)) # Elemwise{switch,no_inplace}.0 print(logpt(x, 5.).eval({})) # -13.418938533204672 print(m.logp({x: 5.})) # -13.418938533204672一般而言如果变量带观测observed参数它是观测随机变量observed RV否则如果带transform属性它是变换随机变量transformed RV否则就是最基础的自由随机变量free RV。注意带观测的随机变量不能再被变换。pm.Model 的磁带机职责pm.Model就像一台磁带机记录被添加到模型中的一切它跟踪随机变量观测或未观测、potential 项需要加入模型 logp 的额外张量以及确定性变换作为记账信息named_varsfree_RVsobserved_RVsdeterministicspotentialsmissing_values模型上下文随后计算一些简单的模型属性构建一个在字典与 numpy/PyTensor ndarray 之间转换的双射映射bijection从而让logp/dlogp函数拥有两种等价形式一种以dict为输入另一种以ndarray为输入。更重要的是pm.Model()包含一些方法来编译以同一模型内初始化的随机变量为输入的 PyTensor 函数with pm.Model() as m: z pm.Normal(z, 0., 10., shape10) x pm.Normal(x, z, 1., shape10) print(m.initial_point) print(m.dict_to_array(m.initial_point)) # m.bijection.map(m.initial_point) print(m.bijection.rmap(np.arange(20))) # {z: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]), x: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])} # [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] # {z: array([10., 11., 12., 13., 14., 15., 16., 17., 18., 19.]), x: array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])}当前仓库中initial_point的实现见 pymc/model/core.py其 docstring 明确说明返回的是变换空间中变换后变量名到数值初值的映射随机种子参数会传给pymc.pytensorf.reseed_rngs。而basic_RVs属性同文件定义为free_RVs observed_RVs即模型真正定义的随机变量项不含确定性变量。Model上以 logp 为名的可用方法非常丰富list(filter(lambda x: logp in x, dir(pm.Model))) # [d2logp, # d2logp_nojac, # datalogpt, # dlogp, # dlogp_array, # dlogp_nojac, # fastd2logp, # fastd2logp_nojac, # fastdlogp, # fastdlogp_nojac, # fastlogp, # fastlogp_nojac, # logp, # logp_array, # logp_dlogp_function, # logp_elemwise, # logp_nojac, # logp_nojact, # logpt, # varlogpt]Logp 与 dlogp 的构建模型收集所有随机变量model.free_RVs与model.observed_RVs中的一切以及 potential 项并求和得到模型 logpproperty def logpt(self): PyTensor scalar of log-probability of the model with self: factors [var.logpt for var in self.basic_RVs] self.potentials logp at.sum([at.sum(factor) for factor in factors]) ... return logp它返回一个 PyTensor 张量其值依赖于模型中的自由参数即 PyTensor 图中的父节点。注意 logp 张量依赖其输入节点因此不能传入新张量来生成 logp 函数。出于同样原因PyMC 大量使用pytensor.clone_replace复制图并替换张量的输入。with pm.Model() as m: z pm.Normal(z, 0., 10., shape10) x pm.Normal(x, z, 1., shape10) y pm.Normal(y, x.sum(), 1., observed2.5) print(m.basic_RVs) # [z, x, y] print(m.free_RVs) # [z, x] type(m.logp) # pytensor.tensor.var.TensorVariable m.logpt.eval(m.initial_point()) # array(-51.25369126)PyMC 随后编译一个带梯度的 logp 函数以model.free_RVs为输入、model.logpt为输出。如果我们想要条件 logp/dlogp 函数输入可以只是model.free_RVs的一个子集def logp_dlogp_function(self, grad_varsNone, **kwargs): if grad_vars is None: grad_vars list(typefilter(self.free_RVs, continuous_types)) else: ... varnames [var.name for var in grad_vars] # 简单情形下只有连续 RV这就是全部 free_RVs extra_vars [var for var in self.free_RVs if var.name not in varnames] return ValueGradFunction(self.logpt, grad_vars, extra_vars, **kwargs)ValueGradFunction是一个可调用类它隔离 PyTensor 图的一部分来编译额外的 PyTensor 函数。PyMC 依赖pytensor.clone_replace复制model.logpt并替换其输入——它不直接编辑或重写图。当前仓库中logp_dlogp_function位于 pymc/model/core.py。高层面上这让我们能轻松构建条件 logp 函数及其梯度。实测示例func m.logp_dlogp_function() func.set_extra_values({}) input_dict m.initial_point() print(input_dict) input_array func.dict_to_array(input_dict) print(input_array) print( ) func(input_array) # {z: array([-0.7202002 , 0.58712205, -1.44120196, -0.53153001, -0.36028732, # -1.49098414, -0.80046792, -0.26351819, 1.91841949, 1.60004128]), x: array([ 0.01490006, 0.60958275, -0.06955203, -0.42430833, -1.43392303, # 1.13713493, 0.31650495, -0.62582879, 0.75642811, 0.50114527])} # [-0.7202002 0.58712205 -1.44120196 -0.53153001 -0.36028732 -1.49098414 # -0.80046792 -0.26351819 1.91841949 1.60004128 0.01490006 0.60958275 # -0.06955203 -0.42430833 -1.43392303 1.13713493 0.31650495 -0.62582879 # 0.75642811 0.50114527] # # (array(-51.0769075), # array([ 0.74230226, 0.01658948, 1.38606194, 0.11253699, -1.07003284, # 2.64302891, 1.12497754, -0.35967542, -1.18117557, -1.11489642, # 0.98281586, 1.69545542, 0.34626619, 1.61069443, 2.79155183, # -0.91020295, 0.60094326, 2.08022672, 2.8799075 , 2.81681213])) irv 1 print(Condition Logp: take %s as input and conditioned on the rest.%(m.free_RVs[irv].name)) func_conditional m.logp_dlogp_function(grad_vars[m.free_RVs[irv]]) func_conditional.set_extra_values(input_dict) input_array2 func_conditional.dict_to_array(input_dict) print(input_array2) print( ) func_conditional(input_array2) # Condition Logp: take x as input and conditioned on the rest. # [ 0.01490006 0.60958275 -0.06955203 -0.42430833 -1.43392303 1.13713493 # 0.31650495 -0.62582879 0.75642811 0.50114527] # # (array(-51.0769075), # array([ 0.98281586, 1.69545542, 0.34626619, 1.61069443, 2.79155183, # -0.91020295, 0.60094326, 2.08022672, 2.8799075 , 2.81681213]))为什么不直接手工编译有人可能会想直接编译一个 logp 函数、自己做记账不就行了例如直接用 PyTensor 构建import pytensor func pytensor.function(m.free_RVs, m.logpt) func(*inputlist) # array(-51.0769075) logpt_grad pytensor.grad(m.logpt, m.free_RVs) func_d pytensor.function(m.free_RVs, logpt_grad) func_d(*inputlist) # [array([ 0.74230226, 0.01658948, 1.38606194, 0.11253699, -1.07003284, # 2.64302891, 1.12497754, -0.35967542, -1.18117557, -1.11489642]), # array([ 0.98281586, 1.69545542, 0.34626619, 1.61069443, 2.79155183, # -0.91020295, 0.60094326, 2.08022672, 2.8799075 , 2.81681213])]类似地构建条件 logpshared pytensor.shared(inputlist[1]) func2 pytensor.function([m.free_RVs[0]], m.logpt, givens[(m.free_RVs[1], shared)]) print(func2(inputlist[0])) # -51.07690750130328 logpt_grad2 pytensor.grad(m.logpt, m.free_RVs[0]) func_d2 pytensor.function([m.free_RVs[0]], logpt_grad2, givens[(m.free_RVs[1], shared)]) print(func_d2(inputlist[0])) # [ 0.74230226 0.01658948 1.38606194 0.11253699 -1.07003284 2.64302891 # 1.12497754 -0.35967542 -1.18117557 -1.11489642]上述方式得到的 logp 与梯度同model.logp_dlogp_function的输出一致。但困难在于把一切编译进单个函数func_logp_and_grad pytensor.function(m.free_RVs, [m.logpt, logpt_grad]) # ERROR我们想要value, grad f(x)这样一次调用同时返回求值及其关于每个输入的梯度但朴素实现做不到。当然可以把两个函数一个 logp、一个 dlogp包起来输出列表但那意味着每次要调用两个函数而且用 Python 逻辑构建条件 logp 时还得自己做记账。使用pytensor.clone_replace后PyTensor 函数的输入始终是一维向量而不是形状各异的 RV 列表因而非常便于做旋转等矩阵运算。补充说明当前的设置相当强大PyTensor 编译函数编译与调用都较快反复调用条件 logp 函数时外部 RV 只需重置一次。但 PyTensor 图与 NumPy 之间传值仍有显著开销——这正是 GPU 上常常看不到优势的原因每次函数调用都要在 GPU 与 CPU 之间拷贝数据对小模型而言 GPU 上的推理反而比 CPU 慢。另外pytensor.clone_replace过于方便PyMC 内部戏称它像毒品一样令人上瘾。如果所有操作都发生在图中包括条件化与取值那么构建模型和运行推理时就没有必要隔离部分图。而如果我们只限定在最擅长解决的问题——所有连续未知参数、可用动态 HMC 采样的模型——甚至更不需要考虑图的克隆/重写。推理MCMC 与 CompoundStep模型实例能够生成条件 logp 与 dlogp 函数这成就了 PyMC 的特色功能之一——CompoundStep复合步进。概念上它是 Metropolis-within-Gibbs 采样器用户可以为不同 RV 指定不同采样器。它也可以看作又一个拦截器pm.sample(...)调用会尝试为不同的free_RVs指派最佳步进方法例如当所有free_RVs都是连续型时使用 NUTS。随后编译条件logp 函数采样器在CompoundStep列表的 for 循环中逐个调用完成一个采样循环。每个采样器实现一个step.step方法执行 MH 更新。每次传入一个字典PyMC 中的point结构与model.initial_point相同输出一个新字典——其中被采样的free_RVs若被接受则带有新值。当前仓库的实现位于 pymc/step_methods/compound.pyCompoundStep.step依次调用列表中每个method.step(point)并汇总统计量def step(self, point) - tuple[PointType, StatsType]: stats [] for method in self.methods: point, sts method.step(point) stats.extend(sts) # model_logp 只能保留最后一个统计量中的值 for sts in stats[:-1]: sts.pop(model_logp, None) return point, stats转移核与 ArrayStep大多数 MCMC 采样器SMC 除外的基类在ArraySteppymc/step_methods/arraystep.py。可以看到step.step()把point映射为数组然后调用self.astep()——这是一个数组进、数组出的函数。PyMC 模型会编译一个条件 logp/dlogp 函数把输入 RV 替换为共享的一维张量原始 RV 的展平与堆叠视图转移核即.astep()接收数组并输出数组例如 Metropolis 采样器。这与 TFP 中的转移核截然不同TFP 的转移核是张量进、张量出且不展平张量例如tfp.mcmc.random_walk_metropolis的new_state_fn接收state parts 列表 seed返回同样类型的 Tensor 列表。动态 HMCPyMC 偏爱 NUTS更准确地说是带复杂停止规则的动态 HMC。这部分完全在 PyTensor 之外完成对 NUTS 而言包括leapfrog 积分、对偶平均、质量矩阵与步长调节、树构建、采样器相关统计量如发散与能量检查。PyMC 曾有一个 PyTensor 版本的 HMC但从未投入使用且已从主仓库移除仍保留在 git 历史中。变分推断VIVI 模块采用了与 MCMC 不同的设计思路——它是函数式的一切都在 PyTensor 内完成包括优化与变分目标构建。变分推断的基类是pymc.variational.Inferencepymc/variational/inference.py它通过如下方式构建目标函数... self.objective op(approx, **kwargs)(tf) ...其中op : Operator 类 approx : Approximation 类或实例 tf : TestFunction 实例 kwargs : 传给 Operator 的关键字参数该设计受 Operator Variational InferencearXiv:1610.09033启发。Inference是 VI 实现中非常高层的一个对象它用 Operator、Approximation、Test function 这些原语组合出单一目标函数。目前 PyMC 不太关心 test function通常不需要也未实现。其余原语在 pymc/variational/opvi.py 中定义为基类通过继承可以轻松实现一大类 VI 方法同时保留充分的扩展灵活性。以 ADVI 为例高层上我们在隐空间中用对角多元高斯近似后验即用高斯逐个近似model.free_RVs中的每个元素。其初始化过程如下def __init__(self, *args, **kwargs): super(ADVI, self).__init__(MeanField(*args, **kwargs)) # 在父类 KLqp 中 super(KLqp, self).__init__(KL, MeanField(*args, **kwargs), None, betabeta) # 在父类 Inference 中 ... self.objective KL(MeanField(*args, **kwargs))(None) ...其中KL是基于 Kullback-Leibler 散度的 Operator不需要 test function... def apply(self, f): return -self.datalogp_norm self.beta * (self.logq_norm - self.varlogp_norm)datalogp_norm、logq_norm、varlogp_norm的定义在 pymc/variational/opvi.py。剥离归一化项后datalogp与varlogp是变分自由 RV 与数据 logp 的期望——PyMC 从模型中克隆 datalogp 与 varlogp把其输入替换为从变分后验采样的 PyTensor 张量。对 ADVI 而言这些样本来自高斯分布。注意来自后验近似的样本通常多一个维度以便计算期望以及期望的梯度通过计算梯度的期望即重参数化技巧。而logq因为是高斯分布求值非常直接。实现 VI 时的一些挑战与洞见基于图的方法很有帮助但 PyTensor 无法直接访问计算图中先前创建的节点。实现中可见大量node_property用法用于缓存节点。TensorFlow 有图工具可能对此有帮助另一方面 TensorFlow 的图管理比预期更棘手高层原因是图是一个只增容器。有不少 bug 起初并不明显。PyTensor 的pytensor.clone_replace在不同层级做大量图替换时需要极其小心。PyMC 内部甚至造了 pytensor.clone_replacecurse克隆替换诅咒这个说法——对它的依赖已经到了无以复加的程度用它向量化模型MCMC 与 VI 都受益加速计算用它为 VI 创建采样图——这是你想要后验预测作为计算图一部分的场景。由于这是 VI 过程的核心PyMC 曾试图在 TF 中复刻该模式。但当pytensor.clone_replace被调用时PyTensor 会创建图的新部分可被垃圾回收而 TF 的图是只增的所以需要以不同方式解决输入替换问题。前向采样如前面所述分布中有方法遍历模型依赖图在 scipy/numpy 中生成前向随机样本。这让我们能用pymc.sampling.sample_prior_predictive做先验预测采样。它是相当快的批量操作但在高维场景下有不少 bug 和边界情况。最大的痛点是自动广播在批量随机生成时我们想生成(n_sample, ) RV.shape的随机样本。在某些情况下广播 RV1 和 RV2 生成多一个批处理维度的 RV3 时会出错更糟的是静默给出错误答案。好消息是这些问题正在被逐步修复。形状处理的挑战与解决方案总结如下with pm.Model() as m: mu pm.Normal(mu, 0., 1., shape(5, 1)) sigma pm.HalfNormal(sigma, 5., shape(1, 10)) pm.Normal(x, mumu, sigmasigma, observednp.random.randn(2, 5, 10)) trace pm.sample_prior_predictive(100) trace[x].shape # 应为 (100, 2, 5, 10) pm.Normal.dist(munp.zeros(2), sigma1).random(size(10, 4))其它与随机样本生成相关的错误也存在例如 Mixture 曾经是坏的。扩展 PyMC开发者文档给出了三类扩展方向均有社区实践案例自定义 Inference 方法例如用 EM 推断线性混合模型、PyMC 中的 Laplace 近似见 junpenglao 的 Planet_Sakaar_Data_Science 仓库中的 notebook。在模型内连接其它库通过创建自定义 PyTensor Op 使用黑盒似然函数或直接使用 emcee。用其它库做推断连接 Julia 求解 ODE其解带梯度可用于 NUTS。这些路径的共同前提是理解 PyMC 的随机变量即符号张量、logp是可编译的图、以及pytensor.clone_replace是替换图中输入的核心工具。我们曾经搞错过什么开发者文档坦诚总结了三个长期痛点理解它们有助于避免踩坑Shape形状形状问题是最常遇到的痛点。TFP 和 Pyro 在这方面的做法当时要严谨得多。PyMC 曾通过 PR 试图修复但很可能是巨大的重构工程。Mixture 分布打了大量补丁但依然不够自然。numpy 中的随机方法从随机变量采样的逻辑非常复杂而且全部在 Python 中实现因此无法进一步变换采样图。PyTensor 没有从各种分布采样的代码而 PyMC 也不愿自己重写一遍。采样器用 Python 编写虽然采样器用 Python 编写带来了极大的灵活性与实验直观性用 PyTensor 写 NUTS 也非常困难但代价是性能损失且使 GPU 上的采样非常低效——因为每次 logp 求值都需要拷贝内存。小结本文以开发者指南为主线梳理了 PyMC 概率编程内核的四块基石分布即算子Distribution→Continuous/Discrete的类层次加上rv_op与dist()API让调用分布类变成创建带分布属性的符号张量logp 即图logp是 PyTensor 计算图中的节点模型 logp 是对所有因子求和编译与梯度均由图引擎完成Model 即磁带机pm.Model上下文管理器收集free_RVs、observed_RVs、potentials、deterministics等并借助 bijection 提供 dict/array 双入口的 logp 家族推理基于条件 logplogp_dlogp_function配合pytensor.clone_replace构建条件 logp/dlogp支撑CompoundStep、动态 HMC 与 VI 的Operator/Approximation抽象。对想深入源码的读者推荐按此路径阅读pymc/distributions/distribution.py分布基类与distAPI、pymc/model/core.pylogp、basic_RVs、initial_point、logp_dlogp_function、pymc/logprob/basic.pylogp/conditional_logp/transformed_conditional_logp、pymc/step_methods/compound.pyCompoundStep以及 pymc/variational/inference.pyInference与目标函数构建。分布实现细节可继续参考 docs/source/contributing/implementing_distribution.md用户向 API 入门可见 docs/source/learn/core_notebooks/pymc_overview.ipynb。【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考