简介这份资源是面向数据分析与机器学习学习者的Python实现PC算法项目源码聚焦因果发现与网络结构学习场景适合具备一定Python与统计基础、希望动手理解条件独立检验与因果图构建的中高级读者。压缩包共13个文件约452KB以4个py源码文件为核心辅以3张png结果示意图、1个csv测试数据集、1个md说明文档及若干配置类文件结构紧凑、便于快速上手。项目围绕部分相关性展开涵盖数据预处理、相关矩阵计算、条件独立测试、定向边剔除与循环迭代等关键环节并给出基于networkx的因果图可视化思路还讨论了非高斯数据与大规模场景下的扩展优化方向。目前已有298人学习读者可借此掌握PC算法的完整实现脉络获得可直接运行的脚本、测试数据与结果图示为因果推断类任务提供可复用的实践参考。1. 从一份 Python 实现 PC 算法源码说起因果结构学习到底能落地在哪很多人第一次听到「Python 实现 PC 算法项目源码」脑子里冒出来的是一堆条件独立检验和邻接矩阵觉得这是论文里的东西离工程很远。但真做过因果推断项目的人会告诉你PC 算法是少数几个「代码量不大、依赖清晰、结果能直接画成图」的因果发现方法特别适合放进数据分析流水线里做变量筛选和结构假设生成。它的核心任务只有一句话给一批观测数据在假设没有隐变量的前提下尽可能还原变量之间的有向无环图。你拿到的源码通常包含三块——条件独立性检验、骨架学习、边定向跑通之后能输出一张可以直接喂给下游建模的因果图。这篇文章面向的是想把这套源码真正用起来的人做数据分析的、搞量化因子研究的、写 Python 但没系统碰过因果推断的都能跟着把环境搭起来、把参数调明白、把坑踩完。2. PC 算法源码拆开看三个核心模块与它们为什么这么设计2.1 条件独立性检验PC 算法的地基PC 算法全称 Peter-Clark 算法它的第一步不是建图而是判断「两个变量在给定某个集合后是否独立」。这个判断做不准后面全是错的。源码里最常见的实现是 Fishers Z 检验适用于连续变量且假设线性高斯如果是离散数据会换成 G² 检验或卡方检验。为什么用条件独立而不是相关系数因为相关系数只能看两两关系而因果结构里真正关键的是「控制掉其他变量后这两个还相不相关」。举个例子A 和 B 看起来强相关但可能都是 C 导致的控制 C 之后相关性消失那 A 和 B 之间就不该有边。源码里这部分通常长这样import numpy as np from scipy.stats import norm def fisher_z_test(data, x, y, cond_set, alpha0.05): 条件独立性检验给定 cond_set判断 x 和 y 是否独立 data: shape (n_samples, n_vars) x, y: 变量索引 cond_set: 条件变量索引列表 alpha: 显著性水平 n data.shape[0] # 计算偏相关系数 if len(cond_set) 0: r np.corrcoef(data[:, x], data[:, y])[0, 1] else: # 用回归残差法计算偏相关 from numpy.linalg import lstsq Z np.column_stack([data[:, c] for c in cond_set] [np.ones(n)]) res_x data[:, x] - Z lstsq(Z, data[:, x], rcondNone)[0] res_y data[:, y] - Z lstsq(Z, data[:, y], rcondNone)[0] r np.corrcoef(res_x, res_y)[0, 1] # Fisher Z 变换 r np.clip(r, -0.999999, 0.999999) z 0.5 * np.log((1 r) / (1 - r)) se 1.0 / np.sqrt(n - len(cond_set) - 3) stat np.abs(z) / se p_value 2 * (1 - norm.cdf(stat)) return p_value alpha # True 表示独立这段代码的逻辑是先用线性回归把条件变量的影响去掉拿残差算偏相关再做 Fisher Z 变换把相关系数转成近似正态统计量最后算 p 值。参数里alpha是最关键的默认 0.05 意味着你允许 5% 的假阳性。调小到 0.01边会更少、图更稀疏但可能漏掉真实边调到 0.1边更多但可能引入伪边。我一般建议先用 0.05 跑一版看结构再根据领域知识微调。cond_set的长度也影响结果条件集越大检验功效越低样本不够时容易把真边判成独立。2.2 骨架学习从完全图到稀疏邻接骨架学习是 PC 算法最耗时的部分。它从一个完全无向图开始逐步删边。具体做法是对每一对相邻节点尝试所有可能的大小为 k 的条件集如果存在某个条件集让它们独立就删掉这条边。k 从 0 开始每轮加 1直到没有边可删或达到最大条件集大小。源码里通常用一个邻接矩阵表示骨架1 表示有边0 表示无边。下面是一个简化实现def learn_skeleton(data, alpha0.05, max_cond_sizeNone): 骨架学习返回邻接矩阵和分离集 n_vars data.shape[1] if max_cond_size is None: max_cond_size n_vars - 2 adj np.ones((n_vars, n_vars)) - np.eye(n_vars) # 完全图 sep_set {} # 记录分离集用于后续定向 for k in range(max_cond_size 1): edges_to_remove [] for i in range(n_vars): for j in range(i 1, n_vars): if adj[i, j] 0: continue # 找 i 和 j 的邻居排除彼此 neighbors [n for n in range(n_vars) if adj[i, n] 1 and n ! j] if len(neighbors) k: continue from itertools import combinations for cond in combinations(neighbors, k): if fisher_z_test(data, i, j, list(cond), alpha): edges_to_remove.append((i, j)) sep_set[(i, j)] list(cond) sep_set[(j, i)] list(cond) break for i, j in edges_to_remove: adj[i, j] 0 adj[j, i] 0 if not edges_to_remove: break return adj, sep_set这里有几个参数值得说清楚。max_cond_size控制条件集的最大大小默认是变量数减 2但实际跑的时候如果变量超过 20 个建议手动设成 3 或 4否则组合爆炸跑一晚上都出不来。alpha和检验函数里保持一致。sep_set是分离集记录哪条边被哪个条件集删掉的后面定向阶段要用。这段代码的时间复杂度是指数级的变量数 10 以内还好20 个以上就要考虑用优化版本或者降维。2.3 边定向从无向骨架到有向因果图骨架只告诉你「谁和谁有关」定向才告诉你「谁导致谁」。PC 算法的定向规则主要有三条碰撞节点识别、避免新碰撞、传播方向。碰撞节点是指 X → Z ← Y 这种结构其中 X 和 Y 不相邻。识别方法是看 Z 是否在 X 和 Y 的分离集中如果不在就是碰撞节点。def orient_edges(adj, sep_set): 边定向返回有向邻接矩阵 n_vars adj.shape[0] directed adj.copy() # 规则1识别碰撞节点 for i in range(n_vars): for j in range(i 1, n_vars): if adj[i, j] 1: continue # 找共同邻居 common [k for k in range(n_vars) if adj[i, k] 1 and adj[j, k] 1] for k in common: if (i, j) not in sep_set or k not in sep_set[(i, j)]: # k 是碰撞节点定向 i - k - j directed[i, k] 1 directed[k, i] 0 directed[j, k] 1 directed[k, j] 0 # 规则2和3避免新碰撞和传播 changed True while changed: changed False for i in range(n_vars): for j in range(n_vars): if directed[i, j] 1 and directed[j, i] 0: # i - j检查 j 的其他邻居 for k in range(n_vars): if k i or directed[j, k] 0: continue if directed[k, j] 1: continue if adj[i, k] 0: # 避免 i - j - k 变成新碰撞 directed[j, k] 1 directed[k, j] 0 changed True return directed定向阶段最容易出问题的是「部分定向」——有些边方向确定不了这是正常的PC 算法只能保证在马尔可夫等价类内唯一。如果你发现输出图里还有双向边或未定向边别慌先检查分离集记录是否正确再确认样本量够不够。样本少于 200 时定向结果波动会很大。3. 把源码跑起来环境配置、数据准备与最小可复现示例3.1 环境配置Python 版本、依赖包与 VSCode 设置这套源码对 Python 版本不挑3.8 到 3.11 都能跑。依赖主要是 numpy、scipy、pandas如果要画图再加 networkx 和 matplotlib。安装命令很直接pip install numpy scipy pandas networkx matplotlib如果你用 VSCode建议装 Python 扩展和 Pylance然后在项目根目录建一个.vscode/settings.json把解释器路径指到你虚拟环境里的 python。虚拟环境用 venv 就行python -m venv venv source venv/bin/activate # Linux/Mac venv\Scripts\activate # Windows pip install -r requirements.txtrequirements.txt 里至少写清楚版本范围避免 numpy 2.0 和旧版 scipy 的兼容问题。我一般会锁 numpy1.21,2.0scipy1.7。如果你在 Linux 上跑系统自带的 python3 可能缺 pip先sudo apt install python3-pip python3-venv再建虚拟环境。3.2 数据准备从 CSV 到算法输入的三个检查点PC 算法对输入数据有隐含假设连续变量、无缺失值、近似线性高斯。拿到一份 CSV 之后先做三件事。第一检查缺失值有缺失就删行或插补别直接扔进去Fisher Z 检验遇到 NaN 会返回 NaN整个骨架学习就废了。第二检查变量类型分类变量要么转成数值编码要么换离散检验方法。第三标准化虽然 Fisher Z 对尺度不敏感但回归残差法在量纲差异大时数值不稳定。import pandas as pd from sklearn.preprocessing import StandardScaler def load_and_prepare(csv_path): df pd.read_csv(csv_path) # 检查缺失 if df.isnull().any().any(): print(存在缺失值删除含缺失的行) df df.dropna() # 只保留数值列 df df.select_dtypes(include[np.number]) # 标准化 scaler StandardScaler() data scaler.fit_transform(df.values) return data, df.columns.tolist()select_dtypes那一步很重要很多人直接把带字符串列的 DataFrame 转 numpy结果报类型错误。标准化用 StandardScaler 就够了不用 MinMax因为 Fisher Z 基于相关系数对线性变换不变。列名要保留后面画图时节点标签用得上。3.3 最小可复现示例20 行代码跑通 PC 算法把前面的模块串起来一个最小可复现的 PC 算法大概是这样import numpy as np from itertools import combinations from scipy.stats import norm def pc_algorithm(data, alpha0.05, max_cond_size3): n_vars data.shape[1] adj np.ones((n_vars, n_vars)) - np.eye(n_vars) sep_set {} # 骨架学习 for k in range(max_cond_size 1): to_remove [] for i in range(n_vars): for j in range(i1, n_vars): if adj[i, j] 0: continue neighbors [n for n in range(n_vars) if adj[i, n] 1 and n ! j] if len(neighbors) k: continue for cond in combinations(neighbors, k): if fisher_z_test(data, i, j, list(cond), alpha): to_remove.append((i, j)) sep_set[(i, j)] list(cond) sep_set[(j, i)] list(cond) break for i, j in to_remove: adj[i, j] adj[j, i] 0 if not to_remove: break # 定向 directed orient_edges(adj, sep_set) return directed, adj, sep_set # 生成模拟数据X - Y - Z np.random.seed(42) n 500 X np.random.randn(n) Y 0.8 * X np.random.randn(n) * 0.5 Z 0.7 * Y np.random.randn(n) * 0.5 data np.column_stack([X, Y, Z]) directed, skeleton, sep pc_algorithm(data, alpha0.05, max_cond_size2) print(有向邻接矩阵) print(directed)跑出来应该能看到 X→Y 和 Y→Z 的边X 和 Z 之间没有直接边。如果结果不对先检查 alpha 是不是太大或太小再检查样本量。500 个样本对 3 个变量来说足够了变量多了要相应增加样本。4. 避坑与排查PC 算法源码落地时最容易翻车的五个地方4.1 样本量不够导致边方向乱跳现象同样的数据换一个随机种子跑定向结果完全不一样。原因Fisher Z 检验在小样本下功效低p 值波动大骨架不稳定定向自然乱。解决样本量至少是变量数的 10 倍最好 20 倍以上。如果数据不够降低 alpha 到 0.01 减少伪边或者先用领域知识固定一部分边。4.2 条件集组合爆炸把内存跑满现象变量数超过 15 个时程序卡死或报 MemoryError。原因骨架学习里combinations(neighbors, k)在 k 增大时组合数指数增长。解决手动设max_cond_size3或 4别用默认值。如果变量确实多先用相关性或 LASSO 做一轮筛选把变量降到 10 个以内再跑 PC。4.3 离散变量直接套 Fisher Z 检验现象输入是 0/1 或分类编码跑出来的图全是边或全没边。原因Fisher Z 假设连续高斯离散数据违反假设p 值不可信。解决离散数据换 G² 检验或卡方检验或者用pgmpy里的chi_square检验函数。如果非要用 Fisher Z先把离散变量做 one-hot 或累积编码但效果一般。4.4 缺失值没处理导致检验结果全为 NaN现象程序不报错但骨架矩阵全是 1一条边都没删。原因数据里有 NaNnp.corrcoef返回 NaNp_value alpha比较返回 False所有边都保留。解决在load_and_prepare里强制dropna()或者用均值/中位数插补。插补后记得检查插补比例超过 10% 就要考虑换数据源。4.5 把 PC 算法输出当成绝对因果结论现象拿着输出图直接写报告说「A 导致 B」。原因PC 算法假设无隐变量、无选择偏差、线性高斯这些在实际数据里很少同时满足。解决把输出图当成「结构假设」而不是「因果结论」用领域知识验证或者用干预数据做进一步确认。如果发现有隐变量嫌疑换 FCI 算法它能处理部分隐变量。5. 进阶技巧用 bootstrap 稳定边、用先验知识约束图、用 networkx 可视化跑通基础版之后真正让 PC 算法源码在生产里可用的是三件事稳定性评估、先验约束、可视化。稳定性评估用 bootstrap对数据有放回抽样 100 次每次跑一遍 PC统计每条边出现的频率频率低于 0.6 的边直接删掉。这个技巧能大幅减少伪边代价是计算时间翻 100 倍变量少的时候值得做。def bootstrap_pc(data, n_bootstrap100, alpha0.05, threshold0.6): n_vars data.shape[1] edge_count np.zeros((n_vars, n_vars)) for b in range(n_bootstrap): idx np.random.choice(data.shape[0], data.shape[0], replaceTrue) sample data[idx] directed, _, _ pc_algorithm(sample, alphaalpha, max_cond_size3) edge_count (directed directed.T 0).astype(int) stable (edge_count / n_bootstrap) threshold return stable.astype(int)先验知识约束是另一个实用技巧。如果你明确知道某些边不可能存在比如时间上后发生的变量不能导致先发生的可以在骨架学习前把对应位置设成 0直接禁止。这能减少搜索空间也能避免算法被数据噪声带偏。可视化用 networkx把有向图转成 DiGraph节点标签用列名边用箭头输出 PNG 或交互式 HTML。import networkx as nx import matplotlib.pyplot as plt def plot_causal_graph(directed, labels): G nx.DiGraph() for i, label in enumerate(labels): G.add_node(i, namelabel) for i in range(directed.shape[0]): for j in range(directed.shape[1]): if directed[i, j] 1: G.add_edge(i, j) pos nx.spring_layout(G, seed42) nx.draw(G, pos, with_labelsTrue, labels{i: labels[i] for i in range(len(labels))}, node_colorlightblue, edge_colorgray, arrowsize20, font_size10) plt.savefig(causal_graph.png, dpi150, bbox_inchestight) plt.close()最后说一个我自己的习惯每次跑完 PC 算法我都会把 alpha、max_cond_size、样本量、边数这四个数记在一个日志文件里。因为因果发现的结果对参数太敏感了不记录的话过两周你根本想不起来当时为什么选了 0.05 而不是 0.01。这个习惯帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取 Python实现PC算法源码解析：因果结构学习实战指南
简介这份资源是面向数据分析与机器学习学习者的Python实现PC算法项目源码聚焦因果发现与网络结构学习场景适合具备一定Python与统计基础、希望动手理解条件独立检验与因果图构建的中高级读者。压缩包共13个文件约452KB以4个py源码文件为核心辅以3张png结果示意图、1个csv测试数据集、1个md说明文档及若干配置类文件结构紧凑、便于快速上手。项目围绕部分相关性展开涵盖数据预处理、相关矩阵计算、条件独立测试、定向边剔除与循环迭代等关键环节并给出基于networkx的因果图可视化思路还讨论了非高斯数据与大规模场景下的扩展优化方向。目前已有298人学习读者可借此掌握PC算法的完整实现脉络获得可直接运行的脚本、测试数据与结果图示为因果推断类任务提供可复用的实践参考。1. 从一份 Python 实现 PC 算法源码说起因果结构学习到底能落地在哪很多人第一次听到「Python 实现 PC 算法项目源码」脑子里冒出来的是一堆条件独立检验和邻接矩阵觉得这是论文里的东西离工程很远。但真做过因果推断项目的人会告诉你PC 算法是少数几个「代码量不大、依赖清晰、结果能直接画成图」的因果发现方法特别适合放进数据分析流水线里做变量筛选和结构假设生成。它的核心任务只有一句话给一批观测数据在假设没有隐变量的前提下尽可能还原变量之间的有向无环图。你拿到的源码通常包含三块——条件独立性检验、骨架学习、边定向跑通之后能输出一张可以直接喂给下游建模的因果图。这篇文章面向的是想把这套源码真正用起来的人做数据分析的、搞量化因子研究的、写 Python 但没系统碰过因果推断的都能跟着把环境搭起来、把参数调明白、把坑踩完。2. PC 算法源码拆开看三个核心模块与它们为什么这么设计2.1 条件独立性检验PC 算法的地基PC 算法全称 Peter-Clark 算法它的第一步不是建图而是判断「两个变量在给定某个集合后是否独立」。这个判断做不准后面全是错的。源码里最常见的实现是 Fishers Z 检验适用于连续变量且假设线性高斯如果是离散数据会换成 G² 检验或卡方检验。为什么用条件独立而不是相关系数因为相关系数只能看两两关系而因果结构里真正关键的是「控制掉其他变量后这两个还相不相关」。举个例子A 和 B 看起来强相关但可能都是 C 导致的控制 C 之后相关性消失那 A 和 B 之间就不该有边。源码里这部分通常长这样import numpy as np from scipy.stats import norm def fisher_z_test(data, x, y, cond_set, alpha0.05): 条件独立性检验给定 cond_set判断 x 和 y 是否独立 data: shape (n_samples, n_vars) x, y: 变量索引 cond_set: 条件变量索引列表 alpha: 显著性水平 n data.shape[0] # 计算偏相关系数 if len(cond_set) 0: r np.corrcoef(data[:, x], data[:, y])[0, 1] else: # 用回归残差法计算偏相关 from numpy.linalg import lstsq Z np.column_stack([data[:, c] for c in cond_set] [np.ones(n)]) res_x data[:, x] - Z lstsq(Z, data[:, x], rcondNone)[0] res_y data[:, y] - Z lstsq(Z, data[:, y], rcondNone)[0] r np.corrcoef(res_x, res_y)[0, 1] # Fisher Z 变换 r np.clip(r, -0.999999, 0.999999) z 0.5 * np.log((1 r) / (1 - r)) se 1.0 / np.sqrt(n - len(cond_set) - 3) stat np.abs(z) / se p_value 2 * (1 - norm.cdf(stat)) return p_value alpha # True 表示独立这段代码的逻辑是先用线性回归把条件变量的影响去掉拿残差算偏相关再做 Fisher Z 变换把相关系数转成近似正态统计量最后算 p 值。参数里alpha是最关键的默认 0.05 意味着你允许 5% 的假阳性。调小到 0.01边会更少、图更稀疏但可能漏掉真实边调到 0.1边更多但可能引入伪边。我一般建议先用 0.05 跑一版看结构再根据领域知识微调。cond_set的长度也影响结果条件集越大检验功效越低样本不够时容易把真边判成独立。2.2 骨架学习从完全图到稀疏邻接骨架学习是 PC 算法最耗时的部分。它从一个完全无向图开始逐步删边。具体做法是对每一对相邻节点尝试所有可能的大小为 k 的条件集如果存在某个条件集让它们独立就删掉这条边。k 从 0 开始每轮加 1直到没有边可删或达到最大条件集大小。源码里通常用一个邻接矩阵表示骨架1 表示有边0 表示无边。下面是一个简化实现def learn_skeleton(data, alpha0.05, max_cond_sizeNone): 骨架学习返回邻接矩阵和分离集 n_vars data.shape[1] if max_cond_size is None: max_cond_size n_vars - 2 adj np.ones((n_vars, n_vars)) - np.eye(n_vars) # 完全图 sep_set {} # 记录分离集用于后续定向 for k in range(max_cond_size 1): edges_to_remove [] for i in range(n_vars): for j in range(i 1, n_vars): if adj[i, j] 0: continue # 找 i 和 j 的邻居排除彼此 neighbors [n for n in range(n_vars) if adj[i, n] 1 and n ! j] if len(neighbors) k: continue from itertools import combinations for cond in combinations(neighbors, k): if fisher_z_test(data, i, j, list(cond), alpha): edges_to_remove.append((i, j)) sep_set[(i, j)] list(cond) sep_set[(j, i)] list(cond) break for i, j in edges_to_remove: adj[i, j] 0 adj[j, i] 0 if not edges_to_remove: break return adj, sep_set这里有几个参数值得说清楚。max_cond_size控制条件集的最大大小默认是变量数减 2但实际跑的时候如果变量超过 20 个建议手动设成 3 或 4否则组合爆炸跑一晚上都出不来。alpha和检验函数里保持一致。sep_set是分离集记录哪条边被哪个条件集删掉的后面定向阶段要用。这段代码的时间复杂度是指数级的变量数 10 以内还好20 个以上就要考虑用优化版本或者降维。2.3 边定向从无向骨架到有向因果图骨架只告诉你「谁和谁有关」定向才告诉你「谁导致谁」。PC 算法的定向规则主要有三条碰撞节点识别、避免新碰撞、传播方向。碰撞节点是指 X → Z ← Y 这种结构其中 X 和 Y 不相邻。识别方法是看 Z 是否在 X 和 Y 的分离集中如果不在就是碰撞节点。def orient_edges(adj, sep_set): 边定向返回有向邻接矩阵 n_vars adj.shape[0] directed adj.copy() # 规则1识别碰撞节点 for i in range(n_vars): for j in range(i 1, n_vars): if adj[i, j] 1: continue # 找共同邻居 common [k for k in range(n_vars) if adj[i, k] 1 and adj[j, k] 1] for k in common: if (i, j) not in sep_set or k not in sep_set[(i, j)]: # k 是碰撞节点定向 i - k - j directed[i, k] 1 directed[k, i] 0 directed[j, k] 1 directed[k, j] 0 # 规则2和3避免新碰撞和传播 changed True while changed: changed False for i in range(n_vars): for j in range(n_vars): if directed[i, j] 1 and directed[j, i] 0: # i - j检查 j 的其他邻居 for k in range(n_vars): if k i or directed[j, k] 0: continue if directed[k, j] 1: continue if adj[i, k] 0: # 避免 i - j - k 变成新碰撞 directed[j, k] 1 directed[k, j] 0 changed True return directed定向阶段最容易出问题的是「部分定向」——有些边方向确定不了这是正常的PC 算法只能保证在马尔可夫等价类内唯一。如果你发现输出图里还有双向边或未定向边别慌先检查分离集记录是否正确再确认样本量够不够。样本少于 200 时定向结果波动会很大。3. 把源码跑起来环境配置、数据准备与最小可复现示例3.1 环境配置Python 版本、依赖包与 VSCode 设置这套源码对 Python 版本不挑3.8 到 3.11 都能跑。依赖主要是 numpy、scipy、pandas如果要画图再加 networkx 和 matplotlib。安装命令很直接pip install numpy scipy pandas networkx matplotlib如果你用 VSCode建议装 Python 扩展和 Pylance然后在项目根目录建一个.vscode/settings.json把解释器路径指到你虚拟环境里的 python。虚拟环境用 venv 就行python -m venv venv source venv/bin/activate # Linux/Mac venv\Scripts\activate # Windows pip install -r requirements.txtrequirements.txt 里至少写清楚版本范围避免 numpy 2.0 和旧版 scipy 的兼容问题。我一般会锁 numpy1.21,2.0scipy1.7。如果你在 Linux 上跑系统自带的 python3 可能缺 pip先sudo apt install python3-pip python3-venv再建虚拟环境。3.2 数据准备从 CSV 到算法输入的三个检查点PC 算法对输入数据有隐含假设连续变量、无缺失值、近似线性高斯。拿到一份 CSV 之后先做三件事。第一检查缺失值有缺失就删行或插补别直接扔进去Fisher Z 检验遇到 NaN 会返回 NaN整个骨架学习就废了。第二检查变量类型分类变量要么转成数值编码要么换离散检验方法。第三标准化虽然 Fisher Z 对尺度不敏感但回归残差法在量纲差异大时数值不稳定。import pandas as pd from sklearn.preprocessing import StandardScaler def load_and_prepare(csv_path): df pd.read_csv(csv_path) # 检查缺失 if df.isnull().any().any(): print(存在缺失值删除含缺失的行) df df.dropna() # 只保留数值列 df df.select_dtypes(include[np.number]) # 标准化 scaler StandardScaler() data scaler.fit_transform(df.values) return data, df.columns.tolist()select_dtypes那一步很重要很多人直接把带字符串列的 DataFrame 转 numpy结果报类型错误。标准化用 StandardScaler 就够了不用 MinMax因为 Fisher Z 基于相关系数对线性变换不变。列名要保留后面画图时节点标签用得上。3.3 最小可复现示例20 行代码跑通 PC 算法把前面的模块串起来一个最小可复现的 PC 算法大概是这样import numpy as np from itertools import combinations from scipy.stats import norm def pc_algorithm(data, alpha0.05, max_cond_size3): n_vars data.shape[1] adj np.ones((n_vars, n_vars)) - np.eye(n_vars) sep_set {} # 骨架学习 for k in range(max_cond_size 1): to_remove [] for i in range(n_vars): for j in range(i1, n_vars): if adj[i, j] 0: continue neighbors [n for n in range(n_vars) if adj[i, n] 1 and n ! j] if len(neighbors) k: continue for cond in combinations(neighbors, k): if fisher_z_test(data, i, j, list(cond), alpha): to_remove.append((i, j)) sep_set[(i, j)] list(cond) sep_set[(j, i)] list(cond) break for i, j in to_remove: adj[i, j] adj[j, i] 0 if not to_remove: break # 定向 directed orient_edges(adj, sep_set) return directed, adj, sep_set # 生成模拟数据X - Y - Z np.random.seed(42) n 500 X np.random.randn(n) Y 0.8 * X np.random.randn(n) * 0.5 Z 0.7 * Y np.random.randn(n) * 0.5 data np.column_stack([X, Y, Z]) directed, skeleton, sep pc_algorithm(data, alpha0.05, max_cond_size2) print(有向邻接矩阵) print(directed)跑出来应该能看到 X→Y 和 Y→Z 的边X 和 Z 之间没有直接边。如果结果不对先检查 alpha 是不是太大或太小再检查样本量。500 个样本对 3 个变量来说足够了变量多了要相应增加样本。4. 避坑与排查PC 算法源码落地时最容易翻车的五个地方4.1 样本量不够导致边方向乱跳现象同样的数据换一个随机种子跑定向结果完全不一样。原因Fisher Z 检验在小样本下功效低p 值波动大骨架不稳定定向自然乱。解决样本量至少是变量数的 10 倍最好 20 倍以上。如果数据不够降低 alpha 到 0.01 减少伪边或者先用领域知识固定一部分边。4.2 条件集组合爆炸把内存跑满现象变量数超过 15 个时程序卡死或报 MemoryError。原因骨架学习里combinations(neighbors, k)在 k 增大时组合数指数增长。解决手动设max_cond_size3或 4别用默认值。如果变量确实多先用相关性或 LASSO 做一轮筛选把变量降到 10 个以内再跑 PC。4.3 离散变量直接套 Fisher Z 检验现象输入是 0/1 或分类编码跑出来的图全是边或全没边。原因Fisher Z 假设连续高斯离散数据违反假设p 值不可信。解决离散数据换 G² 检验或卡方检验或者用pgmpy里的chi_square检验函数。如果非要用 Fisher Z先把离散变量做 one-hot 或累积编码但效果一般。4.4 缺失值没处理导致检验结果全为 NaN现象程序不报错但骨架矩阵全是 1一条边都没删。原因数据里有 NaNnp.corrcoef返回 NaNp_value alpha比较返回 False所有边都保留。解决在load_and_prepare里强制dropna()或者用均值/中位数插补。插补后记得检查插补比例超过 10% 就要考虑换数据源。4.5 把 PC 算法输出当成绝对因果结论现象拿着输出图直接写报告说「A 导致 B」。原因PC 算法假设无隐变量、无选择偏差、线性高斯这些在实际数据里很少同时满足。解决把输出图当成「结构假设」而不是「因果结论」用领域知识验证或者用干预数据做进一步确认。如果发现有隐变量嫌疑换 FCI 算法它能处理部分隐变量。5. 进阶技巧用 bootstrap 稳定边、用先验知识约束图、用 networkx 可视化跑通基础版之后真正让 PC 算法源码在生产里可用的是三件事稳定性评估、先验约束、可视化。稳定性评估用 bootstrap对数据有放回抽样 100 次每次跑一遍 PC统计每条边出现的频率频率低于 0.6 的边直接删掉。这个技巧能大幅减少伪边代价是计算时间翻 100 倍变量少的时候值得做。def bootstrap_pc(data, n_bootstrap100, alpha0.05, threshold0.6): n_vars data.shape[1] edge_count np.zeros((n_vars, n_vars)) for b in range(n_bootstrap): idx np.random.choice(data.shape[0], data.shape[0], replaceTrue) sample data[idx] directed, _, _ pc_algorithm(sample, alphaalpha, max_cond_size3) edge_count (directed directed.T 0).astype(int) stable (edge_count / n_bootstrap) threshold return stable.astype(int)先验知识约束是另一个实用技巧。如果你明确知道某些边不可能存在比如时间上后发生的变量不能导致先发生的可以在骨架学习前把对应位置设成 0直接禁止。这能减少搜索空间也能避免算法被数据噪声带偏。可视化用 networkx把有向图转成 DiGraph节点标签用列名边用箭头输出 PNG 或交互式 HTML。import networkx as nx import matplotlib.pyplot as plt def plot_causal_graph(directed, labels): G nx.DiGraph() for i, label in enumerate(labels): G.add_node(i, namelabel) for i in range(directed.shape[0]): for j in range(directed.shape[1]): if directed[i, j] 1: G.add_edge(i, j) pos nx.spring_layout(G, seed42) nx.draw(G, pos, with_labelsTrue, labels{i: labels[i] for i in range(len(labels))}, node_colorlightblue, edge_colorgray, arrowsize20, font_size10) plt.savefig(causal_graph.png, dpi150, bbox_inchestight) plt.close()最后说一个我自己的习惯每次跑完 PC 算法我都会把 alpha、max_cond_size、样本量、边数这四个数记在一个日志文件里。因为因果发现的结果对参数太敏感了不记录的话过两周你根本想不起来当时为什么选了 0.05 而不是 0.01。这个习惯帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取