拓冰网站定制 
简介这份资源是面向数据分析与因果发现学习者的Python实现PC算法项目源码适合具备一定Python与统计基础、希望理解条件独立检验与因果结构学习的中高级开发者。项目围绕部分相关性思想演示如何从高维数据中剔除间接相关、识别直接因果关系可用于网络结构学习、特征筛选等场景。压缩包共13个文件约452KB以4个py源码文件为核心辅以3张png结果图、1个csv测试数据、1个md说明文档及txt依赖清单等源码、数据与可视化素材齐备。目前已有298人学习下载。读者可获取完整可运行的PC算法实现包括数据预处理、相关矩阵计算、条件独立测试、边剔除与循环迭代等关键环节并借助networkx绘制因果图直观查看结果同时参考README与依赖文件快速复现实验为后续扩展非高斯数据或结合GES、LiNGAM等算法提供实践起点。1. 从一份 Python 实现 PC 算法源码说起因果结构学习到底在算什么很多人第一次拿到「Python 实现 PC 算法项目源码」这类工程时会下意识把它当成一个统计脚本读一份 CSV跑几个相关系数输出一张图。真正跑起来才发现PC 算法要解决的是「在没有先验因果图的前提下从观测数据里恢复变量之间的有向无环图结构」它属于因果发现里的约束型方法核心不是拟合而是不断做条件独立性检验把不可能存在的边一条条删掉再给剩下的边定向。这份源码通常面向三类人做数据分析想从相关性走到因果假设的工程师、做风控或归因建模需要变量结构的人、以及学python基础语法之后想找一个有真实算法含量的项目练手的人。它不依赖深度学习框架纯python代码加 NumPy、pandas、scipy 就能跑对vscode python环境配置或pycharm配置python环境的要求也不高属于「装好python就能读」的那类项目。但它的坑也很集中条件独立性检验的显著性水平怎么设、样本量不够时检验统计量失真、定向阶段遇到冲突边怎么处理。这些问题在源码里往往只体现为几个参数不理解背后的统计含义改一个数结果就完全变样。下面按「原理 → 骨架 → 实现 → 调参排错 → 进阶」的顺序把这份源码拆开讲清楚。2. PC 算法源码背后的条件独立性检验与骨架学习原理2.1 PC 算法为什么先学骨架再定向PC 算法Peter-Clark 算法的名字来自它的两位提出者整个流程分两大阶段骨架学习skeleton discovery和边定向orientation。骨架学习的目标是找出一张无向图图中每条边代表两个变量在给定其他变量子集时仍然统计相关边定向则利用「对撞结构」collider和方向传播规则把无向边变成有向边最终得到 CPDAGcompleted partially directed acyclic graph。之所以先骨架后定向是因为条件独立性检验本身是对称的X 与 Y 在给定 Z 下独立和 Y 与 X 在给定 Z 下独立是一回事检验结果无法直接给出方向。方向信息只能从「X→Z←Y 这种对撞结构里 X 与 Y 边缘独立、但给定 Z 后相关」这种非对称现象中推断。源码里通常把这两步写成两个独立函数骨架阶段用sepset字典记录每条边被哪个条件集删掉定向阶段再回查这个字典。理解这一点很关键如果你只关心变量之间「有没有关系」骨架阶段就够了如果你要的是因果方向定向阶段的正确性完全依赖骨架阶段留下的sepset是否准确。很多源码 bug 都出在sepset记录不全导致定向时找不到对撞结构。2.2 条件独立性检验G² 检验与 Fisher-Z 检验怎么选PC 算法本身不规定用哪种检验源码里常见两种实现基于偏相关的 Fisher-Z 检验和基于离散数据的 G² 检验。选哪种取决于你的数据类型。检验方法适用数据类型统计量源码常见参数Fisher-Z连续、近似正态偏相关系数转 Z 分数alpha、corr_methodG²离散、分类似然比卡方alpha、degree卡方检验离散皮尔逊卡方alphaFisher-Z 的统计量公式是z 0.5 * ln((1r)/(1-r)) * sqrt(n - |S| - 3)其中 r 是给定条件集 S 后的偏相关系数n 是样本量。源码里一般用scipy.stats.norm算双尾 p 值p 值大于alpha就认为条件独立删边。G² 检验的统计量是2 * sum(O * ln(O/E))对每个变量的取值组合求和。离散数据用 G² 更稳因为偏相关对分类变量没有意义。源码里如果同时支持两种通常会有一个indep_test参数让你切换。提示样本量小于 200 时Fisher-Z 的近似会明显偏乐观容易把该删的边留下。这种情况要么增大alpha要么换 G² 并确保每个取值组合有足够样本。2.3 骨架学习的删边顺序与 sepset 记录骨架学习的标准流程是先建完全图然后按条件集大小l从 0 开始递增对每条边 (X, Y)在 X 的邻居集合里取大小为l的子集 S检验 X 与 Y 在给定 S 下是否独立。如果独立删边并把 S 记入sepset[(X,Y)]和sepset[(Y,X)]。这个顺序不能乱l从小到大是为了保证先删掉直接相关的边再处理需要更多条件才能独立的边。如果反过来可能用一个大条件集删掉本该保留的边骨架就错了。源码里通常写成两层循环外层l内层遍历边和邻居子集。# 骨架学习核心循环简化版 for l in range(max_cond_set_size 1): for (x, y) in list(edges): # 取 x 的邻居中排除 y 后大小为 l 的子集 neighbors [n for n in adj[x] if n ! y] if len(neighbors) l: continue for S in combinations(neighbors, l): p indep_test(x, y, list(S), data) if p alpha: # 条件独立删边 edges.discard((x, y)) sepset[(x, y)] set(S) sepset[(y, x)] set(S) break逻辑说明max_cond_set_size一般设成邻居数的上限实际项目里常设 3 到 5再大组合数爆炸且样本不够。indep_test返回 p 值alpha是显著性水平。break很关键一条边只要找到一个条件集使其独立就删不需要继续找。参数说明alpha越小删边越保守图越密alpha越大删边越激进可能误删真边。常见取值 0.01、0.05、0.1源码默认多为 0.05。3. Python 实现 PC 算法源码的模块拆解与最小可跑通代码3.1 源码目录结构与各模块职责一份结构清晰的 PC 算法源码通常长这样pc_algorithm/ ├── main.py # 入口读数据、调算法、输出图 ├── pc.py # 骨架学习 定向主流程 ├── independence.py # 条件独立性检验Fisher-Z / G² ├── orientation.py # 对撞结构识别与方向传播 ├── graph_utils.py # 邻接表、边集、CPDAG 工具函数 └── data/ └── sample.csv # 示例数据independence.py是最该先读的因为它决定了整个算法的统计基础。pc.py是主流程orientation.py最容易出 bug。graph_utils.py里通常有adjacency_to_edges、is_acyclic这类函数调试时很有用。如果你是从python爬虫教程或python爬虫可视化界面转过来的习惯了一个脚本干完所有事读这种分模块的源码会有点不适应。但 PC 算法的定向逻辑复杂拆开写反而好维护也方便你单独替换检验方法。3.2 用 NumPy 和 SciPy 写一个 Fisher-Z 检验函数下面是一个可以直接抄的 Fisher-Z 条件独立性检验实现import numpy as np from scipy.stats import norm def partial_corr(x, y, S, data): 计算给定条件集 S 后 x 与 y 的偏相关系数 if not S: r np.corrcoef(data[:, x], data[:, y])[0, 1] return r # 用回归残差法算偏相关 Z data[:, S] Z np.column_stack([np.ones(Z.shape[0]), Z]) # 加截距 beta_x np.linalg.lstsq(Z, data[:, x], rcondNone)[0] beta_y np.linalg.lstsq(Z, data[:, y], rcondNone)[0] res_x data[:, x] - Z beta_x res_y data[:, y] - Z beta_y return np.corrcoef(res_x, res_y)[0, 1] def fisher_z_test(x, y, S, data, alpha0.05): 返回 p 值p alpha 表示条件独立 n data.shape[0] r partial_corr(x, y, S, data) r np.clip(r, -0.999999, 0.999999) # 防止 log 溢出 z 0.5 * np.log((1 r) / (1 - r)) * np.sqrt(n - len(S) - 3) p 2 * (1 - norm.cdf(abs(z))) return p逻辑说明partial_corr用最小二乘回归残差法算偏相关比直接套公式稳尤其是条件集较大时。np.clip那行是必须的相关系数接近 ±1 时log会溢出源码里如果没这行遇到强相关变量会直接报错或返回 NaN。参数说明alpha是显著性水平n - len(S) - 3是自由度样本量必须大于len(S) 3否则sqrt里是负数。实际项目里如果条件集大小接近样本量检验结果不可信应该限制max_cond_set_size。3.3 骨架学习主循环的最小可跑通版本把检验函数接进骨架学习就是一个能跑的最小版本from itertools import combinations def learn_skeleton(data, alpha0.05, max_cond_set_size3): n_vars data.shape[1] adj {i: set(range(n_vars)) - {i} for i in range(n_vars)} sepset {} for l in range(max_cond_set_size 1): for x in range(n_vars): for y in list(adj[x]): if y not in adj[x]: continue neighbors [n for n in adj[x] if n ! y] if len(neighbors) l: continue for S in combinations(neighbors, l): p fisher_z_test(x, y, list(S), data, alpha) if p alpha: adj[x].discard(y) adj[y].discard(x) sepset[(x, y)] set(S) sepset[(y, x)] set(S) break return adj, sepset逻辑说明adj用字典存邻接表初始是完全图。外层l控制条件集大小内层遍历每条边。if y not in adj[x]是防止边在循环中被删后还继续处理。sepset双向记录定向阶段两个方向都要查。参数说明max_cond_set_size设 3 意味着最多用 3 个变量做条件组合数是 C(k,3)k 是邻居数。邻居超过 10 个时组合数上千跑得慢实际项目里常配合max_cond_set_size2或加剪枝。注意这个最小版本没有做邻居排序优化也没有并行。数据维度超过 30 时建议先做变量筛选或者用python的multiprocessing把不同边的检验并行化。4. 定向阶段、参数调优与常见报错排查4.1 对撞结构识别与方向传播规则骨架学完后定向阶段分三步先找对撞结构unshielded collider再应用 Meek 规则传播方向最后处理冲突边。对撞结构的判定条件是X 和 Y 不相邻但 X—Z—Y 存在且 Z 不在sepset[(X,Y)]里则定向为 X→Z←Y。def orient_colliders(adj, sepset): directed set() nodes list(adj.keys()) for z in nodes: neighbors list(adj[z]) for x, y in combinations(neighbors, 2): if y in adj[x]: # x 和 y 相邻不是对撞 continue if z not in sepset.get((x, y), set()): directed.add((x, z)) directed.add((y, z)) return directed逻辑说明sepset.get((x, y), set())用空集兜底防止 KeyError。z not in sepset是对撞结构的核心判据如果 Z 不在使 X、Y 独立的条件集里说明 X 和 Y 的独立必须经过 Z方向就是 X→Z←Y。参数说明这一步没有可调参数但依赖sepset的完整性。如果骨架阶段用了break提前退出sepset里存的是第一个使边独立的条件集可能不是唯一解定向时会有歧义。有些源码会存所有使边独立的条件集代价是内存和速度。4.2 alpha、max_cond_set_size 与样本量的三角关系这三个参数互相牵制调参时不能单独看参数调大后果调小后果建议alpha删边多图稀疏可能漏真边删边少图密可能留假边从 0.05 起样本大用 0.01max_cond_set_size检验更充分但组合爆炸检验不足骨架偏密邻居数 5 以内设 3否则设 2样本量 n检验更稳偏相关更可信检验失真p 值不可信至少 200最好 1000一个实用经验先固定alpha0.05、max_cond_set_size2跑一遍看骨架边数。如果边数接近完全图说明alpha太小或样本不够如果边数少得可怜说明alpha太大。再根据领域知识微调。4.3 常见报错与排查清单跑 PC 算法源码时报错集中在几类LinAlgError: Singular matrix条件集里的变量线性相关回归矩阵不可逆。排查方法是在partial_corr里加np.linalg.pinv替代lstsq或者先做变量去重。ValueError: math domain errorsqrt(n - len(S) - 3)里是负数说明条件集大小超过样本量。限制max_cond_set_size或增大样本。结果图有环定向阶段 Meek 规则实现有 bug或者骨架阶段sepset记录不全。用is_acyclic函数检查定位到具体边。跑得极慢邻居子集组合数爆炸。加max_cond_set_size上限或者对邻居按相关性排序只取前 k 个做条件集。提示调试时先把max_cond_set_size设成 1跑通全流程再逐步加大。这样能把骨架问题和定向问题分开定位。5. 把 PC 算法源码接进真实数据管道的三个进阶技巧5.1 用先验知识约束骨架减少无谓检验纯数据驱动的 PC 算法在变量多时跑得慢而且容易出假边。实际项目里我一般会加一个blacklist和whitelistblacklist里的边直接不检验whitelist里的边强制保留。比如已知某些变量在时间上有先后或者某些边在业务上不可能存在就可以提前排除。def learn_skeleton_with_prior(data, blacklistset(), whitelistset(), alpha0.05): n_vars data.shape[1] adj {i: set(range(n_vars)) - {i} for i in range(n_vars)} for (x, y) in blacklist: adj[x].discard(y) adj[y].discard(x) # whitelist 的边在删边时跳过 sepset {} for l in range(3): for x in range(n_vars): for y in list(adj[x]): if (x, y) in whitelist or (y, x) in whitelist: continue # ... 同前 return adj, sepset逻辑说明blacklist在初始化后直接删边whitelist在检验循环里跳过。这样既减少了检验次数又保证了业务上必须存在的边不被误删。参数说明blacklist和whitelist都是无向边集合用元组表示。注意whitelist的边在定向阶段仍然要参与只是骨架阶段不删。5.2 用 bootstrap 评估边的稳定性PC 算法对样本扰动敏感同一份数据重采样几次骨架可能不一样。评估稳定性的常见做法是 bootstrap重采样 B 次每次跑一遍骨架学习统计每条边出现的频率。频率低于 0.6 的边可以考虑删掉。import numpy as np def bootstrap_skeleton(data, B100, alpha0.05): n data.shape[0] edge_count {} for b in range(B): idx np.random.choice(n, n, replaceTrue) sample data[idx] adj, _ learn_skeleton(sample, alpha) for x in adj: for y in adj[x]: if x y: edge_count[(x, y)] edge_count.get((x, y), 0) 1 return {e: c / B for e, c in edge_count.items()}逻辑说明每次 bootstrap 重采样后跑骨架统计边出现次数最后除以 B 得到频率。x y是为了避免同一条边统计两次。参数说明B一般取 100 到 500越大越稳但越慢。频率阈值 0.6 到 0.8 之间取决于你对假阳性的容忍度。5.3 输出 CPDAG 并用 networkx 可视化验证最后一步是把结果转成 CPDAG 并画出来。用networkx几行就能画import networkx as nx import matplotlib.pyplot as plt def plot_cpdag(adj, directed): G nx.DiGraph() for x in adj: for y in adj[x]: if (x, y) in directed: G.add_edge(x, y) elif (y, x) not in directed and x y: G.add_edge(x, y, styledashed) # 无向边用虚线 pos nx.spring_layout(G) nx.draw(G, pos, with_labelsTrue, node_colorlightblue) plt.show()逻辑说明有向边画实线无向边画虚线。x y避免无向边画两次。spring_layout是力导向布局变量多时可能重叠可以换kamada_kawai_layout。参数说明directed是定向阶段输出的有向边集合。可视化后重点看有没有明显的对撞结构被漏掉以及无向边是否集中在某个子图里——如果集中说明那部分数据的条件独立性检验不够充分。到这里一份 Python 实现 PC 算法源码从原理到落地的主线就走完了。真正决定结果质量的往往不是代码写得多漂亮而是alpha、max_cond_set_size和样本量这三个数有没有根据你的数据调到位以及sepset有没有记全。本文还有配套的精品资源点击获取 Python实现PC算法源码解析：因果结构学习与条件独立性检验
简介这份资源是面向数据分析与因果发现学习者的Python实现PC算法项目源码适合具备一定Python与统计基础、希望理解条件独立检验与因果结构学习的中高级开发者。项目围绕部分相关性思想演示如何从高维数据中剔除间接相关、识别直接因果关系可用于网络结构学习、特征筛选等场景。压缩包共13个文件约452KB以4个py源码文件为核心辅以3张png结果图、1个csv测试数据、1个md说明文档及txt依赖清单等源码、数据与可视化素材齐备。目前已有298人学习下载。读者可获取完整可运行的PC算法实现包括数据预处理、相关矩阵计算、条件独立测试、边剔除与循环迭代等关键环节并借助networkx绘制因果图直观查看结果同时参考README与依赖文件快速复现实验为后续扩展非高斯数据或结合GES、LiNGAM等算法提供实践起点。1. 从一份 Python 实现 PC 算法源码说起因果结构学习到底在算什么很多人第一次拿到「Python 实现 PC 算法项目源码」这类工程时会下意识把它当成一个统计脚本读一份 CSV跑几个相关系数输出一张图。真正跑起来才发现PC 算法要解决的是「在没有先验因果图的前提下从观测数据里恢复变量之间的有向无环图结构」它属于因果发现里的约束型方法核心不是拟合而是不断做条件独立性检验把不可能存在的边一条条删掉再给剩下的边定向。这份源码通常面向三类人做数据分析想从相关性走到因果假设的工程师、做风控或归因建模需要变量结构的人、以及学python基础语法之后想找一个有真实算法含量的项目练手的人。它不依赖深度学习框架纯python代码加 NumPy、pandas、scipy 就能跑对vscode python环境配置或pycharm配置python环境的要求也不高属于「装好python就能读」的那类项目。但它的坑也很集中条件独立性检验的显著性水平怎么设、样本量不够时检验统计量失真、定向阶段遇到冲突边怎么处理。这些问题在源码里往往只体现为几个参数不理解背后的统计含义改一个数结果就完全变样。下面按「原理 → 骨架 → 实现 → 调参排错 → 进阶」的顺序把这份源码拆开讲清楚。2. PC 算法源码背后的条件独立性检验与骨架学习原理2.1 PC 算法为什么先学骨架再定向PC 算法Peter-Clark 算法的名字来自它的两位提出者整个流程分两大阶段骨架学习skeleton discovery和边定向orientation。骨架学习的目标是找出一张无向图图中每条边代表两个变量在给定其他变量子集时仍然统计相关边定向则利用「对撞结构」collider和方向传播规则把无向边变成有向边最终得到 CPDAGcompleted partially directed acyclic graph。之所以先骨架后定向是因为条件独立性检验本身是对称的X 与 Y 在给定 Z 下独立和 Y 与 X 在给定 Z 下独立是一回事检验结果无法直接给出方向。方向信息只能从「X→Z←Y 这种对撞结构里 X 与 Y 边缘独立、但给定 Z 后相关」这种非对称现象中推断。源码里通常把这两步写成两个独立函数骨架阶段用sepset字典记录每条边被哪个条件集删掉定向阶段再回查这个字典。理解这一点很关键如果你只关心变量之间「有没有关系」骨架阶段就够了如果你要的是因果方向定向阶段的正确性完全依赖骨架阶段留下的sepset是否准确。很多源码 bug 都出在sepset记录不全导致定向时找不到对撞结构。2.2 条件独立性检验G² 检验与 Fisher-Z 检验怎么选PC 算法本身不规定用哪种检验源码里常见两种实现基于偏相关的 Fisher-Z 检验和基于离散数据的 G² 检验。选哪种取决于你的数据类型。检验方法适用数据类型统计量源码常见参数Fisher-Z连续、近似正态偏相关系数转 Z 分数alpha、corr_methodG²离散、分类似然比卡方alpha、degree卡方检验离散皮尔逊卡方alphaFisher-Z 的统计量公式是z 0.5 * ln((1r)/(1-r)) * sqrt(n - |S| - 3)其中 r 是给定条件集 S 后的偏相关系数n 是样本量。源码里一般用scipy.stats.norm算双尾 p 值p 值大于alpha就认为条件独立删边。G² 检验的统计量是2 * sum(O * ln(O/E))对每个变量的取值组合求和。离散数据用 G² 更稳因为偏相关对分类变量没有意义。源码里如果同时支持两种通常会有一个indep_test参数让你切换。提示样本量小于 200 时Fisher-Z 的近似会明显偏乐观容易把该删的边留下。这种情况要么增大alpha要么换 G² 并确保每个取值组合有足够样本。2.3 骨架学习的删边顺序与 sepset 记录骨架学习的标准流程是先建完全图然后按条件集大小l从 0 开始递增对每条边 (X, Y)在 X 的邻居集合里取大小为l的子集 S检验 X 与 Y 在给定 S 下是否独立。如果独立删边并把 S 记入sepset[(X,Y)]和sepset[(Y,X)]。这个顺序不能乱l从小到大是为了保证先删掉直接相关的边再处理需要更多条件才能独立的边。如果反过来可能用一个大条件集删掉本该保留的边骨架就错了。源码里通常写成两层循环外层l内层遍历边和邻居子集。# 骨架学习核心循环简化版 for l in range(max_cond_set_size 1): for (x, y) in list(edges): # 取 x 的邻居中排除 y 后大小为 l 的子集 neighbors [n for n in adj[x] if n ! y] if len(neighbors) l: continue for S in combinations(neighbors, l): p indep_test(x, y, list(S), data) if p alpha: # 条件独立删边 edges.discard((x, y)) sepset[(x, y)] set(S) sepset[(y, x)] set(S) break逻辑说明max_cond_set_size一般设成邻居数的上限实际项目里常设 3 到 5再大组合数爆炸且样本不够。indep_test返回 p 值alpha是显著性水平。break很关键一条边只要找到一个条件集使其独立就删不需要继续找。参数说明alpha越小删边越保守图越密alpha越大删边越激进可能误删真边。常见取值 0.01、0.05、0.1源码默认多为 0.05。3. Python 实现 PC 算法源码的模块拆解与最小可跑通代码3.1 源码目录结构与各模块职责一份结构清晰的 PC 算法源码通常长这样pc_algorithm/ ├── main.py # 入口读数据、调算法、输出图 ├── pc.py # 骨架学习 定向主流程 ├── independence.py # 条件独立性检验Fisher-Z / G² ├── orientation.py # 对撞结构识别与方向传播 ├── graph_utils.py # 邻接表、边集、CPDAG 工具函数 └── data/ └── sample.csv # 示例数据independence.py是最该先读的因为它决定了整个算法的统计基础。pc.py是主流程orientation.py最容易出 bug。graph_utils.py里通常有adjacency_to_edges、is_acyclic这类函数调试时很有用。如果你是从python爬虫教程或python爬虫可视化界面转过来的习惯了一个脚本干完所有事读这种分模块的源码会有点不适应。但 PC 算法的定向逻辑复杂拆开写反而好维护也方便你单独替换检验方法。3.2 用 NumPy 和 SciPy 写一个 Fisher-Z 检验函数下面是一个可以直接抄的 Fisher-Z 条件独立性检验实现import numpy as np from scipy.stats import norm def partial_corr(x, y, S, data): 计算给定条件集 S 后 x 与 y 的偏相关系数 if not S: r np.corrcoef(data[:, x], data[:, y])[0, 1] return r # 用回归残差法算偏相关 Z data[:, S] Z np.column_stack([np.ones(Z.shape[0]), Z]) # 加截距 beta_x np.linalg.lstsq(Z, data[:, x], rcondNone)[0] beta_y np.linalg.lstsq(Z, data[:, y], rcondNone)[0] res_x data[:, x] - Z beta_x res_y data[:, y] - Z beta_y return np.corrcoef(res_x, res_y)[0, 1] def fisher_z_test(x, y, S, data, alpha0.05): 返回 p 值p alpha 表示条件独立 n data.shape[0] r partial_corr(x, y, S, data) r np.clip(r, -0.999999, 0.999999) # 防止 log 溢出 z 0.5 * np.log((1 r) / (1 - r)) * np.sqrt(n - len(S) - 3) p 2 * (1 - norm.cdf(abs(z))) return p逻辑说明partial_corr用最小二乘回归残差法算偏相关比直接套公式稳尤其是条件集较大时。np.clip那行是必须的相关系数接近 ±1 时log会溢出源码里如果没这行遇到强相关变量会直接报错或返回 NaN。参数说明alpha是显著性水平n - len(S) - 3是自由度样本量必须大于len(S) 3否则sqrt里是负数。实际项目里如果条件集大小接近样本量检验结果不可信应该限制max_cond_set_size。3.3 骨架学习主循环的最小可跑通版本把检验函数接进骨架学习就是一个能跑的最小版本from itertools import combinations def learn_skeleton(data, alpha0.05, max_cond_set_size3): n_vars data.shape[1] adj {i: set(range(n_vars)) - {i} for i in range(n_vars)} sepset {} for l in range(max_cond_set_size 1): for x in range(n_vars): for y in list(adj[x]): if y not in adj[x]: continue neighbors [n for n in adj[x] if n ! y] if len(neighbors) l: continue for S in combinations(neighbors, l): p fisher_z_test(x, y, list(S), data, alpha) if p alpha: adj[x].discard(y) adj[y].discard(x) sepset[(x, y)] set(S) sepset[(y, x)] set(S) break return adj, sepset逻辑说明adj用字典存邻接表初始是完全图。外层l控制条件集大小内层遍历每条边。if y not in adj[x]是防止边在循环中被删后还继续处理。sepset双向记录定向阶段两个方向都要查。参数说明max_cond_set_size设 3 意味着最多用 3 个变量做条件组合数是 C(k,3)k 是邻居数。邻居超过 10 个时组合数上千跑得慢实际项目里常配合max_cond_set_size2或加剪枝。注意这个最小版本没有做邻居排序优化也没有并行。数据维度超过 30 时建议先做变量筛选或者用python的multiprocessing把不同边的检验并行化。4. 定向阶段、参数调优与常见报错排查4.1 对撞结构识别与方向传播规则骨架学完后定向阶段分三步先找对撞结构unshielded collider再应用 Meek 规则传播方向最后处理冲突边。对撞结构的判定条件是X 和 Y 不相邻但 X—Z—Y 存在且 Z 不在sepset[(X,Y)]里则定向为 X→Z←Y。def orient_colliders(adj, sepset): directed set() nodes list(adj.keys()) for z in nodes: neighbors list(adj[z]) for x, y in combinations(neighbors, 2): if y in adj[x]: # x 和 y 相邻不是对撞 continue if z not in sepset.get((x, y), set()): directed.add((x, z)) directed.add((y, z)) return directed逻辑说明sepset.get((x, y), set())用空集兜底防止 KeyError。z not in sepset是对撞结构的核心判据如果 Z 不在使 X、Y 独立的条件集里说明 X 和 Y 的独立必须经过 Z方向就是 X→Z←Y。参数说明这一步没有可调参数但依赖sepset的完整性。如果骨架阶段用了break提前退出sepset里存的是第一个使边独立的条件集可能不是唯一解定向时会有歧义。有些源码会存所有使边独立的条件集代价是内存和速度。4.2 alpha、max_cond_set_size 与样本量的三角关系这三个参数互相牵制调参时不能单独看参数调大后果调小后果建议alpha删边多图稀疏可能漏真边删边少图密可能留假边从 0.05 起样本大用 0.01max_cond_set_size检验更充分但组合爆炸检验不足骨架偏密邻居数 5 以内设 3否则设 2样本量 n检验更稳偏相关更可信检验失真p 值不可信至少 200最好 1000一个实用经验先固定alpha0.05、max_cond_set_size2跑一遍看骨架边数。如果边数接近完全图说明alpha太小或样本不够如果边数少得可怜说明alpha太大。再根据领域知识微调。4.3 常见报错与排查清单跑 PC 算法源码时报错集中在几类LinAlgError: Singular matrix条件集里的变量线性相关回归矩阵不可逆。排查方法是在partial_corr里加np.linalg.pinv替代lstsq或者先做变量去重。ValueError: math domain errorsqrt(n - len(S) - 3)里是负数说明条件集大小超过样本量。限制max_cond_set_size或增大样本。结果图有环定向阶段 Meek 规则实现有 bug或者骨架阶段sepset记录不全。用is_acyclic函数检查定位到具体边。跑得极慢邻居子集组合数爆炸。加max_cond_set_size上限或者对邻居按相关性排序只取前 k 个做条件集。提示调试时先把max_cond_set_size设成 1跑通全流程再逐步加大。这样能把骨架问题和定向问题分开定位。5. 把 PC 算法源码接进真实数据管道的三个进阶技巧5.1 用先验知识约束骨架减少无谓检验纯数据驱动的 PC 算法在变量多时跑得慢而且容易出假边。实际项目里我一般会加一个blacklist和whitelistblacklist里的边直接不检验whitelist里的边强制保留。比如已知某些变量在时间上有先后或者某些边在业务上不可能存在就可以提前排除。def learn_skeleton_with_prior(data, blacklistset(), whitelistset(), alpha0.05): n_vars data.shape[1] adj {i: set(range(n_vars)) - {i} for i in range(n_vars)} for (x, y) in blacklist: adj[x].discard(y) adj[y].discard(x) # whitelist 的边在删边时跳过 sepset {} for l in range(3): for x in range(n_vars): for y in list(adj[x]): if (x, y) in whitelist or (y, x) in whitelist: continue # ... 同前 return adj, sepset逻辑说明blacklist在初始化后直接删边whitelist在检验循环里跳过。这样既减少了检验次数又保证了业务上必须存在的边不被误删。参数说明blacklist和whitelist都是无向边集合用元组表示。注意whitelist的边在定向阶段仍然要参与只是骨架阶段不删。5.2 用 bootstrap 评估边的稳定性PC 算法对样本扰动敏感同一份数据重采样几次骨架可能不一样。评估稳定性的常见做法是 bootstrap重采样 B 次每次跑一遍骨架学习统计每条边出现的频率。频率低于 0.6 的边可以考虑删掉。import numpy as np def bootstrap_skeleton(data, B100, alpha0.05): n data.shape[0] edge_count {} for b in range(B): idx np.random.choice(n, n, replaceTrue) sample data[idx] adj, _ learn_skeleton(sample, alpha) for x in adj: for y in adj[x]: if x y: edge_count[(x, y)] edge_count.get((x, y), 0) 1 return {e: c / B for e, c in edge_count.items()}逻辑说明每次 bootstrap 重采样后跑骨架统计边出现次数最后除以 B 得到频率。x y是为了避免同一条边统计两次。参数说明B一般取 100 到 500越大越稳但越慢。频率阈值 0.6 到 0.8 之间取决于你对假阳性的容忍度。5.3 输出 CPDAG 并用 networkx 可视化验证最后一步是把结果转成 CPDAG 并画出来。用networkx几行就能画import networkx as nx import matplotlib.pyplot as plt def plot_cpdag(adj, directed): G nx.DiGraph() for x in adj: for y in adj[x]: if (x, y) in directed: G.add_edge(x, y) elif (y, x) not in directed and x y: G.add_edge(x, y, styledashed) # 无向边用虚线 pos nx.spring_layout(G) nx.draw(G, pos, with_labelsTrue, node_colorlightblue) plt.show()逻辑说明有向边画实线无向边画虚线。x y避免无向边画两次。spring_layout是力导向布局变量多时可能重叠可以换kamada_kawai_layout。参数说明directed是定向阶段输出的有向边集合。可视化后重点看有没有明显的对撞结构被漏掉以及无向边是否集中在某个子图里——如果集中说明那部分数据的条件独立性检验不够充分。到这里一份 Python 实现 PC 算法源码从原理到落地的主线就走完了。真正决定结果质量的往往不是代码写得多漂亮而是alpha、max_cond_set_size和样本量这三个数有没有根据你的数据调到位以及sepset有没有记全。本文还有配套的精品资源点击获取