简介一份基于 Python 的 PC 算法实现源码面向从事因果发现、网络结构学习或高维数据条件独立性分析的研究者与开发者。项目围绕偏相关这一核心概念实现从数据预处理、相关矩阵计算、条件独立检验到边剔除与结果可视化的完整流程适配典型的表格式数据便于快速复现并理解基础版本。压缩包共 13 个文件以 Python 脚本为主辅以 CSV 样本数据、实验结果示意图、依赖清单和说明文档整体仅 452KB结构轻量适合直接阅读与二次开发。目前已有 298 人学习下载。通过阅读核心源码和注释可掌握条件独立检验、显著性阈值设置和有向图输出等关键细节还可以了解如何借助 networkx 绘制因果网络项目中还给出了向非高斯数据、并行化计算等方向扩展优化的思路便于进一步改造和落地。1. PC算法在因果发现里处在什么位置为什么你最终会回来手写它你手里有一批观测数据变量十几个想从里面找因果关系而不是相关关系——搜索“Python实现PC算法项目源码”的人大多走到了这一步传统回归已经在业务上解释不了“为什么两个指标总是同涨同跌”而实验又做不了。PC算法Peter-Clark算法就是在这种情况下被反复提起的名字。它是从观测数据推断有向无环图最经典的入门算法绝大多数因果发现项目都以它为起点。但很多人第一步就踩进黑匣子网上教程全是“调用现成库跑通”遇到输出异常时既不知道调alpha还是换检验方式更判断不了结果是真结构还是假设被违背的假象。这篇我按自己实际落地的方案讲先把PC算法的三个假设拆开再给出最小可运行的手写代码最后讲参数参考和五个血泪坑。适合正在用Python做因果发现的从业者尤其是数据来自观测、没有受控实验可做的人。相信我你迟早会从“调库跑通”回到“手写一遍”的。2. 先立住理论地基PC算法的三个假设与两阶段流程2.1 三个假设忠实性、因果充分性与正态性少一个结果都不可信PC算法输出一张有向无环图的勇气完全建立在三个假设上。第一个是因果充分性causal sufficiency所有变量都被观测到不存在未观测的共同原因。这在现实里是最难满足的。比如你观测了点击量和订单量却没记下促销活动这个上游变量两个变量会被一个隐藏共因推着走PC算法会把虚假的边误认为直接因果。数据侧的信号是回归残差之间仍明显相关或者把数据按时间分段后PC结果剧烈变化。第二个是忠实性faithfulness真实图中存在的每条因果路径都会在概率上留下可观测的条件依赖反过来数据里出现的独立关系必须来自图结构本身不是概率巧合。这句话很绕但破坏它的例子很常见——抵消效应X→Y→Z 和 X→Z 两条路同时存在如果两条路径的系数恰好大小相等符号相反X和Z在给定Y时仍然独立PC算法就把那条真实的边X→Z剪掉了。忠实性假设不允许这种“数值巧合”存在但真实数据天天发生。第三个是正态性或至少可交换性偏相关检验和Fisher z变换的分布假设建立在数据近似高斯上这是一条实用性假设。遇到厚尾、多模态或者离散变量偏相关检验会成片误判。离散数据要换成卡方或互信息检验这是我强烈建议一开始就决定的——它直接影响你待会儿代码里调用哪个函数。三个假设不是学术摆设你拿着PC算法的结果给业务方看的时候能说得清“我们的结论在什么前提下成立”才不会被一句“是不是两个变量碰巧相关”问倒。2.2 偏相关检验PC算法的条件独立引擎PC算法判断“X和Y是否直接相连”不靠普通相关靠条件独立检验。给定一组变量S如果X和Y在控制S后仍然相关说明S没有完全解释二者之间的依赖X和Y之间应该有一条边至少在骨架层保留。判断方法在高斯数据下就是偏相关。偏相关用精度矩阵算非常简洁协方差矩阵的逆precision matrix归一化后就是偏相关矩阵。具体公式是ρ_xy|S -prec_xy / sqrt(prec_xx × prec_yy)当S为空集时这个值退化成普通相关系数。用精度矩阵的好处是一次求逆可以反复取任意条件集的偏相关不用每次都重新回归坏处是当样本量不足或者变量高度共线时协方差矩阵接近奇异求逆就是第一步翻车点。显著性检验用Fisher z变换z 0.5 × ln((1r)/(1-r))再乘以 sqrt(n − |S| − 3)近似标准正态p值小于alpha就认为不独立。这套检验有个隐含要求样本量至少比条件集大小多4否则方差估计是负的。后面代码里我会对条件集大小做上限保护。2.3 两阶段流程骨架学习与边的定向PC算法分两步。第一步叫骨架学习skeleton从完全无向图开始对每一对相邻节点用它们的共同邻居做候选条件集从条件集大小为0开始逐层增大一旦找到某个条件集使两个节点条件独立就把这条边删掉并把该条件集记入d分离集d-sep set。关键点是删掉一条边后邻居集合立刻变化后续检验必须基于更新后的图——这也是它在效率和精度上区别于朴素条件独立检验的地方。第二步叫定向orientation。先找v结构对任意链 i−j−k如果i和k不直接相邻且j不在分离i和k的d分离集中那么j是碰撞器collider方向为 i→j←k。紧接着用Meek规则把方向向整个无向骨架传播。做完这些后图中仍会留下一些无法确定方向的边PC算法输出的其实是CPDAG完全部分有向无环图代表一个等价类而不是唯一的一张DAG。这个认知特别重要——很多人在复现后看到“还有边没定向”就以为代码错了其实这才是PC算法的真实输出。明确这一点再动手后面每一步都不会慌。3. 从零手写PC算法三段代码与参数详解3.1 数据准备与偏相关检验先算对独立性我习惯先把数据标准化再走精度矩阵路线算偏相关。标准化不是必须的但数值稳定性好很多尤其在变量量纲差异大时。下面这段是偏相关矩阵和Fisher z检验的最小实现我常用它替代任何统计库。import numpy as np import pandas as pd def partial_corr_matrix(data: np.ndarray) - np.ndarray: 计算偏相关矩阵。 data: shape (n_samples, n_features)行是样本列是变量。 返回偏相关矩阵对角线为1。 # 标准化到零均值单位方差数值更稳 data (data - data.mean(axis0)) / data.std(axis0, ddof1) corr np.corrcoef(data, rowvarFalse) # 用伪逆而不是逆协方差接近奇异时伪逆不会直接崩 prec np.linalg.pinv(corr) diag np.sqrt(np.abs(np.diag(prec))) pcorr -prec / np.outer(diag, diag) np.fill_diagonal(pcorr, 1.0) return pcorr def fisher_z_test(r: float, n: int, cond_size: int, alpha: float 0.05) - bool: 检验偏相关 r 在给定 n 个样本、条件集大小 cond_size 下是否显著。 返回 True 表示“不独立”拒绝原假设。 z 0.5 * np.log((1 r) / max(1 - r, 1e-12)) stat np.sqrt(max(n - cond_size - 3, 1)) * z # |stat| 大于标准正态分布 alpha/2 分位数则拒绝独立 return abs(stat) 1.96 # alpha0.05 的分位数可换成 scipy.stats.norm.ppf逻辑说明我先把数据标准化保证协方差矩阵的量纲统一然后求相关矩阵的伪逆pinv它和真实逆在满秩时等价在秩亏时给出最小二乘意义下的解不会抛奇异矩阵异常。偏相关由精度矩阵归一化得出等你后面看其他实现时会发现也有人用“协方差矩阵求逆后除以对角元勾股”的写法本质一样。Fisher z 检验里我给了一个 1.96 的硬编码这是 alpha0.05 的常用分位数如果你需要换显著性水平建议换成 scipy.stats.norm.ppf(1 - alpha / 2)别到处改魔法数字。在 PC 算法里我很少直接用整个偏相关矩阵而是对每一对节点调用一个“给定条件集后计算偏相关”的函数。上面的矩阵版用来快速看数据概况逐对计算版则嵌入骨架学习循环中。你写代码时建议先跑一遍矩阵版确认数据里高分列有没有异常——如果偏相关矩阵对角线外全接近 1 或全接近 0先回去检查数据清洗别急着跑算法。3.2 骨架学习用邻接矩阵做逐层剪边骨架学习的核心是维护一个邻接矩阵 G矩阵里 1 代表当前还有边。对每一对相邻节点在所有共同邻居中暴力枚举条件集。代码我按最常见的实现方式写优先保证可读性def pc_skeleton(data: np.ndarray, alpha: float 0.05, max_cond_size: int 5): n_features data.shape[1] n_samples data.shape[0] # 初始化为完全图对角线为0自己连自己没有意义 G np.ones((n_features, n_features)) np.fill_diagonal(G, 0) # dsep[i][j] 记录使 i,j 条件独立的那个条件集列表 dsep {} def neighbors(v): return {u for u in range(n_features) if G[v][u] 1 and u ! v} for l in range(max_cond_size 1): changed False for i in range(n_features): for j in range(n_features): if G[i][j] 0 or i j: continue nbrs neighbors(i) - {j} if len(nbrs) l: continue # 在邻居中取大小为 l 的组合做候选条件集 for cond in itertools.combinations(nbrs, l): r partial_corr_between(data, i, j, list(cond)) if not fisher_z_test(r, n_samples, len(cond), alpha): G[i][j] G[j][i] 0 dsep[(i, j)] list(cond) dsep[(j, i)] list(cond) changed True break # 找到一组独立就剪边不必枚举完 if not changed: break # 这一层没有边被剪掉可以提前退出 return G, dsep逻辑说明最外层循环控制条件集大小 l从 0 开始逐步加深。内层对每对仍然相邻的节点做检验候选条件集只从 i 的邻居中取并且要排除 j 本身这是 PC 算法的标准做法——条件集里不需要包含它想解释的那对节点。注意我把断边后的邻居集合实时更新所以每层都在“当前骨架”上做检验不是原始完全图。dsep字典在定向阶段是命脉务必保留。参数说明alpha是独立性检验的显著性水平调大比如 0.1会让检验更容易判定为不独立结果图更稠密调小比如 0.01则剪边更狠图更稀疏。max_cond_size是条件集大小的上限PC 算法理论上是穷尽所有层的但暴力枚举组合的开销是组合数级的实际数据里条件集超过 4~5 就会很慢而且样本量也撑不住。我一般设 5数据集变量少且样本充足时可以提到 6 或 7。itertools.combinations是最直观的枚举方式但如果邻居很多建议换成按“当前层 l 生成组合”的迭代器写法否则内存会先爆。3.3 定向碰撞器识别与方向规则传播骨架学习把图剪成无向图后定向分两步。先找所有 v 结构再用 Meek 规则把方向传播出去。v 结构是整篇 PC 算法里唯一能从数据里直接“看出”方向的局部结构其余方向靠不产生环和碰撞器的传播规则推导。def orient_graph(G, dsep, alpha, data): n G.shape[0] # 先复制骨架 dir_graph G.copy().astype(int) # 1 表示无向边2 表示有向边用方向信息存储 # 用一个 dict 记录已定向的边: (from, to) - True oriented {} # 第一步v结构识别。对任意 i-j-ki 与 k 不相邻且 j 不在 dsep(i,k) 中 for i in range(n): for j in range(n): if G[i][j] 0 or i j: continue for k in range(n): if k i or k j: continue if G[i][k] 0 or G[j][k] 0: continue if G[i][k] 1 and k in dsep.get((i, k), []): continue # j 不是分离 i,k 的条件则 j 是碰撞器 if j not in dsep.get((i, k), []): oriented[(i, j)] True oriented[(k, j)] True # 第二步Meek规则传播这里实现最简单的一条 # 如果 i - j 且 j - k无向且 i 与 k 不相邻则定向 j - k changed True while changed: changed False for j in range(n): for i in range(n): if (i, j) in oriented: for k in range(n): if G[j][k] 0 or j k: continue if (j, k) in oriented or (k, j) in oriented: continue if G[i][k] 0: # i 与 k 不相邻 oriented[(j, k)] True changed True return oriented逻辑说明v 结构识别的核心不是看图而是查 dsep 集。条件 i 与 k 不相邻然后问 j 是否在分离 i 和 k 的条件集里。如果 j 不在说明 j 不能被其他变量解释成中间变量最合理的解读是 i 和 k 都指向 j也就是碰撞器。这里最容易错的地方是把“不相邻”误解成“没有路径”实际是骨架层已经删掉了直接边才轮到 v 结构出场。Meek 规则这一段只是四条规则中的第一条完整实现要处理环和更多方向传播情况我在小图够用图大了建议直接看 causal-learn 源码里的meek_rules实现。定向结果里一定有边保持无向这是 CPDAG 的正常样子别强加规则把所有边都定向——那反而会违背 PC 的保守前提。在代码里我刻意返回oriented字典而不是直接改dir_graph这样方便复现时检查每一步的定向逻辑。你会看到有一类输出极其常见骨架正确v 结构正确但方向传播完之后局部出现“有向环”。这不是 PC 理论错了多半是 Meek 规则实现漏了一条通常是规则 2 和规则 3涉及碰撞器保护和路径检查。如果你也遇到去对照规则全集。4. 用现成库快速落地causal-learn 的最小复现路径4.1 为什么有现成库还要自己手写一遍手写实现的代码价值在于当你拿着一个结果去和别人讨论时你能说出“这条边为什么留下了”“那个 v 结构为什么被定向”而不是“库跑出来就是这样”。现成库的价值则是交付和生产尤其 causal-learn 这类库把 PC 算法在性能上做了大量优化包括本地化的条件集剪枝、优先级的节点选择策略这些我手写的版本都没有。我的习惯是手写版本跑一遍理解行为库跑一遍确认结果没差如果手写和库结果明显不同那问题几乎一定出在偏相关检验的实现细节上比如是用了相关矩阵还是协方差矩阵求逆、条件集排除方式是不是把目标节点算进了邻居、alpha 被重复使用等。4.2 最小再现代码与参数配置causal-learn 是目前 Python 生态里最接近 PC 算法论文的库我用它做验证已经两年了。最小代码如下跑通后再逐步调参# 安装pip install causal-learn import pandas as pd from causallearn.search.ConstraintBased.PC import pc from causallearn.utils.GraphUtils import GraphUtils df pd.read_csv(your_data.csv) # 列名建议给中文加拼音库对列名没有要求 data df.values # 核心调用alpha 是最关键的参数 cg pc(data, alpha0.05, indep_testfisherz, depth3) # 结果可视化保存为图片容易复盘 GraphUtils.plot(cg.G, pc_result.png)参数说明alpha和手写版一致控制独立性检验的阈值indep_test是独立性检验方式连续数据用fisherz离散数据可以换chisq或gsq这与你最早选定的统计检验对应depth相当于我手写版里的max_cond_size它传 3 表示条件集最大大小 3传 -1 表示无限制但样本量不足时无限制会让运行时长涨得离谱。pc()返回的cg.G就是学习到的图对象cg.G.graph里存邻接矩阵方向箭头在cg.G.directed_edges()里可以拿出来做业务展示。注意不同版本的 causal-learn 对pc()的函数签名的写法有细微差异我建议跑通之后再print(cg.G.graph)看一眼矩阵核对是不是 n×n 的 0/1/2 矩阵0 无边1 无向2 有向。我见过直接把cg.G存下来当 pandas DataFrame 用的人结果列名和顺序全乱了可视化也错位。5. PC算法落地避坑5个我实际踩过的数据与参数问题5.1 数据侧的三连坑矩阵奇异、类别错配、样本量虚胖坑1协方差矩阵奇异导致偏相关全是假的。现象代码在求逆那一步直接报 LinAlgError或者没报错但偏相关矩阵里有大量 NaN。原因变量数量接近甚至超过样本量或者有两个变量完全线性相关比如百分比和它对应的小数。解决先np.linalg.matrix_rank查秩删掉冗余列再把inv换成pinv最后实在不行降低max_cond_size让条件集更小减少对高维精度矩阵的需求。坑2离散数据硬上 fisherz 检验。现象同一份数据用fisherz跑出 8 条边换成chisq跑出 4 条而且业务经验明显支持后者。原因偏相关的高斯假设在离散变量上系统性失效尤其当变量是 0/1 编码、取值稀疏时独立性检验的方差被低估。解决先看数据类型的分布离散列超过 2 个就换chisq或gsq如果数据里混装连续和离散最好做分箱后的互信息检验而不是在一个模型里用两种检验。坑3样本量“看着多”但条件集撑不住。现象骨架阶段 l 增到 3 就开始剧烈变慢增到 4 直接卡死中途中断后重新跑结果不一样。原因样本量虽然上万但某个变量的邻居数很多组合枚举量在 l3 时就已经超过百万规模另外 PC 算法的顺序依懒性导致中断后状态不一致。解决把max_cond_size设成 3 跑通基线再用random.sample做多次 bootstrap 重采样看骨架稳定性而不是一次跑完再调参。5.2 算法侧的翻车点剪边顺序、定向循环与“确定性幻觉”坑4剪边顺序和条件集枚举顺序影响了最终结果。现象我把itertools.combinations换成按变量索引顺序枚举后结果从 5 条边变成 6 条另一条边时有时无。原因PC 算法在理论上顺序无关但实际偏相关估计和检验有数值误差先剪哪条边会改变后续邻居集合放大这个误差。解决手写版接受“顺序略有影响”是正常的要稳定结果就固定随机种子、固定变量名顺序并在报告中注明采用“按变量顺序枚举”的策略真对比方法优劣时用因果结构恢复的 Hamming 距离不要比单次出行。坑5定向后输出里出现了有向环。现象oriented字典里同时存在 (A,B) 和 (B,A)可视化时直观就是两根箭头和 DAG 前提直接冲突。原因我在 v 结构定向后Meek 规则没有做碰撞器保护规则 2 和规则 3 的本职就是保证“不引入环、不覆盖已有碰撞器”我只实现了规则 1 就会漏。解决要么补齐四条 Meek 规则推荐看 Koller 教材里的表述要么在传播循环后加一个校验函数检测环然后直接删除那两条冲突边。别装作没看见——业务方看到环时会立刻对算法失去信任。最后还有一个心态上的坑PC 算法输出里留有大量无向边不是 bug是数据信息量不足的真实反映。很多人看到 10 条边有 4 条没定向就开始堆规则强行定向这是制造“确定性幻觉”。正确做法是把无向边原样交付在报告里标注“此处需要实验或先验知识干预”。这个习惯帮你至少省掉三次业务方追问。6. 最后一步用稳定性检验与可视化验收 PC 算法结果PC 算法跑完只是起点交付前我至少做两件事稳定性检验和方向可信度标注。稳定性检验用 bootstrap 最简单——对原始数据有放回地抽样 100 次每次只重算骨架不跑定向统计每对变量之间出现边的频率输出边的频次矩阵。频率低于 0.6 的边我会在最终图里画成虚线并在说明里写清楚“这条边只在 60% 的样本里出现置信度偏弱”。这比单次跑一版图后使劲解释要专业得多。方向可信度标注同样重要。对 v 结构定向出来的边我在图里用实心箭头对 Meek 规则传播出来的边我用空心箭头代表“结构推断非直接检验”。两者混在一张图里业务方一眼能看到哪些方向是硬证据哪些是软推断。我还常用一个小技巧把数据里最强的共因变量先手动拿出来单独跑一次 PC再把它的预测值作为回归残差灌回去跑第二次两次结果对比能明显暴露因果充分性假设被破坏的位置。最近一个项目里我靠这个技巧成功找到一个隐藏的批次效应变量——原始数据跑出来 A 与 B 高度相关但残差版跑出来两者完全独立。教训是PC 算法不是“喂数据出真相”的机器它是“在假设前提下的结构探测器”。你越早承认它的边界越能把它用出价值。希望帮到你至少少走我当年翻车的那几步。本文还有配套的精品资源点击获取