简介：这份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 的方向可以确定，否则会形成一个新的对撞结构。

第三条，避免环。定向过程中不能产生有向环，如果某条边的方向会导致环，就反向或者保持无向。

这三条规则反复应用直到没有新的边可以被定向。最终输出的图叫 CPDAG（Completed 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, rcond=None) coef_y, _, _, _ = lstsq(Z, data[y].values, rcond=None) 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, alpha=0.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, alpha=0.05, max_cond_size=None): """ 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, alpha=0.05, max_cond_size=None): """ 完整的 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, alpha=0.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 算法只完成了第一步。希望这些经验能帮到你。

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