简介面向时间序列因果推断领域的科研人员和数据科学从业者这份压缩包针对单一方法难以应对复杂因果关系的问题聚焦基于噪声与基于约束两类算法的混合策略完整呈现NBCB与CBNB两条技术路线的复现流程。资源为1个docx文档约25KB虽体积精巧但内容密度高涵盖环境配置与工具包安装命令、线性动态结构因果模型SCM的数据生成、VarLiNGAM在因果顺序发现中的应用以及PCMI在边修剪中的条件独立性检验实现每一步都配有可运行代码与解释。文档还讨论了复杂条件独立性检测、并行化加速等优化手段并延伸至非线性因果处理和隐混淆变量应对方法。对希望从原理走向编程落地、理解混合因果发现内在逻辑的中高级读者尤其友好目前已有71人学习。借助文档中的分步推导、示例数据与结果对照读者能快速搭建自己的实验框架并复现核心算法。1. 为什么同时用噪声和约束混合因果发现的两条路线拿到一段时间序列直接算相关矩阵或跑LSTM得到的只是关联而不是因果。基于约束的方法比如PC、PCMCI这类擅长用条件独立性检验把伪相关边剪掉却经常在方向判定上翻车基于噪声的方法比如LiNGAM家族依赖非高斯噪声的结构强行识别方向对条件集的选择和样本量又极其敏感。NBCBNoise-Based-then-Constraint-Based与CBNBConstraint-Based-then-Noise-Based正是两种把这两派串成流水线的混合算法前者先用噪声方法定因果顺序再用约束方法删冗余边后者反过来先剪骨架后定向。这套实现覆盖了从数据生成、VarLiNGAM、PCMCI到非线性扩展和隐藏混杂处理的完整代码适合做论文基线对比、真实场景前的仿真验证以及想摸清两个流派边界和坑位的工程师。2. 从滞后线性SCM开始时间序列数据生成与参数设计2.1 结构因果模型如何落到差分方程时间序列里的结构因果模型SCM通常会写成一个带滞后的线性系统当前时刻的每个变量由过去若干个时刻所有变量的线性组合加上一个噪声项决定。这个设定既符合大多数时序因果发现论文的默认假设也让后面的VarLiNGAM和PCMCI能直接套用线性回归与相关检验。生成代码里的系数张量是coefficients[i, j, lag - 1]含义是变量 j 在 t-lag 时刻的值对变量 i 在 t 时刻的贡献。也就是说第一维是目标变量 i、第二维是源变量 j、第三维是滞后阶。这个维度约定后面会直接影响因果图矩阵的读写方向建议在读代码时先固定下来不然后面NBCB的方向修正容易绕晕。需要先消歧两个词这里的“约束”指统计学里的条件独立性约束和FPGA时序约束sdc/xdc里那些时钟mux约束、IO约束完全是两回事这里的“混合算法”也跟混合整数线性规划MILP无关。搜索时经常撞名但技术栈不在一个领域。2.2 generate_time_series 与滞后对齐import numpy as np def generate_time_series(n_samples, n_variables, max_lag, noise_std0.1): # 系数在[-1,1]均匀采样弱边直接置零得到一个稀疏的真实因果结构 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循环从t max_lag开始是因为前 max_lag 个时间点没有完整的滞后历史可供计算只能作为边界丢在一边。三重循环的写法在 n_samples1000、n_variables5、max_lag2 时只有五万次乘加跑起来秒级完成但变量数涨到 20 以上时三层循环加滞后会明显变慢这时候要么把系数矩阵乘法向量化要么用 numba 编译内层循环。coefficients[np.abs(coefficients) 0.1] 0是一个 L0 风格的正则化约束目的是让生成数据背后的因果图保持稀疏。注意这一步也会把一部分本来就不算强的真实边直接抹掉所以生成数据的“真实结构”要以置零后的系数为准后面算 SHD 评估时用的 ground truth 也要从这里取。2.3 参数表与实验配置参数取值作用调参建议n_samples1000样本量决定回归与检验稳定性低于 200 时 Pearson 检验方差变大边剪不稳n_variables5变量数决定因果图规模超过 10 时朴素 PCMCI 的检验次数指数膨胀max_lag2最大滞后期用 BIC/AIC 或交叉验证选参考 5.3noise_std0.1加性噪声标准差控制信噪比噪声过大时真边容易被误删稀疏阈值0.1弱系数直接置零调大图更稀疏调小更稠密生成结束后真实因果图可以通过np.abs(coefficients).sum(axis2) 0恢复出来每个有非零系数的变量对就是一条真实因果边。注意这里得到的是跨滞后汇总的图也就是把 max_lag2 时刻的两条滞后边合成了一条变量级边后续算法输出的也是这种汇总图评估时两类图要对齐语义。3. VarLiNGAM与PCMCI两个核心模块的工程实现3.1 VarLiNGAM的残差峰度迭代LiNGAM这一支方法的核心假设是数据由线性结构方程生成噪声是非高斯的且因果图无环。在这个假设下果变量对因变量做回归后残差中如果还残留其他变量的信号残差的非高斯性会变强反之真正外生变量的残差只包含它自己的噪声。因此迭代选出“残差峰度最小”的变量作为当前最外生变量再逐步从剩余变量中剥离就能恢复因果顺序。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().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这段代码有两点值得细看。第一峰度计算只取residuals[max_lag:, i]是因为前 max_lag 行残差是无人认领的 0会把峰度直接拉向负值这个切片是所有残差类方法都要做的对齐操作。第二np.argmin选的是峰度最小的变量而这只是论文复现里的一种写法。多数 LiNGAM 开源实现比如 lingam 库里的 DirectLiNGAM倾向于用np.argmax(np.abs(kurtosis))也就是优先提取非高斯性最强的分量因为 LiNGAM 的理论起点恰恰是“噪声非高斯”。建议两种判据都跑一遍用 5.3 的 SHD 去对照真实图哪个中位数低就固定哪个。提示scipy.stats.kurtosis默认返回的是超额峰度正态分布对应 0。如果数据噪声形态接近高斯峰度判据会非常不稳定这时可以换成scipy.stats.normaltest的 p 值作为选择依据逻辑是从“峰度最小”改成“正态性最显著”。3.2 自回归特征矩阵的构造细节3.1 的代码里有一个关键操作容易被一眼带过X np.hstack([data[max_lag - lag:-lag, :] for lag in range(1, max_lag 1)])。这个拼接构造出的是自回归特征矩阵每一行对应一个时间点 t包含该时刻所有变量的前 max_lag 个滞后观测。X_lag1 data[max_lag - 1:-1, :] # 每个 t 对应的 t-1 时刻 X_lag2 data[max_lag - 2:-2, :] # 每个 t 对应的 t-2 时刻 X np.hstack([X_lag1, X_lag2]) # 行数 n_samples - max_lag行对齐关系是这样的X[0]对应 t max_lag此时X_lag1[0]恰好是 data[max_lag - 1] 即 t-1X_lag2[0]恰好是 data[max_lag - 2] 即 t-2。所以不需要按时间索引循环对齐行序天然对应。这个特征矩阵的维度是(n_samples - max_lag, n_variables * max_lag)和data[max_lag:, i]的长度完全匹配。工程上这里有个小优化原始代码把 X 的构造放在每个变量 i 的循环内部等于同样的矩阵拼了 n_variables 次。数据量小无所谓但变量数一涨建议把 X 提出来只算一次。另外如果后续要做滞后阶选择这个 X 的列顺序是“先 lag1 的全部变量再 lag2 的全部变量”不是按变量分组的回归系数的解释要按这个顺序切分。3.3 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) np.fill_diagonal(graph, False) 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完整的 PCMCI 分两个阶段先用 PC1 在各滞后时刻做条件独立筛选找出每个变量的候选父母集再用 MCIMomentary Conditional Independence做最终检验以候选父母集作为条件集来控制高维下的假阳性。这份复现代码为了可读性退化成“所有变量对做无条件 Pearson 检验”p 值大于 alpha 就把边删掉实际上只保留了 PCMCI 第二阶段的简化形式。alpha0.05 在 1000 样本下通常会留下不少弱边alpha 降到 0.01 又会误删真边真实 PCMCI 通过 MCI 的父集条件缓解的正是这个矛盾。初始化部分也可以简化graph np.ones(...)之后只需要把对角线置 False原代码里那段 for 循环的 else 分支每次都在赋 True属于冗余操作等价于什么都没做。3.4 条件独立 vs 无条件独立什么时候该换检验检验方式条件集能剔除的伪相关计算成本适用场景无条件 Pearson空集线性相关性弱的伪边O(V^2)快速扫描、基线对照条件 PearsonPC 风格单个或一组变量公共原因、中间变量导致的伪相关O(V^3 * 条件集大小)标准因果发现高斯过程残差检验连续多个变量非线性伪相关最高涉及多次 GP 拟合非线性数据见第 5 章无条件检验的问题在于它分不清“X 与 Y 相关”是直接因果还是经由第三个变量传递。比如 X 影响 Z、Z 影响 Y无条件相关会把 X 和 Y 也连上条件独立检验加入 Z 作为条件集后就能把这条伪边拆掉。PCMCI 的 MCI 阶段本质上就是在做这件事所以后续要上真实数据第一个升级动作就是把这个无条件 Pearson 替换成偏相关或条件独立检验。4. NBCB与CBNB的组合逻辑、方向修正与输出解读4.1 两条流水线的分工NBCB 的流程是先用 VarLiNGAM 算出全变量因果顺序再用 PCMCI 删掉条件独立意义下的冗余边最后按因果顺序统一方向。约束方法在这里只承担剪枝不承担定向方向完全由噪声方法给出的全序决定。CBNB 恰好反过来先让 PCMCI 剪出一个无向骨架然后只在骨架保留的候选边上用 VarLiNGAM 给出的顺序去定向。噪声方法的全局分布假设即使局部失效骨架还在兜底不会把整张图带偏。维度NBCBCBNB第一步VarLiNGAM 定全序PCMCI 剪无向骨架第二步PCMCI 删冗余边VarLiNGAM 在剩余边内定向主要风险排序错了后面无法纠正骨架漏边则定向无从谈起适合场景非高斯性显著、变量数少弱相关多、需要保守剪枝两种方法的共同点是“两个模块的结果按流水线叠加”不是各出一张图再合并。这也是混本文还有配套的精品资源点击获取 
长河君 · 2026/9/19 15:27:54 · 阅读 1.2k 噪声与约束融合的混合因果发现：NBCB与CBNB实现解析
简介面向时间序列因果推断领域的科研人员和数据科学从业者这份压缩包针对单一方法难以应对复杂因果关系的问题聚焦基于噪声与基于约束两类算法的混合策略完整呈现NBCB与CBNB两条技术路线的复现流程。资源为1个docx文档约25KB虽体积精巧但内容密度高涵盖环境配置与工具包安装命令、线性动态结构因果模型SCM的数据生成、VarLiNGAM在因果顺序发现中的应用以及PCMI在边修剪中的条件独立性检验实现每一步都配有可运行代码与解释。文档还讨论了复杂条件独立性检测、并行化加速等优化手段并延伸至非线性因果处理和隐混淆变量应对方法。对希望从原理走向编程落地、理解混合因果发现内在逻辑的中高级读者尤其友好目前已有71人学习。借助文档中的分步推导、示例数据与结果对照读者能快速搭建自己的实验框架并复现核心算法。1. 为什么同时用噪声和约束混合因果发现的两条路线拿到一段时间序列直接算相关矩阵或跑LSTM得到的只是关联而不是因果。基于约束的方法比如PC、PCMCI这类擅长用条件独立性检验把伪相关边剪掉却经常在方向判定上翻车基于噪声的方法比如LiNGAM家族依赖非高斯噪声的结构强行识别方向对条件集的选择和样本量又极其敏感。NBCBNoise-Based-then-Constraint-Based与CBNBConstraint-Based-then-Noise-Based正是两种把这两派串成流水线的混合算法前者先用噪声方法定因果顺序再用约束方法删冗余边后者反过来先剪骨架后定向。这套实现覆盖了从数据生成、VarLiNGAM、PCMCI到非线性扩展和隐藏混杂处理的完整代码适合做论文基线对比、真实场景前的仿真验证以及想摸清两个流派边界和坑位的工程师。2. 从滞后线性SCM开始时间序列数据生成与参数设计2.1 结构因果模型如何落到差分方程时间序列里的结构因果模型SCM通常会写成一个带滞后的线性系统当前时刻的每个变量由过去若干个时刻所有变量的线性组合加上一个噪声项决定。这个设定既符合大多数时序因果发现论文的默认假设也让后面的VarLiNGAM和PCMCI能直接套用线性回归与相关检验。生成代码里的系数张量是coefficients[i, j, lag - 1]含义是变量 j 在 t-lag 时刻的值对变量 i 在 t 时刻的贡献。也就是说第一维是目标变量 i、第二维是源变量 j、第三维是滞后阶。这个维度约定后面会直接影响因果图矩阵的读写方向建议在读代码时先固定下来不然后面NBCB的方向修正容易绕晕。需要先消歧两个词这里的“约束”指统计学里的条件独立性约束和FPGA时序约束sdc/xdc里那些时钟mux约束、IO约束完全是两回事这里的“混合算法”也跟混合整数线性规划MILP无关。搜索时经常撞名但技术栈不在一个领域。2.2 generate_time_series 与滞后对齐import numpy as np def generate_time_series(n_samples, n_variables, max_lag, noise_std0.1): # 系数在[-1,1]均匀采样弱边直接置零得到一个稀疏的真实因果结构 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循环从t max_lag开始是因为前 max_lag 个时间点没有完整的滞后历史可供计算只能作为边界丢在一边。三重循环的写法在 n_samples1000、n_variables5、max_lag2 时只有五万次乘加跑起来秒级完成但变量数涨到 20 以上时三层循环加滞后会明显变慢这时候要么把系数矩阵乘法向量化要么用 numba 编译内层循环。coefficients[np.abs(coefficients) 0.1] 0是一个 L0 风格的正则化约束目的是让生成数据背后的因果图保持稀疏。注意这一步也会把一部分本来就不算强的真实边直接抹掉所以生成数据的“真实结构”要以置零后的系数为准后面算 SHD 评估时用的 ground truth 也要从这里取。2.3 参数表与实验配置参数取值作用调参建议n_samples1000样本量决定回归与检验稳定性低于 200 时 Pearson 检验方差变大边剪不稳n_variables5变量数决定因果图规模超过 10 时朴素 PCMCI 的检验次数指数膨胀max_lag2最大滞后期用 BIC/AIC 或交叉验证选参考 5.3noise_std0.1加性噪声标准差控制信噪比噪声过大时真边容易被误删稀疏阈值0.1弱系数直接置零调大图更稀疏调小更稠密生成结束后真实因果图可以通过np.abs(coefficients).sum(axis2) 0恢复出来每个有非零系数的变量对就是一条真实因果边。注意这里得到的是跨滞后汇总的图也就是把 max_lag2 时刻的两条滞后边合成了一条变量级边后续算法输出的也是这种汇总图评估时两类图要对齐语义。3. VarLiNGAM与PCMCI两个核心模块的工程实现3.1 VarLiNGAM的残差峰度迭代LiNGAM这一支方法的核心假设是数据由线性结构方程生成噪声是非高斯的且因果图无环。在这个假设下果变量对因变量做回归后残差中如果还残留其他变量的信号残差的非高斯性会变强反之真正外生变量的残差只包含它自己的噪声。因此迭代选出“残差峰度最小”的变量作为当前最外生变量再逐步从剩余变量中剥离就能恢复因果顺序。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().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这段代码有两点值得细看。第一峰度计算只取residuals[max_lag:, i]是因为前 max_lag 行残差是无人认领的 0会把峰度直接拉向负值这个切片是所有残差类方法都要做的对齐操作。第二np.argmin选的是峰度最小的变量而这只是论文复现里的一种写法。多数 LiNGAM 开源实现比如 lingam 库里的 DirectLiNGAM倾向于用np.argmax(np.abs(kurtosis))也就是优先提取非高斯性最强的分量因为 LiNGAM 的理论起点恰恰是“噪声非高斯”。建议两种判据都跑一遍用 5.3 的 SHD 去对照真实图哪个中位数低就固定哪个。提示scipy.stats.kurtosis默认返回的是超额峰度正态分布对应 0。如果数据噪声形态接近高斯峰度判据会非常不稳定这时可以换成scipy.stats.normaltest的 p 值作为选择依据逻辑是从“峰度最小”改成“正态性最显著”。3.2 自回归特征矩阵的构造细节3.1 的代码里有一个关键操作容易被一眼带过X np.hstack([data[max_lag - lag:-lag, :] for lag in range(1, max_lag 1)])。这个拼接构造出的是自回归特征矩阵每一行对应一个时间点 t包含该时刻所有变量的前 max_lag 个滞后观测。X_lag1 data[max_lag - 1:-1, :] # 每个 t 对应的 t-1 时刻 X_lag2 data[max_lag - 2:-2, :] # 每个 t 对应的 t-2 时刻 X np.hstack([X_lag1, X_lag2]) # 行数 n_samples - max_lag行对齐关系是这样的X[0]对应 t max_lag此时X_lag1[0]恰好是 data[max_lag - 1] 即 t-1X_lag2[0]恰好是 data[max_lag - 2] 即 t-2。所以不需要按时间索引循环对齐行序天然对应。这个特征矩阵的维度是(n_samples - max_lag, n_variables * max_lag)和data[max_lag:, i]的长度完全匹配。工程上这里有个小优化原始代码把 X 的构造放在每个变量 i 的循环内部等于同样的矩阵拼了 n_variables 次。数据量小无所谓但变量数一涨建议把 X 提出来只算一次。另外如果后续要做滞后阶选择这个 X 的列顺序是“先 lag1 的全部变量再 lag2 的全部变量”不是按变量分组的回归系数的解释要按这个顺序切分。3.3 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) np.fill_diagonal(graph, False) 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完整的 PCMCI 分两个阶段先用 PC1 在各滞后时刻做条件独立筛选找出每个变量的候选父母集再用 MCIMomentary Conditional Independence做最终检验以候选父母集作为条件集来控制高维下的假阳性。这份复现代码为了可读性退化成“所有变量对做无条件 Pearson 检验”p 值大于 alpha 就把边删掉实际上只保留了 PCMCI 第二阶段的简化形式。alpha0.05 在 1000 样本下通常会留下不少弱边alpha 降到 0.01 又会误删真边真实 PCMCI 通过 MCI 的父集条件缓解的正是这个矛盾。初始化部分也可以简化graph np.ones(...)之后只需要把对角线置 False原代码里那段 for 循环的 else 分支每次都在赋 True属于冗余操作等价于什么都没做。3.4 条件独立 vs 无条件独立什么时候该换检验检验方式条件集能剔除的伪相关计算成本适用场景无条件 Pearson空集线性相关性弱的伪边O(V^2)快速扫描、基线对照条件 PearsonPC 风格单个或一组变量公共原因、中间变量导致的伪相关O(V^3 * 条件集大小)标准因果发现高斯过程残差检验连续多个变量非线性伪相关最高涉及多次 GP 拟合非线性数据见第 5 章无条件检验的问题在于它分不清“X 与 Y 相关”是直接因果还是经由第三个变量传递。比如 X 影响 Z、Z 影响 Y无条件相关会把 X 和 Y 也连上条件独立检验加入 Z 作为条件集后就能把这条伪边拆掉。PCMCI 的 MCI 阶段本质上就是在做这件事所以后续要上真实数据第一个升级动作就是把这个无条件 Pearson 替换成偏相关或条件独立检验。4. NBCB与CBNB的组合逻辑、方向修正与输出解读4.1 两条流水线的分工NBCB 的流程是先用 VarLiNGAM 算出全变量因果顺序再用 PCMCI 删掉条件独立意义下的冗余边最后按因果顺序统一方向。约束方法在这里只承担剪枝不承担定向方向完全由噪声方法给出的全序决定。CBNB 恰好反过来先让 PCMCI 剪出一个无向骨架然后只在骨架保留的候选边上用 VarLiNGAM 给出的顺序去定向。噪声方法的全局分布假设即使局部失效骨架还在兜底不会把整张图带偏。维度NBCBCBNB第一步VarLiNGAM 定全序PCMCI 剪无向骨架第二步PCMCI 删冗余边VarLiNGAM 在剩余边内定向主要风险排序错了后面无法纠正骨架漏边则定向无从谈起适合场景非高斯性显著、变量数少弱相关多、需要保守剪枝两种方法的共同点是“两个模块的结果按流水线叠加”不是各出一张图再合并。这也是混本文还有配套的精品资源点击获取