PyMC Simulator 分布详解用似然无关的随机模拟器驱动 SMC-ABC 近似贝叶斯推断【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc导读本文围绕 PyMC 中的Simulator分布展开讲解如何把任意的 Python 随机模拟函数包装成一个概率分布并借助pm.sample_smc()以近似贝叶斯计算ABC的方式完成参数推断。读完本文你将掌握Simulator的全部构造参数fn、params、distance、sum_stat、epsilon等、内置与自定义距离/汇总统计量的用法以及它背后的实现原理与测试验证方式。本文以 API 参考文档 docs/source/api/distributions/simulator.rst 为骨架结合 pymc/distributions/simulator.py 源码与其测试 tests/distributions/test_simulator.py 展开。一、Simulator是什么为没有显式似然的模型提供推断通道Simulator是 PyMC 中一个特殊的分布类官方用途非常明确用于配合pm.sample_smc()进行近似贝叶斯推断Approximate Bayesian Inference。许多科学模型只能正向模拟——给定参数就能生成数据却写不出解析的似然函数。Simulator正是为这类模型设计的它由用户提供一个随机模拟函数fn并定义一个随机伪对数似然stochastic pseudo-loglikelihood其取值由观测数据与模拟数据之间的某种距离度量决定再用超参数epsilon对距离进行缩放。这样PyMC 的 SMC 采样器就可以在没有显式对数似然的情况下依据模拟结果与观测的吻合程度来筛选、加权粒子最终得到参数的后验分布。从源码结构看pymc/distributions/simulator.pySimulator继承自pymc.distributions.distribution.Distribution并配套定义了SimulatorRV继承 PyTensor 的RandomVariable作为其随机变量类型。这意味着 Simulator 变量与普通 PyMC 随机变量一样可以进入模型图参与logp计算、先验/后验预测采样。二、核心参数详解完整参数表Simulator的构造参数在源码 docstring 与dist方法签名中均有完整定义pymc/distributions/simulator.py如下表所示参数类型/默认值说明fncallablePython 随机模拟函数签名必须为(rng, arg1, arg2, ..., argN, size)其中rng是numpy.random.Generatorsize定义期望样本大小*unnamed_paramslist of TensorVariable按位置顺序传给fn的参数例如param1, param2, ..., paramNparamslist of TensorVariableunnamed_params的关键字形式两者只能二选一同时传入会抛出ValueErrordistancePyTensor Op / callable / str默认gaussian距离函数可选gaussian、laplace、kullback_leibler或用户自定义函数/Op自定义函数需接收epsilon、观测数据汇总统计、模拟数据汇总统计三个输入sum_statPyTensor Op / callable / str默认identity汇总统计量函数可选identity、sort、mean、median自定义函数或 Op应返回一维 numpy 数组或 PyTensor 向量epsilontensor_like of float默认1.0距离函数的缩放超参数可以是标量也可以是和sum_stat输出等长的数组ndim_suppint默认0SimulatorRV 的支撑维度0 表示标量1 表示向量等ndims_paramslist of int可选每个参数的最小维度数例如两个标量输入对应[0, 0]默认是长度为参数个数的全 0 列表class_namestr可选包装 Simulator 方法的 RandomVariable 类的后缀名2.1 关于fn的签名约定fn是 Simulator 的灵魂。它必须接受一个numpy.random.Generator作为第一个参数并在内部调用该 rng 的随机方法后续参数与size按位置传入。官方文档示例def simulator_fn(rng, loc, scale, size): return rng.normal(loc, scale, sizesize)测试中的等价写法tests/distributions/test_simulator.py也印证了这一约定staticmethod def normal_sim(rng, a, b, size): return rng.normal(a, b, sizesize)2.2 参数传递规则位置参数与params互斥源码 pymc/distributions/simulator.py 明确规定了参数传递逻辑if params is None: params unnamed_params else: if unnamed_params: raise ValueError(Cannot pass both unnamed parameters and params)对应测试test_simulator_error_msg验证了三种错误用法会抛出异常tests/distributions/test_simulator.py无效的distance名、无效的sum_stat名、以及同时传入位置参数与params。这提醒我们使用 Simulator 时参数传递方式必须二选一。三、距离函数与汇总统计量伪对数似然的三个旋钮Simulator 的伪对数似然由三部分组成epsilon缩放、sum_stat汇总、distance距离度量。三者共同决定了模拟数据离观测数据多近才算好。3.1 内置距离函数及其数学形式gaussian高斯核-0.5 * ((obs_data - sim_data) / epsilon) ** 2laplace拉普拉斯核-abs((obs_data - sim_data) / epsilon)kullback_leibler近似 KL 散度(d/n) * (1/epsilon) * sum(-log(nu_d / rho_d)) log_r出自 Pérez-Cruz, F. (2008) 对连续分布 KL 散度的估计方法其中gaussian、laplace在源码中直接以 Python 函数形式实现pymc/distributions/simulator.pyKullbackLeibler则基于scipy.spatial.cKDTree的近邻距离估计实现pymc/distributions/simulator.py。值得注意的两个组合等价关系官方文档明确给出distancegaussiansum_statsort等价于一维 2-Wasserstein 距离distancelaplacesum_statsort等价于一维 1-Wasserstein 距离。这意味着通过组合内置选项就能得到具有明确统计意义的距离度量无需自行实现。注意在当前仓库中kullback_leibler分支仍未完成重构源码会在选择该选项时直接抛出NotImplementedErrorpymc/distributions/simulator.py对应测试test_automatic_use_of_sort也标记为xfail。因此目前可稳定使用的是gaussian与laplace。3.2 内置汇总统计量identity原样返回即不做汇总默认sort对数据排序映射到 PyTensor 的pt.sortmean均值pt.meanmedian中位数——PyTensor 本身缺失该算子源码通过create_sum_stat_op_from_fn(np.median)包装实现pymc/distributions/simulator.py。3.3 自定义距离与汇总统计量自动包装机制distance与sum_stat都接受自定义 callable源码提供了两套自动包装逻辑create_sum_stat_op_from_fnpymc/distributions/simulator.py先用虚拟输入调用该函数若返回TensorVariable则直接作为图内函数使用否则包装为 PyTensorOp在perform中用np.atleast_1d(fn(x))计算结果。create_distance_op_from_fnpymc/distributions/simulator.py同样的逻辑但虚拟输入为(epsilon, obs_data, sim_data)三个参数。测试 tests/distributions/test_simulator.py 中的自定义示例可直接复用staticmethod def abs_diff(eps, obs_data, sim_data): return np.mean(np.abs((obs_data - sim_data) / eps)) staticmethod def quantiles(x): return np.quantile(x, [0.25, 0.5, 0.75])test_custom_dist_sum_stat与test_custom_dist_sum_stat_scalar验证了自定义函数在float32/float64两种精度下、以及标量观测场景下的可用性tests/distributions/test_simulator.py。四、从入门到实战完整可运行示例官方文档给出了最简示例用一个正态模拟器推断位置参数loc与尺度参数scale的后验。import pymc as pm def simulator_fn(rng, loc, scale, size): return rng.normal(loc, scale, sizesize) with pm.Model() as m: loc pm.Normal(loc, 0, 1) scale pm.HalfNormal(scale, 1) simulator pm.Simulator(simulator, simulator_fn, loc, scale, observeddata) idata pm.sample_smc()4.1 结合sample_smc的推断流程Simulator与pm.sample_smc()是固定搭配。sample_smc的完整签名位于 pymc/smc/sampling.py常用参数包括draws每阶段粒子数默认为 2000kernelSMC 内核默认IMH独立 Metropolis-Hastingschains/cores链数与并行核数默认取 CPU 数与 2 中的较大者thresholdIMH 内核参数决定各阶段间逆温度beta的步长介于 0 与 1 之间越大阶段越多默认 0.5correlation_thresholdIMH 内核参数控制自动计算的 MCMC 步数越小步数越多默认 0.01。SMC 通过一系列逆温度beta从 0 到 1逼近目标分布其中beta0时退化为先验、beta1时为后验采样过程在粒子权重重采样与 MCMC 链的交替中推进。源码 pymc/smc/sampling.py 给出了完整的十步算法摘要可据此理解为什么 Simulator 的随机伪对数似然能够被 SMC 的加权机制平滑处理。4.2 多种距离/汇总统计量的组合实战测试test_multiple_simulators展示了同一模型内并存两个不同配置的 Simulator并验证各自的 logp 图确实使用了对应的算子tests/distributions/test_simulator.pysim1 pm.Simulator( sim1, normal_sim, a, 0.1, distancegaussian, sum_statsort, observeddata1, ) sim2 pm.Simulator( sim2, normal_sim, b, 0.1, distancelaplace, sum_statmean, epsilon0.1, observeddata2, )测试还会检查编译后的 logp 图中是否包含SortOp节点sum_statsort的 sim1 应包含sum_statmean的 sim2 不应包含这为调试我的距离/统计量到底有没有生效提供了可验证的思路。4.3 嵌套 Simulator模拟结果作为另一个模拟器的输入test_nested_simulators展示了更高级的用法——把未观测的sim1无observed作为sim2的参数从而建模模拟套模拟的层级结构tests/distributions/test_simulator.pysim1 pm.Simulator( sim1, normal_sim, params(0, 4), distancegaussian, sum_statidentity, ) sim2 pm.Simulator( sim2, normal_sim, params(sim1, 0.1), distancegaussian, sum_statmean, epsilon0.1, observeddata, )五、底层原理Simulator 如何被编译进 PyTensor 计算图5.1SimulatorRV动态生成的 RandomVariable 子类SimulatorRV是所有 Simulator 随机变量的基类pymc/distributions/simulator.py它把用户传入的fn、_distance、_sum_stat、epsilon挂到类属性上并提供统一的rng_fn、distance、sum_stat三个类方法。Simulator.dist在内部完成字符串选项到可调用对象的解析后会调用rv_op通过type()动态创建一个以class_name命名、继承SimulatorRV的新类并把fn、distance、sum_stat、epsilon、签名等全部注入其中最终以该 Op 构造随机变量pymc/distributions/simulator.py。Simulator.__new__会自动把class_name默认设置为fSimulator_{name}即与变量名绑定pymc/distributions/simulator.py。5.2 伪对数似然的计算路径关键实现在simulator_logp注册于SimulatorRV的_logprob见 pymc/distributions/simulator.py其核心逻辑为创建新的共享随机数生成器simulator_rng避免并行采样时随机性退化用与原 RV 相同的输入重新构造一个sim_value模拟输出返回op.distance(op.epsilon, op.sum_stat(value), op.sum_stat(sim_value))——即观测值汇总统计与模拟值汇总统计之间的距离。这就是文档所说由观测与模拟数据之间的距离度量定义、并由epsilon调节的随机伪对数似然的具体实现。也正因如此Simulator 的 logp 求值天然带有随机性——测试test_upstream_rngs_not_in_compiled_logp验证了对同一输入反复求值结果不同且编译图中只保留一个共享 RNGtests/distributions/test_simulator.py。5.3 支持点support point用 10 次模拟的均值初始化Simulator还注册了_support_pointpymc/distributions/simulator.py通过额外 10 次模拟取平均作为参数初始化的支持点。test_simulator_support_point用 z 检验验证了该支持点与 10 次采样均值的期望一致tests/distributions/test_simulator.py支持标量、向量参数与不同size的组合。5.4 采样器兼容性虽然官方推荐搭配sample_smc但 Simulator 并不排斥传统 MCMCtest_simulator_metropolis_mcmc用pm.Metropolis对包含 Simulator 的模型采样并恢复了真实参数tests/distributions/test_simulator.py说明在似然可计算的前提下Simulator 的伪对数似然也可以被步进方法直接使用。六、常见错误与规避综合源码与测试以下用法会直接报错无效的距离名称distancenot_real抛ValueError: The distance metric not_real is not implemented无效的汇总统计名称sum_statnot_real抛ValueError: The summary statistic not_real is not implemented位置参数与params混用抛ValueError: Cannot pass both unnamed parameters andparams选择kullback_leibler当前版本抛NotImplementedError待重构。此外两个实战建议其一epsilon应依据sum_stat输出的尺度调整必要时可传入与汇总统计等长的数组逐元素缩放其二distance/sum_stat的自定义函数若不能天然映射为 PyTensor 张量运算会被自动包装成 Op 执行此时建议用测试中的模式先以虚拟输入探测是否返回TensorVariable验证其图内行为。七、延伸阅读API 入口与同目录其他分布页面docs/source/api/distributions.rstSMC 采样器完整说明与内核参数docs/source/api/smc.rst源码实现pymc/distributions/simulator.py、pymc/smc/sampling.py测试用例含多模拟器、嵌套模拟器、自定义距离等场景tests/distributions/test_simulator.py本文所依据的 API 文档页面由 autoclass 模板生成docs/source/_templates/distribution.rst【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考