# PyMC 建模使用 PyMC v5+ 进行现代贝叶斯建模。主要默认设置：nutpie 采样器（快 2-5 倍）、分层模型的非中心参数化、精确 GP 上的 HSGP、可读推理数据的坐标/尺寸，以及提前保存工作流程以防止后期崩溃导致数据丢失。 **建模策略**：迭代地构建模型 - 从简单开始，检查先前的预测，拟合和诊断，检查后验的预测，一次扩展一个部分。请参阅 [references/workflow.md]() 了解完整的工作流程。 **笔记本偏好**：使用 marimo 进行交互式建模，除非项目已经使用 Jupyter。 ## 型号规格 ### 基本结构将 pymc 导入为 pm 将 arviz 导入为 az 以 pm.Model(coords=coords) 作为模型： # 数据容器（用于样本外预测） x = pm.Data("x", x_obs, dims="obs") # 先验 beta = pm.Normal("beta", mu=0, sigma=1, dims="features") 西格玛 = pm.HalfNormal("西格玛", 西格玛=1) # 可能性 mu = pm.math.dot(x, beta) y = pm.Normal("y", mu=mu, sigma=sigma, 观察=y_obs, dims="obs") # 推论 idata = pm.sample(nuts_sampler="nutpie", random_seed=42) ### 坐标和尺寸当模型具有有意义的结构时，使用 coords/dims 进行可解释的推理数据：坐标={ “obs”：np.arange（n_obs）， “特征”：[“截距”，“年龄”，“收入”]， “组”：组标签， } 跳过开销超过收益的简单模型。 ### 参数化对于数据较弱的分层模型，首选非中心参数化： # 非中心（更适合发散）偏移 = pm.Normal("偏移", 0, 1, dims="组") alpha = mu_alpha + sigma_alpha * 偏移量 # 居中（有强大的数据效果更好） alpha = pm.Normal("alpha", mu_alpha, sigma_alpha, dims="group") ## 推论 ### 默认采样（nutpie 首选）型号： idata = pm.sample( 绘制=1000，调整=1000，链=4， nut_sampler =“坚果派”，随机种子=42，） idata.to_netcdf("results.nc") # 采样后立即保存 **重要**：nutpie 不会自动存储 log_likelihood （它默默地忽略 `idata_kwargs={"log_likelihood": True}`）。如果需要 LOO-CV 或模型比较，请在采样后计算： pm.compute_log_likelihood（idata，模型=模型） ### 何时使用 PyMC 的默认 NUTS nutpie 无法处理离散参数或某些转换（例如，使用“OrderedLogistic”/“OrderedProbit”进行“ordered”转换）。对于这些模型，省略 `nuts_sampler="nutpie"`： idata = pm.sample(绘制=1000，调整=1000，链=4，random_seed=42) 切勿更改模型规范来解决采样器限制。如果未安装 nutpie，请安装它（`pip install nutpie`）或回退到 `nuts_sampler="numpyro"`。 ### 替代 MCMC 后端请参阅 [references/inference.md]() 了解： * **NumPyro/JAX** ：GPU 加速，矢量化链 ### 近似推理对于快速（但不精确）的后验近似： * **ADVI/DADVI** ：采用高斯近似的变分推理 * **Pathfinder** ：用于初始化或筛选的拟牛顿优化 ## 诊断和 ArviZ 工作流程 **最小工作流程清单** - 每个模型脚本应包括： 1. 事前预测检查（`pm.sample_prior_predictive`） 2. 采样后立即保存结果(`idata.to_netcdf(...)`) 3. 散度计数+r_hat+ESS检查 4.后验预测检查（`pm.sample_posterior_predictive`）每次采样运行后请遵循此系统工作流程： ### 第 1 阶段：立即检查（必需） # 1.检查分歧（必须为0或接近0） n_div = idata.sample_stats["发散"].sum().item() print(f"分歧：{n_div}") # 2. 收敛诊断摘要 Summary = az.summary(idata, var_names=["~offset"]) # 排除辅助 print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "ess_bulk", "ess_tail", "r_hat"]]) # 3.视觉收敛检查 az.plot_trace(idata, 紧凑=True) az.plot_rank(idata, var_names=["beta", "sigma"]) **通过标准**（必须全部通过才能继续）： * 零散度（或 < 0.1% 且随机分散） * 所有参数的“r_hat < 1.01” * `ess_bulk > 400` 和 `ess_tail > 400` * 迹线图显示良好的混合（重叠密度、模糊毛毛虫） ### 第 2 阶段：深度融合（如果第 1 阶段处于边缘状态） # ESS 演化（应该线性增长） az.plot_ess(idata, kind="进化") # 能量诊断（HMC 健康） az.plot_energy(idata) # 自相关（应该快速衰减） az.plot_autocorr(idata, var_names=["beta"]) ### 第 3 阶段：模型批评（必填） # 生成后验预测型号： pm.sample_posterior_predictive（idata，extend_inferencedata = True） # 模型是否捕获数据？ az.plot_ppc(idata, kind="累积") # 校准检查 az.plot_loo_pit(idata, y="y") **关键规则**：在阶段 1-3 通过之前切勿解释参数。 ### 第 4 阶段：参数解释 # 后文总结 az.plot_posterior(idata, var_names=["beta"], ref_val=0) # 分层参数的森林图 az.plot_forest(idata, var_names=["alpha"], 组合=True) # 参数相关性（识别不可识别性） az.plot_pair(idata, var_names=["alpha", "beta", "sigma"]) 请参阅 [references/arviz.md]() 了解 ArviZ 的全面用法。请参阅 [references/diagnostics.md]() 进行故障排除。 ## 事前和事后预测检查 ### 先验预测（拟合前）在安装之前一定要检查先前的影响：型号： Prior_pred = pm.sample_prior_predictive(绘制=500) az.plot_ppc(prior_pred, group="prior", kind="cumulative") Prior_y = Prior_pred.prior_predictive["y"].values.flatten() print(f"先验预测范围: [{prior_y.min():.1f}, {prior_y.max():.1f}]") **规则**：在任何新模型上的“pm.sample()”之前运行预先的预测检查。如果范围不可信（负计数，概率 > 1），请在继续之前调整先验。 ### 后验预测（拟合后）型号： pm.sample_posterior_predictive（idata，extend_inferencedata = True） az.plot_ppc(idata, kind="累积") az.plot_loo_pit(idata, y="y") 观察到的数据（黑线）应属于后验预测分布。有关详细解释，请参阅 [references/arviz.md]()。 ## 模型调试在采样之前，使用 model.debug() 和 model.point_logps() 验证模型。使用“print(model)”作为结构，使用“pm.model_to_graphviz(model)”作为 DAG 可视化。 ### 常见问题症状|可能的原因 |修复 ---|---|--- `ValueError：形状不匹配` |参数与观察维度|使用索引向量：`alpha[group_idx]` `初步评估失败` |数据外部分发支持|检查边界；使用`init =“adapt_diag”` `质量矩阵包含零` |无尺度预测变量或平坦先验 |标准化功能；使用信息较弱的先验高发散计数 |漏斗几何形状|非中心参数化对数概率中的“NaN” |无效的参数组合 |检查参数约束，添加边界 `-inf` 对数概率 |可能性支持之外的观察|验证数据与分发域匹配慢速离散采样| NUTS 与离散 | 不兼容边缘化离散变量请参阅 [references/troubleshooting.md]() 获取全面的问题解决指南。要调试分歧，请使用 az.plot_pair(idata, divergences=True) 来定位集群。请参阅 [references/diagnostics.md]() § 分歧故障排除。有关分析慢速模型的信息，请参阅 [references/troubleshooting.md]() § 性能问题。 ## 型号对比 ### LOO-CV（首选） # 使用逐点诊断计算 LOO loo = az.loo(idata, pointwise=True) print(f"ELPD: {loo.elpd_loo:.1f} ± {loo.se:.1f}") # 检查 Pareto k 值（必须 < 0.7 才能获得可靠的 LOO） print(f"坏 k (>0.7): {(loo.pareto_k > 0.7).sum().item()}") az.plot_khat(idata) ### 比较模型 # 如果使用 nutpie，首先计算对数似然（nutpie 不会自动存储它） pm.compute_log_likelihood(idata_a, model=model_a) pm.compute_log_likelihood(idata_b, model=model_b) 比较 = az.compare({ “模型_a”：idata_a， “模型_b”：idata_b， }, ic="loo") 打印（比较[[“排名”，“elpd_loo”，“elpd_diff”，“权重”]]） az.plot_compare（比较） **决策规则**：如果两个模型具有相似的堆叠权重，则它们实际上是等效的。有关详细的模型比较工作流程，请参阅 [references/arviz.md]()。 ### 迭代模型构建逐步构建复杂性：首先拟合最简单的合理模型，对其进行诊断，检查后验预测，然后一次添加一个复杂性。通过 LOO 比较每个扩展。如果堆叠权重相似，则模型实际上是等效的。有关完整的迭代工作流程，请参阅 [references/workflow.md]()。 ## 保存和加载结果 ### 推理数据持久化保存采样结果以供以后分析或共享： # 保存为NetCDF（推荐格式） idata.to_netcdf(“结果/model_v1.nc”) # 加载 idata = az.from_netcdf("结果/model_v1.nc") 有关大型 InferenceData 对象的压缩存储，请参阅 [references/workflow.md]()。 **关键**：采样后立即保存 - 后期崩溃会破坏有效结果：型号： idata = pm.sample(nuts_sampler="nutpie") idata.to_netcdf("results.nc") # 在任何后处理之前保存！型号： pm.sample_posterior_predictive（idata，extend_inferencedata = True） idata.to_netcdf("results.nc") # 使用后验预测进行更新 ## 事先选择请参阅 [references/priors.md]() 了解： * 按分布类型划分的默认信息较弱 * 事前预测检查工作流程 * 特定领域的建议 ## 常见模式 ### 分层/多级将 pm.Model(coords={"group": groups, "obs": obs_idx}) 作为分层： # 超先验 mu_alpha = pm.Normal("mu_alpha", 0, 1) sigma_alpha = pm.HalfNormal("sigma_alpha", 1) # 组级别（非中心） alpha_offset = pm.Normal("alpha_offset", 0, 1, dims="group") alpha = pm.Deterministic("alpha", mu_alpha + sigma_alpha * alpha_offset, dims="group") # 可能性 y = pm.Normal("y", alpha[group_idx], sigma, 观察=y_obs, dims="obs") ### GLM # 逻辑回归使用 pm.Model() 作为逻辑： alpha = pm.Normal("alpha", 0, 2.5) beta = pm.Normal("beta", 0, 2.5, dims="features") p = pm.math.sigmoid(alpha + pm.math.dot(X, beta)) y = pm.Bernoulli("y", p=p, 观察=y_obs) # 泊松回归将 pm.Model() 作为泊松： beta = pm.Normal("beta", 0, 1, dims="features") y = pm.Poisson("y", mu=pm.math.exp(pm.math.dot(X, beta)), 观察=y_obs) ### 高斯过程对于 1-3D 输入的 GP 问题，**始终首选 HSGP**。它是 O(nm) 而不是 O(n³)，即使在 n=200 时，精确的 GP (`pm.gp.Marginal`) 对于 MCMC 来说也慢得令人望而却步：将 pm.Model() 作为 gp_model： # 超参数 ell = pm.InverseGamma("ell", alpha=5, beta=5) eta = pm.HalfNormal("eta", sigma=2) 西格玛 = pm.HalfNormal("西格玛", 西格玛=0.5) # 协方差函数（推荐Matern52） cov = eta**2 * pm.gp.cov.Matern52(1, ls=ell) # HSGP 近似 gp = pm.gp.HSGP(m=[20], c=1.5, cov_func=cov) f = gp.prior("f", X=X[:, None]) # X 必须是 2D # 可能性 y = pm.Normal("y", mu=f, sigma=sigma, 观察=y_obs) 对于周期性模式，请使用“pm.gp.HSGPPeriodic”。仅对特别需要精确推理的非常小的数据集 (n < ~50) 使用 `pm.gp.Marginal` 或 `pm.gp.Latent`。有关 HSGP 参数选择（m、c）、HSGPPeriodic、协方差函数和常见模式，请参阅 [references/gp.md]()。 ### 时间序列使用 pm.Model(coords={"time": range(T)}) 作为 ar_model： rho = pm.Uniform("rho", -1, 1) 西格玛 = pm.HalfNormal("西格玛", 西格玛=1) y = pm.AR("y", rho=[rho], sigma=sigma, 常量=True, 观察=y_obs，暗淡=“时间”）有关 AR/ARMA、随机游走、结构时间序列、状态空间模型和预测模式，请参阅 [references/timeseries.md]()。 ### BART（贝叶斯加性回归树）将 pymc_bart 导入为 pmb 将 pm.Model() 作为 bart_model： mu = pmb.BART("mu", X=X, Y=y, m=50) 西格玛 = pm.HalfNormal("西格玛", 1) y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, 观察=y) 有关回归/分类、变量重要性和配置，请参阅 [references/bart.md]()。 ### 混合模型将 numpy 导入为 np 坐标= {“组件”：范围（K）} 将 pm.Model(coords=coords) 作为 gmm： # 混合权重 w = pm.Dirichlet("w", a=np.ones(K), dims="component") # 组件参数（按顺序排列以避免标签切换） mu = pm.Normal("mu", mu=0, sigma=10, dims="组件", 变换=pm.distributions.transforms.ordered, initval=np.linspace(y_obs.min(), y_obs.max(), K)) sigma = pm.HalfNormal("sigma", sigma=2, dims="组件") # 混合可能性 y = pm.NormalMixture("y", w=w, mu=mu, sigma=sigma, 观察=y_obs) **重要**：混合模型通常需要 `target_accept=0.9` 或更高，以避免与多模态几何结构的分歧。始终在有序方式上提供“initval” - 如果没有它，组件可能会开始重叠，并且采样器很难将它们分开。有关标签切换解决方案、边缘化混合物和混合物诊断，请参阅 [references/mixtures.md]()。 ### 稀疏回归/马蹄形使用正则化（芬兰）马蹄先验进行具有预期稀疏性的高维回归。马蹄先验创建双漏斗几何——使用“target_accept=0.95”或更高。请参阅 [references/priors.md]() 了解完整的正则化马蹄代码、拉普拉斯、R2D2 和 spike-and-slab 替代方案。 ### 专门的可能性 # 零膨胀泊松（多余的零）将 pm.Model() 作为 zip_model： psi = pm.Beta("psi", alpha=2, beta=2) # P(结构零) mu = pm.Exponential("mu", lam=1) y = pm.ZeroInflatedPoisson("y", psi=psi, mu=mu, 观察=y_obs) # 删失数据（例如，右删失生存率）将 pm.Model() 作为 censored_model： mu = pm.Normal("mu", mu=0, sigma=10) 西格玛 = pm.HalfNormal("西格玛", 西格玛=5) y = pm.Censored("y", dist=pm.Normal.dist(mu=mu, sigma=sigma), 下=无，上=审查时间，观察=y_obs） # 序数回归以 pm.Model() 作为序数： beta = pm.Normal("beta", mu=0, sigma=2, dims="features") 切点 = pm.Normal("切点", mu=0, sigma=2, 变换=pm.distributions.transforms.ordered, 形状=n_类别 - 1) y = pm.OrderedLogistic("y", eta=pm.math.dot(X, beta), 切点=切点，观察=y_obs） **注意**：不要对变量和维度使用相同的名称。例如，如果您有一个名为“cutpoints”的维度，请勿同时命名变量“cutpoints”——这会导致形状错误。请参阅 [references/specialized_likelihoods.md]() 了解零膨胀、跨栏、审查/截断、序数和稳健回归模型。 ## 常见陷阱请参阅 [references/troubleshooting.md]() 获取全面的问题解决指南，其中包括： * 形状和尺寸错误、初始化失败、质量矩阵问题 * 分歧和几何问题（中心与非中心） * PyMC API 问题（变量命名、不推荐使用的参数） * 性能问题（GP、大型确定性、重新编译） * 可识别性、多重共线性、先验数据冲突 * 离散变量挑战、数据容器、预测 ## 因果推理操作 PyMC 支持因果查询的 do-calculus： # pm.do — 干预（破坏传入的边缘）使用 pm.do(causal_model, {"x": 2}) 作为干预模型： idata = pm.sample_prior_predictive() # P(y, z | do(x=2)) # pm.observe — 条件（保留因果结构）使用 pm.observe(causal_model, {"y": 1}) 作为 conditioned_model： idata = pm.sample(nuts_sampler="nutpie") # P(x, z | y=1) # 组合：P(y | do(x=2), z=0) 将 pm.do(causal_model, {"x": 2}) 作为 m1：将 pm.observe(m1, {"z": 0}) 作为 m2： idata = pm.sample(nuts_sampler="nutpie") 有关详细的因果推理模式，请参阅 [references/causal.md]()。 ## pymc-extras 通过“import pymc_extras as pmx”进行关键扩展： * `pmx.marginalize(model, ["discrete_var"])` — 边缘化 NUTS 的离散参数 * `pmx.R2D2M2CP(...)` — R2D2 先验回归（参见 [references/priors.md]()） * `pmx.fit_laplace(model)` — 用于快速推理的拉普拉斯近似 ## 自定义分布和模型组件 # 通过势的软约束将 pytensor.tensor 导入为 pt pm.Potential("sum_to_zero", -100 * pt.sqr(alpha.sum())) 请参阅 [references/custom_models.md]() 了解“pm.DensityDist”、“pm.Potential”、“pm.Simulator”和“pm.CustomDist”。

Pymc 建模

包含技能