企业官网建设流程全解析

简介：概率图模型（比如贝叶斯网络）是处理不确定性和因果推断的重要工具，其核心环节是结构学习。结构学习解决如何从数据中自动发现变量依赖关系的问题，常见的有基于条件独立性检验的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(columns=drop_cols, inplace=True) # 用众数填补离散变量的缺失值 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], q=4, labels=[0, 1, 2, 3], duplicates="drop") df.drop(columns=[col], inplace=True) 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_size=0.3, random_state=42, stratify=df["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(threshold=0.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(data=train_data) dag_pc = pc.estimate( variant="stable", ci_test="chi_square", significance_level=0.05, max_cond_vars=5 ) end_pc = time.time() print("PC算法耗时:", round(end_pc - start_pc, 4), "秒") print("PC学习到的边数:", len(dag_pc.edges()))

variant="stable"是PC的稳定版本，它改变了冲突边的删除顺序，结果不再受变量输入顺序影响，这在工程上非常重要。ci_test="chi_square"适用于离散数据，max_cond_vars=5限制了条件集的维度，防止高维条件下的检验因样本稀疏而失效。如果数据量较小，这个值建议设为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, errors="coerce").values.astype(np.float64) X_train = np.nan_to_num(X_train, nan=0.0) start_notears = time.time() # 使用默认lambda参数，lambda越大图越稀疏 W_est = notears.linear_model(X_train, lambda1=0.05, loss="l2") 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用于把优化结果中接近零的权重清零，因为增广拉格朗日法收敛后会有数值噪声，不会精确等于零。loss="l2"表示使用最小二乘损失，适用连续数据。需要说明的是，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, estimator=MaximumLikelihoodEstimator) # 用Notears学到的结构构建模型 bayes_notears = BayesianNetwork(dag_notears.edges()) bayes_notears.fit(train_data, estimator=MaximumLikelihoodEstimator) # 查看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, estimator=BayesianEstimator, prior_type="BDeu", equivalent_sample_size=10)

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, target="Depression"): 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], evidence=evidence) # 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_label="1") recall = recall_score(y_true, y_pred, pos_label="1") f1 = f1_score(y_true, y_pred, pos_label="1") cm = confusion_matrix(y_true, y_pred) print(f"===== {model_name} =====") print(f"准确率: {acc:.4f}") print(f"精确率: {precision:.4f}") print(f"召回率: {recall:.4f}") print(f"F1-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"], evidence=evidence) 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, target="Depression", pos_class="1"): 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], evidence=evidence) 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_label="1") fpr_nt, tpr_nt, _ = roc_curve(y_true, proba_notears, pos_label="1") 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, estimator=MaximumLikelihoodEstimator) 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_label="1") 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三个维度的数据汇总成一张对比表，就是整个贝叶斯网络建模项目完整的交付物。

本文还有配套的精品资源，点击获取