原作者：Erdogan T

贝叶斯统计拥有极其广泛的应用领域。然而，在刚开始接触时，人们往往很难理解不同技术之间如何对应不同的问题与应用场景。

在这篇文章中，我将带你整体了解贝叶斯方法的应用版图，并说明各种应用如何对应不同的因果发现方法。

我们正逐渐进入因果发现研究的一个关键时期：我们拥有海量数据；拥有可扩展的计算能力；同时也拥有强大的贝叶斯分析框架。

现在正是深入学习贝叶斯统计的最佳时机。

换句话说，我们希望回答以下问题：

• 如何利用离散数据或连续数据构建一个因果网络（Directed Acyclic Graph，DAG，有向无环图）？
• 在存在或不存在响应变量（response）或处理变量（treatment）的情况下，是否能够确定因果结构？
• 应该如何选择不同的搜索算法，例如 PC 算法、Hill Climb Search（爬山搜索） 等？

在阅读完本文之后，你将能够理解：

如何为自己的具体应用场景选择合适的贝叶斯因果发现技术，并迈出实践的第一步。

请放慢节奏，泡一杯茶，享受阅读过程。

如果你觉得这篇文章对你有所帮助，那将是我的荣幸。

我也强烈建议你亲自尝试文中的实践示例。这将帮助你：

• 学得更快，
• 理解得更深入，
• 记忆得更持久。

准备好你的茶，开始探索吧！

简要历史

在数据科学领域中，贝叶斯方法正逐渐成为“新晋热门技术”。

为了理解贝叶斯统计未来的发展方向，了解其历史背景是非常有帮助的。

事实上，贝叶斯统计已经存在了相当长的时间。

故事要从 托马斯·贝叶斯（Thomas Bayes） 说起。他在一篇论文中提出了后来被称为**贝叶斯定理（Bayes’ Theorem）**的思想，并于 1763 年发表该成果 [1]。

这项工作奠定了现代贝叶斯统计的基础。

然而，在随后的很长一段时间里，贝叶斯方法几乎被学术界忽视。直到 20 世纪 50 至 60 年代，马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC） 被提出。

MCMC 的出现具有决定性意义，因为它提供了一种近似计算方法，使我们能够在以下情况下进行概率估计：

当精确求解复杂概率分布变得困难甚至不可能时。

时间飞逝。即使在贝叶斯最初工作提出 263 年之后，我们仍然在其理论基础之上不断发展新的方法与技术。

什么是贝叶斯定理？

贝叶斯定理（Bayes’ Theorem） 是贝叶斯网络（Bayesian Networks）的理论基础。

贝叶斯公式本身用于更新模型中的信息或信念，其基本形式如下：

该公式包含四个重要的概率部分。下面我们逐一进行说明。

1. 先验概率（Prior Probability）

最常被使用的概率是 先验概率（prior），也称为信念概率（belief probability）。

这是人们在日常生活中经常无意识使用的一种概率概念，它表示：

在观察到证据之前，对某个假设的初始判断。

换句话说，它代表我们基于经验、直觉或历史信息所形成的认知。

例如：

• 医生在尚未获得实验室检测结果之前，可能根据患者的年龄和症状判断其患某种疾病的概率为 20%；
• 又或者在阴天时，你凭经验认为有 80% 的概率会下雨，于是决定带上雨伞。

2. 条件概率 / 似然（Likelihood）

公式中的另一部分是 条件概率（conditional probability），也称为似然（likelihood）。

它表示：

在假设成立的情况下，观察到当前证据的概率。

这一概率通常来源于数据本身。

3. 边际概率（Marginal Probability）

接下来是 边际概率（marginal probability）。

它描述的是：

在所有可能假设下，新证据出现的总体概率。

这一项通常需要通过对所有可能情况进行计算或求和得到。

4. 后验概率（Posterior Probability）

最后得到的是 后验概率（posterior probability）：P(Z∣X)

即：

在观察到证据 X 之后，事件 Z 发生的概率。

简单来说：后验=先验×数据证据÷整体可能性

走向流行的漫长道路

“贝叶斯还是非贝叶斯（频率学派）？”

这个问题至今仍在统计学界持续争论。

究竟哪一种方法更好？

首先需要明确的是：

两种方法都能够从数据中发现模式。

但它们之间存在一个关键差异：只有贝叶斯统计能够显式地纳入先验知识。

这意味着：研究者可以在分析开始之前加入合理假设，从而：

• 获得更好的起点；
• 提高分析的稳定性与可靠性。

因此，贝叶斯方法在某些领域尤为强大，例如：医学领域。

在医学研究中，经常存在：

• 数据缺失，
• 疾病信息不完整，
• 治疗数据有限。

但医生通常拥有丰富的领域经验，而这些经验可以自然地作为先验信息纳入模型。

相比之下：频率学派方法（Frequentist approach） 通常要求在实验开始之前就固定实验设计，而贝叶斯分析则更加灵活，不需要提前完全确定分析结构。

为什么频率学派长期占主导？

在很长一段时间内，频率学派方法占据主导地位，主要原因是：

它们在大样本情况下计算推断更加容易。

原作者本人在 2010 年代早期从事全基因组分子数据研究时也经历过这一问题。

当时的研究目标是：

• 建模成千上万个特征（基因表达数据），
• 分析治疗相关的复杂因果关系。

然而，这类问题通常需要：对高维概率分布进行复杂积分计算

在当时：

• 数据规模巨大，
• 计算能力有限。

因此，更现实的选择是采用具有**解析解（closed-form solutions）**的统计方法。

换句话说：在当时的计算条件下，一些基因组层面的因果问题几乎无法通过贝叶斯方法解决。至少，当时的计算能力尚不足以支持。

贝叶斯方法的重新兴起

如今，贝叶斯技术正在重新获得关注。

但我们仍未完全解决所有问题。

1995–2018 年期间使用贝叶斯统计的医学论文数量增长趋势

图片来源：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6406060/

近年来，计算能力迅速提升，这可能推动贝叶斯方法得到更全面的应用。

然而挑战依然存在。除了计算成本较高之外，下一个关键问题是：

可解释性（Interpretability）

一般而言：无论一种技术、多复杂的方法或统计模型，如果我们无法正确理解、应用或解释它，那么它本身就没有实际意义。

因此，尽管计算能力已经显著提高，贝叶斯方法当前面临的新挑战在于：

如何有效地应用并正确解释贝叶斯统计结果。

贝叶斯思维的兴起

贝叶斯思维是一种非常自然的方式，用于不断更新我们对某一情境的理解。随着我们收集到越来越多的数据，我们的结论也会变得更加准确。

贝叶斯方法的许多核心基础是在过去几个世纪中逐步发展起来的，而且主要诞生于学术研究环境之中。然而，如今的发展已经不仅局限于学术界：

• 开源社区正在持续推动相关技术的发展；
• 大型科技公司也在做出越来越多的贡献。

因此，除了理论基础之外，我们现在正看到越来越多可用于构建贝叶斯模型的框架与工具出现。我们正在迈向因果发现领域的一个“完美风暴”时期：

• 海量数据的出现，
• 可扩展计算能力的提升，
• 强大的贝叶斯分析框架。

这些条件正在同时成熟。

近年来，许多行业一直致力于构建数据驱动（data-driven）解决方案，其中机器学习方法发挥了关键作用。而如今，这一趋势正在进一步发展为：

数据驱动的决策制定（data-driven decision-making）。

数据科学家逐渐认识到：

在模型中引入先验知识能够提升预测能力，并产生更优的结果。

尤其是在可以进行**干预（intervention）**以进一步优化系统的场景中，这一点尤为重要。

一些广为人知的贝叶斯技术包括：

• 统计推断（Statistical Inference）
• 贝叶斯网络用于复杂决策建模
• 贝叶斯优化（Bayesian Optimization）

在接下来的部分中，作者将介绍一些可以通过贝叶斯方法建模的真实世界应用。

因果发现的整体版图

贝叶斯统计的整体框架之所以显得复杂，是因为：

“贝叶斯统计”实际上是一个总称。

它涵盖了大量基于贝叶斯定理构建的统计方法与分析途径。

为了更清晰地组织这些方法，作者采用了数据科学领域中常见的分类方式：

• 监督（Supervised）
• 无监督（Unsupervised）

需要注意的是，在贝叶斯统计领域中，这种分类方式并不总是被正式使用。但作者认为这种划分方式非常直观，因为在现实问题中：

并非总能进行带有目标变量的受控实验（即监督学习）。

从**非受控或观测数据（observational data）**中发现因果关系的方法（即无监督方法）具有重要价值，但它们通常采用与监督方法不同的策略。

因此，每一类问题通常对应：

• 不同的方法，
• 不同的建模思路，
• 不同类型的输入数据。

在本节中，作者将介绍一些能够通过贝叶斯建模解决的现实问题。

从整体上看，贝叶斯策略大致可以分为以下几类：

1. 优化（Optimization）
2. 时间序列（Time Series）
3. 因果发现（监督式）
4. 因果发现（无监督式）

贝叶斯应用与工具库的四大类别：时间序列、优化、监督式因果发现与无监督式因果发现

贝叶斯优化

优化方法在实际应用中具有基础性的重要作用，因为它们能够：

• 显著加快计算速度，
• 降低计算成本。

一个典型例子是 Bayesian Optimization Library，其目标是在尽可能少的迭代次数中找到未知函数的最大值。

另一个重要应用场景是：

机器学习模型的超参数搜索（Hyperparameter Tuning）。

例如在基于树的算法中：

• XGBoost
• LightGBM
• CatBoost

传统的 网格搜索（Grid Search） 方法通常需要大量计算资源与内存，因为它依赖暴力枚举所有参数组合。

而使用贝叶斯优化时，可以更加高效地对超参数进行调优。

一个广泛使用的相关工具库是 Hyperopt。

该库通过贝叶斯优化改进传统网格搜索，并包含三种优化算法：

• Random Search（随机搜索）
• Tree-Structured Parzen Estimator（TPE）
• Adaptive TPE

不过，这些方法通常是**与具体模型无关（method-agnostic）**的，因此仍需要用户手动将其集成到目标模型中。

为了解决这一问题，Hgboost 将 Hyperopt 的优化策略与多种树模型进行了整合，包括：

• XGBoost
• LightGBM
• CatBoost

从而能够更加高效地对这些模型进行自动化参数优化。

贝叶斯时间序列

贝叶斯建模的第二个重要应用领域，是随时间估计某项干预措施所产生的因果效应。

设想这样一个场景：

你是一位企业经营者，刚刚推出了一项市场营销活动，随后你发现销售额出现了明显增长。

这听起来当然是好消息。但真正的问题是：

销售增长究竟是由于营销活动带来的，还是本来就会发生的？

例如，这种增长可能来源于：

• 季节性波动（比如冬天的厚衣服卖得好）
• 天气因素
• 市场周期
• 或仅仅是巧合

这正是贝叶斯时间序列分析发挥作用的地方。

它能够帮助我们：

将真实的干预效果与数据中的随机波动区分开来。

在上述例子中，该方法可以估计：

营销活动的真实影响，以及该活动是否真正导致了销售增长这一结果。

目前已有多种库可以利用贝叶斯统计对时间趋势进行建模，其中两个较为知名的工具是：

• PyMC
• CausalImpact

CausalImpact

CausalImpact 库基于 贝叶斯结构化时间序列（Bayesian Structural Time Series） 方法。

它通过比较：

• 预期时间序列（expected time series）
• 实际观测时间序列（observed time series）

之间的差异，从而分析某项干预措施所产生的影响。

PyMC

相比之下，PyMC 的应用范围更加广泛，可以用于：

• 因果推断（Causal Inference）
• 高斯过程（Gaussian Processes）
• 空间分析（Spatial Analysis）
• 时间序列分析（Time Series Analysis）

贝叶斯时间序列的案例:

这些应用有一个共同点：

在一个已有系统或流程中引入了某种干预（intervention），并观察其随时间产生的变化。

核心问题在于：如何量化该干预的因果影响。

然而，与所有统计方法一样：

强结论依赖于强假设。

因此，并非所有应用场景都适合进行精确量化。

例如，对于以下类型的结果变量：

• “生产效率（productivity）”
• “员工士气（employee morale）”

这类主观或难以测量的指标往往较难进行严格的因果量化。

原因包括：

• 混杂变量（confounding variables）的存在；
• 难以将干预效果与其他影响因素完全分离。

因果发现（监督式）

在数据集中确定变量之间的因果关系，大致可以分为两种形式。

类似于机器学习方法：

• 监督式方法（Supervised）
• 无监督式方法（Unsupervised）

在监督式因果发现中，数据集中通常包含：

• 一个目标变量（target / outcome variable）
• 以及（在某些场景下）一个处理变量（treatment variable）

在这种贝叶斯建模情境下，可以形象地理解为：

你进入犯罪现场时，已经有一个明确的嫌疑人。

你的任务是收集证据，以确认或否定这一假设。

因此，这是一种受控方式来研究因果关系。同时，当引入特定处理变量时，你甚至可以主动影响结果。

一个典型示例来自 DoWhy 项目网站：

研究目标是分析酒店订单取消行为的因果关系。

其中：

• 目标变量（Outcome）：订单取消数量
• 处理变量（Treatment）：“是否分配不同房间（different room assigned）”

研究问题是：

改变房间分配策略是否能够降低订单取消率？

能够处理此类监督式因果发现问题的常见库包括：

• PyMC
• DoWhy

这些工具能够在给定：结果变量（outcome），处理变量（treatment）的情况下估计因果效应。下表列出了一些包含目标变量与处理变量的应用示例。

监督式因果发现示例：

在下一部分中，作者将进一步深入介绍：

如何使用无监督方法进行因果发现。

因果发现（无监督）

寻找因果关系的第四类方法是无监督贝叶斯方法。与监督式方法相比，这类方法通常更加具有挑战性。

我们再次使用之前的“犯罪现场”类比来说明因果发现的过程。

在无监督方法中，你进入犯罪现场时：

并没有任何明确的嫌疑人。

不过，现场遗留下了一些线索（变量），这些线索可能最终指向真正的嫌疑人。

这意味着：你需要分析不同线索之间的组合关系，从而推断出因果关系，并最终找到真正的原因。

于是问题变成：

面对这些留下来的线索，你应该如何进行调查？

你可以采取不同策略：

• 逐个独立检查每条线索，并排除无效信息；
• 或者尝试所有可能的线索组合，并以系统化方式逐一分析。

当线索数量较少时，这项任务相对简单，因为你可以检查所有可能情况，并判断它们是否指向某个嫌疑对象。

然而，当线索数量达到：数百个，甚至数千个，情况就会变得极其耗时，甚至几乎无法完成。

在后续章节中作者将说明：

可组合的数量会随着变量数量呈指数级增长。

这使得在变量过多时寻找真正原因变得非常困难。

幸运的是，我们仍然有解决方案。

关键在于：

需要一种高效的搜索策略（search strategy），能够以智能方式遍历全部可能的线索组合。

无监督贝叶斯方法的目标是：

在没有预先指导信息的情况下，从大量变量中发现隐藏的因果关系。

核心挑战在于：

如何在数量巨大的变量组合空间中，高效识别真正的因果连接。

因此，当以下情况出现时，无监督因果发现尤其适用：

• 无法进行受控实验；
• 只能依赖观测数据（observational data）；
• 仍希望发现潜在的因果结构。

通过这种分析，我们随后可以判断：是否能够对某些因果因素进行干预（intervention）。

然而，如前所述，无监督分析本身具有较高难度，原因包括：

• 数据中不存在明确的响应变量（response variable）或目标变量；
• 数据类型可能多种多样，例如： 连续数据（continuous） 离散数据（discrete） 混合数据（hybrid）

这些都会进一步增加建模复杂度。

常用的无监督因果发现工具

目前较为知名的库包括：

Pgmpy

提供底层概率图模型与统计计算功能。

Bnlearn

提供可直接使用的因果发现流程，包括：

• 结构学习（Structure Learning）
• 参数学习（Parameter Learning）
• 概率推断（Inference）

为了降低理解难度，作者将后续内容分为三种主要方法：

1. 基于约束的方法（Constraint-based）
2. 基于评分的方法（Score-based）
3. 混合结构学习方法（Hybrid methods）

在这些方法中，会频繁出现以下三个核心概念：

结构学习（Structure Learning）

估计最符合数据的有向无环图（Directed Acyclic Graph, DAG），从而刻画变量之间的依赖关系。

参数学习（Parameter Learning）

在已知数据与 DAG 结构的情况下，估计各变量的（条件）概率分布。

概率推断（Probabilistic Inference）

在给定 DAG 结构与模型参数后，计算特定查询问题对应的概率结果。

无监督因果发现的策略

The Strategies for Unsupervised Causal Discovery

在无监督因果发现中，通常存在三种主要策略：

1. 基于约束的方法（Constraint-based）
2. 基于评分的方法（Score-based）
3. 混合策略（Hybrid Strategy）

下图展示了这些策略及其对应的统计方法概览。

每一种策略都具有特定的统计特性，并最终能够构建出一个因果网络（causal network）。

然而，就像机器学习中的无监督方法一样：

因果发现结果的质量，很大程度上取决于所选择的统计方法与数据集特性之间的匹配程度。

同时，也必须始终考虑实际应用目标。

这意味着在现实应用中：

往往需要在方法、计算成本与结果解释性之间做出权衡。

为了决定哪一种统计方法最适合你的数据与应用场景，我们需要理解这三种策略之间的差异与特点。

需要注意的是，下图展示的方法属于 Bnlearn 库的一部分。该框架提供了多种策略、搜索方法以及评分函数，帮助用户为具体问题选择最合适的方法。

因果发现网络的质量取决于统计方法与数据集特性之间的匹配程度。

Bnlearn 库提供现成流程，用于通过结构学习、参数学习与概率推断执行无监督因果发现。

基于约束的方法（Constraint-based，PC算法）

在基于约束的策略中，例如 PC（Peter-Clark）算法，贝叶斯网络的结构是通过检验变量之间的**条件独立性（conditional independence）**来学习得到的。

从直觉上讲，这是一种较为经典的方法：通过逻辑推理与排除法逐步缩小可能性空间。

例如以“喷水器（sprinkler）数据集”为例：

如果草坪一整周都是湿的，但没有证据表明这一周一直在下雨，那么我们就可以合理地排除“降雨”这一变量作为原因。

该方法依赖统计检验来判断：

在给定其他变量条件下，两变量是否相互独立。

PC 算法的基本流程为：

1. 初始时，将所有变量连接成一个无向图；
2. 通过条件独立性检验，逐步删除变量之间不必要的连接；
3. 当所有条件独立关系确定后；
4. 根据一组方向规则，将边定向，最终形成有向无环图（DAG）。

PC 算法的特点是：

• 不依赖评分函数；
• 重点关注数据中的条件独立关系。

基于评分的方法（Score-based）

这种方法主要包含两个核心组成部分：

• 搜索策略（Search Strategy）
• 评分函数（Scoring Function）

在设定搜索策略与评分函数之后，算法会在所有可能的 DAG 空间中进行搜索。

具体过程是：

• 构建多个候选 DAG；
• 使用评分函数对每一个 DAG 进行评价；
• 找到评分最高的网络结构。

目标是找到：

在给定数据条件下最可能的因果结构。

换句话说，该 DAG 表示变量之间最合理的依赖关系。

在 DAG 中：

• 每个节点代表一个变量；
• 每条有向边表示变量之间的依赖关系（父节点 → 子节点）。

混合结构学习（Hybrid Structure Learning）

混合方法结合了：

• 基于约束的方法
• 基于评分的方法

其基本思想是：

1. 首先使用基于约束的方法识别部分条件独立关系，从而缩小搜索空间；
2. 然后使用基于评分的方法，在缩小后的空间中寻找最佳网络结构。

这种方法试图同时利用两类方法的优势：

• 基于约束的方法能够快速减少可能的 DAG 数量；
• 基于评分的方法能够进一步优化网络结构，使其更符合数据。

因此，在包含数百个变量的大规模数据集中，这种方法通常更加：

• 高效，
• 且准确。

它同时利用：

• 条件独立信息，
• 以及评分函数评估机制。

在下一部分中，作者将通过一些示例数据集（toy examples）演示：

• 如何选择合适的因果发现方法，
• 以及如何在具体应用场景中进行实践。

Python 中的 Bnlearn 库

下面简要介绍本文所有分析中使用的 bnlearn 库。

bnlearn 为 Python 提供了一套用于：

• 创建
• 分析
• 可视化

贝叶斯网络（Bayesian Networks） 的工具。

该库支持：

• 离散数据
• 连续数据
• 混合数据

并且设计目标是：

在保持易用性的同时，提供完整的贝叶斯因果学习流程。

使用 bnlearn，你可以执行：

• 结构学习
• 参数估计
• 概率推断
• 生成合成数据

因此，它既适用于探索性分析，也适用于高级因果发现研究。

通过简单设置参数即可应用多种统计检验方法。

此外，还提供辅助函数用于：

• 数据转换
• 计算图的拓扑排序
• 比较不同网络结构
• 生成多种可视化图表

核心功能包括

主要流程（Key Pipelines）

Structure Learning（结构学习）
从数据中学习网络结构，或结合专家知识构建网络。

Parameter Learning（参数学习）
根据观测数据估计条件概率分布。

Inference（概率推断）
对网络进行查询，计算概率、模拟干预并测试因果效应。

Synthetic Data（合成数据生成）
从学习得到的贝叶斯网络生成新数据，用于模拟实验、基准测试或数据增强。

Discretize Data（数据离散化）
利用内置方法将连续变量转换为离散状态，以便建模与推断。

Model Evaluation（模型评估）
通过评分函数、统计检验与性能指标比较不同网络。

Visualization（可视化）
提供交互式图结构绘制、概率热图以及结构对比展示。

bnlearn 相比其他贝叶斯实现的优势

• 包含最常用的贝叶斯分析流程；
• 使用方式简单直观；
• 开源（MIT License）；
• 提供完善的文档与博客支持；
• GitHub 超过 600+ Star；
• 下载量超过 70 万次。

安装方式：

pip install bnlearn

使用基于评分方法进行无监督因果发现的直观实践示例

为了直观地理解因果发现（causal discovery）的概念，我将从 Sprinkler（洒水器）数据集 开始，并演示如何学习数据结构。我们首先会生成所有可能的解（不使用搜索策略），然后对每一个结构进行评分。

示例

假设你有一个花园，并希望弄清楚为什么草地会变湿。在过去 1000 天 中，你测量了四个变量：

• Rain（降雨）
• Cloudy（多云）
• Sprinkler system（洒水系统）
• Wet grass（草地是否湿润）

每个变量都有两个状态：Yes 或 No

如果我们能够弄清楚这四个变量之间的连接关系，就可以通过干预措施来避免草地变湿。

在接下来的章节中，我将通过实际操作示例演示：如何根据数据集确定因果网络。我们将从仅包含四个变量的简单 Sprinkler 数据集开始，随后也会分析包含更多变量以及不同数据类型（如分类变量与连续变量）的数据集。

数据在 1000 天内对四个变量进行了测量：
Sprinkler system、Rain、Cloudy 和 Wet grass。

# Load library
import bnlearn as bn

# Load sprinkler dataset
df = bn.import_example('sprinkler')

# Print to screen for illustration
print(df)

+----+----------+-------------+--------+-------------+
|    |   Cloudy |   Sprinkler |   Rain |   Wet_Grass |
+====+==========+=============+========+=============+
|  0 |        0 |           0 |      0 |           0 |
+----+----------+-------------+--------+-------------+
|  1 |        1 |           0 |      1 |           1 |
+----+----------+-------------+--------+-------------+
|  2 |        0 |           1 |      0 |           1 |
+----+----------+-------------+--------+-------------+
| .. |        1 |           1 |      1 |           1 |
+----+----------+-------------+--------+-------------+
|999 |        1 |           1 |      1 |           1 |
+----+----------+-------------+--------+-------------+

如何找到最优因果网络？

最优因果网络，换句话说，即 有向无环图（Directed Acyclic Graph，DAG），表示能够最好描述数据集中变量关系的结构。

在本演示中，我们使用仅包含 4 个变量 的 Sprinkler 数据集。一般来说，由于可能组合数量巨大，寻找因果 DAG 是一项具有挑战性的任务。

在这个例子中，仅四个变量就已经存在 543 个可能的 DAG。这意味着可以构建 543 个不同的网络结构，并且每一个都需要根据其对数据集的拟合程度进行评估。

虽然看起来 4 个节点似乎不应该产生如此多的图结构，但复杂性来源于节点之间有向边的排列方式（同时必须保证不存在环）。

每一对节点之间可能存在：

• 没有边
• 从一个节点指向另一个节点的有向边
• 相反方向的有向边（但不能形成循环）

此外，每个变量具有两个状态（例如 yes / no）。如果状态数量增加，复杂度还会进一步提升。

具有 4 个节点且每个节点包含 2 个状态时，唯一 DAG 的总数为 543。

DAG 的总数量可以通过如下公式计算，其中总数 (G) 是节点数量 (n) 的函数 (G(n))，并且随 (n) 呈超指数增长。

用于计算给定节点数量 (n) 时可能 DAG 数量的公式

如果我们计算不同节点数量（每个节点具有两个状态）下的组合总数，可以发现 DAG 数量会随着节点数呈指数级增长。例如：

• 4 个节点 → 543 个组合
• 5 个节点 → 29,281 个组合
• 更多节点 → 数量迅速爆炸

你可以在下面的代码示例中自行尝试。我创建了一个运行示例，你会发现这是一个计算开销非常大的任务。

随着节点数量增加，可能 DAG 数量呈指数增长。
请尝试使用下面代码重复该实验。

import itertools

# Function to check if a graph is acyclic.
def is_acyclic(graph, n):
    visited = [False] * n
    rec_stack = [False] * n
    def visit(v):
        if rec_stack[v]:
            return False  # Found a cycle
        if visited[v]:
            return True  # Already visited
        visited[v] = True
        rec_stack[v] = True
        for u in range(n):
            if graph[v][u] and not visit(u):
                return False
        rec_stack[v] = False
        return True
    return all(visit(v) for v in range(n) if not visited[v])

# Number of variables (nodes)
n = 4

# Generate all possible directed edges between the 4 variables
nodes = range(n)
edges = list(itertools.permutations(nodes, 2))
total_graphs = 0

# Check all subsets of edges
for subset in itertools.product([0, 1], repeat=len(edges)):
    graph = [[0] * n for _ in range(n)]
    for i, (u, v) in enumerate(edges):
        if subset[i]:
            graph[u][v] = 1
    if is_acyclic(graph, n):
        total_graphs += 1

print(f"Total number of possible DAGs for {n} variables: {total_graphs}")

为了确定最佳因果 DAG，我们需要根据数据集对每一个 DAG 进行评分。

一种最直接的方法是：

对所有可能的 DAG 进行穷举评分（exhaustive search），然后按照得分排序，并选择得分最高的 DAG。

这种方法并不少见，被称为 穷举搜索。

优点

该方法会评估所有 DAG，因此可以确定找到全局最优解。

缺点

当变量数量增加时，所需计算资源极其庞大。

因此，穷举搜索通常只推荐用于变量数量较少的数据集。

下面的代码示例中，我们将加载 Sprinkler 数据集，使用穷举搜索方法对每个 DAG 进行评分，然后绘制最佳 DAG。

pip install bnlearn

import matplotlib.pyplot as plt
import bnlearn as bn

df = bn.import_example('sprinkler')

model = bn.structure_learning.fit(
    df,
    methodtype='exhaustivesearch',
    scoretype='bic',
    return_all_dags=True
)

print(model['adjmat'])

target     Cloudy   Rain  Sprinkler  Wet_Grass
source
Cloudy      False   True       True      False
Rain        False  False      False       True
Sprinkler   False  False      False       True
Wet_Grass   False  False      False      False

G = bn.plot(model, interactive=False)
dot = bn.plot_graphviz(model)

Sprinkler 系统的最佳 DAG 示例。

最佳评分 DAG 的逻辑如下：

• 草地变湿的概率取决于 Sprinkler 和 Rain
• 洒水系统是否开启取决于 Cloudy
• 是否下雨取决于 Cloudy

在上图中，我们绘制了每个 DAG 的 BIC 分数。

• 最佳 DAG 得分约为 -1953
• 最差 DAG 得分约为 -2700

这些分数仅在当前实验内部具有相对意义，不能与其他实验直接比较。

一个重要概念是 Markov 等价性（Markov equivalence）。

这意味着：

多个不同的网络结构，可能在仅使用观测数据的情况下，同样能够很好地解释变量之间的关系。

因此，如果没有额外信息（例如干预数据或领域知识），识别真实的底层因果结构会更加困难。

在下一节中，我将演示如何处理包含更多变量的数据集。请继续阅读！

在包含大量变量的数据集中寻找最优因果网络

当变量数量较少时，通过穷举测试所有可能组合来寻找 DAG 是一种很好的方法。然而，在真实世界场景中，我们通常会遇到 超过 10 个变量，并且每个节点可能具有多个状态。这意味着由于计算负担过大，我们已经无法再对所有可能的 DAG 进行穷举测试。

因此，我们需要一种搜索策略（search strategy），能够在不测试所有 DAG 的情况下，高效地遍历整个 DAG 搜索空间，从而找到最优 DAG。目前已经存在多种搜索策略，每一种方法都有其自身的特性。

搜索策略可以在无需逐一测试每个 DAG 的情况下，高效地遍历整个 DAG 搜索空间并找到最优结构。

让我们在一个中等规模到大型的数据集上进行实验。我们将加载 Alarm 监控系统数据集 [7]，该数据集包含 37 个离散变量，并且每个变量可能具有 2 个或更多状态。因此，所有可能唯一 DAG 所形成的搜索空间是极其巨大的。

此时已经无法再使用穷举搜索，因此我们需要一种更加智能的搜索策略。在这种情况下，我们需要一种：

• 运行速度快
• 适用于中小规模数据集

的搜索方法。

一个典型例子是 HillClimbSearch（爬山搜索）。这是一种启发式搜索方法（heuristic search），采用贪心局部搜索策略：

• 从一个完全无连接的 DAG 开始；
• 通过不断进行单条边的修改；
• 每一步选择能够最大化评分提升的操作；
• 当达到局部最优（local maximum）时停止搜索。

# Load libraries
import matplotlib.pyplot as plt
import bnlearn as bn

# Load Alarm dataset
df = bn.import_example('alarm')

# Learn the DAG using hillclimbsearch and BIC
model = bn.structure_learning.fit(
    df,
    methodtype='hillclimbsearch',
    scoretype='bic'
)

# Compute edge weights using ChiSquare independence test.
model = bn.independence_test(model, df, test='chi_square', prune=False)

# Plot the best DAG
bn.plot(
    model,
    edge_labels='pvalue',
    params_static={
        'maxscale': 4,
        'figsize': (15, 15),
        'font_size': 14,
        'arrowsize': 10
    }
)

# Plot using graphviz
dot = bn.plot_graphviz(model)

# Store to pdf
dot.view(filename='bnlearn_alarm')

基于 HillClimbSearch 搜索策略 与 BIC 评分函数 得到的 Alarm 数据集最佳 DAG。
边上的数值表示通过 卡方检验（Chi-square test） 得到的 。
图片由作者使用 Bnlearn 库生成。

到目前为止，我们已经看到并实际体验到：

当变量数量较多时，搜索策略并不是可选项，而是必须使用的工具。

因为它能够高效地遍历所有可能 DAG 的巨大搜索空间。

此外，每一个 DAG 仍然需要通过**评分函数（scoring function）**来评估贝叶斯网络对数据集的拟合程度。

整个流程被称为：基于评分的策略（Score-based Strategy）

在下一节中，我们将比较不同搜索策略与评分函数之间的差异。

当使用不同搜索策略与评分函数时，最终结果会发生变化

在使用基于评分的方法进行因果发现时：

搜索策略 + 评分函数 的组合共同决定最终得到的 DAG。

这本质上是在 速度（speed） 与 准确性（accuracy） 之间取得平衡。

因此，需要确保所选择的方法能够匹配：

• 具体应用场景
• 数据集结构

搜索策略

搜索策略是一种算法，用于高效遍历整个 DAG 搜索空间，从而在不测试所有结构的情况下找到最优 DAG。

不同方法在以下方面存在权衡：

• 计算复杂度
• 搜索效率
• 对变量依赖关系的刻画能力

方法选择取决于：

• 数据规模
• 数据结构
• 具体任务目标

穷举搜索

顾名思义，该方法会评估所有可能的 DAG，并选择评分最高的结构。

优点：

• 能保证找到最优结构（全局最优）。

缺点：

• 计算成本极高；
• 随 DAG 数量指数增长迅速失效。

通常仅适用于：

• 最大约 5–10 个节点 的数据集（前提是计算资源足够）。

爬山搜索

一种启发式搜索策略（Bnlearn 默认方法），适用于较大的网络。

工作方式：

• 从无连接 DAG 开始；
• 逐步进行单边修改；
• 每一步选择评分提升最大的操作；
• 在达到局部最优时停止。

缺点：

• 可能陷入局部最优；
• 多次运行结果可能不同。

Chow-Liu 算法

该算法专门用于学习树结构网络。

特点：

• 构建最大似然树；
• 每个节点最多只有一个父节点；
• 计算效率非常高。

此外：

• 根节点参数可选；
• 可控制 DAG 的起始节点。

适用于较简单的依赖模型。

TAN (树增强朴素贝叶斯）

TAN 是树结构方法的扩展版本，可表示更复杂关系。

特点：

• 在 Naive Bayes 基础上增加树结构依赖；
• 能捕获更多变量间关系；
• 在大规模数据集中表现良好。

要求：

• 必须指定根节点参数。

但在某些应用中，这种限制可能并不理想。

朴素贝叶斯

最简单的概率模型形式。

核心假设：

在给定类别标签条件下，各变量之间条件独立。

虽然这一假设在现实数据中很少完全成立，但在以下情况下表现良好：

• 大规模分类任务
• 特征依赖较弱
• 需要快速模型

缺点：

• 无法表示变量之间的相互依赖关系；
• 灵活性较低。

同样需要指定根节点。

Direct-LiNGAM

该方法专门用于：

• 连续数据
• 混合数据

采用半参数建模方法，并假设：

• 变量之间存在线性关系；
• 误差项服从非高斯分布；
• 图结构必须保持无环。

ICA-LiNGAM

用于从数据集中估计因果图结构的方法。

基本假设：

• 数据满足线性模型；
• 噪声为非高斯分布；
• 噪声在变量之间相互独立；
• 不存在未观测混杂因素（unobserved confounders）。

汇总

评分函数（Scoring Function）

评分函数用于量化某一个特定 DAG 对观测数据的解释程度。下面列出的评分函数都是常见且被广泛使用的方法。理解它们各自的特性，有助于在实际应用中选择最合适的搜索策略与评分函数组合。

Bayesian Information Criterion（BIC，贝叶斯信息准则）

BIC 在模型拟合优度与模型复杂度之间进行权衡。
它会对包含更多边的 DAG 施加惩罚，从而避免模型出现过拟合。

Bayesian Dirichlet Equivalent Uniform（BDeu）评分

该评分基于贝叶斯统计方法，使用 Dirichlet 分布 作为概率表的先验分布。

其中：

• 使用多项分布（multinomial distribution）来建模离散类别结果的概率；
• 在观察数据之后，后验分布仍然保持 Dirichlet 分布形式。

这一性质使得在包含离散变量的贝叶斯网络中：

• 信念更新
• 后验概率计算

过程更加简化。

K2 评分（K2 Score）

K2 评分假设：

• 所有可能的网络结构具有均匀先验分布（uniform prior）。

该方法通过变量之间的父子关系来估计观测数据的似然，从而评估网络结构的拟合程度。

与使用 Dirichlet 等特定先验的方法不同，K2 方法通过假设父节点之间相互独立来简化计算，使概率计算更加可行。

这种方法能够在不依赖复杂先验分布的情况下，高效地寻找较好的网络结构。

BDS 评分（Bayesian Dirichlet with Structure）

BDS 评分是在 BDeu 的基础上进行扩展的方法，它引入了网络结构的先验知识。

其特点包括：

• 同样假设 Dirichlet 参数先验；
• 同时对可能的 DAG 结构 引入先验分布。

当你拥有：

• 领域知识，
• 或对变量关系存在结构约束时，

该方法尤其有用。

BDS 在以下因素之间取得平衡：

• 数据拟合程度
• 网络复杂度
• 结构先验强度

通过同时引入参数先验与结构先验，BDS 往往能够得到更加准确且更符合实际问题的模型，尤其适用于小到中等规模数据集。

AIC 评分（Akaike Information Criterion）

AIC 是一种频率学派方法（frequentist method），用于在模型拟合效果与模型复杂度之间进行权衡。

其计算基于：

• 模型的最大似然估计（Maximum Likelihood Estimation）
• 对模型自由参数数量的惩罚项

模型越复杂，惩罚越大。

与 BIC 不同的是：

• BIC 在大数据集情况下更倾向选择简单模型；
• AIC 对复杂度的惩罚较弱，因此更容易选择包含更多参数的模型。

AIC 通常在以下情况下更受欢迎：

• 数据集较大；
• 数据为连续变量；
• 模型简洁性不是首要目标。

评分函数对比表

连续数据与混合数据的因果发现

到目前为止，我们主要讨论的是适用于离散数据集的基于评分方法。然而，在真实世界中，数据往往包含连续变量。

当数据集中包含连续变量或混合类型变量时，大致存在三种学习 DAG 的方法，每种方法都有其自身的优缺点：

• 利用领域知识，手动将连续数据离散化，然后再应用贝叶斯方法。
• 以自动化方式对连续数据进行离散化处理，随后使用贝叶斯方法建模。
• 直接使用贝叶斯方法对连续数据或混合数据进行建模。

通过手动离散化连续变量进行因果发现

基于领域知识对连续变量进行手动离散化是最直接的方法。然而，这通常也是最耗费精力的方法，因为它需要根据对数据背景的理解，将连续信号转换为一组离散区间。此外，还可能需要考虑变量之间的关系，这会进一步增加复杂性。

不过，如果离散化过程执行得当，就能够形成具有现实意义的数据分组，从而反映真实世界属性，并提高模型的可解释性。

为了演示该方法，我将加载 auto_mpg 数据集 [8]，该数据集来源于 UCI Machine Learning Repository。例如，我们可以根据现实意义定义汽车马力（horsepower）的类别：

• Low：马力小于 100（小型、省油汽车）
• Medium：马力在 100–150 之间（中等性能汽车）
• High：马力大于 150（高性能车辆）

该数据集包含 392 个样本 和 6 个变量，如下代码所示。

# Libraries
import pandas as pd
import datazets as ds
import bnlearn as bn

# Download dataset
df = pd.read_csv(
    'http://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data-original',
    delim_whitespace=True,
    header=None,
    names=['mpg','cylinders','displacement','horsepower','weight','acceleration','model year','origin','car name']
)

# Cleaning
df.dropna(inplace=True)
df.drop(['model year','origin','car name'], axis=1, inplace=True)

# Show dataset
df.head()

#     mpg  cylinders  displacement  horsepower  weight  acceleration
# 0  18.0        8.0         307.0       130.0  3504.0          12.0
# 1  15.0        8.0         350.0       165.0  3693.0          11.5
# 2  18.0        8.0         318.0       150.0  3436.0          11.0
# 3  16.0        8.0         304.0       150.0  3433.0          12.0
# 4  17.0        8.0         302.0       140.0  3449.0          10.5
#
# [392 rows x 6 columns]

根据领域知识定义马力区间：

# Define horsepower bins based on domain knowledge
bins = [0, 100, 150, df['horsepower'].max()]
labels = ['low', 'medium', 'high']

# Discretize horsepower
df['horsepower_category'] = pd.cut(
    df['horsepower'],
    bins=bins,
    labels=labels,
    include_lowest=True
)

print(df[['horsepower','horsepower_category']].head())

#    horsepower horsepower_category
# 0       130.0              medium
# 1       165.0                high
# 2       150.0              medium
# 3       150.0              medium
# 4       140.0              medium

有时也可以直接将连续值转换为整数类别，例如气缸数量：

# Set the cylinder to integers
df['cylinders'] = df['cylinders'].astype(int)

当所有连续变量都完成分类之后，就可以应用标准的结构学习流程。

# Structure learning
model = bn.structure_learning.fit(df, methodtype='hc')

# Compute edge strength
model = bn.independence_test(model, df)

# Plot DAG
bn.plot(model, edge_labels='pvalue')
dotgraph = bn.plot_graphviz(model, edge_labels='pvalue')
dotgraph

通过自动离散化连续变量进行因果发现

与手动离散化不同，我们也可以自动确定每个变量的最佳分箱方式。然而，与人工分箱相比，这类方法需要更加谨慎，因为自动生成的区间可能并不一定符合领域知识中的真实边界。

为了获得比简单的：

• 等宽分箱（equal-width）
• 等频分箱（equal-frequency）

更有意义的分类方式，我们可以先确定最符合数据的概率分布，然后利用 95% 置信区间 来划分：

• low
• medium
• high

三个类别。

此外，建议通过分布图进行可视化检查，以判断划分是否合理。

例如，对 acceleration（加速时间） 自动分箱后可以发现：

• 8–11 秒 → 快速车辆（fast）
• 20–24 秒 → 慢速车辆（slow）
• 中间区间 → 正常车辆（normal）

这一划分在现实中是合理的，因此可以继续使用这些类别。

如何寻找最适合数据的理论分布

# Install library
pip install distfit

from distfit import distfit
import matplotlib.pyplot as plt

# Initialize and set 95% confidence interval
dist = distfit(alpha=0.05)

# Fit distribution
dist.fit_transform(df['acceleration'])

# Plot distribution
dist.plot()
plt.show()

bins = [
    df['acceleration'].min(),
    dist.model['CII_min_alpha'],
    dist.model['CII_max_alpha'],
    df['acceleration'].max()
]

# Discretize acceleration
df['acceleration_category'] = pd.cut(
    df['acceleration'],
    bins=bins,
    labels=['fast','normal','slow'],
    include_lowest=True
)

del df['acceleration']

使用 distfit 进行概率密度拟合。

如果存在多个连续变量，也可以自动执行分布拟合：

cols = ['mpg','displacement','weight']

for col in cols:
    dist = distfit(alpha=0.05)
    dist.fit_transform(df[col])

    dist.plot()
    plt.show()

    bins = [
        df[col].min(),
        dist.model['CII_min_alpha'],
        dist.model['CII_max_alpha'],
        df[col].max()
    ]

    df[col + '_category'] = pd.cut(
        df[col],
        bins=bins,
        labels=['low','medium','high'],
        include_lowest=True
    )

    del df[col]

当所有连续变量完成分类后，即可再次应用因果发现中的结构学习方法。

# Structure learning
model = bn.structure_learning.fit(df, methodtype='hc')

# Compute edge strength
model = bn.independence_test(model, df)

# Plot DAG
bn.plot(model, edge_labels='pvalue')
dotgraph = bn.plot_graphviz(model, edge_labels='pvalue')
dotgraph

基于 HillClimbSearch 搜索策略 与 BIC 评分函数 得到的 Auto MPG 数据集最佳 DAG。
所有连续变量均已自动离散化。
使用 distfit 确定 95% 置信区间，通过分布图验证分类合理性。
图中边的数值为通过卡方检验计算得到的 (-\log_{10}(P\text{-value}))。
该图由 Bnlearn 生成。

通过直接建模连续变量进行因果发现（概念验证）

使用离散化方法具有重要意义，因为它能够开启大量可应用的贝叶斯统计方法。然而，在实际情况中，并不总是可以或希望对变量进行离散化。因此，也存在可以直接建模连续数据集的方法，例如 LiNGAM 方法，该方法同样已经在 Bnlearn 中实现。

Direct-LiNGAM 方法是一种半参数方法，它假设观测变量之间存在线性关系，同时要求误差项服从非高斯分布，并且保证所得图结构保持无环（acyclic）。

该方法通过反复执行：

• 回归分析
• 独立性检验

来完成结构学习，其中使用的是基于最小二乘法的线性回归。

在每一次回归中：

• 一个变量作为因变量（outcome）
• 另一个变量作为自变量（predictor）

该过程会对所有变量重复执行。

换句话说，Direct-LiNGAM 方法允许对连续数据和混合数据进行建模。

一个缺点是：
当使用该方法时，结构学习（causal discovery）本身就是最终结果，因此无法进一步执行：

• 参数学习（parameter learning）
• 概率推断（probabilistic inference）

为了演示其工作方式，我们将创建一个包含 6 个变量 的玩具数据集。该数据集的目标是展示不同变量的贡献及其对因变量的因果影响。

样本数量设为：n = 1000

连续变量服从均匀分布。

随后我们人为引入变量之间的依赖关系，并使用模型去恢复这些原始关系。

数据生成步骤

Step 1:
x3 为根节点，并初始化为均匀分布。

Step 2:
x0 与 x2 通过乘以 x3 的值得到，因此依赖于 x3。

Step 3:
x5 通过乘以 x0 得到，因此依赖于 x0。

Step 4:
x1 与 x4 通过乘以 x0 得到，因此依赖于 x0。

实验 DAG 及其边上的因果系数如下（使用 Bnlearn 绘制）。

import numpy as np
import pandas as pd
from lingam.utils import make_dot

# Number of samples
n = 1000

# step 1
x3 = np.random.uniform(size=n)

# step 2
x0 = 3.0 * x3 + np.random.uniform(size=n)
x2 = 6.0 * x3 + np.random.uniform(size=n)

# step 3
x5 = 4.0 * x0 + np.random.uniform(size=n)

# step 4
x1 = 3.0 * x0 + 2.0 * x2 + np.random.uniform(size=n)
x4 = 8.0 * x0 - 1.0 * x2 + np.random.uniform(size=n)

df = pd.DataFrame(
    np.array([x0, x1, x2, x3, x4, x5]).T,
    columns=['x0', 'x1', 'x2', 'x3', 'x4', 'x5']
)

print(df.head())

邻接矩阵 / 系数矩阵：

m = np.array([
    [0.0, 0.0, 0.0, 3.0, 0.0, 0.0],
    [3.0, 0.0, 2.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 6.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [8.0, 0.0, -1.0, 0.0, 0.0, 0.0],
    [4.0, 0.0, 0.0, 0.0, 0.0, 0.0],
])

dot = make_dot(m)
dot

结构学习可以通过 direct-lingam 方法执行。

# Load library
import bnlearn as bn

# Structure learning
model = bn.structure_learning.fit(df, methodtype='direct-lingam')

当我们查看输出结果时，可以发现模型很好地恢复了变量之间的依赖关系：

print(model['adjmat'])

# target        x0        x1       x2   x3        x4       x5
# source
# x0      0.000000  2.987320  0.00000  0.0  8.057757  3.99624
# x1      0.000000  0.000000  0.00000  0.0  0.000000  0.00000
# x2      0.000000  2.010043  0.00000  0.0 -0.915306  0.00000
# x3      2.971198  0.000000  5.98564  0.0 -0.704964  0.00000
# x4      0.000000  0.000000  0.00000  0.0  0.000000  0.00000
# x5      0.000000  0.000000  0.00000  0.0  0.000000  0.00000

bn.plot(model)

计算边强度（使用卡方检验统计量）：

model = bn.independence_test(model, df, prune=False)
print(model['adjmat'])

因果顺序（causal ordering）可以通过如下方式获得：

print(model['causal_order'])
# ['x3', 'x0', 'x5', 'x2', 'x1', 'x4']

我们也可以绘制因果图：

bn.plot(model)

最终结果如下 DAG 所示，可以看到：

• 因果网络结构被正确恢复
• 我们最初设定的因果系数同样被成功估计

该因果图由 Bnlearn 生成。

使用 Auto MPG 数据集进行连续数据的因果发现

作为最后一个示例，我们再次加载 auto mpg 数据集，但这一次不进行离散化，而是直接对连续变量建模。

在该方法中，我们将：

• 方法设置为 direct-lingam
• 使用 bnlearn 执行结构学习

import bnlearn as bn

# Load example mixed dataset
df = bn.import_example(data='auto_mpg')
del df['origin']

# Structure learning / causal discovery
model = bn.structure_learning.fit(
    df,
    methodtype='direct-lingam',
    params_lingam={'random_state': 2}
)

绘制结果：

bn.plot(model)

计算边强度：

model = bn.independence_test(model, df, prune=True)

Plotting

bn.plot(model)
bn.plot(model, edge_labels='pvalue')

# Interactive plot
bn.plot(model, interactive=True)

# Graphviz plot
dotgraph = bn.plot_graphviz(model)
dotgraph

# Export pdf
dotgraph.view(filename='auto_mpg_bnlearn.pdf')

# Graphviz with p-values
dotgraph2 = bn.plot_graphviz(model, edge_labels='pvalue')
dotgraph2

因果发现（结构学习）的最终结果如下图所示。

使用 Direct-LiNGAM 方法 对连续 Auto MPG 数据集 进行因果发现：

• 图 A：表示边上的因果权重
• 图 B：表示变量之间统计显著性（−Log10(P-value)）

图像由 Bnlearn 生成。

当你只有专家知识或少量样本时生成合成数据

在许多领域，例如 医疗、金融、网络安全以及自动驾驶系统 中，真实世界的数据往往具有以下特点：

• 数据敏感
• 获取成本高
• 数据分布不平衡
• 难以收集（尤其是罕见或极端场景）

此时，合成数据（Synthetic Data） 就成为一种非常有力的替代方案。

从整体上看，生成合成数据大致可以分为两大类：

• 概率式方法（Probabilistic）
• 生成式方法（Generative）

如果你希望进一步了解，那请期待后续的文章。

贝叶斯统计很强大，但是……

到目前为止，我对贝叶斯方法一直保持相当积极的态度。但同样重要的是理解它的一些局限性。

首先，并非所有复杂系统都可以表示为有向无环图（DAG）。
某些系统本身是循环系统（cyclical systems），而这类结构无法通过标准贝叶斯网络建模。

此外，因果模型往往包含随机性效应（stochastic effects）。这意味着：

即使使用相同的贝叶斯方法，多次实验也可能得到不同结果（例如落入不同的局部最优）。

马尔可夫等价性（Markov Equivalence）

另一个重要问题是 Markov 等价性。

这意味着：

多个不同的网络结构，仅依靠观测数据，都可能同样良好地解释变量之间的关系。

因此，如果没有额外信息，例如：

• 干预数据（interventional data）
• 领域知识（domain knowledge）

就很难确定真正的底层因果结构。

模型选择的重要性

还需要注意的一点是：

贝叶斯模型的结果会随着所使用模型的不同而变化。

在 Auto MPG 数据集 的示例中已经展示：

不同建模方法可能产生不同的 DAG。

例如：

**BIC（贝叶斯信息准则）**通常假设：

更简单的结构更优。

但在复杂系统中，这一假设并不总是成立。

因此，应始终根据以下因素选择方法：

• 搜索策略（search method）
• 评分函数（scoring function）
• 数据集特性
• 实际应用场景

库中的默认参数并不一定适用于你的问题。

数据不完整问题

现实中的数据集通常是不完整的，例如：

• 缺少变量
• 样本数量不足

这些问题同样会影响模型结果。

如果希望基于贝叶斯方法得出强结论，必须认识到：

这通常依赖于较强的建模假设。

换句话说，在复杂真实世界事件中构建应用可能会变得困难，因此往往需要：

独立验证实验（independent validation tests）

例如：

在疾病发病机制建模等高风险场景中，可以通过实验室实验来验证模型结果。

最后的话

到这里，你已经学习了如何在以下情况下应用因果发现方法：

• 小规模与大规模数据集
• 离散变量
• 连续变量
• 混合类型变量

我们实验了多种流行的贝叶斯因果建模库，包括：

监督式因果方法

用于包含处理变量与干预变量的因果发现：

• DoWhy
• PyMC

无监督因果发现方法

用于从观测数据中自动发现关系：

• Pgmpy（提供底层统计函数）
• Bnlearn（提供完整流程管道）

Bnlearn 使得在小型到大型数据集中检测因果关系变得更加直接。

希望这篇入门指南能够帮助你理解如何应用因果发现方法。

成功的关键在于：

理解哪种统计方法与工具最适合你的问题。

因此，请持续：

• 尝试不同方法
• 比较实验结果
• 深化对因果建模的理解

最后，我强烈建议你亲自运行本文中的实践示例，这将帮助你：

• 学得更快
• 理解更深
• 记忆更久

来杯咖啡，好好探索吧 ☕。

References

Bayesian Statistics, Wikipedia

Taskesen E. Chat with Your Dataset using Bayesian Inferences. Data Science Collective (DSC), Nov. 2025

Taskesen E. Human-Machine Collaboration with Bayesian Modeling: Learn To Combine Knowledge With Data, Data Science Collective (DSC), Okt. 2025

Taskesen E. The Starters Guide to Causal Structure Learning with Bayesian Methods in Python, Data Science Collective (DSC), Sep. 2025

Kay H. et al. Inferring causal impact using Bayesian structural time-series models, 2015, Annals of Applied Statistics (247–274, vol9)

I. A. Beinlich et al. The ALARM Monitoring System: A Case Study with Two Probabilistic Inference Techniques for Belief Networks, European Conference on Artificial Intelligence in Medicine, 1989.

UCI Machine Learning Repository — Auto MPG dataset (CC BY 4.0)

Taskesen, E. (2020). Learning Bayesian Networks with the bnlearn Python Package (Version 0.3.22)

Taskesen, E. Six Causal Libraries Compared: Which Bayesian Approach Finds Hidden Causes in Your Data?, Data Science Collective (DSC), Okt. 2025

Kay H. et al. Inferring causal impact using Bayesian structural time-series models, 2015, Annals of Applied Statistics (247–274, vol9)