PyMC 贝叶斯建模工作流与常用模式实践指南从先验检查到模型比较的完整技术栈【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills本指南基于 scientific-agent-skills 仓库中 PyMC 技能模块的核心参考文档 workflows.md 编写系统讲解使用 PyMC6.x构建、验证与分析贝叶斯模型的完整工作流涵盖标准九步工作流模板、线性/逻辑/分层/泊松/时间序列/混合六类模型构建模式、数据准备与先验选择的最佳实践以及基于 PSIS-LOO 的模型比较与采样诊断排障方案。读者学完后可直接复刻一套可运行的端到端贝叶斯分析流程并结合仓库提供的 model_diagnostics.py、model_comparison.py 脚本和两份可定制模板完成实际项目落地。一、标准贝叶斯工作流九步端到端模板PyMC 技能的核心纪律是先检查后采样Never sample first and check later。SKILL.md 将其归纳为八个阶段而 workflows.md 给出了每一步的完整可运行代码。以下模板综合两者可作为任何贝叶斯回归分析的基本骨架import pymc as pm import arviz as az import numpy as np import matplotlib.pyplot as plt # 1. PREPARE DATA # X ... # Predictor variables y ... # Observed outcomes # Standardize predictors for better sampling X_scaled (X - X.mean(axis0)) / X.std(axis0) # 2. BUILD MODEL # coords { predictors: [var1, var2, var3], obs_id: np.arange(len(y)) } with pm.Model(coordscoords) as model: X_data pm.Data(X_scaled, X_scaled, dims(obs_id, predictors)) # Priors alpha pm.Normal(alpha, mu0, sigma1) beta pm.Normal(beta, mu0, sigma1, dimspredictors) sigma pm.HalfNormal(sigma, sigma1) # Linear predictor mu alpha pm.math.dot(X_data, beta) # Tie observed shape to X_data so out-of-sample prediction can resize it y_obs pm.Normal(y_obs, mumu, sigmasigma, observedy, shapeX_data.shape[0], dimsobs_id) # 3. PRIOR PREDICTIVE CHECK # with model: prior_pred pm.sample_prior_predictive(draws1000, random_seed42) az.plot_ppc(prior_pred, groupprior, num_pp_samples100) plt.title(Prior Predictive Check) plt.show() # 4. FIT MODEL # with model: # Quick VI exploration (optional) approx pm.fit(n20000, random_seed42) # Full MCMC inference idata pm.sample( draws2000, tune1000, chains4, target_accept0.9, random_seed42, idata_kwargs{log_likelihood: True} # For model comparison ) # 5. CHECK DIAGNOSTICS # print(az.summary(idata, var_names[alpha, beta, sigma])) summary az.summary(idata) if (summary[r_hat] 1.01).any(): print(WARNING: Some R-hat values 1.01, chains may not have converged) if (summary[ess_bulk] 400).any(): print(WARNING: Some ESS values 400, consider more samples) divergences idata.sample_stats.diverging.sum().item() print(fNumber of divergences: {divergences}) az.plot_trace_dist(idata, var_names[alpha, beta, sigma]) plt.tight_layout() plt.show() # 6. POSTERIOR PREDICTIVE CHECK # with model: pm.sample_posterior_predictive(idata, extend_inferencedataTrue, random_seed42) az.plot_ppc(idata, num_pp_samples100) plt.title(Posterior Predictive Check) plt.show() # 7. ANALYZE RESULTS # az.plot_posterior(idata, var_names[alpha, beta, sigma]) plt.tight_layout() plt.show() az.plot_forest(idata, var_names[beta], combinedTrue) plt.title(Coefficient Estimates) plt.show() # 8. PREDICTIONS FOR NEW DATA # X_new ... # New predictor values X_new_scaled (X_new - X.mean(axis0)) / X.std(axis0) with model: pm.set_data({X_scaled: X_new_scaled}, coords{obs_id: np.arange(len(X_new_scaled))}) post_pred pm.sample_posterior_predictive( idata, var_names[y_obs], predictionsTrue, random_seed42 ) y_pred_mean post_pred.predictions[y_obs].mean(dim[chain, draw]) y_pred_hdi az.hdi(post_pred.predictions, var_names[y_obs]) # 9. SAVE RESULTS # idata.to_netcdf(model_results.nc) # Save for later工作流关键设计点标准化预测变量对连续型预测变量做 (X - mean) / std 标准化让先验具有可解释性例如beta ~ Normal(0,1)表示每变化一个标准差结果变化约 1 个单位同时显著改善采样几何形态。用pm.Data包数据模型中的X_data是可变的 data container后续通过pm.set_data()替换新数据即可做样本外预测y_obs的shape显式绑定到X_data.shape[0]这样预测时新样本的行数可以不同。先验预测检查必须在拟合之前先验不合理会直接污染后验先确认先验能产生合理的数据范围再进入采样。idata_kwargs{log_likelihood: True}采样时即计算逐点对数似然这是后续模型比较LOO/WAIC的前提。PyMC 6 版本注意避免使用已弃用的nuts_sampler_kwargs需要透传采样器专属参数时用显式的采样器关键字字典如nuts{target_accept: 0.9}。仓库模板佐证仓库 assets/linear_regression_template.py 将上述流程落成了可直接运行的完整脚本含合成数据、诊断输出、预测区间打印与结果保存而 assets/hierarchical_model_template.py 则提供了分层模型的完整版本两者均以# TODO标注需要按业务自定义的位置。二、模型构建模式六类常用模型代码骨架workflows.md 整理了六类可复用的模型结构与其姊妹文档 model_patterns.md 相互印证。1. 线性回归连续型结果with pm.Model() as linear_model: # Priors alpha pm.Normal(alpha, mu0, sigma10) beta pm.Normal(beta, mu0, sigma10, shapen_predictors) sigma pm.HalfNormal(sigma, sigma1) # Linear predictor mu alpha pm.math.dot(X, beta) # Likelihood y pm.Normal(y, mumu, sigmasigma, observedy_obs)2. 逻辑回归二分类结果with pm.Model() as logistic_model: alpha pm.Normal(alpha, mu0, sigma10) beta pm.Normal(beta, mu0, sigma10, shapen_predictors) logit_p alpha pm.math.dot(X, beta) # Likelihood y pm.Bernoulli(y, logit_plogit_p, observedy_obs)注意pm.Bernoulli同时支持pp与logit_plogit_p两种接口后者让线性预测器直接作用于 logit 尺度数值上更稳定。3. 分层 / 多水平模型分组数据with pm.Model(coords{group: group_names, obs: np.arange(n_obs)}) as hierarchical_model: # Hyperpriors mu_alpha pm.Normal(mu_alpha, mu0, sigma10) sigma_alpha pm.HalfNormal(sigma_alpha, sigma1) mu_beta pm.Normal(mu_beta, mu0, sigma10) sigma_beta pm.HalfNormal(sigma_beta, sigma1) # Group-level parameters (non-centered) alpha_offset pm.Normal(alpha_offset, mu0, sigma1, dimsgroup) alpha pm.Deterministic(alpha, mu_alpha sigma_alpha * alpha_offset, dimsgroup) beta_offset pm.Normal(beta_offset, mu0, sigma1, dimsgroup) beta pm.Deterministic(beta, mu_beta sigma_beta * beta_offset, dimsgroup) # Observation-level model mu alpha[group_idx] beta[group_idx] * X sigma pm.HalfNormal(sigma, sigma1) y pm.Normal(y, mumu, sigmasigma, observedy_obs, dimsobs)关键纪律分层模型必须使用非中心化参数化non-centered parameterization——先让offset服从标准正态Normal(0,1)再通过alpha mu_alpha sigma_alpha * alpha_offset组合回真实参数。这与 sampling_inference.md 中中心化参数化容易引发散度的警告一致也是避免 NUTS 在漏斗形后验上产生散度的核心手段。4. 泊松回归计数数据with pm.Model() as poisson_model: alpha pm.Normal(alpha, mu0, sigma10) beta pm.Normal(beta, mu0, sigma10, shapen_predictors) # Linear predictor on log scale log_lambda alpha pm.math.dot(X, beta) # Likelihood y pm.Poisson(y, mupm.math.exp(log_lambda), observedy_obs)若计数数据存在过度离散方差大于均值model_patterns.md 建议改用NegativeBinomialSKILL.md 还补充了零膨胀场景下的ZeroInflatedPoisson与兼具零膨胀与过度离散的HurdleNegativeBinomial。5. 时间序列自回归过程with pm.Model() as ar_model: # Innovation standard deviation sigma pm.HalfNormal(sigma, sigma1) # AR coefficients rho pm.Normal(rho, mu0, sigma0.5, shapear_order) # Initial distribution init_dist pm.Normal.dist(mu0, sigmasigma) # AR process y pm.AR(y, rhorho, sigmasigma, init_distinit_dist, observedy_obs)自回归系数先验Normal(0, 0.5)比默认的宽先验更紧是为了让过程保持平稳性。6. 混合模型多组分分布with pm.Model() as mixture_model: # Component weights w pm.Dirichlet(w, anp.ones(n_components)) # Component parameters mu pm.Normal(mu, mu0, sigma10, shapen_components) sigma pm.HalfNormal(sigma, sigma1, shapen_components) # Mixture components [pm.Normal.dist(mumu[i], sigmasigma[i]) for i in range(n_components)] y pm.Mixture(y, ww, comp_distscomponents, observedy_obs)权重w用Dirichlet保证和为 1各组分的分布通过pm.Normal.dist()构造后传入pm.Mixture。三、数据准备最佳实践workflows.md 给出三类高频数据处理模式1. 标准化与系数回变换# Standardize X_mean X.mean(axis0) X_std X.std(axis0) X_scaled (X - X_mean) / X_std # Model with scaled data with pm.Model() as model: beta_scaled pm.Normal(beta_scaled, 0, 1) # ... rest of model ... # Transform back to original scale beta_original beta_scaled / X_std alpha_original alpha - (beta_scaled * X_mean / X_std).sum()标准化后先验Normal(0, 1)对每个系数而言都是弱信息先验展示结果时再通过公式回变换到原始尺度保持解释性。2. 缺失数据处理作为参数插补missing_idx np.isnan(X) X_observed np.where(missing_idx, 0, X) # Placeholder with pm.Model() as model: # Prior for missing values X_missing pm.Normal(X_missing, mu0, sigma1, shapemissing_idx.sum()) # Combine observed and imputed X_complete pm.math.switch(missing_idx.flatten(), X_missing, X_observed.flatten()) # ... rest of model using X_complete ...贝叶斯方法的优势在于缺失值不再被填充一次而是作为待估参数其后验自动传播不确定性。用pm.math.switch在观测值与缺失值之间按掩码切换。3. 中心化X_centered X - X.mean(axis0) y_centered y - y.mean() with pm.Model() as model: alpha pm.Normal(alpha, mu0, sigma1) # Intercept near 0 when centered beta pm.Normal(beta, mu0, sigma1, shapen_predictors) mu alpha pm.math.dot(X_centered, beta) sigma pm.HalfNormal(sigma, sigma1) y_obs pm.Normal(y_obs, mumu, sigmasigma, observedy_centered)结果中心化后截距天然接近 0可以用Normal(0, 1)的紧先验降低参数之间的相关性。四、先验选择指南弱信息先验信息有限时的默认选择# For standardized predictors beta pm.Normal(beta, mu0, sigma1) # For scale parameters sigma pm.HalfNormal(sigma, sigma1) # For probabilities p pm.Beta(p, alpha2, beta2) # Slight preference for middle values信息性先验利用领域知识# Effect size from literature: Cohens d ≈ 0.3 beta pm.Normal(beta, mu0.3, sigma0.1) # Physical constraint: probability between 0.7-0.9 p pm.Beta(p, alpha8, beta2) # Check with prior predictive!信息性先验必须用先验预测检查验证其隐含的数据分布是否与领域认知一致——Beta(8, 2)的均值约 0.8符合概率在 0.7~0.9的物理约束。先验预测检查必做with model: prior_pred pm.sample_prior_predictive(draws1000) # Check if predictions are reasonable print(fPrior predictive range: {prior_pred.prior_predictive[y].min():.2f} to {prior_pred.prior_predictive[y].max():.2f}) print(fObserved range: {y_obs.min():.2f} to {y_obs.max():.2f}) # Visualize az.plot_ppc(prior_pred, groupprior)判断标准先验预测的数据范围是否覆盖观测范围、极值是否合理。若先验生成的模拟数据与观测数据量级完全不符说明先验设定有问题必须调整后重新检查。分布选择快速参考来自 SKILL.md尺度参数σ, τHalfNormal(sigma1)默认、Exponential(lam1)备选、Gamma(alpha2, beta1)信息更强。无界参数Normal(mu0, sigma1)标准化数据、StudentT(nu3, mu0, sigma1)对异常值稳健。正数参数LogNormal(mu0, sigma1)、Gamma(alpha2, beta1)。概率Beta(2, 2)弱信息、Uniform(0,1)无信息谨慎使用。相关矩阵LKJCholeskyCov首选协方差先验、LKJCorr(n, eta2)仅相关先验eta1 为均匀eta1 偏向单位矩阵。计数似然Poisson等离散、NegativeBinomial过度离散、ZeroInflatedPoisson多余零、HurdleNegativeBinomial多余零过度离散。五、模型比较工作流多模型 LOO 比较import arviz as az models {} idatas {} # 每个模型采样时都必须带上 log_likelihood with pm.Model() as models[linear]: # ... define model ... idatas[linear] pm.sample(idata_kwargs{log_likelihood: True}) with pm.Model() as models[interaction]: # ... define model ... idatas[interaction] pm.sample(idata_kwargs{log_likelihood: True}) with pm.Model() as models[hierarchical]: # ... define model ... idatas[hierarchical] pm.sample(idata_kwargs{log_likelihood: True}) # Compare using LOO comparison az.compare(idatas, icloo) print(comparison) # Visualize comparison az.plot_compare(comparison) plt.show() # Check LOO reliability for name, idata in idatas.items(): loo az.loo(idata, pointwiseTrue) high_pareto_k (loo.pareto_k 0.7).sum().item() if high_pareto_k 0: print(fWarning: {name} has {high_pareto_k} observations with high Pareto-k)解读 ELPD 差异model_patterns.md 明确指出ArviZ 1.x 的compare()只基于 PSIS-LOO ELPD 排序报告尺度为 ELPD越高越好最优模型elpd_diff为 0其余为负。解读准则|elpd_diff| 4模型表现相似选择更简单的模型|elpd_diff| 4 但小于 2 倍dse对更优模型只有中等证据|elpd_diff| 4 且超过 2 倍dse对更优模型有强证据。Pareto-k 可靠性阈值k 0.7 时 LOO 可靠k 0.7 时建议改用 WAIC 或 k 折交叉验证。仓库 scripts/model_comparison.py 中compare_models()在ic参数上只接受loo或elpdArviZ 1.x 无ic开关、无偏差尺度WAIC 需通过az.waic()逐模型单独计算check_loo_reliability()则自动按 k0.5好、0.5≤k0.7可接受、0.7≤k1.0差、k≥1.0很差四个区间统计并给出处理建议。模型权重与加权平均weights comparison[weight].values print(Model probabilities:) for name, weight in zip(comparison.index, weights): print(f {name}: {weight:.2%}) # Model averaging (weighted predictions) def weighted_predictions(idatas, weights): preds [] for (name, idata), weight in zip(idatas.items(), weights): group idata.predictions if hasattr(idata, predictions) else idata.posterior_predictive pred group[y_obs].mean(dim[chain, draw]) preds.append(weight * pred) return sum(preds) averaged_pred weighted_predictions(idatas, weights)weight列来自 stacking 权重可视为各模型的后验概率。当模型间 ELPD 差异不显著时直接使用scripts/model_comparison.py中的model_averaging()做加权平均预测比硬选单一模型更稳健该函数在未提供权重时会自动调用compare_models()取 stacking 权重并归一化。当 LOO 因高 Pareto-k 不可靠时同脚本的cross_validation_comparison()提供了 k 折交叉验证的完整实施指南对每个 fold 重拟合、pm.compute_log_likelihood计算留出集对数似然并累加比较。六、诊断与排障综合诊断函数workflows.md 提供了一个完整的diagnose_sampling()函数模板顺序检查四大类指标仓库 scripts/model_diagnostics.py 将同类逻辑工程化为check_diagnostics()与create_diagnostic_report()from scripts.model_diagnostics import check_diagnostics, create_diagnostic_report # 快速检查R-hat、ESS、散度、树深 results check_diagnostics(idata, var_names[alpha, beta, sigma]) # 完整报告trace/rank/自相关/能量/局部ESS图 汇总CSV create_diagnostic_report(idata, var_names[alpha, beta, sigma], output_dirdiagnostics/)源码中的几个关键实现细节值得注意check_diagnostics()调用az.summary(idata, round_tonone)——ArviZ 1.x 默认会为显示格式化汇总返回字符串导致所有数值比较抛TypeErrorround_tonone是必需的见 model_diagnostics.py 注释。散度检查会计算散度占比divergences / total_samples * 100树深检查以默认最大树深 10 为阈值统计命中比例model_diagnostics.py。create_diagnostic_report()通过plot_collection.savefig(...)保存图片因为 ArviZ 1.x 的绘图函数返回PlotCollection而非 pyplot 当前图直接用plt.savefig()会写出空白图model_diagnostics.py。脚本还提供compare_prior_posterior()通过把先验蓝色与后验绿色分布叠加在同一坐标轴上可视化先验被数据更新的程度。常见问题与修复速查表问题解决方案散度Divergences提高target_accept0.95改用非中心化参数化有效样本量低Low ESS增加 draws重参数化以降低相关性R-hat 高运行更长的链检查多模态改进初始化采样慢用 ADVI 初始化重参数化降低模型复杂度后验有偏检查先验预测确认似然设定正确SKILL.md 的补充建议相关性强的回归预测变量可用 QR 分解重参数化分层模型散度优先重参数化而非仅提高target_accept可并行增加cores8, chains8加速。七、命名维度dims的正确用法使用 dims 的收益代码更可读维度名即语义更容易子集与后续分析按坐标名筛选与 xarray 生态深度集成采样结果自动保留维度结构。coords { predictors: [age, income, education], groups: [A, B, C], time: pd.date_range(2020-01-01, periods100, freqD) } with pm.Model(coordscoords) as model: # Use dims instead of shape beta pm.Normal(beta, mu0, sigma1, dimspredictors) alpha pm.Normal(alpha, mu0, sigma1, dimsgroups) y pm.Normal(y, mu0, sigma1, dims[groups, time], observeddata) # After sampling, dimensions are preserved idata pm.sample() # Easy subsetting beta_age idata.posterior[beta].sel(predictorsage) group_A idata.posterior[alpha].sel(groupsA)dims与coords组合后后验数据以 xarrayDataTree组织sel()、mean(dim...)等操作直接按坐标名进行无需记忆索引位置。PyMC 6 / ArviZ 1 在保留.posterior、.posterior_predictive等熟悉分组的前提下全面转向DataTree结构。八、保存与加载结果# Save posterior data tree idata.to_netcdf(results.nc) # Load saved posterior data loaded_idata az.from_netcdf(results.nc) # Save model for later predictions # Only unpickle model files you created and trust; prefer NetCDF for sampled results. import pickle with open(model.pkl, wb) as f: pickle.dump({model: model, idata: idata}, f) # Load model with open(model.pkl, rb) as f: saved pickle.load(f) model saved[model] idata saved[idata]官方推荐以 NetCDF 为主格式保存采样结果跨环境可移植、无需反序列化代码pickle 保存整个模型对象虽便于继续做pm.set_data 预测但存在反序列化安全隐患只能加载自己创建并信任的模型文件。仓库模板还会额外把az.summary输出导出为 CSV如 linear_regression_template.py 中的summary.to_csv(model_summary.csv)便于归档与报告。九、配套资源导航本技能模块在仓库中的完整配套包括参考文档standard_workflow.md八步工作流详解、workflows.md本文核心、常见模式 cookbook、model_patterns.md模型模式与比较解读、sampling_inference.mdNUTS/Metropolis/SMC、ADVI/SVGD 与重参数化技巧、distributions.md按类别组织的分布目录。可执行脚本model_diagnostics.py诊断检查与报告生成、model_comparison.pyLOO 比较、可靠性检查、模型平均、k 折 CV 指南。完整模板linear_regression_template.py线性回归全流程、hierarchical_model_template.py分层模型全流程含组间收缩可视化与新组预测。运行环境SKILL.md 要求 Python 3.12 与 PyMC 6.0.1 兼容依赖推荐uv pip install pymc[nutpie]6.0.1nutpieextra 提供更快的 Rust/Numba NUTS 实现如需 NumPyro 或 BlackJAX 采样器需在锁文件中单独固定 JAX 兼容依赖。十、结语把工作流变成习惯贯穿本指南始终的核心理念是 SKILL.md 中反复强调的先检查后采样与逐步迭代先验预测检查 → 拟合 → 诊断 → 后验预测检查任何一步发现问题都应回到模型规范进行调整而不是提高target_accept后寄希望于运气。将 workflows.md 的九步模板、六类模型骨架与仓库提供的脚本、模板组合使用即可在科学计算项目中建立一套规范、可复现、可比较的贝叶斯分析流水线。【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
PyMC 贝叶斯建模工作流与常用模式实践指南：从先验检查到模型比较的完整技术栈
PyMC 贝叶斯建模工作流与常用模式实践指南从先验检查到模型比较的完整技术栈【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills本指南基于 scientific-agent-skills 仓库中 PyMC 技能模块的核心参考文档 workflows.md 编写系统讲解使用 PyMC6.x构建、验证与分析贝叶斯模型的完整工作流涵盖标准九步工作流模板、线性/逻辑/分层/泊松/时间序列/混合六类模型构建模式、数据准备与先验选择的最佳实践以及基于 PSIS-LOO 的模型比较与采样诊断排障方案。读者学完后可直接复刻一套可运行的端到端贝叶斯分析流程并结合仓库提供的 model_diagnostics.py、model_comparison.py 脚本和两份可定制模板完成实际项目落地。一、标准贝叶斯工作流九步端到端模板PyMC 技能的核心纪律是先检查后采样Never sample first and check later。SKILL.md 将其归纳为八个阶段而 workflows.md 给出了每一步的完整可运行代码。以下模板综合两者可作为任何贝叶斯回归分析的基本骨架import pymc as pm import arviz as az import numpy as np import matplotlib.pyplot as plt # 1. PREPARE DATA # X ... # Predictor variables y ... # Observed outcomes # Standardize predictors for better sampling X_scaled (X - X.mean(axis0)) / X.std(axis0) # 2. BUILD MODEL # coords { predictors: [var1, var2, var3], obs_id: np.arange(len(y)) } with pm.Model(coordscoords) as model: X_data pm.Data(X_scaled, X_scaled, dims(obs_id, predictors)) # Priors alpha pm.Normal(alpha, mu0, sigma1) beta pm.Normal(beta, mu0, sigma1, dimspredictors) sigma pm.HalfNormal(sigma, sigma1) # Linear predictor mu alpha pm.math.dot(X_data, beta) # Tie observed shape to X_data so out-of-sample prediction can resize it y_obs pm.Normal(y_obs, mumu, sigmasigma, observedy, shapeX_data.shape[0], dimsobs_id) # 3. PRIOR PREDICTIVE CHECK # with model: prior_pred pm.sample_prior_predictive(draws1000, random_seed42) az.plot_ppc(prior_pred, groupprior, num_pp_samples100) plt.title(Prior Predictive Check) plt.show() # 4. FIT MODEL # with model: # Quick VI exploration (optional) approx pm.fit(n20000, random_seed42) # Full MCMC inference idata pm.sample( draws2000, tune1000, chains4, target_accept0.9, random_seed42, idata_kwargs{log_likelihood: True} # For model comparison ) # 5. CHECK DIAGNOSTICS # print(az.summary(idata, var_names[alpha, beta, sigma])) summary az.summary(idata) if (summary[r_hat] 1.01).any(): print(WARNING: Some R-hat values 1.01, chains may not have converged) if (summary[ess_bulk] 400).any(): print(WARNING: Some ESS values 400, consider more samples) divergences idata.sample_stats.diverging.sum().item() print(fNumber of divergences: {divergences}) az.plot_trace_dist(idata, var_names[alpha, beta, sigma]) plt.tight_layout() plt.show() # 6. POSTERIOR PREDICTIVE CHECK # with model: pm.sample_posterior_predictive(idata, extend_inferencedataTrue, random_seed42) az.plot_ppc(idata, num_pp_samples100) plt.title(Posterior Predictive Check) plt.show() # 7. ANALYZE RESULTS # az.plot_posterior(idata, var_names[alpha, beta, sigma]) plt.tight_layout() plt.show() az.plot_forest(idata, var_names[beta], combinedTrue) plt.title(Coefficient Estimates) plt.show() # 8. PREDICTIONS FOR NEW DATA # X_new ... # New predictor values X_new_scaled (X_new - X.mean(axis0)) / X.std(axis0) with model: pm.set_data({X_scaled: X_new_scaled}, coords{obs_id: np.arange(len(X_new_scaled))}) post_pred pm.sample_posterior_predictive( idata, var_names[y_obs], predictionsTrue, random_seed42 ) y_pred_mean post_pred.predictions[y_obs].mean(dim[chain, draw]) y_pred_hdi az.hdi(post_pred.predictions, var_names[y_obs]) # 9. SAVE RESULTS # idata.to_netcdf(model_results.nc) # Save for later工作流关键设计点标准化预测变量对连续型预测变量做 (X - mean) / std 标准化让先验具有可解释性例如beta ~ Normal(0,1)表示每变化一个标准差结果变化约 1 个单位同时显著改善采样几何形态。用pm.Data包数据模型中的X_data是可变的 data container后续通过pm.set_data()替换新数据即可做样本外预测y_obs的shape显式绑定到X_data.shape[0]这样预测时新样本的行数可以不同。先验预测检查必须在拟合之前先验不合理会直接污染后验先确认先验能产生合理的数据范围再进入采样。idata_kwargs{log_likelihood: True}采样时即计算逐点对数似然这是后续模型比较LOO/WAIC的前提。PyMC 6 版本注意避免使用已弃用的nuts_sampler_kwargs需要透传采样器专属参数时用显式的采样器关键字字典如nuts{target_accept: 0.9}。仓库模板佐证仓库 assets/linear_regression_template.py 将上述流程落成了可直接运行的完整脚本含合成数据、诊断输出、预测区间打印与结果保存而 assets/hierarchical_model_template.py 则提供了分层模型的完整版本两者均以# TODO标注需要按业务自定义的位置。二、模型构建模式六类常用模型代码骨架workflows.md 整理了六类可复用的模型结构与其姊妹文档 model_patterns.md 相互印证。1. 线性回归连续型结果with pm.Model() as linear_model: # Priors alpha pm.Normal(alpha, mu0, sigma10) beta pm.Normal(beta, mu0, sigma10, shapen_predictors) sigma pm.HalfNormal(sigma, sigma1) # Linear predictor mu alpha pm.math.dot(X, beta) # Likelihood y pm.Normal(y, mumu, sigmasigma, observedy_obs)2. 逻辑回归二分类结果with pm.Model() as logistic_model: alpha pm.Normal(alpha, mu0, sigma10) beta pm.Normal(beta, mu0, sigma10, shapen_predictors) logit_p alpha pm.math.dot(X, beta) # Likelihood y pm.Bernoulli(y, logit_plogit_p, observedy_obs)注意pm.Bernoulli同时支持pp与logit_plogit_p两种接口后者让线性预测器直接作用于 logit 尺度数值上更稳定。3. 分层 / 多水平模型分组数据with pm.Model(coords{group: group_names, obs: np.arange(n_obs)}) as hierarchical_model: # Hyperpriors mu_alpha pm.Normal(mu_alpha, mu0, sigma10) sigma_alpha pm.HalfNormal(sigma_alpha, sigma1) mu_beta pm.Normal(mu_beta, mu0, sigma10) sigma_beta pm.HalfNormal(sigma_beta, sigma1) # Group-level parameters (non-centered) alpha_offset pm.Normal(alpha_offset, mu0, sigma1, dimsgroup) alpha pm.Deterministic(alpha, mu_alpha sigma_alpha * alpha_offset, dimsgroup) beta_offset pm.Normal(beta_offset, mu0, sigma1, dimsgroup) beta pm.Deterministic(beta, mu_beta sigma_beta * beta_offset, dimsgroup) # Observation-level model mu alpha[group_idx] beta[group_idx] * X sigma pm.HalfNormal(sigma, sigma1) y pm.Normal(y, mumu, sigmasigma, observedy_obs, dimsobs)关键纪律分层模型必须使用非中心化参数化non-centered parameterization——先让offset服从标准正态Normal(0,1)再通过alpha mu_alpha sigma_alpha * alpha_offset组合回真实参数。这与 sampling_inference.md 中中心化参数化容易引发散度的警告一致也是避免 NUTS 在漏斗形后验上产生散度的核心手段。4. 泊松回归计数数据with pm.Model() as poisson_model: alpha pm.Normal(alpha, mu0, sigma10) beta pm.Normal(beta, mu0, sigma10, shapen_predictors) # Linear predictor on log scale log_lambda alpha pm.math.dot(X, beta) # Likelihood y pm.Poisson(y, mupm.math.exp(log_lambda), observedy_obs)若计数数据存在过度离散方差大于均值model_patterns.md 建议改用NegativeBinomialSKILL.md 还补充了零膨胀场景下的ZeroInflatedPoisson与兼具零膨胀与过度离散的HurdleNegativeBinomial。5. 时间序列自回归过程with pm.Model() as ar_model: # Innovation standard deviation sigma pm.HalfNormal(sigma, sigma1) # AR coefficients rho pm.Normal(rho, mu0, sigma0.5, shapear_order) # Initial distribution init_dist pm.Normal.dist(mu0, sigmasigma) # AR process y pm.AR(y, rhorho, sigmasigma, init_distinit_dist, observedy_obs)自回归系数先验Normal(0, 0.5)比默认的宽先验更紧是为了让过程保持平稳性。6. 混合模型多组分分布with pm.Model() as mixture_model: # Component weights w pm.Dirichlet(w, anp.ones(n_components)) # Component parameters mu pm.Normal(mu, mu0, sigma10, shapen_components) sigma pm.HalfNormal(sigma, sigma1, shapen_components) # Mixture components [pm.Normal.dist(mumu[i], sigmasigma[i]) for i in range(n_components)] y pm.Mixture(y, ww, comp_distscomponents, observedy_obs)权重w用Dirichlet保证和为 1各组分的分布通过pm.Normal.dist()构造后传入pm.Mixture。三、数据准备最佳实践workflows.md 给出三类高频数据处理模式1. 标准化与系数回变换# Standardize X_mean X.mean(axis0) X_std X.std(axis0) X_scaled (X - X_mean) / X_std # Model with scaled data with pm.Model() as model: beta_scaled pm.Normal(beta_scaled, 0, 1) # ... rest of model ... # Transform back to original scale beta_original beta_scaled / X_std alpha_original alpha - (beta_scaled * X_mean / X_std).sum()标准化后先验Normal(0, 1)对每个系数而言都是弱信息先验展示结果时再通过公式回变换到原始尺度保持解释性。2. 缺失数据处理作为参数插补missing_idx np.isnan(X) X_observed np.where(missing_idx, 0, X) # Placeholder with pm.Model() as model: # Prior for missing values X_missing pm.Normal(X_missing, mu0, sigma1, shapemissing_idx.sum()) # Combine observed and imputed X_complete pm.math.switch(missing_idx.flatten(), X_missing, X_observed.flatten()) # ... rest of model using X_complete ...贝叶斯方法的优势在于缺失值不再被填充一次而是作为待估参数其后验自动传播不确定性。用pm.math.switch在观测值与缺失值之间按掩码切换。3. 中心化X_centered X - X.mean(axis0) y_centered y - y.mean() with pm.Model() as model: alpha pm.Normal(alpha, mu0, sigma1) # Intercept near 0 when centered beta pm.Normal(beta, mu0, sigma1, shapen_predictors) mu alpha pm.math.dot(X_centered, beta) sigma pm.HalfNormal(sigma, sigma1) y_obs pm.Normal(y_obs, mumu, sigmasigma, observedy_centered)结果中心化后截距天然接近 0可以用Normal(0, 1)的紧先验降低参数之间的相关性。四、先验选择指南弱信息先验信息有限时的默认选择# For standardized predictors beta pm.Normal(beta, mu0, sigma1) # For scale parameters sigma pm.HalfNormal(sigma, sigma1) # For probabilities p pm.Beta(p, alpha2, beta2) # Slight preference for middle values信息性先验利用领域知识# Effect size from literature: Cohens d ≈ 0.3 beta pm.Normal(beta, mu0.3, sigma0.1) # Physical constraint: probability between 0.7-0.9 p pm.Beta(p, alpha8, beta2) # Check with prior predictive!信息性先验必须用先验预测检查验证其隐含的数据分布是否与领域认知一致——Beta(8, 2)的均值约 0.8符合概率在 0.7~0.9的物理约束。先验预测检查必做with model: prior_pred pm.sample_prior_predictive(draws1000) # Check if predictions are reasonable print(fPrior predictive range: {prior_pred.prior_predictive[y].min():.2f} to {prior_pred.prior_predictive[y].max():.2f}) print(fObserved range: {y_obs.min():.2f} to {y_obs.max():.2f}) # Visualize az.plot_ppc(prior_pred, groupprior)判断标准先验预测的数据范围是否覆盖观测范围、极值是否合理。若先验生成的模拟数据与观测数据量级完全不符说明先验设定有问题必须调整后重新检查。分布选择快速参考来自 SKILL.md尺度参数σ, τHalfNormal(sigma1)默认、Exponential(lam1)备选、Gamma(alpha2, beta1)信息更强。无界参数Normal(mu0, sigma1)标准化数据、StudentT(nu3, mu0, sigma1)对异常值稳健。正数参数LogNormal(mu0, sigma1)、Gamma(alpha2, beta1)。概率Beta(2, 2)弱信息、Uniform(0,1)无信息谨慎使用。相关矩阵LKJCholeskyCov首选协方差先验、LKJCorr(n, eta2)仅相关先验eta1 为均匀eta1 偏向单位矩阵。计数似然Poisson等离散、NegativeBinomial过度离散、ZeroInflatedPoisson多余零、HurdleNegativeBinomial多余零过度离散。五、模型比较工作流多模型 LOO 比较import arviz as az models {} idatas {} # 每个模型采样时都必须带上 log_likelihood with pm.Model() as models[linear]: # ... define model ... idatas[linear] pm.sample(idata_kwargs{log_likelihood: True}) with pm.Model() as models[interaction]: # ... define model ... idatas[interaction] pm.sample(idata_kwargs{log_likelihood: True}) with pm.Model() as models[hierarchical]: # ... define model ... idatas[hierarchical] pm.sample(idata_kwargs{log_likelihood: True}) # Compare using LOO comparison az.compare(idatas, icloo) print(comparison) # Visualize comparison az.plot_compare(comparison) plt.show() # Check LOO reliability for name, idata in idatas.items(): loo az.loo(idata, pointwiseTrue) high_pareto_k (loo.pareto_k 0.7).sum().item() if high_pareto_k 0: print(fWarning: {name} has {high_pareto_k} observations with high Pareto-k)解读 ELPD 差异model_patterns.md 明确指出ArviZ 1.x 的compare()只基于 PSIS-LOO ELPD 排序报告尺度为 ELPD越高越好最优模型elpd_diff为 0其余为负。解读准则|elpd_diff| 4模型表现相似选择更简单的模型|elpd_diff| 4 但小于 2 倍dse对更优模型只有中等证据|elpd_diff| 4 且超过 2 倍dse对更优模型有强证据。Pareto-k 可靠性阈值k 0.7 时 LOO 可靠k 0.7 时建议改用 WAIC 或 k 折交叉验证。仓库 scripts/model_comparison.py 中compare_models()在ic参数上只接受loo或elpdArviZ 1.x 无ic开关、无偏差尺度WAIC 需通过az.waic()逐模型单独计算check_loo_reliability()则自动按 k0.5好、0.5≤k0.7可接受、0.7≤k1.0差、k≥1.0很差四个区间统计并给出处理建议。模型权重与加权平均weights comparison[weight].values print(Model probabilities:) for name, weight in zip(comparison.index, weights): print(f {name}: {weight:.2%}) # Model averaging (weighted predictions) def weighted_predictions(idatas, weights): preds [] for (name, idata), weight in zip(idatas.items(), weights): group idata.predictions if hasattr(idata, predictions) else idata.posterior_predictive pred group[y_obs].mean(dim[chain, draw]) preds.append(weight * pred) return sum(preds) averaged_pred weighted_predictions(idatas, weights)weight列来自 stacking 权重可视为各模型的后验概率。当模型间 ELPD 差异不显著时直接使用scripts/model_comparison.py中的model_averaging()做加权平均预测比硬选单一模型更稳健该函数在未提供权重时会自动调用compare_models()取 stacking 权重并归一化。当 LOO 因高 Pareto-k 不可靠时同脚本的cross_validation_comparison()提供了 k 折交叉验证的完整实施指南对每个 fold 重拟合、pm.compute_log_likelihood计算留出集对数似然并累加比较。六、诊断与排障综合诊断函数workflows.md 提供了一个完整的diagnose_sampling()函数模板顺序检查四大类指标仓库 scripts/model_diagnostics.py 将同类逻辑工程化为check_diagnostics()与create_diagnostic_report()from scripts.model_diagnostics import check_diagnostics, create_diagnostic_report # 快速检查R-hat、ESS、散度、树深 results check_diagnostics(idata, var_names[alpha, beta, sigma]) # 完整报告trace/rank/自相关/能量/局部ESS图 汇总CSV create_diagnostic_report(idata, var_names[alpha, beta, sigma], output_dirdiagnostics/)源码中的几个关键实现细节值得注意check_diagnostics()调用az.summary(idata, round_tonone)——ArviZ 1.x 默认会为显示格式化汇总返回字符串导致所有数值比较抛TypeErrorround_tonone是必需的见 model_diagnostics.py 注释。散度检查会计算散度占比divergences / total_samples * 100树深检查以默认最大树深 10 为阈值统计命中比例model_diagnostics.py。create_diagnostic_report()通过plot_collection.savefig(...)保存图片因为 ArviZ 1.x 的绘图函数返回PlotCollection而非 pyplot 当前图直接用plt.savefig()会写出空白图model_diagnostics.py。脚本还提供compare_prior_posterior()通过把先验蓝色与后验绿色分布叠加在同一坐标轴上可视化先验被数据更新的程度。常见问题与修复速查表问题解决方案散度Divergences提高target_accept0.95改用非中心化参数化有效样本量低Low ESS增加 draws重参数化以降低相关性R-hat 高运行更长的链检查多模态改进初始化采样慢用 ADVI 初始化重参数化降低模型复杂度后验有偏检查先验预测确认似然设定正确SKILL.md 的补充建议相关性强的回归预测变量可用 QR 分解重参数化分层模型散度优先重参数化而非仅提高target_accept可并行增加cores8, chains8加速。七、命名维度dims的正确用法使用 dims 的收益代码更可读维度名即语义更容易子集与后续分析按坐标名筛选与 xarray 生态深度集成采样结果自动保留维度结构。coords { predictors: [age, income, education], groups: [A, B, C], time: pd.date_range(2020-01-01, periods100, freqD) } with pm.Model(coordscoords) as model: # Use dims instead of shape beta pm.Normal(beta, mu0, sigma1, dimspredictors) alpha pm.Normal(alpha, mu0, sigma1, dimsgroups) y pm.Normal(y, mu0, sigma1, dims[groups, time], observeddata) # After sampling, dimensions are preserved idata pm.sample() # Easy subsetting beta_age idata.posterior[beta].sel(predictorsage) group_A idata.posterior[alpha].sel(groupsA)dims与coords组合后后验数据以 xarrayDataTree组织sel()、mean(dim...)等操作直接按坐标名进行无需记忆索引位置。PyMC 6 / ArviZ 1 在保留.posterior、.posterior_predictive等熟悉分组的前提下全面转向DataTree结构。八、保存与加载结果# Save posterior data tree idata.to_netcdf(results.nc) # Load saved posterior data loaded_idata az.from_netcdf(results.nc) # Save model for later predictions # Only unpickle model files you created and trust; prefer NetCDF for sampled results. import pickle with open(model.pkl, wb) as f: pickle.dump({model: model, idata: idata}, f) # Load model with open(model.pkl, rb) as f: saved pickle.load(f) model saved[model] idata saved[idata]官方推荐以 NetCDF 为主格式保存采样结果跨环境可移植、无需反序列化代码pickle 保存整个模型对象虽便于继续做pm.set_data 预测但存在反序列化安全隐患只能加载自己创建并信任的模型文件。仓库模板还会额外把az.summary输出导出为 CSV如 linear_regression_template.py 中的summary.to_csv(model_summary.csv)便于归档与报告。九、配套资源导航本技能模块在仓库中的完整配套包括参考文档standard_workflow.md八步工作流详解、workflows.md本文核心、常见模式 cookbook、model_patterns.md模型模式与比较解读、sampling_inference.mdNUTS/Metropolis/SMC、ADVI/SVGD 与重参数化技巧、distributions.md按类别组织的分布目录。可执行脚本model_diagnostics.py诊断检查与报告生成、model_comparison.pyLOO 比较、可靠性检查、模型平均、k 折 CV 指南。完整模板linear_regression_template.py线性回归全流程、hierarchical_model_template.py分层模型全流程含组间收缩可视化与新组预测。运行环境SKILL.md 要求 Python 3.12 与 PyMC 6.0.1 兼容依赖推荐uv pip install pymc[nutpie]6.0.1nutpieextra 提供更快的 Rust/Numba NUTS 实现如需 NumPyro 或 BlackJAX 采样器需在锁文件中单独固定 JAX 兼容依赖。十、结语把工作流变成习惯贯穿本指南始终的核心理念是 SKILL.md 中反复强调的先检查后采样与逐步迭代先验预测检查 → 拟合 → 诊断 → 后验预测检查任何一步发现问题都应回到模型规范进行调整而不是提高target_accept后寄希望于运气。将 workflows.md 的九步模板、六类模型骨架与仓库提供的脚本、模板组合使用即可在科学计算项目中建立一套规范、可复现、可比较的贝叶斯分析流水线。【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
