简介一套基于Python实现PC算法的项目源码面向需要探索变量间条件独立关系开展因果发现与网络结构学习的数据分析或机器学习人员。资源将PC算法原理转化为可运行的工程代码包含数据预处理、相关矩阵计算、条件独立检验、定向边剔除等核心环节并提供测试数据与可视化支持。包体共13个文件涵盖Python源码算法实现与运行脚本、CSV测试数据、PNG结果示意图、README说明文档及若干工程配置文件总体积约452KB目前已有298人学习下载。通过阅读项目中的核心模块和注释可快速掌握PC算法落地思路基于清晰的目录结构能方便地替换自有数据集或结合networkx、matplotlib等库进一步优化和二次开发。项目还提供扩展方向参考适合用作因果学习相关课程设计或论文实验的基础框架。1. PC算法在解决什么问题从观测数据反推因果骨架给你一张没有时间戳、不加干预的观测数据表你有多大把握说清楚“是A推动了B还是B推动了A”营销增量测算、推荐归因和医疗数据分析里每天都在问这类问题而PC算法是回答它们的经典起点。PC算法全称Peter-Clark算法由Peter Spirtes和Clark Glymour在1993年提出注意不是主成分分析PCA。流程分两步先把所有变量放进一张完全无向图用条件独立性检验逐条剪掉不相关的边再用V结构和方向规则把剩余边的方向定出来。整个过程不需要干预实验不要求节点顺序输入只有数据矩阵和一个显著性水平alpha。下面从理论、Python源码实现、参数整定到验证排错逐个讲清。2. 理论立住再写码条件独立性、V结构与定向规则2.1 骨架学习完全图靠条件独立性剪枝PC算法第一步输出无向骨架。初始化时p个变量构成完全图任意两点都有边。算法对每条边X-Y反复提问在给定变量集Z的条件下X和Y仍然相关吗如果检验结论是条件独立就删除这条边并把Z记录为分离集。对于近似多元高斯数据条件独立等价于偏相关系数为0。给定条件集Z时偏相关系数可以直接从精度矩阵读出以X、Y和Z的联合相关矩阵Σ求逆得到PX与Y的偏相关系数为 -P_ij / sqrt(P_ii·P_jj)。把偏相关系数做Fisher Z变换z arctanh(r)检验统计量 sqrt(n - |Z| - 3) · |z|统计量近似标准正态双尾p值大于alpha就接受独立、删边。条件集深度必须从0开始递增因为分离集要尽可能小同一对变量可能被多个集合分开先取最小集合V结构判定才可靠。如果从大条件集开始记录到的分离集可能不是最小的方向判定会出错。2.2 方向判定先找V结构再用Meek规则收尾骨架只说明变量是否相邻。PC算法用两类信息定向。第一类是V结构若三元组X-Z-Y中X与Y不相邻且Z不在sep_set(X,Y)里说明X和Y在给定Z的条件下反而相关只有X→Z←Y这种碰撞子结构可以解释于是把两条边向外锁定。找到全部V结构后仍有未定向边。第二类信息是Meek规则它保证在既不产生新V结构也不产生环的前提下尽量定向剩余边。最常用的三条R1若a→b且b—c相连、a与c不相邻则定向为b→cR2若a→b→c且a—c相连则定向为a→cR3若a—c→b且a—d→b且a与b不相邻、c与d不相邻则定向为a→b规则循环应用直到没有新定向边。完整PC实现还会加一条R4它对理论完全性有意义对骨架和绝大多数边方向影响不大这里的源码实现保留这三条。2.3 为什么不直接用相关矩阵或回归PC算法全程都在做“相关性检测”那为什么不直接跑互信息或回归差别在方向。互信息是无向的回归系数虽然带方向但依赖你预先指定的变量顺序而这个顺序正是因果发现要回答的问题。PC算法能整图输出因果方向在忠实性、充分样本和无隐变量假设下理论上有可证明的正确性条件。causal-learn、gCastle等库的PC都沿这套流程自己实现一遍是把alpha和max_depth如何影响结果看懂的最直接方式而不是把包当黑盒。三者对比看下表方法输出方向依赖父变量顺序整图输出互信息无向不依赖不能回归有向强依赖不能PC算法部分或全部不依赖可以3. Python源码实现骨架搜索、V结构定向与主流程3.1 数据结构与独立性检验先明确三个对象graph是n×n布尔矩阵True表示骨架保留边dir_matrix是n×n整数矩阵1表示i→j-1表示j→i0表示未定向sep_set是字典键为(i,j)值为条件变量列表存双向方便V结构查询。独立性检验用Fisher Z变换这里是完整函数import itertools import numpy as np from scipy import stats def fisher_z_test(data, x, y, cond_set): 偏相关Fisher Z检验返回双尾p值。 p alpha 时认为x与y在给定cond_set下条件独立。 n data.shape[0] # 收集本次检验涉及的全部变量列号并按序去重 var_idx list(dict.fromkeys([x, y] [c for c in cond_set])) if len(var_idx) n: # 样本数不够估计相关矩阵保守返回不独立 return 0.0 sub data[:, var_idx] corr np.corrcoef(sub, rowvarFalse) try: prec np.linalg.inv(corr) except np.linalg.LinAlgError: # 变量共线导致矩阵奇异保守保留边 return 0.0 xi var_idx.index(x) yi var_idx.index(y) pcorr -prec[xi, yi] / np.sqrt(prec[xi, xi] * prec[yi, yi]) pcorr np.clip(pcorr, -0.999, 0.999) # arctanh即为Fisher Z变换 z_val np.arctanh(pcorr) stat np.sqrt(n - len(var_idx) - 3) * abs(z_val) return 2 * (1 - stats.norm.cdf(stat))三个关键决策解释一下。用精度矩阵求偏相关而不是按条件集大小展开递推式好处是cond_set任意长度都只需要一次矩阵求逆。var_idx用dict.fromkeys去重并保持顺序避免同一列重复进入相关矩阵。奇异矩阵返回0.0意味着“不独立”骨架不会因为一次数值问题被错误删边代价是可能多保留几条假边属于保守策略。3.2 骨架搜索主循环骨架循环按条件集深度从0递增对每条现存边枚举候选条件集组合命中一个分离集就删除def pc_skeleton(data, alpha0.05, max_depthNone): PC算法第一步从完全图剪到无向骨架。 参数 ---------- data : np.ndarray, shape (n_samples, n_vars) alpha : float, 显著性水平 max_depth : int | None, 最大条件集大小, 默认n_vars-2 返回 ---------- graph : np.ndarray, 骨架邻接矩阵 sep_set : dict, 分离集字典 n_vars data.shape[1] graph np.ones((n_vars, n_vars), dtypebool) np.fill_diagonal(graph, False) sep_set {} if max_depth is None: max_depth n_vars - 2 depth 0 while depth max_depth: changed False for i in range(n_vars): for j in range(i 1, n_vars): if not graph[i, j]: continue # 当前邻居快照里去掉j作为候选条件集 nbrs_i [k for k in range(n_vars) if graph[i, k] and k ! j] if len(nbrs_i) depth: continue for cond in itertools.combinations(nbrs_i, depth): p fisher_z_test(data, i, j, list(cond)) if p alpha: graph[i, j] graph[j, i] False sep_set[(i, j)] list(cond) sep_set[(j, i)] list(cond) changed True break # 找到一组分离条件就删 if not changed: break depth 1 return graph, sep_set易错点在于nbrs_i必须在每对变量处理时重新取。如果在外层把邻居快照保存删边后的循环会继续用到旧邻居既增加无用组合又可能把条件集扩大导致本可删的边没有删。标准PC在“每个depth下逐条删边”就是这种即时更新语义。提示sep_set会同时存(i,j)和(j,i)后续V结构查询更安全。查找时用sep_set.get((i, j), [])避免KeyError误判成“不在分离集”而制造出假的V结构。3.3 V结构定向识别V结构扫描的是不相邻的变量对找它们的共同邻居def orient_v_structures(graph, sep_set): 识别V结构并返回方向矩阵。 dir_matrix[i][j]1 表示 i-j-1 表示 j-i0 表示未定向。 n graph.shape[0] dir_matrix np.zeros((n, n), dtypeint) for i in range(n): for j in range(i 1, n): if graph[i, j]: continue # 只处理不相邻的变量对 # 找共同邻居 for k in range(n): if not (graph[i, k] and graph[j, k]): continue if k not in sep_set.get((i, j), []): dir_matrix[i, k] 1 # i - k dir_matrix[k, i] -1 dir_matrix[j, k] 1 # j - k dir_matrix[k, j] -1 return dir_matrix判断依据是k不在i和j的分离集里意味着i与j在给定k后变得相关k只能是碰撞子。矩阵的对称写法要成对维护后面Meek规则判断dir时依赖两条记录同时存在。3.4 Meek规则完成剩余定向用一个while循环反复应用R1到R3直到没有新方向出现def apply_meek_rules(graph, dir_matrix): 用Meek规则做方向传播原位修改dir_matrix。 n graph.shape[0] changed True while changed: changed False # R1: a-b 且 b-c 存在, a-c 不相邻 b-c for a in range(n): for b in range(n): if dir_matrix[a, b] ! 1: continue for c in range(n): if not graph[b, c]: continue if dir_matrix[b, c] ! 0: continue if graph[a, c]: continue dir_matrix[b, c] 1 dir_matrix[c, b] -1 changed True # R2: a-b-c 且 a-c 存在未定向 a-c for a in range(n): for b in range(n): if dir_matrix[a, b] ! 1: continue for c in range(n): if dir_matrix[b, c] ! 1: continue if not graph[a, c]: continue if dir_matrix[a, c] ! 0: continue dir_matrix[a, c] 1 dir_matrix[c, a] -1 changed True # R3: 有 a-c-b 和 a-d-b, # 且 a 与 b 不相邻, c 与 d 不相邻 a-b for b in range(n): incoming [c for c in range(n) if dir_matrix[c, b] 1] if len(incoming) 2: continue for idx1 in range(len(incoming)): for idx2 in range(idx1 1, len(incoming)): c, d incoming[idx1], incoming[idx2] if graph[c, d]: continue for a in range(n): if a b or a c or a d: continue if graph[a, b] or not graph[a, c] or not graph[a, d]: continue if dir_matrix[a, c] ! 0 or dir_matrix[a, d] ! 0: continue dir_matrix[a, b] 1 dir_matrix[b, a] -1 changed True return dir_matrix主流程只有三行顺序不要颠倒骨架先行V结构定向必须在Meek之前因为Meek规则会基于已有方向继续传播。如果想支持非高斯数据只需要替换fisher_z_test的返回逻辑骨架和方向部分不用动。def pc(data, alpha0.05, max_depthNone): graph, sep_set pc_skeleton(data, alpha, max_depth) dir_matrix orient_v_structures(graph, sep_set) dir_matrix apply_meek_rules(graph, dir_matrix) return graph, dir_matrix4. 合成数据实验用已知DAG验证实现并调参4.1 生成已知结构的线性高斯数据先构造一个知道标准答案的DAG再用线性结构方程生成观测数据def simulate_data(n_vars6, n_samples1000, prob_edge0.3, seed42): 生成随机DAG和线性高斯观测数据。 变量按索引从小到大连接天然保证无环。 rng np.random.default_rng(seed) dag np.zeros((n_vars, n_vars), dtypebool) for i in range(n_vars): for j in range(i 1, n_vars): if rng.random() prob_edge: dag[i, j] True # i - j data np.zeros((n_samples, n_vars)) for j in range(n_vars): parents np.where(dag[:, j])[0] if len(parents) 0: data[:, j] 0.6 * data[:, parents].sum(axis1) data[:, j] rng.normal(0, 1, n_samples) return dag, data每个变量只与索引更小的变量相连所以生成的DAG天然无环。边数由prob_edge控制0.3在6变量下平均带来约4到5条边属于中等稀疏场景。跑算法时要先把dag转成无向骨架再对比PC算法只承诺恢复骨架和部分方向def evaluate(dag, graph, dir_matrix): n dag.shape[0] true_skel dag | dag.T pred_skel graph | graph.T # 骨架的精确率、召回率 tp int((true_skel pred_skel).sum()) fp int((~true_skel pred_skel).sum()) fn int((true_skel ~pred_skel).sum()) precision tp / (tp fp) if tp fp else 0 recall tp / (tp fn) if tp fn else 0 # 方向正确数只统计算法输出且真实存在的有向边 correct_dir sum( 1 for i in range(n) for j in range(n) if dir_matrix[i, j] 1 and dag[i, j] ) total_dir int((dir_matrix 1).sum()) return precision, recall, correct_dir, total_dir dag, data simulate_data() graph, dir_matrix pc(data, alpha0.05) precision, recall, cd, td evaluate(dag, graph, dir_matrix) print(fprecision{precision:.2f}, recall{recall:.2f}, directed{cd}/{td})一次运行会打印类似precision0.90, recall0.85, directed5/10这样的输出。注意directed分数里分母是算法判断为有向的边数不是真实DAG的边数当PC只恢复出部分方向时未定向边不计入分子分母所以这个指标衡量的是“已经定向的边里猜对的比例”。4.2 三个关键参数alpha、max_depth、样本量直接影响PC结果的参数集中在下面这张表参数位置/来源默认值作用调整建议alphapc()第二参数0.05条件独立性检验阈值palpha才接受独立并删边小样本可以调到0.1变量多时避免过小否则假阳性边会变多max_depthpc()第三参数n_vars-2最大条件集大小稀疏图设2到4即可设太大会让组合数指数级上升样本量数据本身无固定出现在Fisher Z统计量中影响检验效力最好大于最大条件集再加3500为经验下限alpha的影响要看骨架收缩率alpha从0.01调到0.1删掉的边越多骨架越稀疏。如果两张图在高alpha下还残留大量短路径说明数据接近忠实性违背需要回到数据质量本身。max_depth则是一把双刃剑设小了高次条件独立关系发现不了设大了在30个变量以上时组合爆炸非常快运行时间按组合数增长。4.3 典型失败模式和如何区分最常见失败是缺失边。样本量不足时边真实的因果效应很弱Fisher Z检验没有足够功效拒绝独立性边就被删了。识别方法是看数据生成时的系数把simulate_data中的0.6改成0.1大多数边都会消失这是检验功效问题而不是实现bug。第二种是假阳性边多出现在alpha过小或数据存在强相关但非因果关系时比如两个变量同时受第三个隐性因素影响。第三种是方向错误常见于样本不足导致V结构判定错位。遇到这类情况先检查sep_set是不是空集——空分离集会大量出现在小样本中说明骨架搜索提前终止于过低的depth。5. 排错与验证给你的PC实现做一次快速体检5.1 三个低成本检查先查分离集对称性。sep_set里(i,j)与(j,i)应完全一致。不一致说明同一条件集下两次检验结论不同多半是数据非高斯或样本太少。再查p值分布。抽取fisher_z_test返回的所有p值做直方图假设成立时p值近似均匀大量p值压在0附近骨架会混入伪相关边大量压到1附近则检验功效不足真边保不住。最后看骨架密度。真实稀疏因果图的骨架密度应显著低于0.5如果几乎全连通先怀疑alpha过大或者max_depth过小。5.2 与参考实现做交叉验证不引入重依赖的验证方法是把graph矩阵导出成三元组列表与另一个PC实现比如causal-learn中的同名函数输出的骨架列表做集合比对。相差1到2条边属于正常差太多检查是不是把邻居快照取在了外层以及条件集深度的递增顺序是否一致。实现细节上的剪枝顺序可能不同但分离集应当一致。5.3 独立性检验替换接口数据明显非高斯时把Fisher Z换掉即可骨架搜索和方向规则不用动。替换时保持函数签名一致def my_test(data, x, y, cond_set): # 返回 p 值p alpha 表示条件独立 ...然后在pc_skeleton里只改一行p my_test(data, i, j, list(cond)) # 原来是 fisher_z_test(...)换成距离相关的置换检验后单次检验会从微秒级升到毫秒级max_depth必须同步调小。把alpha和max_depth设成匹配自身数据规模的取值再配合前面的体检方法PC算法输出的方向边在进入后续干预分析前才算真正可信。本文还有配套的精品资源点击获取