简介这份Python实现PC算法的项目源码面向具备一定统计与编程基础、希望深入因果发现与网络结构学习的数据分析者和机器学习学习者。PC算法通过部分相关性检验变量间的条件独立关系可用于高维数据中剔除间接相关、揭示直接因果结构。资源包共13个文件约452KB以4个py源码文件为核心辅以3张png结果示意图、1个csv测试数据集、1个md说明文档及若干配置文件结构紧凑、便于快速上手。项目围绕数据预处理、相关矩阵计算、条件独立测试、定向边剔除与循环迭代等核心步骤展开并给出基于networkx的图结构可视化思路以及面向非高斯数据和大规模数据的扩展优化方向。已有298人学习适合想理解PC算法原理并动手复现因果发现流程的读者参考。1. PC算法到底在算什么从条件独立到因果骨架的那条线PC 算法Peter-Clark Algorithm是因果发现领域最经典的约束型算法之一它要解决的问题很具体手上只有一堆观测数据没有任何先验的因果方向信息能不能把变量之间的因果骨架图给还原出来。很多人第一次听到「Python实现PC算法项目源码」这个标题脑子里浮现的是某个现成的 pip 包装完调个函数就出图。实际情况是PC 算法本身不复杂但它的实现细节里藏着大量统计学上的坑条件独立检验选什么、显著性水平怎么定、样本量够不够每一个都会直接改变最终输出的图结构。这篇文章面向两类人一类是想把 PC 算法真正跑起来、拿到可解释因果图的工程师和数据分析师另一类是想读懂 PC 算法源码、搞清楚每一步在做什么、方便自己改造成业务版本的人。我会从算法骨架讲起然后给出一份可以直接复现的 Python 实现再逐层拆解参数、踩坑点和验证方法。读完你应该能自己写出一份不依赖第三方因果库的 PC 算法代码并且知道在什么数据规模下它还能用、什么时候该换算法。2. PC算法的骨架拆解为什么先做骨架再做定向2.1 从完全图到稀疏骨架条件独立检验在做什么PC 算法的核心思想可以用一句话概括如果两个变量 X 和 Y 在给定某个变量集合 S 的条件下独立那么 X 和 Y 之间就不应该有边。算法从一个完全无向图出发对每一对变量尝试寻找一个条件集使得它们在给定该条件集时条件独立。如果找到了就把这条边删掉。这里的关键在于「条件集从哪来」。PC 算法的做法是对于当前还和 X 相邻的节点按邻接集大小从 0 开始逐层搜索。第 0 层就是看 X 和 Y 的边缘独立性第 1 层看给定一个邻居时是否独立第 2 层看给定两个邻居时是否独立以此类推。这个逐层扩展的过程保证了算法不会漏掉低阶的条件独立关系同时把搜索空间控制在可接受的范围内。条件独立检验通常用偏相关系数或者 Fisher Z 检验。对于线性高斯数据偏相关系数为零等价于条件独立所以用 Fisher Z 变换把偏相关系数转成近似正态统计量再做显著性检验这是最标准的做法。对于离散数据一般用 G² 检验或卡方检验。选哪种检验取决于你的数据类型选错了后面全错。2.2 定向规则把无向骨架变成部分有向图骨架建好之后PC 算法用一组定向规则把无向边变成有向边。最核心的是三条第一条对撞结构。如果 X 和 Z 不相邻但 X→Y←Z 这种结构存在也就是 Y 同时和 X、Z 相邻而 X 和 Z 之间没有边并且 Y 不在 X 和 Z 的条件集里那么 X→Y 和 Z→Y 的方向就确定了。这条规则是 PC 算法能定向的根本原因。第二条避免新对撞。如果已经确定了一条有向边 X→Y而 Y 和 Z 相邻X 和 Z 不相邻那么 Y→Z 的方向可以确定否则会形成一个新的对撞结构。第三条避免环。定向过程中不能产生有向环如果某条边的方向会导致环就反向或者保持无向。这三条规则反复应用直到没有新的边可以被定向。最终输出的图叫 CPDAGCompleted Partially Directed Acyclic Graph里面有些边是有向的有些还是无向的无向边表示在当前数据下无法确定方向。2.3 为什么样本量和检验阈值会直接改变图结构PC 算法对样本量非常敏感。条件独立检验的统计功效随样本量增加而提高样本量不够的时候本来应该被删掉的边删不掉图会偏密样本量太大的时候微弱的依赖关系也会被检验出来图可能偏密。更麻烦的是PC 算法是逐层搜索的一旦某一层删错了一条边后面的定向规则会基于错误的骨架继续推错误会累积。显著性水平 α 的选择也是玄学。α 设大了条件独立被误判为依赖边删不掉α 设小了依赖被误判为独立边被误删。实践中一般从 0.01 到 0.05 之间试但真正靠谱的做法是做敏感性分析看不同 α 下图结构的变化。如果 α 从 0.01 变到 0.05 图结构剧烈变化说明数据量不够或者变量间关系太弱这时候 PC 算法的输出不可信。3. 用Python从零实现PC算法核心代码与参数说明3.1 环境准备与依赖选择实现 PC 算法不需要太重的依赖。核心就是 numpy 做矩阵运算scipy 做统计检验networkx 做图结构管理。不建议一上来就用 causal-learn 或者 pgmpy那些库封装太厚出了问题不好排查。自己写一遍后面调参和改造都方便。pip install numpy scipy networkx pandas matplotlib这几个包都是常规科学计算栈版本没有特别要求numpy 1.20 以上、scipy 1.6 以上就行。pandas 用来读数据matplotlib 用来画图不是必须的但调试的时候很有用。3.2 条件独立检验偏相关系数与Fisher Z变换偏相关系数的计算是 PC 算法的计算核心。给定变量集合 S要算 X 和 Y 在给定 S 下的偏相关系数标准做法是用回归残差的相关性来算。import numpy as np from scipy import stats def partial_corr(x, y, S, data): 计算 x 和 y 在给定 S 下的偏相关系数 x, y: 变量名或列索引 S: 条件集列表形式 data: pandas DataFrame if len(S) 0: r np.corrcoef(data[x], data[y])[0, 1] return r # 用线性回归去掉 S 的影响 from numpy.linalg import lstsq Z data[S].values Z np.column_stack([np.ones(len(Z)), Z]) # 加截距项 # 对 x 和 y 分别回归取残差 coef_x, _, _, _ lstsq(Z, data[x].values, rcondNone) coef_y, _, _, _ lstsq(Z, data[y].values, rcondNone) res_x data[x].values - Z coef_x res_y data[y].values - Z coef_y r np.corrcoef(res_x, res_y)[0, 1] return r def fisher_z_test(x, y, S, data, alpha0.05): Fisher Z 检验返回是否独立 n len(data) r partial_corr(x, y, S, data) # Fisher Z 变换 z 0.5 * np.log((1 r) / (1 - r)) # 标准差 se 1.0 / np.sqrt(n - len(S) - 3) # 检验统计量 stat abs(z) / se # 双尾检验 p_value 2 * (1 - stats.norm.cdf(stat)) return p_value alpha, p_value这段代码里有两个关键点。第一偏相关系数的计算用的是回归残差法这是最直观也最稳定的做法比直接套公式算协方差矩阵的逆要稳。第二Fisher Z 变换的自由度是 n - len(S) - 3这个 3 是固定的来自变换本身的方差近似。如果样本量 n 小于 len(S) 3这个检验就没法做了实践中要保证 n 至少是条件集大小的 5 到 10 倍。参数 alpha 是显著性水平默认 0.05。返回的 p_value 可以用来做敏感性分析看不同阈值下哪些边会被删掉。3.3 骨架搜索逐层条件集与邻接表更新骨架搜索是 PC 算法最耗时的部分。核心逻辑是对每一对相邻节点从空集开始逐步扩大条件集直到找到一组条件使得它们独立或者条件集大小超过当前邻接集。def skeleton_discovery(data, alpha0.05, max_cond_sizeNone): PC 算法骨架搜索 返回无向图邻接表和分离集 variables list(data.columns) n_vars len(variables) # 初始化完全图 adj {v: set(variables) - {v} for v in variables} sep_set {} # 记录分离集用于后续定向 if max_cond_size is None: max_cond_size n_vars - 2 for cond_size in range(max_cond_size 1): # 收集所有需要检验的边 edges_to_check [] for x in variables: for y in adj[x]: if x y: # 避免重复 edges_to_check.append((x, y)) for x, y in edges_to_check: if y not in adj[x]: continue # 边已经被删了 # 候选条件集从 x 的邻居中选排除 y neighbors adj[x] - {y} if len(neighbors) cond_size: continue from itertools import combinations for S in combinations(neighbors, cond_size): independent, p_val fisher_z_test(x, y, list(S), data, alpha) if independent: # 删除边 adj[x].discard(y) adj[y].discard(x) sep_set[(x, y)] list(S) sep_set[(y, x)] list(S) break return adj, sep_set这段代码有几个工程上的细节值得说。第一条件集是从 x 的邻居里选的不是从所有变量里选这是 PC 算法的标准做法能大幅减少检验次数。第二用 combinations 生成条件集当邻居数量多的时候组合数会爆炸所以实践中要限制 max_cond_size一般不超过 3 到 4。第三sep_set 记录了每对变量是在哪个条件集下被判定独立的这个信息在定向阶段要用。参数 max_cond_size 控制搜索深度。设太小会漏掉高阶条件独立关系图偏密设太大会导致计算量指数增长。经验值是变量数的三分之一到一半但不超过 5。3.4 定向规则实现对撞结构与避免新对撞骨架建好之后定向阶段要把无向边变成有向边。核心是识别对撞结构然后传播方向。def orient_edges(adj, sep_set, variables): 定向阶段识别对撞结构并传播方向 directed set() # 有向边 (x, y) 表示 x - y undirected set() # 收集所有无向边 for x in variables: for y in adj[x]: if x y: undirected.add((x, y)) # 规则1识别对撞结构 X - Y - Z for y in variables: neighbors list(adj[y]) for i in range(len(neighbors)): for j in range(i 1, len(neighbors)): x, z neighbors[i], neighbors[j] if z in adj[x]: continue # x 和 z 相邻不是对撞 # 检查 y 是否在 sep_set[(x, z)] 中 sep sep_set.get((x, z), []) if y not in sep: # 对撞结构成立 directed.add((x, y)) directed.add((z, y)) undirected.discard((min(x, y), max(x, y))) undirected.discard((min(z, y), max(z, y))) # 规则2避免新对撞 changed True while changed: changed False for x, y in list(undirected): # 如果 x - y 已经确定检查 y 的其他邻居 if (x, y) in directed: for z in adj[y]: if z x: continue if (min(y, z), max(y, z)) in undirected: if z not in adj[x]: # y - z 方向确定 directed.add((y, z)) undirected.discard((min(y, z), max(y, z))) changed True return directed, undirected定向规则里最容易出错的是对撞结构的判断条件。必须同时满足三个条件X 和 Z 不相邻、Y 同时和 X 和 Z 相邻、Y 不在 X 和 Z 的分离集中。第三个条件最容易被忽略如果 Y 在分离集中说明 X 和 Z 的独立性是 Y 导致的这时候不能定向为对撞。规则2 的实现是一个迭代过程因为定向一条边可能会触发新的定向。循环直到没有新的边可以被定向为止。3.5 完整流程串联与输出解读把上面的模块串起来就是一个完整的 PC 算法实现。def pc_algorithm(data, alpha0.05, max_cond_sizeNone): 完整的 PC 算法 variables list(data.columns) # 第一步骨架搜索 adj, sep_set skeleton_discovery(data, alpha, max_cond_size) # 第二步定向 directed, undirected orient_edges(adj, sep_set, variables) return { adjacency: adj, sep_set: sep_set, directed: directed, undirected: undirected } # 使用示例 import pandas as pd import numpy as np # 生成模拟数据X - Y - Z, X - Z np.random.seed(42) n 1000 X np.random.randn(n) Y 0.8 * X np.random.randn(n) * 0.5 Z 0.6 * X 0.7 * Y np.random.randn(n) * 0.5 data pd.DataFrame({X: X, Y: Y, Z: Z}) result pc_algorithm(data, alpha0.01) print(有向边, result[directed]) print(无向边, result[undirected])输出解读的时候要注意PC 算法输出的是 CPDAG不是唯一的因果图。有向边表示在所有马尔可夫等价类中方向一致无向边表示方向不确定。如果你看到 X→Y 和 Y→Z 都是有向的但 X 和 Z 之间没有边这不一定意味着 X 和 Z 独立可能只是条件独立检验没找到它们之间的直接依赖。4. PC算法落地时的避坑清单从数据预处理到结果验证4.1 数据预处理没做好后面全白搭PC 算法对数据的假设是连续变量、线性关系、高斯噪声。如果你的数据里有分类变量直接扔进去算偏相关系数会得到完全错误的结果。分类变量要么先做独热编码然后当连续变量处理效果一般要么换用基于互信息的条件独立检验。缺失值也是大问题PC 算法没有内置的缺失值处理机制要么删样本要么插补插补方法的选择会直接影响条件独立检验的结果。另一个容易被忽略的是变量尺度。偏相关系数本身对尺度不敏感但如果你在预处理阶段做了标准化要注意标准化是在全量数据上做的还是分训练测试集做的。PC 算法一般不需要划分训练测试集但如果你要做交叉验证来选 α标准化必须在每折内部独立做否则会信息泄露。4.2 条件独立检验选错类型图结构完全变样连续高斯数据用 Fisher Z 检验离散数据用 G² 检验混合数据用条件互信息或者基于核的方法。选错了检验类型不是精度下降的问题是根本性错误。比如离散数据用 Fisher Z偏相关系数根本没有意义算出来的 p 值也是错的。还有一个隐蔽的坑Fisher Z 检验假设变量间是线性关系。如果真实关系是非线性的比如 Y X²偏相关系数可能接近零检验会判定独立但实际上 X 和 Y 有强依赖。这种情况下需要换用基于互信息或距离相关的方法但那些方法的计算量会大很多。4.3 样本量不够时PC算法输出的图不可信PC 算法的最低样本量要求没有严格公式但经验上每个变量至少需要 10 到 20 个样本条件集大小每增加 1样本量需求大概翻倍。如果你有 20 个变量、200 个样本跑 PC 算法大概率会得到一张乱七八糟的图。判断样本量够不够的一个实用方法是跑多次 Bootstrap看边出现的频率。如果某条边在 80% 以上的 Bootstrap 样本中都出现可以认为它是稳定的如果只有 50% 左右说明这条边不可靠。这个做法比单纯看 p 值要靠谱得多。4.4 定向规则实现中的边界情况对撞结构识别的时候如果 X 和 Z 之间本来有边但后来被删了sep_set 里可能没有记录。这时候要检查 sep_set 的默认值处理不能直接假设 sep_set[(x, z)] 存在。另外如果条件集为空sep_set 里记录的是空列表判断 y not in sep 的时候空列表会让所有 y 都满足条件这会导致过度定向。规则2 的迭代终止条件也要注意。如果实现不当可能会在两条边之间反复定向形成死循环。加一个最大迭代次数或者用集合记录已处理的边可以避免这个问题。4.5 结果验证不要只看图要看边稳定性PC 算法的输出是一张图但图上的每条边可信度是不一样的。除了 Bootstrap 频率还可以做以下验证第一用不同的 α 值跑多次看哪些边在所有 α 下都稳定存在。第二用不同的条件独立检验方法跑看结果是否一致。第三如果有领域知识检查输出的边是否符合已知的因果关系。如果一条边在数据上显著但领域上不可能大概率是混杂因素没控制好。还有一个实用的技巧把 PC 算法的输出和基于评分的算法比如 GES、NOTEARS的输出做对比。如果两种方法得到的图高度一致可信度就高如果差异很大说明数据本身的信息不足以确定因果结构这时候任何算法的输出都要谨慎对待。5. 让PC算法跑得更稳几个我反复用到的调参习惯第一个习惯是先用小样本快速试跑。拿 100 到 200 个样本、5 到 8 个变量先跑一遍看骨架搜索和定向阶段有没有报错图结构是不是合理。这一步主要是验证代码逻辑不是验证结果。代码没问题了再上全量数据。第二个习惯是固定随机种子。PC 算法本身是确定性的但如果你在预处理阶段用了任何随机方法比如插补、Bootstrap固定种子能保证结果可复现。我一般会在脚本开头写np.random.seed(42)然后所有随机操作都基于这个种子。第三个习惯是记录每次运行的参数和输出。PC 算法的结果对 α 和 max_cond_size 很敏感不记录参数的话过两天回头看图都不知道是怎么跑出来的。我一般会在输出目录里存一个 params.json记录 alpha、max_cond_size、样本量、变量列表以及每条边的 p 值和分离集。第四个习惯是对输出的边做排序。按 p 值从小到大排p 值最小的边最可信。如果时间有限优先验证排在前面的边。这个排序在调试阶段特别有用能快速定位哪些边是噪声。最后一个习惯是不要迷信 PC 算法的输出。PC 算法是探索性工具不是确认性工具。它给出的图是一个假设需要后续用干预实验或者领域知识来验证。我见过太多人把 PC 算法的输出直接当成因果结论写进报告这是很危险的。因果发现的第一步是发现可能的因果结构第二步是验证PC 算法只完成了第一步。希望这些经验能帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/3 3:53:57 · 阅读 1.2k Python从零实现PC算法：因果发现中的条件独立检验与定向规则
简介这份Python实现PC算法的项目源码面向具备一定统计与编程基础、希望深入因果发现与网络结构学习的数据分析者和机器学习学习者。PC算法通过部分相关性检验变量间的条件独立关系可用于高维数据中剔除间接相关、揭示直接因果结构。资源包共13个文件约452KB以4个py源码文件为核心辅以3张png结果示意图、1个csv测试数据集、1个md说明文档及若干配置文件结构紧凑、便于快速上手。项目围绕数据预处理、相关矩阵计算、条件独立测试、定向边剔除与循环迭代等核心步骤展开并给出基于networkx的图结构可视化思路以及面向非高斯数据和大规模数据的扩展优化方向。已有298人学习适合想理解PC算法原理并动手复现因果发现流程的读者参考。1. PC算法到底在算什么从条件独立到因果骨架的那条线PC 算法Peter-Clark Algorithm是因果发现领域最经典的约束型算法之一它要解决的问题很具体手上只有一堆观测数据没有任何先验的因果方向信息能不能把变量之间的因果骨架图给还原出来。很多人第一次听到「Python实现PC算法项目源码」这个标题脑子里浮现的是某个现成的 pip 包装完调个函数就出图。实际情况是PC 算法本身不复杂但它的实现细节里藏着大量统计学上的坑条件独立检验选什么、显著性水平怎么定、样本量够不够每一个都会直接改变最终输出的图结构。这篇文章面向两类人一类是想把 PC 算法真正跑起来、拿到可解释因果图的工程师和数据分析师另一类是想读懂 PC 算法源码、搞清楚每一步在做什么、方便自己改造成业务版本的人。我会从算法骨架讲起然后给出一份可以直接复现的 Python 实现再逐层拆解参数、踩坑点和验证方法。读完你应该能自己写出一份不依赖第三方因果库的 PC 算法代码并且知道在什么数据规模下它还能用、什么时候该换算法。2. PC算法的骨架拆解为什么先做骨架再做定向2.1 从完全图到稀疏骨架条件独立检验在做什么PC 算法的核心思想可以用一句话概括如果两个变量 X 和 Y 在给定某个变量集合 S 的条件下独立那么 X 和 Y 之间就不应该有边。算法从一个完全无向图出发对每一对变量尝试寻找一个条件集使得它们在给定该条件集时条件独立。如果找到了就把这条边删掉。这里的关键在于「条件集从哪来」。PC 算法的做法是对于当前还和 X 相邻的节点按邻接集大小从 0 开始逐层搜索。第 0 层就是看 X 和 Y 的边缘独立性第 1 层看给定一个邻居时是否独立第 2 层看给定两个邻居时是否独立以此类推。这个逐层扩展的过程保证了算法不会漏掉低阶的条件独立关系同时把搜索空间控制在可接受的范围内。条件独立检验通常用偏相关系数或者 Fisher Z 检验。对于线性高斯数据偏相关系数为零等价于条件独立所以用 Fisher Z 变换把偏相关系数转成近似正态统计量再做显著性检验这是最标准的做法。对于离散数据一般用 G² 检验或卡方检验。选哪种检验取决于你的数据类型选错了后面全错。2.2 定向规则把无向骨架变成部分有向图骨架建好之后PC 算法用一组定向规则把无向边变成有向边。最核心的是三条第一条对撞结构。如果 X 和 Z 不相邻但 X→Y←Z 这种结构存在也就是 Y 同时和 X、Z 相邻而 X 和 Z 之间没有边并且 Y 不在 X 和 Z 的条件集里那么 X→Y 和 Z→Y 的方向就确定了。这条规则是 PC 算法能定向的根本原因。第二条避免新对撞。如果已经确定了一条有向边 X→Y而 Y 和 Z 相邻X 和 Z 不相邻那么 Y→Z 的方向可以确定否则会形成一个新的对撞结构。第三条避免环。定向过程中不能产生有向环如果某条边的方向会导致环就反向或者保持无向。这三条规则反复应用直到没有新的边可以被定向。最终输出的图叫 CPDAGCompleted Partially Directed Acyclic Graph里面有些边是有向的有些还是无向的无向边表示在当前数据下无法确定方向。2.3 为什么样本量和检验阈值会直接改变图结构PC 算法对样本量非常敏感。条件独立检验的统计功效随样本量增加而提高样本量不够的时候本来应该被删掉的边删不掉图会偏密样本量太大的时候微弱的依赖关系也会被检验出来图可能偏密。更麻烦的是PC 算法是逐层搜索的一旦某一层删错了一条边后面的定向规则会基于错误的骨架继续推错误会累积。显著性水平 α 的选择也是玄学。α 设大了条件独立被误判为依赖边删不掉α 设小了依赖被误判为独立边被误删。实践中一般从 0.01 到 0.05 之间试但真正靠谱的做法是做敏感性分析看不同 α 下图结构的变化。如果 α 从 0.01 变到 0.05 图结构剧烈变化说明数据量不够或者变量间关系太弱这时候 PC 算法的输出不可信。3. 用Python从零实现PC算法核心代码与参数说明3.1 环境准备与依赖选择实现 PC 算法不需要太重的依赖。核心就是 numpy 做矩阵运算scipy 做统计检验networkx 做图结构管理。不建议一上来就用 causal-learn 或者 pgmpy那些库封装太厚出了问题不好排查。自己写一遍后面调参和改造都方便。pip install numpy scipy networkx pandas matplotlib这几个包都是常规科学计算栈版本没有特别要求numpy 1.20 以上、scipy 1.6 以上就行。pandas 用来读数据matplotlib 用来画图不是必须的但调试的时候很有用。3.2 条件独立检验偏相关系数与Fisher Z变换偏相关系数的计算是 PC 算法的计算核心。给定变量集合 S要算 X 和 Y 在给定 S 下的偏相关系数标准做法是用回归残差的相关性来算。import numpy as np from scipy import stats def partial_corr(x, y, S, data): 计算 x 和 y 在给定 S 下的偏相关系数 x, y: 变量名或列索引 S: 条件集列表形式 data: pandas DataFrame if len(S) 0: r np.corrcoef(data[x], data[y])[0, 1] return r # 用线性回归去掉 S 的影响 from numpy.linalg import lstsq Z data[S].values Z np.column_stack([np.ones(len(Z)), Z]) # 加截距项 # 对 x 和 y 分别回归取残差 coef_x, _, _, _ lstsq(Z, data[x].values, rcondNone) coef_y, _, _, _ lstsq(Z, data[y].values, rcondNone) res_x data[x].values - Z coef_x res_y data[y].values - Z coef_y r np.corrcoef(res_x, res_y)[0, 1] return r def fisher_z_test(x, y, S, data, alpha0.05): Fisher Z 检验返回是否独立 n len(data) r partial_corr(x, y, S, data) # Fisher Z 变换 z 0.5 * np.log((1 r) / (1 - r)) # 标准差 se 1.0 / np.sqrt(n - len(S) - 3) # 检验统计量 stat abs(z) / se # 双尾检验 p_value 2 * (1 - stats.norm.cdf(stat)) return p_value alpha, p_value这段代码里有两个关键点。第一偏相关系数的计算用的是回归残差法这是最直观也最稳定的做法比直接套公式算协方差矩阵的逆要稳。第二Fisher Z 变换的自由度是 n - len(S) - 3这个 3 是固定的来自变换本身的方差近似。如果样本量 n 小于 len(S) 3这个检验就没法做了实践中要保证 n 至少是条件集大小的 5 到 10 倍。参数 alpha 是显著性水平默认 0.05。返回的 p_value 可以用来做敏感性分析看不同阈值下哪些边会被删掉。3.3 骨架搜索逐层条件集与邻接表更新骨架搜索是 PC 算法最耗时的部分。核心逻辑是对每一对相邻节点从空集开始逐步扩大条件集直到找到一组条件使得它们独立或者条件集大小超过当前邻接集。def skeleton_discovery(data, alpha0.05, max_cond_sizeNone): PC 算法骨架搜索 返回无向图邻接表和分离集 variables list(data.columns) n_vars len(variables) # 初始化完全图 adj {v: set(variables) - {v} for v in variables} sep_set {} # 记录分离集用于后续定向 if max_cond_size is None: max_cond_size n_vars - 2 for cond_size in range(max_cond_size 1): # 收集所有需要检验的边 edges_to_check [] for x in variables: for y in adj[x]: if x y: # 避免重复 edges_to_check.append((x, y)) for x, y in edges_to_check: if y not in adj[x]: continue # 边已经被删了 # 候选条件集从 x 的邻居中选排除 y neighbors adj[x] - {y} if len(neighbors) cond_size: continue from itertools import combinations for S in combinations(neighbors, cond_size): independent, p_val fisher_z_test(x, y, list(S), data, alpha) if independent: # 删除边 adj[x].discard(y) adj[y].discard(x) sep_set[(x, y)] list(S) sep_set[(y, x)] list(S) break return adj, sep_set这段代码有几个工程上的细节值得说。第一条件集是从 x 的邻居里选的不是从所有变量里选这是 PC 算法的标准做法能大幅减少检验次数。第二用 combinations 生成条件集当邻居数量多的时候组合数会爆炸所以实践中要限制 max_cond_size一般不超过 3 到 4。第三sep_set 记录了每对变量是在哪个条件集下被判定独立的这个信息在定向阶段要用。参数 max_cond_size 控制搜索深度。设太小会漏掉高阶条件独立关系图偏密设太大会导致计算量指数增长。经验值是变量数的三分之一到一半但不超过 5。3.4 定向规则实现对撞结构与避免新对撞骨架建好之后定向阶段要把无向边变成有向边。核心是识别对撞结构然后传播方向。def orient_edges(adj, sep_set, variables): 定向阶段识别对撞结构并传播方向 directed set() # 有向边 (x, y) 表示 x - y undirected set() # 收集所有无向边 for x in variables: for y in adj[x]: if x y: undirected.add((x, y)) # 规则1识别对撞结构 X - Y - Z for y in variables: neighbors list(adj[y]) for i in range(len(neighbors)): for j in range(i 1, len(neighbors)): x, z neighbors[i], neighbors[j] if z in adj[x]: continue # x 和 z 相邻不是对撞 # 检查 y 是否在 sep_set[(x, z)] 中 sep sep_set.get((x, z), []) if y not in sep: # 对撞结构成立 directed.add((x, y)) directed.add((z, y)) undirected.discard((min(x, y), max(x, y))) undirected.discard((min(z, y), max(z, y))) # 规则2避免新对撞 changed True while changed: changed False for x, y in list(undirected): # 如果 x - y 已经确定检查 y 的其他邻居 if (x, y) in directed: for z in adj[y]: if z x: continue if (min(y, z), max(y, z)) in undirected: if z not in adj[x]: # y - z 方向确定 directed.add((y, z)) undirected.discard((min(y, z), max(y, z))) changed True return directed, undirected定向规则里最容易出错的是对撞结构的判断条件。必须同时满足三个条件X 和 Z 不相邻、Y 同时和 X 和 Z 相邻、Y 不在 X 和 Z 的分离集中。第三个条件最容易被忽略如果 Y 在分离集中说明 X 和 Z 的独立性是 Y 导致的这时候不能定向为对撞。规则2 的实现是一个迭代过程因为定向一条边可能会触发新的定向。循环直到没有新的边可以被定向为止。3.5 完整流程串联与输出解读把上面的模块串起来就是一个完整的 PC 算法实现。def pc_algorithm(data, alpha0.05, max_cond_sizeNone): 完整的 PC 算法 variables list(data.columns) # 第一步骨架搜索 adj, sep_set skeleton_discovery(data, alpha, max_cond_size) # 第二步定向 directed, undirected orient_edges(adj, sep_set, variables) return { adjacency: adj, sep_set: sep_set, directed: directed, undirected: undirected } # 使用示例 import pandas as pd import numpy as np # 生成模拟数据X - Y - Z, X - Z np.random.seed(42) n 1000 X np.random.randn(n) Y 0.8 * X np.random.randn(n) * 0.5 Z 0.6 * X 0.7 * Y np.random.randn(n) * 0.5 data pd.DataFrame({X: X, Y: Y, Z: Z}) result pc_algorithm(data, alpha0.01) print(有向边, result[directed]) print(无向边, result[undirected])输出解读的时候要注意PC 算法输出的是 CPDAG不是唯一的因果图。有向边表示在所有马尔可夫等价类中方向一致无向边表示方向不确定。如果你看到 X→Y 和 Y→Z 都是有向的但 X 和 Z 之间没有边这不一定意味着 X 和 Z 独立可能只是条件独立检验没找到它们之间的直接依赖。4. PC算法落地时的避坑清单从数据预处理到结果验证4.1 数据预处理没做好后面全白搭PC 算法对数据的假设是连续变量、线性关系、高斯噪声。如果你的数据里有分类变量直接扔进去算偏相关系数会得到完全错误的结果。分类变量要么先做独热编码然后当连续变量处理效果一般要么换用基于互信息的条件独立检验。缺失值也是大问题PC 算法没有内置的缺失值处理机制要么删样本要么插补插补方法的选择会直接影响条件独立检验的结果。另一个容易被忽略的是变量尺度。偏相关系数本身对尺度不敏感但如果你在预处理阶段做了标准化要注意标准化是在全量数据上做的还是分训练测试集做的。PC 算法一般不需要划分训练测试集但如果你要做交叉验证来选 α标准化必须在每折内部独立做否则会信息泄露。4.2 条件独立检验选错类型图结构完全变样连续高斯数据用 Fisher Z 检验离散数据用 G² 检验混合数据用条件互信息或者基于核的方法。选错了检验类型不是精度下降的问题是根本性错误。比如离散数据用 Fisher Z偏相关系数根本没有意义算出来的 p 值也是错的。还有一个隐蔽的坑Fisher Z 检验假设变量间是线性关系。如果真实关系是非线性的比如 Y X²偏相关系数可能接近零检验会判定独立但实际上 X 和 Y 有强依赖。这种情况下需要换用基于互信息或距离相关的方法但那些方法的计算量会大很多。4.3 样本量不够时PC算法输出的图不可信PC 算法的最低样本量要求没有严格公式但经验上每个变量至少需要 10 到 20 个样本条件集大小每增加 1样本量需求大概翻倍。如果你有 20 个变量、200 个样本跑 PC 算法大概率会得到一张乱七八糟的图。判断样本量够不够的一个实用方法是跑多次 Bootstrap看边出现的频率。如果某条边在 80% 以上的 Bootstrap 样本中都出现可以认为它是稳定的如果只有 50% 左右说明这条边不可靠。这个做法比单纯看 p 值要靠谱得多。4.4 定向规则实现中的边界情况对撞结构识别的时候如果 X 和 Z 之间本来有边但后来被删了sep_set 里可能没有记录。这时候要检查 sep_set 的默认值处理不能直接假设 sep_set[(x, z)] 存在。另外如果条件集为空sep_set 里记录的是空列表判断 y not in sep 的时候空列表会让所有 y 都满足条件这会导致过度定向。规则2 的迭代终止条件也要注意。如果实现不当可能会在两条边之间反复定向形成死循环。加一个最大迭代次数或者用集合记录已处理的边可以避免这个问题。4.5 结果验证不要只看图要看边稳定性PC 算法的输出是一张图但图上的每条边可信度是不一样的。除了 Bootstrap 频率还可以做以下验证第一用不同的 α 值跑多次看哪些边在所有 α 下都稳定存在。第二用不同的条件独立检验方法跑看结果是否一致。第三如果有领域知识检查输出的边是否符合已知的因果关系。如果一条边在数据上显著但领域上不可能大概率是混杂因素没控制好。还有一个实用的技巧把 PC 算法的输出和基于评分的算法比如 GES、NOTEARS的输出做对比。如果两种方法得到的图高度一致可信度就高如果差异很大说明数据本身的信息不足以确定因果结构这时候任何算法的输出都要谨慎对待。5. 让PC算法跑得更稳几个我反复用到的调参习惯第一个习惯是先用小样本快速试跑。拿 100 到 200 个样本、5 到 8 个变量先跑一遍看骨架搜索和定向阶段有没有报错图结构是不是合理。这一步主要是验证代码逻辑不是验证结果。代码没问题了再上全量数据。第二个习惯是固定随机种子。PC 算法本身是确定性的但如果你在预处理阶段用了任何随机方法比如插补、Bootstrap固定种子能保证结果可复现。我一般会在脚本开头写np.random.seed(42)然后所有随机操作都基于这个种子。第三个习惯是记录每次运行的参数和输出。PC 算法的结果对 α 和 max_cond_size 很敏感不记录参数的话过两天回头看图都不知道是怎么跑出来的。我一般会在输出目录里存一个 params.json记录 alpha、max_cond_size、样本量、变量列表以及每条边的 p 值和分离集。第四个习惯是对输出的边做排序。按 p 值从小到大排p 值最小的边最可信。如果时间有限优先验证排在前面的边。这个排序在调试阶段特别有用能快速定位哪些边是噪声。最后一个习惯是不要迷信 PC 算法的输出。PC 算法是探索性工具不是确认性工具。它给出的图是一个假设需要后续用干预实验或者领域知识来验证。我见过太多人把 PC 算法的输出直接当成因果结论写进报告这是很危险的。因果发现的第一步是发现可能的因果结构第二步是验证PC 算法只完成了第一步。希望这些经验能帮到你。本文还有配套的精品资源点击获取