简介概率图模型比如贝叶斯网络是处理不确定性和因果推断的重要工具其核心环节是结构学习。结构学习解决如何从数据中自动发现变量依赖关系的问题常见的有基于条件独立性检验的PC算法以及将问题转化为连续优化求解的Notears算法。PC算法逐层剪枝、逻辑直观Notears利用DAG约束的连续化处理可直接用梯度下降求解。两方法在计算效率和结构精度上各有取舍直接影响后续MLE参数学习和模型推理效果。在实际工业项目中从数据清洗、离散化到变量筛选再到结构学习、参数学习与测试集评估如AUC、logp需要一套完整的建模流程。本文基于抑郁症调查数据对比这两种结构学习算法在真实场景中的表现为贝叶斯网络建模选型提供实践参考。1. 结构学习算法之争从Notears与PC的对比看贝叶斯网络落地选型贝叶斯网络的结构学习一直是因果推断和概率图模型落地中最耗时的环节。PC算法基于条件独立性检验逐层剪枝候选边逻辑直观但面对高维数据时容易因检验次数过多而效率低下Notears则将结构学习转化为连续优化问题把离散的DAG约束改写为光滑的等式约束可以直接用梯度下降求解。这两条路线一个偏统计检验、一个偏数值优化在抑郁症相关的调查数据上表现差异显著。对于需要处理几十个变量、又要做预测评估的工程师来说选错算法可能意味着几小时的等待和精度损失。这篇文章围绕一个实际项目完整走一遍从数据预处理到结构学习、参数学习、模型评估的流程并给出两个算法的时间与性能对比适合正在做贝叶斯网络建模、但不想只停留在理论上的从业者参考。2. 数据预处理与数据划分让结构学习算法拿到干净输入2.1 数据清洗与类型转换的常见做法贝叶斯网络的结构学习对数据质量非常敏感尤其是PC算法它依赖条件独立性检验的p值来判断边的去留而p值计算的前提是数据分布符合算法的假设。pgmpy内置的PC实现基于卡方检验或G检验要求输入变量为离散类别型Notears虽然可以处理连续数据但在混合类型数据上也需要统一编码。因此第一步是先做数据清洗处理缺失值、剔除常数列、统一变量类型。import pandas as pd import numpy as np df pd.read_csv(depression_data.csv) print(原始数据形状:, df.shape) print(缺失值统计:\n, df.isnull().sum()) # 删除缺失比例超过40%的列 miss_ratio df.isnull().mean() drop_cols miss_ratio[miss_ratio 0.4].index.tolist() df.drop(columnsdrop_cols, inplaceTrue) # 用众数填补离散变量的缺失值 cat_cols df.select_dtypes(include[object]).columns for col in cat_cols: df[col] df[col].fillna(df[col].mode()[0]) # 数值列用中位数填补对离群值更稳健 num_cols df.select_dtypes(include[np.number]).columns for col in num_cols: df[col] df[col].fillna(df[col].median())这里的处理逻辑是分层补齐先删高缺失列再用众数/中位数分别填补。之所以离散变量用众数、连续变量用中位数是因为贝叶斯网络参数学习阶段需要统计频次众数不会引入额外分布偏移中位数对极端值不敏感避免个别离群样本把数据分布拉偏。缺失值处理完后还需要把连续变量离散化。2.2 连续变量离散化与数据划分Notears的官方源码设计目标是连续数据场景但本项目最终要使用MLE参数学习和pgmpy的推理接口这些组件都要求离散输入。建议做一个对照实验一组数据分箱后同时跑PC和Notears保证算法在相同数据尺度下比较。常见的分箱方法有等宽分箱、等频分箱和基于决策树的熵分箱等频分箱在这个场景下更合适因为贝叶斯网络的CPD表格需要每个取值组合有足够的样本支撑。from sklearn.model_selection import train_test_split # 等频分箱分为4个区间用qcut处理偏态分布 for col in num_cols: try: df[col _bin] pd.qcut(df[col], q4, labels[0, 1, 2, 3], duplicatesdrop) df.drop(columns[col], inplaceTrue) except ValueError: # 取值种类过少时直接保留原值 df[col] df[col].astype(category).cat.codes # 所有列转成字符串类别型pgmpy要求显式声明状态 for col in df.columns: df[col] df[col].astype(str) train_data, test_data train_test_split(df, test_size0.3, random_state42, stratifydf[Depression]) print(训练集样本数:, len(train_data), 测试集样本数:, len(test_data))qcut在这里做的是等频分箱每个箱子样本数接近避免等宽分箱下数据集中在某个区间导致CPD中的某些条件概率为零。stratify参数按Depression变量分层抽样保证训练集和测试集中正负样本比例一致后续计算AUC-ROC时不会因为样本不均衡而虚高。注意最后把所有列转成字符串这是pgmpy中BayesianNetwork对象的硬性约束。2.3 变量筛选与共线性处理数据预处理阶段还有一个容易被忽略的步骤变量筛选。如果数据中包含ID列、时间戳列或与目标变量存在机械相关性的字段结构学习算法会把它们识别为强相关节点产生误导性的有向边。另外高度共线的变量会让PC算法的条件独立性检验频繁出现数值不稳定。# 去除低方差列 from sklearn.feature_selection import VarianceThreshold selector VarianceThreshold(threshold0.05) X df.drop(columns[Depression]) selector.fit(X) keep_cols X.columns[selector.get_support()].tolist() print(保留变量:, keep_cols)低方差过滤的阈值设为0.05意味着某列如果95%以上的样本取值相同就不保留。这类列对结构学习没有贡献还可能让检验统计量的自由度计算出错。经过这一步数据集中的特征会收敛在10~20个真正有区分度的变量上PC算法的检验次数也随之减少效率能提升不少。3. PC与Notears结构学习两种算法在同一数据上的实现3.1 PC算法在pgmpy中的调用与参数说明PC算法的核心逻辑是从完全无向图出发逐层增加条件集大小检验两节点是否在给定条件集下独立不独立则保留边最后通过v-structure定向和Meek规则确定边的方向。pgmpy库封装了这一流程使用门槛低但参数选择直接影响结果质量。from pgmpy.estimators import PC from pgmpy.base import DAG import time start_pc time.time() pc PC(datatrain_data) dag_pc pc.estimate( variantstable, ci_testchi_square, significance_level0.05, max_cond_vars5 ) end_pc time.time() print(PC算法耗时:, round(end_pc - start_pc, 4), 秒) print(PC学习到的边数:, len(dag_pc.edges()))variantstable是PC的稳定版本它改变了冲突边的删除顺序结果不再受变量输入顺序影响这在工程上非常重要。ci_testchi_square适用于离散数据max_cond_vars5限制了条件集的维度防止高维条件下的检验因样本稀疏而失效。如果数据量较小这个值建议设为3样本量上万时5也比较保守可以按变量总数取log。耗时的打印是必要的因为Notears的时间对比是项目的核心输出之一。3.2 Notears算法的连续优化求解原理与实现Notears的全称是Non-combinatorial Optimization via Trace Exponential and Augmented lagRangian for Structure learning它的关键思想是把DAG约束写成h(W) tr(e^{W \circ W}) - d 0其中W是加权邻接矩阵d是节点数。这个等式约束使得DAG问题可以被标准的增广拉格朗日方法求解梯度下降每步更新W直到收敛。相比PC的逐对检验Notears的复杂度与变量数直接相关在密集图上通常更快。import numpy as np import notears # 将训练数据转为数值矩阵注意Notears要求连续值输入 X_train train_data.apply(pd.to_numeric, errorscoerce).values.astype(np.float64) X_train np.nan_to_num(X_train, nan0.0) start_notears time.time() # 使用默认lambda参数lambda越大图越稀疏 W_est notears.linear_model(X_train, lambda10.05, lossl2) end_notears time.time() print(Notears算法耗时:, round(end_notears - start_notears, 4), 秒) # 根据阈值剪掉弱边得到DAG W_thresh np.where(np.abs(W_est) 0.3, W_est, 0) print(Notears稀疏化后非零边数:, np.count_nonzero(W_thresh))这里的lambda1是L1正则系数控制稀疏程度值越大得到的边越少反之保留的候选边越多。阈值0.3用于把优化结果中接近零的权重清零因为增广拉格朗日法收敛后会有数值噪声不会精确等于零。lossl2表示使用最小二乘损失适用连续数据。需要说明的是Notears原始实现针对连续变量设计本项目在分箱离散化后仍然可以用但更严谨的做法是在离散化之前对连续特征单独跑一次Notears作为对照对比两种数据形态下学习到的结构差异。3.3 DAG合法性检查与问题定位结构学习完成后不能直接拿去参数学习必须做合法性校验。两个算法都可能产出带环结构尤其当数据噪声较大时PC的定向规则可能保持部分边为无向状态而Notears在阈值处理后也可能出现残余的有向环。from pgmpy.base import DAG # 检查PC结果是否为有向无环图 print(PC结果有效DAG:, dag_pc.is_dag()) # 将Notears的邻接矩阵转为DAG对象 nodes list(train_data.columns) edges_notears [] for i in range(len(nodes)): for j in range(len(nodes)): if W_thresh[i, j] ! 0 and i ! j: edges_notears.append((nodes[i], nodes[j])) dag_notears DAG() dag_notears.add_nodes_from(nodes) try: dag_notears.add_edges_from(edges_notears) print(Notears结果有效DAG:, dag_notears.is_dag()) except ValueError as e: print(Notears结果存在环需要处理:, e)如果Notears的结果有环常见做法是增加lambda1或提高阈值让图变得更稀疏。环的本质是反馈回路通常来自权重接近阈值边界的那几条边把它们剪掉即可。PC结果如果有无向边pgmpy的estimate在某些情况下会保留部分无向边这时需要手动检查变量间的实际语义关系来定向或者直接调用pdag_to_dag这类工具做进一步处理。4. MLE参数学习与预测评估从DAG到可推理的贝叶斯网络4.1 MLE参数学习方法与CPD生成结构学习得到DAG后参数学习要解决的问题是给定图结构每个节点的条件概率分布应该是什么。最大似然估计在这里就是在数据上统计每个父节点取值组合下子节点取值的频次归一化后得到条件概率表。逻辑上和朴素贝叶斯的参数估计一致区别是这里的父节点集合来自前面学到的结构而不是假设所有特征都独立。from pgmpy.models import BayesianNetwork from pgmpy.estimators import MaximumLikelihoodEstimator # 用PC学到的结构构建模型 bayes_pc BayesianNetwork(dag_pc.edges()) bayes_pc.fit(train_data, estimatorMaximumLikelihoodEstimator) # 用Notears学到的结构构建模型 bayes_notears BayesianNetwork(dag_notears.edges()) bayes_notears.fit(train_data, estimatorMaximumLikelihoodEstimator) # 查看Depression节点的CPD print(bayes_pc.get_cpds(Depression)) print(bayes_notears.get_cpds(Depression))fit会对每个节点单独计算CPD写入模型对象。值得注意的一个坑是如果前面赋给BayesianNetwork的边包含未在训练数据中出现的节点组合fit会直接报错所以需要先检查DAG的节点集合是否和数据列完全一致。另外当某个父变量组合在训练集中没有样本时MLE会得到概率为零的条目这在后续推理中会导致预测结果出现0概率需要用拉普拉斯平滑处理。# 手工检查CPD中是否存在零概率条目 cpd bayes_pc.get_cpds(Depression) zero_cells (cpd.values 0).sum() print(PC模型CPD零概率单元格数:, zero_cells) if zero_cells 0: # 改用贝叶斯估计加入伪计数做平滑 from pgmpy.estimators import BayesianEstimator bayes_pc BayesianNetwork(dag_pc.edges()) bayes_pc.fit(train_data, estimatorBayesianEstimator, prior_typeBDeu, equivalent_sample_size10)BDeu先验等价于在统计频次时额外注入等效样本等效样本数设10是一个常见的保守选择。样本量较大时这个先验的影响可以忽略但它能保证CPD中没有绝对的零概率让后面的推理计算稳定。4.2 测试集评估借助VariableElimination进行预测模型评估阶段需要用测试集数据结合学到的网络推理Depression变量的取值。这里用变量消除推理引擎对每个测试样本给定其他变量的观测值计算出Depression取各状态的后验概率取最大概率作为预测类别。from pgmpy.inference import VariableElimination infer_pc VariableElimination(bayes_pc) infer_notears VariableElimination(bayes_notears) def predict_bn(infer, test_df, targetDepression): preds [] for _, row in test_df.iterrows(): evidence {col: str(row[col]) for col in test_df.columns if col ! target} try: result infer.query([target], evidenceevidence) # result.values保存各状态的概率 probs result.values pred result.state_names[target][int(np.argmax(probs))] preds.append(pred) except Exception: preds.append(test_df[target].mode()[0]) return preds y_true test_data[Depression].values y_pred_pc predict_bn(infer_pc, test_data) y_pred_notears predict_bn(infer_notears, test_data)这里有个工程上的细节query方法中evidence的键是变量名值必须是字符串因为之前数据都转成了str类型。异常捕获用于兜底推理过程中可能出现的概率计算不收敛问题此时用众数预测作为回退。这种逐行循环的方式在几百条测试集上可以接受如果测试集上万条建议分块传入并缓存已经算过的evidence组合。4.3 混淆矩阵、AUC-ROC与logp对比拿到预测结果后sklearn提供了完整的评估工具。逻辑回归等模型返回连续的预测分数而贝叶斯网络推理返回的是后验概率分布两者格式略有差异。这里不仅计算标准分类指标还要额外加上logp在模型下观测到整个测试集的联合对数似然。logp高说明网络结构的拟合度高预测稳定性更好。from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score from sklearn.metrics import confusion_matrix, roc_curve, auc import matplotlib.pyplot as plt def evaluate_model(y_true, y_pred, model_name): acc accuracy_score(y_true, y_pred) precision precision_score(y_true, y_pred, pos_label1) recall recall_score(y_true, y_pred, pos_label1) f1 f1_score(y_true, y_pred, pos_label1) cm confusion_matrix(y_true, y_pred) print(f {model_name} ) print(f准确率: {acc:.4f}) print(f精确率: {precision:.4f}) print(f召回率: {recall:.4f}) print(fF1-score: {f1:.4f}) print(混淆矩阵:\n, cm) return acc, precision, recall, f1, cm # 补充logp的简单近似的计算 def calculate_logp(model, test_df): import numpy as np from pgmpy.inference import VariableElimination infer VariableElimination(model) logp_sum 0 for _, row in test_df.iterrows(): evidence {col: str(row[col]) for col in test_df.columns if col ! Depression} q infer.query([Depression], evidenceevidence) prob q.values[list(q.state_names[Depression]).index(str(row[Depression]))] logp_sum np.log(prob if prob 1e-12 else 1e-12) return logp_sum logp_pc calculate_logp(bayes_pc, test_data.head(200)) logp_notears calculate_logp(bayes_notears, test_data.head(200)) print(PC logp(前200样本):, logp_pc) print(Notears logp(前200样本):, logp_notears)logp计算时取了前200条样本因为每条样本都要做一次完整推理测试集大时比较耗时。1e-12的截断是为了防止log(0)出现负无穷工程上常见的处理。这里统计的联合对数似然其实是一个近似实际计算全联合分布需要对所有变量做消元代价更大用Depression的条件概率替代可以作为相对对比的参考指标。5. 模型评价指标的工程化应用阈值调整与结构稀疏性的联动5.1 从分类概率到决策阈值前面用np.argmax取后验概率最大的类别这是一种默认决策策略但它未必是最优的。当数据集中正负样本不均衡时后验概率本身会偏斜需要额外看AUC-ROC来排除阈值的影响。AUC-ROC描述的是模型对正负样本排序能力的强弱和阈值无关。可以在推理时保留所有测试样本的Depression后验概率然后动态调整阈值。from sklearn.metrics import roc_curve, auc def predict_proba_bn(infer, test_df, targetDepression, pos_class1): prob_pos [] for _, row in test_df.iterrows(): evidence {col: str(row[col]) for col in test_df.columns if col ! target} q infer.query([target], evidenceevidence) idx list(q.state_names[target]).index(pos_class) prob_pos.append(q.values[idx]) return np.array(prob_pos) proba_pc predict_proba_bn(infer_pc, test_data) proba_notears predict_proba_bn(infer_notears, test_data) fpr_pc, tpr_pc, _ roc_curve(y_true, proba_pc, pos_label1) fpr_nt, tpr_nt, _ roc_curve(y_true, proba_notears, pos_label1) print(PC AUC:, auc(fpr_pc, tpr_pc)) print(Notears AUC:, auc(fpr_nt, tpr_nt))AUC的解读要结合业务背景。抑郁症预测场景中召回率通常比精确率更重要因为漏报的代价更高。这时可以优先选择让tpr更高的阈值区间而不是默认的0.5。这里给出的概率输出接口也为后续绘制ROC曲线提供了数据直接plot即可。5.2 结构稀疏性与预测性能的联动分析从实际项目结果来看一个值得注意的现象是Notears学出的结构往往比PC更稠密。原因在于PC的独立性检验会剪掉条件独立的大量弱边而Notears的L1正则虽然约束了整体稀疏性但边的保留策略并不同——它倾向于保留权重较大的边集以保证全局损失最小。稠密结构在训练集上logp通常更高但泛化表现未必更好。做一个简单的验证实验把Notears的阈值从0.3提高到0.5再看AUC和F1的变化。阈值越高稀疏化后的边数越少结构越接近PC的结果。这个实验能直观展示「结构复杂度」和「预测性能」之间的权衡关系。for thresh in [0.2, 0.3, 0.4, 0.5]: W_tmp np.where(np.abs(W_est) thresh, W_est, 0) n_edges np.count_nonzero(W_tmp) # 构建临时bayes_net并评估 edges_tmp [] for i in range(len(nodes)): for j in range(len(nodes)): if W_tmp[i, j] ! 0 and i ! j: edges_tmp.append((nodes[i], nodes[j])) bn_tmp BayesianNetwork(edges_tmp) bn_tmp.fit(train_data, estimatorMaximumLikelihoodEstimator) infer_tmp VariableElimination(bn_tmp) proba_tmp predict_proba_bn(infer_tmp, test_data) fpr_t, tpr_t, _ roc_curve(y_true, proba_tmp, pos_label1) print(f阈值{thresh}: 边数{n_edges}条, AUC{auc(fpr_t, tpr_t):.4f})这个循环实际上在做一个结构级的超参数扫描每次迭代都完成一次完整的参数学习和推理评估。项目里可以将这个结果画成一条曲线横轴是阈值纵轴是AUC曲线峰值对应的就是当前数据下最优的稀疏化配置。做这类对比的另一个参考维度是计算时间在同样数据上PC的检验次数随条件集大小指数增长Notears的优化轮数相对稳定但每轮涉及矩阵指数运算在变量数超过30时单轮开销明显上涨。从工程角度PC算法在变量少于15个、且样本量在几千级别时表现稳定结构解释性强Notears更适合变量数中等、需要快速得到可微结构的场景。最后建议把两份DAG都导出成BIF格式或dot格式便于查看算法之间的结构差异。from pgmpy.readwrite import BIFWriter BIFWriter(bayes_pc).write_bif(bayes_pc.bif) BIFWriter(bayes_notears).write_bif(bayes_notears.bif) print(模型已导出为bif文件可用可视化工具打开检查)这里的BIF格式是贝叶斯网络的通用交换格式保存了节点、边和CPD信息。导出后用第三方软件浏览结构可以快速定位两个算法产生分歧的边结合领域知识判断哪边的因果方向更合理这是量化指标之外最重要的定性校验。将时间、AUC、logp三个维度的数据汇总成一张对比表就是整个贝叶斯网络建模项目完整的交付物。本文还有配套的精品资源点击获取