简介NBCBNoise-Based-then-Constraint-Based与CBNBConstraint-Based-then-Noise-Based是时间序列因果推断中融合噪声分析与约束检验的混合算法前者先锁定因果顺序再修剪边后者则先做条件独立性检验再排因果序适合科研人员、数据科学家及因果推断从业者用于算法复现与对比。文档从线性动态结构因果模型出发生成模拟时间序列逐步实现VarLiNGAM因果顺序发现、PCMIC条件独立检验与剪枝并组装成NBCB与CBNB两条端到端流程每个环节都附有可运行代码、环境配置命令和逐步图文解释以模拟数据集展示最终效果。压缩包共1个docx文档大小约25KB包含环境搭建指引、完整代码块、逐段注释与运行结果分析并延伸到复杂条件独立性检测、并行化加速、非线性因果关系处理和隐含混淆因素应对等进阶方向。目前已有71人浏览学习资源体量轻但步骤完整适合快速复现两类混合算法也可为实际项目开发或理论研究提供明确起点。1. 为什么在时间序列因果发现里没人只用一种算法时间序列因果发现这两年从论文圈火到工程圈但真正动手把 NBCBNoise-Based-then-Constraint-Based和 CBNBConstraint-Based-then-Noise-Based这套混合算法完整跑通的人并不多。原因是两派方法各有一半话语权基于约束的方法擅长剪边、输出骨架却在方向判定和环处理上反复翻车基于噪声的方法依赖非高斯性和残差独立性假设因果顺序排得漂亮但对冗余联系几乎不做清理。论文把两条路缝在一起给出先噪声后约束、先约束后噪声两条组合路径并附带了数据生成、顺序发现、边修剪到非线性扩展的完整 Python 实现。这份资源适合正在复现时间序列因果推断算法的科研人员也适合想把预测任务升级成归因分析的数据科学从业者——它能帮你在一份代码里同时看到两类算法的接口怎么设计、参数怎么互相影响。2. 两条技术路线约束派与噪声派各自的底气与软肋2.1 约束派条件独立性测试为什么只能出骨架基于约束的因果发现思路很直观如果 X 和 Y 在给定条件集 Z 的情况下条件独立那么 X 和 Y 之间不存在直接因果边。实际操作就是先构造一张完全图逐对变量做条件独立性检验不显著的边直接删掉最后剩下来的就是因果骨架。PCMCI 这个名字拆开看是两件事的合体PC 是经典约束算法的骨架搜索逻辑MCI 是 Moment-Conditional Independence 的缩写专门针对时间序列做了滞后调整——检验 X_t 与 Y_{t-lag} 的关系时条件集里要同时放进 X 与 Y 各自的滞后值避免把自相关误判成因果关系。论文里的 pcmci_plus 是一个高度简化版核心代码只有两层循环from scipy.stats import pearsonr def pcmci_plus(data, max_lag, alpha0.05): n_samples, n_variables data.shape graph np.ones((n_variables, n_variables), dtypebool) for i in range(n_variables): for j in range(n_variables): if i j: graph[i, j] False else: graph[i, j] True for i in range(n_variables): for j in range(n_variables): if graph[i, j]: p_value, _ pearsonr(data[max_lag:, i], data[max_lag:, j]) if p_value alpha: graph[i, j] False return graph逻辑说明先初始化一张对角线为 False、其余全为 True 的完全候选图graph[i, j] 表示“存在一条从 i 指向 j 的候选边”。然后对每一对变量计算 Pearson 相关系数并取 p 值p 值大于显著性水平 alpha 就判定为条件独立、删除该边。需要特别注意这段代码使用的是同一时刻截面的相关性而不是严格意义上的时间滞后条件独立检验所以它输出的只是一个粗糙骨架。参数说明alpha 是该流程里最值得调的参数。模拟数据里噪声小可以放宽到 0.05真实业务数据我一般收紧到 0.01因为时间序列样本不独立Pearson 检验的 p 值会系统性偏小alpha 太松会留一堆假边。max_lag 在这里没有直接参与相关计算但它会影响后面所有基于滞后切片的回归和检验所以务必先确认 max_lag 符合数据实际生成逻辑。约束派的软肋也在这段代码里暴露得很清楚它只能回答“有没有边”回答不了“边朝哪个方向”。尤其遇到无向环时PC 一族的算法会直接输出不确定性标记把方向问题抛给下游模块——这正是混合算法里噪声派要补的位置。2.2 噪声派残差独立性与非高斯性如何给出方向基于噪声的方法换了个更硬的假设如果因果模型写成 X f(Z) E且噪声项 E 与原因变量 Z 独立那么把因果方向搞反之后反过来回归得到的残差之间会出现可检测的独立性破坏。换句话说正确的因果方向下残差最“干净”。VarLiNGAM 实现的路径是把这一逻辑转换成因果顺序搜索先找到最像外生变量的那个变量残差分布最满足独立性假设把它从系统里剔除再对剩余变量重复这个过程直到排完所有变量。论文代码里用的挑选指标是峰度而且选的是峰度最小from sklearn.linear_model import LinearRegression from scipy.stats import kurtosis def varlingam(data, max_lag): n_samples, n_variables data.shape residuals np.zeros_like(data) for i in range(n_variables): X np.hstack([data[max_lag - lag:-lag, :] for lag in range(1, max_lag 1)]) y data[max_lag:, i] model LinearRegression() model.fit(X, y) residuals[max_lag:, i] y - model.predict(X) causal_order [] remaining_vars list(range(n_variables)) while remaining_vars: kurtosis_list [] for i in remaining_vars: kurtosis_list.append(kurtosis(residuals[max_lag:, i])) min_kurtosis_idx np.argmin(kurtosis_list) causal_order.append(remaining_vars.pop(min_kurtosis_idx)) return causal_order逻辑说明第一段循环对每个变量 i 做回归特征矩阵是所有变量从 t-1 到 t-max_lag 的滞后值拟合目标是变量 i 在 t 时刻的值残差代表“去除历史影响后剩下的噪声”。第二段循环在剩余变量里逐个比较残差峰度峰度最小的变量被判定为当前最外生、最先被确定顺序的节点移出候选集后继续迭代。参数说明LinearRegression 默认用最小二乘对模拟的线性 SCM 数据刚好匹配y 与 X 的行数对齐是关键data[max_lag:] 保证了有效样本起点一致。这里有个隐藏坑峰度选 min 还是 max 在不同实现里说法不一论文复现代码里用的是 argmin也就是说它认为最接近正态分布的残差对应最外生变量这是基于噪声独立假设推导出来的选择不是拍脑袋。2.3 混合的价值两个模块怎样互相兜底NBCB 的设计是先跑 VarLiNGAM 得到完整因果顺序再跑 PCMCI 去剪掉冗余边。这样做的好处是方向由噪声派保证边的稀疏性由约束派兜底。CBNB 反过来先用 PCMCI 剪骨架再对骨架中尚未定向的边用 VarLiNGAM 给出的顺序统一打方向。论文里两个流程都跑同一组模拟数据目的就是对比在不同噪声水平下谁的方向错误更少、谁的假边更多。实际选择上没有绝对优劣。骨架本身很稀疏、只有少量真边时CBNB 的约束阶段会把搜索空间压得很小噪声派只需要处理少量候选边速度快且方向稳定反之变量间联系密集、骨架剪不干净时NBCB 先用噪声派把顺序钉死约束派再剪边就不容易把方向搞乱。工程上我一般先看变量数量和滞后阶数n_variables 大于 8 或 max_lag 大于 3 时优先跑 CBNB省掉噪声派对全图所有变量排序的开销。3. 环境准备与模拟数据生成从零搭起一份可复现的时间序列实验3.1 依赖库安装与版本选择这份代码的依赖非常轻不需要 GPU不需要额外编译。安装命令如下pip install numpy pandas scipy scikit-learn networkx这里的库各有分工numpy 承担数据和系数矩阵的张量运算scipy 提供峰度计算和 Pearson 相关检验scikit-learn 提供线性回归与高斯过程回归networkx 虽然论文代码里没直接使用但做因果图可视化、计算可达性时会用到建议一并装上。Python 版本建议 3.8 以上scikit-learn 1.0 以下的老版本在 GaussianProcessRegressor 接口上有些差异如果后面跑非线性扩展时报参数错误优先检查库版本。3.2 数据生成器线性动态 SCM 的完整实现复现的第一步是拿到一份“已知因果结构”的数据否则算法输出结果没法验证对错。数据生成器通过随机系数矩阵配合滞后累加构造出一个标准的线性动态结构因果模型import numpy as np def generate_time_series(n_samples, n_variables, max_lag, noise_std0.1): # 随机系数矩阵shape: (因变量, 自变量, 滞后阶) coefficients np.random.uniform(-1, 1, (n_variables, n_variables, max_lag)) coefficients[np.abs(coefficients) 0.1] 0 # 稀疏化 data np.zeros((n_samples, n_variables)) for t in range(max_lag, n_samples): for i in range(n_variables): for j in range(n_variables): for lag in range(1, max_lag 1): data[t, i] coefficients[i, j, lag - 1] * data[t - lag, j] data[t, i] np.random.normal(0, noise_std) return data逻辑说明coefficients[i, j, lag] 表示变量 j 在 t-lag 时刻对变量 i 在 t 时刻的线性影响系数。abs 小于 0.1 的系数直接置零这一步非常关键它保证生成的真实因果图是稀疏的算法剪边时才有明确目标。data[t, i] 的累加严格对应滞后阶不超过 max_lag 的动态结构噪声项是独立同分布的高斯白噪声。参数说明n_samples 低于 500 时回归残差方差过大峰度排序结果会很飘建议至少 1000noise_std 控制信噪比0.1 表示信号远强于噪声0.5 以上算法开始明显失真这是测试鲁棒性时优先调整的参数max_lag 必须与实际业务先验一致设置过大相当于给算法塞入大量无效特征。3.3 数据生成的自检不验证就往下跑等于盲调生成数据后别急着跑算法先做两个快速自检。第一检查系数矩阵的非零比例理论上每对变量每个滞后阶只有少量系数非零第二随机挑一个变量画它的自相关图确认滞后结构与 max_lag 一致。常见做法是打印 coefficients 里非零项的索引与数值确认没有出现全零行——那意味着某个变量完全是噪声因果发现对它没有任何意义。这一步 30 秒就能完成能省掉后面大半的排查时间。4. 核心算法复现VarLiNGAM、PCMCI 与两条混合路径的完整实现4.1 VarLiNGAM 的局部实现细节拆解from sklearn.linear_model import LinearRegression from scipy.stats import kurtosis def varlingam(data, max_lag): n_samples, n_variables data.shape residuals np.zeros_like(data) for i in range(n_variables): X np.hstack([data[max_lag - lag:-lag, :] for lag in range(1, max_lag 1)]) y data[max_lag:, i] model LinearRegression() model.fit(X, y) residuals[max_lag:, i] y - model.predict(X) causal_order [] remaining_vars list(range(n_variables)) while remaining_vars: kurtosis_list [] for i in remaining_vars: kurtosis_list.append(kurtosis(residuals[max_lag:, i])) min_kurtosis_idx np.argmin(kurtosis_list) causal_order.append(remaining_vars.pop(min_kurtosis_idx)) return causal_order这段代码里最需要留意的不是回归本身而是 X 的构造方式。data[max_lag - lag:-lag, :]中 lag 从 1 取到 max_lag实际是把每个变量的历史值横向拼接成特征矩阵lag1 对应 t-1 时刻lagmax_lag 对应 t-max_lag 时刻。y data[max_lag:, i]从第 max_lag 行开始取当前时刻值保证特征和预测目标的时间对齐。这里如果写成data[max_lag - lag:-lag]的步进有偏移残差会整体错位峰度排序结果会变得完全没有意义。另一个细节是峰度计算只用residuals[max_lag:, i]排除前 max_lag 行的无效残差。这行切片在原文里反复出现所有基于残差的统计量都必须保持同样的起点否则前段零值会污染分布形态。4.2 PCMCI 剪边函数的实现与局限PCMCI 在这份代码里的角色是“边修剪器”。它接收数据、滞后阶数和显著性水平返回一张布尔邻接矩阵表示变量间候选边的去留from scipy.stats import pearsonr def pcmci_plus(data, max_lag, alpha0.05): n_samples, n_variables data.shape graph np.ones((n_variables, n_variables), dtypebool) for i in range(n_variables): for j in range(n_variables): if i j: graph[i, j] False else: graph[i, j] True for i in range(n_variables): for j in range(n_variables): if graph[i, j]: p_value, _ pearsonr(data[max_lag:, i], data[max_lag:, j]) if p_value alpha: graph[i, j] False return graph逻辑说明graph[i, j] 为 True 表示保留从 i 到 j 的边为 False 表示删除。Pearson 相关在这里同时检验 i 和 j 的同周期线性关联p 值大于 alpha 说明“不能拒绝无关假设”于是删边。这个实现是高度压缩版的 PCMCI真正的 PCMCI 需要做多重条件集搜索和滞后条件检验工程上可以替换成 tigramite 库里的原版函数但论文复现场景下这个简化版本已经足够体现混合算法的骨架逻辑。4.3 NBCB 与 CBNB 的完整组合流程NBCB 的组合方式是“先排顺序、再剪边”。Causal order 确定后通过双层循环把排在后面的变量指向前面变量的边全部删除确保图中不出现反向边def nbcb(data, max_lag, alpha0.05): # 先用噪声派确定因果顺序 causal_order varlingam(data, max_lag) # 再用约束派剪掉不显著的边 graph pcmci_plus(data, max_lag, alpha) # 根据因果顺序修剪边 for i in range(len(causal_order)): for j in range(i 1, len(causal_order)): graph[causal_order[i], causal_order[j]] False return graph逻辑说明causal_order 里越靠前的变量越早被确定为外生变量即“原因优先出现”。因此对于任意一对变量 i j正确的方向应该是causal_order[i]指向causal_order[j]反向边应当被删除。代码里把前者的边置为 False 是一种保守做法——在骨架不确定时宁可不输出方向也不输出错误方向。我在这份代码的复现里做了一个小改动按论文原始逻辑NBCB 最后这段方向修剪存在歧义有人解读为“仅保留从因到果的边”有人解读为“删掉全部跨顺序边”。我采用的策略是前者因为它和 CBNB 的方向逻辑正好互补对比实验结果更干净。CBNB 的组合流程反过来先剪边再定向def cbnb(data, max_lag, alpha0.05): # 先用约束派得到骨架 graph pcmci_plus(data, max_lag, alpha) # 再用噪声派确定因果顺序 causal_order varlingam(data, max_lag) # 根据因果顺序定向边 for i in range(len(causal_order)): for j in range(i 1, len(causal_order)): if graph[causal_order[i], causal_order[j]]: graph[causal_order[i], causal_order[j]] True graph[causal_order[j], causal_order[i]] False return graph逻辑说明骨架阶段保留下来的边可能是无向的这一步只对骨架中确实存在的边做定向。graph[causal_order[i], causal_order[j]]置 True 表示确认因到果方向反向位置置 False。如果原本该边因不显著已被删除则跳过不做任何操作。参数说明两个组合函数共享 alpha 和 max_lagalpha 控制约束派剪边的宽松程度max_lag 同时影响回归特征矩阵的构造和残差的有效样本起点。两个参数一动后面所有输出都会变。4.4 跑通全流程的测试入口# 生成数据 n_samples 1000 n_variables 5 max_lag 2 data generate_time_series(n_samples, n_variables, max_lag) # 运行NBCB算法 nbcb_graph nbcb(data, max_lag) # 运行CBNB算法 cbnb_graph cbnb(data, max_lag) print(NBCB Graph:) print(nbcb_graph) print(CBNB Graph:) print(cbnb_graph)逻辑说明这里固定 n_samples1000、n_variables5、max_lag2是把复现成本压到最低的配置。n_variables5 意味着因果顺序搜索空间只有 5 个变量峰度排序循环的次数很少任何一台笔记本都能秒级跑完。输出两张 5x5 的布尔矩阵行表示原因列表示结果。我建议拿到输出后先对比两张图的差异集中在哪些变量对上。如果 NBCB 和 CBNB 对同一对变量的方向判断相反说明这两个变量的信噪比过低或者存在隐藏的第三变量共同驱动。这种情况下不是算法 bug而是数据本身就不适合做因果定向需要回到数据生成器增大信号强度或增加样本量。5. 避坑这版代码跑通后最该检查的五个细节5.1 方向完全反了NBCB 的修剪逻辑和直觉不一样现象NBCB 输出结果里因果方向看起来和 CBNB 完全相反甚至出现“结果指向原因”的边。原因原始代码里graph[causal_order[i], causal_order[j]] False把因到果的边删掉了保留的是果到因的反向候选边。这是论文代码里最容易误读的一行。不同复现者对“顺序排列方向”的理解不同有人把 causal_order[0] 当作最早的原因有人把它当作最下游的结果两种解读下同一行代码的效果截然相反。解决先打印causal_order再用模拟数据里已知的真实系数方向对照。我的习惯是直接用 CBNB 的方向作为基准——CBNB 里graph[causal_order[i], causal_order[j]] True是明确保留因到果的边语义清晰。对齐 NBCB 时把修剪循环改成只删反向边即可。5.2 峰度判据选 argmin 还是 argmax结果天差地别现象换了数据生成器的随机种子后因果顺序频繁跳动甚至直接反转。原因当残差分布接近高斯时峰度值非常接近argmin 和 argmax 的区分度都很差。论文代码用 argmin即挑选峰度最小的变量作为最先确定的顺序节点这在误差项为超高斯分布时合理但如果噪声实际是亚高斯峰度小于 3排序结果会变得很脆。解决先对残差分布做一次峰度统计观察数值范围。如果所有峰度值都在 2.8 到 3.2 之间说明数据接近高斯噪声派方法本身就不适用。这时要么增大噪声的非高斯性生成数据时用指数分布或 t 分布噪声要么换用基于熵的独立性度量替代峰度。5.3 max_lag 设错所有残差和相关性检验全部失真现象算法运行不报错但输出图与真实结构严重不符而且无论怎么调 alpha 都没有改善。原因max_lag 决定回归的滞后特征数。设为 0 时模型变成纯截面回归完全丢掉了时序信息设为过大时特征矩阵膨胀线性回归在小样本上过拟合残差被过度解释。解决用数据生成器里的真实滞后阶数作为基准。工程数据没有先验时常见做法是用 Ljung-Box 检验或偏自相关函数判断最大滞后阶也可以直接跑一个简单的时序交叉验证观察不同 max_lag 下预测误差最低的值。论文复现场景下我建议固定 max_lag2这是模拟数据 SCM 的标准配置。5.4 PCMCI 剪边不彻底假边残留现象输出骨架里始终存在一对变量在所有条件下都显著相关但真实因果图里它们并没有直接边。原因简化版 pcmci_plus 只做双变量 Pearson 检验没有做条件独立性检验。两个变量中间隔着一个共同原因变量时它们的无条件相关系数必然显著但给定共同原因后边应该消失。简化实现没有这一步假边自然残留。解决如果只是复现论文可以在分析时忽略这类残留如果做严肃实验需要把 pcmci_plus 替换成真正的多条件检验。最少工作量的方案是引入 tigramite 库的 ParCorr 条件独立性测试它原生支持滞后条件集效果和原版 PCMCI 接近。5.5 数据生成器跑出来的结果不可复现现象每次运行生成的数据不一样算法输出也不一样无法和论文结果对比。原因np.random.uniform和np.random.normal都没有固定随机种子。模拟实验需要可复现性随机种子不固定等于每次都在新数据集上测试。解决在脚本开头加np.random.seed(42)。更严谨的做法是直接在 generate_time_series 函数里加一个 seed 参数默认设为固定值这样在不同机器上跑同一份代码能得到完全一致的系数矩阵和噪声序列。6. 进阶用法非线性条件独立检验与隐藏混杂因子处理6.1 用高斯过程回归把 CI 测试扩展到非线性场景线性 SCM 生成的模拟数据里所有关系都是一阶线性滞后Pearson 相关够用。真实业务数据的因果关系往往是非线性的——比如气温对负荷的影响存在饱和效应线性相关无法捕捉这种局部依赖。一个可操作的升级方式是用高斯过程回归替代线性回归把残差中的非线性依赖剥出来再检验独立性from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, ConstantKernel from scipy.stats import pearsonr def nonlinear_conditional_independence_test(X, Y, Z, alpha0.05): # 高斯过程回归核函数常数核乘RBF kernel ConstantKernel(1.0) * RBF(length_scale1.0) gp GaussianProcessRegressor(kernelkernel) # 用条件集Z分别回归Y和X取残差 gp.fit(Z, Y) residuals_Y Y - gp.predict(Z) gp.fit(Z, X) residuals_X X - gp.predict(Z) # 残差相关性显著则拒绝条件独立 p_value, _ pearsonr(residuals_X, residuals_Y) return p_value alpha # 返回True表示条件独立逻辑说明这个函数的核心逻辑是“用条件集 Z 把 X 和 Y 中能被 Z 解释的部分都去掉剩下的残差如果还显著相关说明 X 和 Y 之间存在 Z 解释不了的关系”。高斯过程回归在这里的作用是拟合非线性映射比线性回归更能捕捉 Z 对 X、Y 的复杂作用。需要说明的是严格的条件独立性检验还需要考虑残差的独立性和分布假设这里用 Pearson 相关是一种工程近似。参数说明ConstantKernel(1.0) * RBF(length_scale1.0)是高斯过程的默认核配置。length_scale 控制拟合的平滑程度值越小拟合越细但过小会过拟合噪声。真实数据上可以调成 0.1 到 10 之间观察输出稳定性但复现论文时保持默认即可。6.2 FCI 的简化实现与隐藏混杂因子边界from itertools import combinations def fci(data, max_lag, alpha0.05): n_samples, n_variables data.shape graph np.ones((n_variables, n_variables), dtypebool) for i in range(n_variables): for j in range(n_variables): if i j: graph[i, j] False else: graph[i, j] True # 在候选边存在时遍历第三个变量作为条件集 for i in range(n_variables): for j in range(n_variables): if graph[i, j]: for k in range(n_variables): if k ! i and k ! j: p_value, _ pearsonr(data[max_lag:, i], data[max_lag:, j]) if p_value alpha: graph[i, j] False break return graph逻辑说明这段代码是标准的 FCI 简化框架。真正的 FCI 算法会系统性地枚举所有可能的条件集判断 X 与 Y 是否在某个条件下独立输出部分有向无环图。这里的简化版本只遍历单变量条件集好处是计算量从指数级降到 O(n³)坏处是包含两个以上混杂因子时无法正确剪边。论文最后给出的hybrid_causal_discovery函数是把 6.1 节的条件独立测试和 FCI 的骨架搜索结合到一起替换掉 Pearson 相关就可以在同一框架下处理非线性关系和隐藏混杂因子。我实际跑下来最明显的感受是函数本身不复杂复杂的是判断“什么时候该相信骨架、什么时候该相信顺序”——两条路径输出一致时结果可信度高输出矛盾时大概率是数据质量或假设条件出了问题。从那以后我每次跑混合因果发现都会强制走一遍固定随机种子、打印因果顺序、对比 NBCB 与 CBNB 方向差集这三个步骤已经形成肌肉记忆。这个习惯帮我避开过至少三次“看似跑通、实则方向完全反了”的假成功希望帮到你。本文还有配套的精品资源点击获取