简介：一套基于K2算法的贝叶斯网络结构学习实现，面向机器学习、生物信息等需要从观测数据中推断网络结构的研究者与开发者。资源聚焦K2评分搜索策略，通过MATLAB脚本与C源码配合，演示在节点顺序约束下贪婪搜索网络结构的过程，可帮助解决中等规模变量的因果图构建问题。压缩包共含7个文件，以4个m文件为主，另含1个C源文件、1个txt说明文档和1个mat样例数据，整体仅10KB。m文件覆盖主控函数、辅助函数与对象构建等模块，C文件提供底层计算逻辑，txt记录运行要点，mat数据便于直接复现示例。目前已有2167人学习下载，属于结构学习方向较为轻量的小型资源。读者可从中获得可运行的K2算法原型、带注释的MATLAB函数及配套样例数据，深入理解贝叶斯网络结构学习中评分与搜索的协作机制；同时C源码为二次开发或将算法嵌入到自己的实验框架提供了参考基础。
1. 为什么结构学习绕不开 K2 算法：一个直观场景
你手上有一批业务系统的离散观测数据，变量有十来个，你知道其中几个明显有先后关系，却不清楚完整因果结构。直接用贝叶斯网络结构学习去搜一张有向无环图，候选数量会随变量数爆炸式增长；而 K2 算法用一个很朴素的思路避开这个坎：先认定一个节点顺序，再逐节点挑选父节点，边挑边看评分涨不涨。这篇文章把 K2 算法从“概率图模型里的名词”拆成能立刻落地的流程：评分函数为什么有效、节点顺序怎么给、父节点上限怎么调，以及你会遇到的真实翻车现场。适合手里有离散数据、想快速得到一张可解释 DAG 的工程师或分析人员；如果你是第一次接触贝叶斯网络，也能从这份实现里看懂结构学习的大致套路。
2. 贝叶斯网络结构学习的评分思路：K2 算法在找什么
2.1 结构学习不是黑匣子：评分函数回答“哪个网络更好”
贝叶斯网络结构由一组有向边和条件独立假设组成。结构学习的目标是从数据中找一张 DAG，使得网络能较好地解释观测数据，又不至于复杂到没法用。完整搜索是不可行的：n 个节点的 DAG 数量超指数，所以常见做法是把问题拆成“评分函数 + 搜索策略”两个模块，K2 算法正是在这个框架下工作的。
评分函数给每个候选父节点集打一个分。对某个节点来说，父节点选得越好，它自身的条件分布就越能解释数据。所有节点的分数叠加，就是整张网络的结构分。搜索策略负责在巨大的父集组合中寻找分数高的结构，而不是把每一条可能边都试一遍。K2 的搜索策略很朴素：给定一个节点顺序，对每个节点只从排在前面的节点中挑选父节点；无环约束被节点顺序自动保证，这是它的第一个关键简化。
2.2 K2 的局部评分：父节点如何决定一个变量的条件分布
假设当前要学的是节点 Xi，候选父集合是 Pa。在离散数据上，父节点的每种取值组合给出一个条件概率表，表里是 Xi 各状态的分布。拟合数据越好，似然越高；但只追求似然会鼓励把越多父节点塞进来，所以评分函数还要加一个复杂度惩罚。
K2 原始评分背后的想法是：把参数先验当成伪计数，计算在给定父集条件下 Xi 的边际似然；BDeu 是带 Dirichlet 先验的扩展，BIC 则是它的渐近近似。用式子看更清楚：
BIC 局部分 = log P(D | Pa) - 0.5 * d * log N
其中 N 是样本量，d 是父节点组合下的自由参数数量。这个惩罚项会让收益不大的弱边被砍掉。如果你自己写实现，最省事的是直接算 BIC；如果样本量小或者某个父配置下样本极度稀疏，BDeu 的平滑先验更稳妥。我一般会在代码里把伪计数 alpha 作为参数留着，数据量特别大时 alpha 的影响会变小，数据量小时它能救命。
这里维度 d 的算法是：父节点组合数乘以 (当前节点状态数 - 1)。比如 Xi 有 3 个状态，一个父节点有 2 个状态，那么 d = 2 * (3 - 1) = 4。如果父节点组合里许多格子的样本数量为 0，普通的最大似然会给零概率，评分直接变成负无穷；所以平滑不是可选项，是必选项。alpha 设 1 时等价于拉普拉斯平滑，设更小会让模型对稀疏组合更敏感，设更大则会把弱依赖也抹平。实际调参时要结合变量状态数和样本量一起看。
2.3 贪心搜索策略：按节点顺序逐个“养父节点”
K2 的输入除了数据，还有一个节点顺序。顺序的含义是：排在后面的节点可以依赖排在前面的节点，反过来不行。这既是限制，也是优势。因为边只能从前面指向后面，学习过程天然不会成环，搜索时不用每加一条边就做一次环路检测。
算法对每个节点 Xi 这样走：先假设它没有父节点，算一个初始分数；然后在所有候选父节点里逐个尝试加入一个父节点，看局部评分提升多少；挑提升最大的那个加入。重复这个过程，直到没有父节点能继续提升评分，或者父节点数量达到你设定的上限 max_parents 为止。
我用大白话把执行步骤拆成五步：
- 把所有变量按顺序排好，前面的变量才能当后面的父节点。
- 从第二个节点开始处理，第一个节点没有候选父节点，父集固定为空。
- 对当前节点，把前面所有节点当作候选父节点，逐个算一遍加入前后的局部分差。
- 选择分差最大的候选加入；如果最大分差小于等于 0，停止这个节点。
- 如果达到 max_parents，也停止，换下一个节点。
这里的贪心是关键：每次只看“下一个父节点”能不能提升，选了 A 就不会再考虑用 B 换掉 A。因此 K2 不能保证全局最优，它给的是局部最优解。实际工程中这通常够用，尤其是变量顺序有一部分来自业务时间线时，结果的可解释性比其他黑匣子优化算法好得多。运行复杂度也可以估算：每个节点最多尝试 max_parents 轮，每轮遍历候选父节点集合，最坏情况是 O(n^2 * max_parents * 评分计算)，几十个变量内完全跑得动。
举个例子，如果真实因果是 A 和 B 共同影响 C，而顺序写成了 A、C、B，那么 C 无法把 B 当父节点，K2 只能在 C 后面再加一条 C -> B 或 A -> C 之类的边来补偿。这个补偿看似提升了分数，但结构已经偏离真实机制。所以每跑完一组顺序，我都会做一件事：把结果里所有边的方向和业务时间线逐条对一遍，对不上的边单独标出来。这不是可选项，而是 K2 使用的基本纪律。
2.4 和 PC、MMHC 相比，K2 什么时候值得用
PC 算法走的是条件独立检验路线：先检验无条件独立的变量对，再逐步引入条件集，最后定向边。它对样本量、显著水平非常敏感，变量一多，条件集组合也会爆炸。MMHC 是约束 + 评分混合方法，先用条件独立检验筛掉很多边，再做贝叶斯评分搜索，实现复杂度明显更高。K2 是三者里最容易实现、最容易解释的。
| 方法 | 是否要节点顺序 | 是否检查无环 | 主要瓶颈 | 适用场景 |
|---|---|---|---|---|
| K2 | 是 | 不需要，顺序保证 | 节点顺序准确性 | 变量几十个以内、有先验顺序的离散数据 |
| PC | 否 | 定向时处理 | 条件独立检验对样本量敏感 | 变量少、样本足、希望弱假设 |
| MMHC | 否 | 是 | 候选筛选和搜索复杂 | 变量规模较大、追求稳定性 |
什么情况不要优先选 K2？如果你的业务变量里有一个明显的隐藏共因，它不在数据列里，又总是排在结果前面或后面，K2 学出来的边就会是错的，因为算法只能从观测变量里选父节点。这时候你应该先补变量，或者改用能容纳隐变量的方法。如果变量之间主要是强相关但没有方向先验，K2 也容易把相关错当成因果，此时要用多种顺序做敏感性验证，后文我会专门谈这个坑。因此，先把 K2 当作“快速得到一张可作为 baseline 的 DAG”，再考虑更复杂的结构学习算法。
判断 K2 是否适用，可以问三个问题：变量是否都离散或可分箱？是否至少能给出一个粗略顺序？样本量是否足以支撑条件概率表？如果你三个答案都是否，K2 大概率会让你很失望。工程上我不会在没有顺序的项目里硬上 K2，宁可先用相关性网络做探索，再回头补领域顺序。
3. 从零开始跑通 K2：最小代码与逐步拆解
3.1 先有一份离散化数据：状态计数决定一切
K2 要求数据是离散的，连续变量必须先分箱。为了让你能直接复现，我构造一份 1000 条样本的小数据：rain 是下雨状态，sprinkler 是洒水器开关，mud 是地面泥泞。真实关系定为 rain 影响 sprinkler，rain 和 sprinkler 共同影响 mud。这样后面能清楚看出 K2 是否找回依赖方向。
import numpy as np import pandas as pd rng = np.random.default_rng(42) rain = rng.binomial(1, 0.4, 1000) sprinkler = np.where(rain == 1, rng.binomial(1, 0.05, 1000), rng.binomial(1, 0.35, 1000)) mud = ((rain == 1) | (sprinkler == 1)).astype(int) data = pd.DataFrame({'rain': rain, 'sprinkler': sprinkler, 'mud': mud})这里rng.binomial(1, p, 1000)生成 0/1 状态的二项分布样本；np.where让 sprinkler 的条件概率依赖 rain。mud 的生成用了“或”逻辑，让 rain 和 sprinkler 都成为父节点。变量状态数越小，K2 的条件概率表越密，学出来的边越稳。真实业务里如果变量是连续值，一般用分位数或业务阈值先离散化成 3 到 5 档，不要直接喂原始数值给 K2。
3.2 手写 K2：一个小而完整的可运行实现
下面这份代码是 K2 最核心的部分。它计算每个节点的局部 BIC 分数，并按顺序贪心添加父节点。我没有用现成库，是为了让你看清每一步在算什么。
import math from collections import Counter def local_score(data, node, parents, alpha=1.0): states = sorted(data[node].dropna().unique()) r_i = len(states) n_total = len(data) log_lik = 0.0 if not parents: counts = data[node].value_counts(dropna=False) for s in states: n_ijk = counts.get(s, 0) p = (n_ijk + alpha) / (n_total + alpha * r_i) log_lik += n_ijk * math.log(p) dim = r_i - 1 else: for _, group in data.groupby(parents, observed=True): n_ij = len(group) counts = group[node].value_counts(dropna=False) for s in states: n_ijk = counts.get(s, 0) p = (n_ijk + alpha) / (n_ij + alpha * r_i) log_lik += n_ijk * math.log(p) dim = 1 for pa in parents: dim *= data[pa].nunique(dropna=False) dim *= (r_i - 1) bic = log_lik - 0.5 * dim * math.log(n_total) return bic def k2_learn(data, node_order, max_parents=3, alpha=1.0): parents = {node: [] for node in node_order} for i, node in enumerate(node_order): candidates = [n for n in node_order[:i]] current_score = local_score(data, node, [], alpha) while candidates and len(parents[node]) < max_parents: best_gain = 0.0 best_parent = None for pa in candidates: if pa in parents[node]: continue score = local_score(data, node, parents[node] + [pa], alpha) gain = score - current_score if gain > best_gain: best_gain = gain best_parent = pa if best_parent is None: break parents[node].append(best_parent) candidates.remove(best_parent) current_score += best_gain return parentslocal_score的核心是统计父节点组合下的条件计数，再用平滑后的概率算对数似然，最后减去dim * 0.5 * log(N)的复杂度惩罚。alpha是伪计数，默认 1，对应拉普拉斯平滑；如果你的数据里某个父配置的样本非常少，可以把 alpha 调大一点，但要清楚这会弱化真实差异。
k2_learn里有个细节：候选父节点只取node_order[:i]。这让无环约束被顺序天然保证，不需要像其他搜索算法那样反复做拓扑检查。max_parents控制树的宽度，默认 3；超过 3 的父节点组合会让条件概率表碎掉，除非样本量非常大，否则不建议一上来就调高。
这份代码故意不用成熟库，原因有二：一是 K2 在主流 Python 库里往往只是评分函数，结构搜索还要配合 HillClimbSearch 等通用贪心算法，容易让你误以为已经用了真正的 K2；二是真实业务里你需要控制顺序、max_parents、平滑参数，手写这个小循环反而更容易理解边界。如果你想快速验证，可以在 pgmpy 里用 K2Score 搭配 HillClimbSearch 跑一遍，但它不能强制使用你给的节点顺序，只能当作参考。
3.3 跑一次并对比错误节点顺序：顺序是 K2 的磨心
现在用两套顺序跑同一个数据。第一套是真实生成顺序：rain 在前，sprinkler 在中间，mud 最后。第二套故意把 mud 放到最前面，模拟顺序完全反过来。
order_true = ['rain', 'sprinkler', 'mud'] order_reverse = ['mud', 'rain', 'sprinkler'] print(k2_learn(data, order_true)) print(k2_learn(data, order_reverse))输出结果会明显不同。正确顺序下，K2 大概率找到rain -> sprinkler和rain -> mud、sprinkler -> mud；错误顺序下，mud 被当成最上游，它只能学出从 mud 出发的边，方向完全逆转。这个实验说明 K2 不是“从数据里凭空发现因果”，它是在你给定的顺序约束下寻找最合理的结构。
如果你手上没有可靠顺序，我建议多跑几组候选顺序，把结果里的边放进一个集合，看看哪些边在不同顺序下都能出现。只有那些稳定的边才值得写进报告。顺序对 K2 的影响排在所有参数的第一位，比 max_parents 和评分函数更值得先花时间。
4. 调参是结构学习的一半：节点顺序、父节点上限、评分函数
4.1 节点顺序怎么定：从领域因果链到数据校验
节点顺序在 K2 里不是可有可无的输入，它决定了边方向的搜索空间。最常见的定顺序方式是业务时间线：上游流程的变量放前面，下游结果放后面。比如在客服工单数据里，渠道来源、客户分级、订单金额都发生在“是否投诉”之前，那它们就都应该排在投诉前面。医学里更明显：暴露因素在中间变量之前，中间变量在结局之前。
如果业务上拿不出完整顺序，也有补救办法。你可以先用逻辑回归或互信息做一个初步关联排序，把和目标变量关联弱的排后面；也可以把多个顺序跑出来的结构做对比，挑数据评分更高的一组。注意，评分高不一定代表因果正确，但至少说明在同样的评分函数下它更匹配数据。顺序是 K2 的“半监督先验”，你对业务秩序了解越深，K2 的结果越可信。反过来说，完全随机排序时，K2 学出来的方向基本可以视为随机，别把它当成因果证据。
4.2 max_parents：宁缺毋滥，网络稀疏化的第一道闸门
max_parents 是每个节点最多能连上的父节点数。K2 是贪心添加父节点，如果没有上限，评分函数会一直尝试加边，最后形成一张又密又难解释的网。限制到 2 到 4，网络会稀疏很多，也更容易在一个节点上做人工校验。
| 参数取值 | 节点条件表大小 | 常见效果 |
|---|---|---|
| max_parents=2 | 最多 2 个父节点组合 | 结构稀疏，边容易解释，适合 500 样本以内的数据 |
| max_parents=3 | 3 个父节点组合 | 能表达常见交互，适合 1000 样本以上 |
| max_parents=5 | 5 个父节点组合 | 只在样本量很大时才敢用，否则参数估计方差很大 |
我实际调参时不会只看分数。我会从 max_parents=2 开始跑，逐步升到 3、4，观察新增的边是不是业务上能解释的。如果新增边让某个节点的条件概率表出现大量“父组合无样本”，那就降回来。另一种验证是留一法或交叉验证：把数据切成训练集和验证集，分别比较学到的结构在验证集上的平均对数似然，分数下降时就说明 max_parents 加过头了。
4.3 BIC 与 BDeu 怎么选：先验强度和样本量决定评分走向
K2 原始评分和 BDeu 都源自贝叶斯边际似然，区别在于先验假设。K2 原始形式假设参数先验与样本无关，BDeu 用等价样本量参数 es 控制先验强度。BIC 则是一个更简单的渐近近似，它的惩罚项只和样本量、自由参数数有关，和先验无关。
如果你在工程里没有特殊理由，我会建议默认 BIC 或 BDeu 都行，但先验参数要小心。BDeu 的等价样本量 es 越大，说明越相信参数接近均匀分布，学出来的结构倾向更稀疏；es 越小，越相信数据本身，结构倾向更复杂。es 从 1 改到 10，结构变化是正常的。样本量超过几千时，先验影响变小，BIC 和 BDeu 的结果会趋于一致；样本量只有几百时，BDeu 的平滑能让评分稳定很多。
| 评分 | 是否需要先验参数 | 主要特点 | 推荐场景 |
|---|---|---|---|
| BIC | 否 | 简洁、易实现，惩罚固定 | 样本量较大的 baseline |
| BDeu | es | 平滑效果好，对稀疏表友好 | 小样本、状态数多 |
| K2 原始 | 伪计数 | 计算直观 | 教学和二次开发 |
建议把评分函数当成可替换模块，不要写死在代码里。上面那份手写代码里，我用的 BIC 简洁但吃样本，如果切换到 BDeu，只需要把伪计数 alpha 和维度惩罚的公式替换成 BDeu 的边际似然，其余结构不变。
4.4 离散化：分箱数对结构学习的直接影响
K2 处理不了连续变量。常见的做法是等宽分箱或分位数分箱，但分箱数的选择不是随意的。分箱太少，比如两箱，会把非线性关系压成简单的“高/低”，遗漏条件依赖；分箱太多，比如十箱，父节点组合数爆炸，许多组合只剩几条样本，评分容易把弱相关也当成稳定边。
我通常先用分位数分成三档：低、中、高。跑通 K2 后，再看关键节点在父节点不同取值下的条件概率是否单调或有序。如果发现某个变量在分箱后几乎没有区分度，就把分箱数降到两档再跑一遍；如果核心边的方向没有变化，说明离散化对结论影响不大。你在论文或技术报告里要写明分箱方法，否则别人复现不了。
这里还有一个容易忽略的点：离散化要在结构学习之前完成，且不能把未来时间的信息带进去，尤其是做动态贝叶斯网络时，每个时间片要按同一套分箱边界处理。后续第 6 章我会讲到怎么往动态方向试探，离散化的一致性就是那一步的地基。
5. 避坑清单：K2 结构学习最常见的 5 个踩坑点
5.1 现象：换一个节点顺序，学出完全不同的 DAG
这是新手遇到最多的翻车现场。同一个数据，顺序 A 学到X -> Y，顺序 B 学到Y -> X，两套网络在验证集上的评分还差不多。
原因在于 K2 本身不识别因果方向，它只保证“边方向符合给定顺序”。如果 X 和 Y 只是强相关，设置顺序时谁在前面，K2 就倾向让谁当父节点。解决问题要用业务顺序做约束，并用多组顺序做敏感性测试。对每一个可能方向，记录它出现在多少个顺序的结果里；只有在绝大多数顺序下方向都一致时，才能把它当作可用结论。
5.2 现象：结果全是单方向链，缺少业务里的共因变量
有些数据源自同一个上游因子，比如“季节”同时影响“温度”“销量”“客流”。如果季节变量不在数据集里，K2 找不到共同父节点，它会用一条链去近似，比如让温度当销量父节点，客流当温度父节点。业务上这些边很容易被误读。
原因是 K2 只能从已给变量里选父节点，隐藏共因不会凭空出现。解决方法是先检查变量清单，把已知的季节、地区、渠道等共因变量放进顺序靠前的位置；如果确实采集不到，就不要把 K2 输出当成因果结构，最多当相关性网络。另一个思路是换成同时考虑隐变量的结构学习算法，但工程复杂度会高很多。
5.3 现象：评分一直上升，最后每个节点挂满父节点
如果 max_parents 设得很高，比如 6，BIC 的惩罚项也没拦住，你可能会得到一张完全连通的网络。原因通常是样本量足够大，模型可以在条件概率表里塞下任何弱相关，惩罚项相对变弱；也可能是离散化太细，每个父组合下样本数很少，似然对稀疏组合异常敏感。
先降 max_parents，一般不超过 4。然后做交叉验证：分别比较 max_parents=2、3、4、5 时结构在验证集上的对数似然，如果验证集分数峰值出现在某个较小值附近，说明大父集只是在拟合训练集噪声。若确认离散化过细，把分箱数降到 3 到 5 档，再观察边的数量变化。
5.4 现象：日志里出现 NaN 或报错 math domain error
报错通常出现在局部计数为 0 时。比如某个父节点组合下，数据里完全没有当前节点的某个状态，未经平滑的log(0)直接让程序崩溃，或者评分变成负无穷后所有候选都“无法改进”。
原因很简单：代码没有做拉普拉斯平滑，或者平滑用的 alpha 相对状态数太小。解决方法是把alpha=1.0写成默认参数，并确认states来自全量唯一值而不是当前分组里的唯一值；后者会漏掉零计数导致维度计算错误。这个坑我在自己写的脚本里踩过多次，所以上面的代码坚持用全量states列表做循环。
5.5 现象：用 K2 学完，放进动态贝叶斯网络做推理，概率严重偏移
这里涉及到动态贝叶斯网络。静态 K2 学的是同一个时间截面里的依赖，它不知道变量在 t 和 t+1 时间片之间的关系。你直接把静态边平移到动态网络里，等于假设“影响是瞬时的”，推理结果当然会偏离真实时序。
正确做法是先把同一时刻的变量放进一个时间片，用 K2 学片内边；再对跨时间片候选边做同样的评分扩展，比如rain_t -> sprinkler_{t+1}。要保证训练样本的时间对齐，不能把不同时刻的数据混成一行。动态贝叶斯网络下的 K2 变体仍然用局部评分，但候选父节点既包括同片前序节点，也包括上一片的节点，搜索逻辑和静态 K2 一致。
6. 用 bootstrap 给每条边算置信度，顺便接住动态贝叶斯网络
静态 K2 给的是唯一一张 DAG，但工程上更关心这张图稳不稳。我现在的习惯是学完结构后必做 bootstrap：对原始数据有放回抽样，跑 20 到 50 轮 K2，统计每条边出现的频率。频率超过 0.7 的边才算稳定，低于 0.3 的直接丢进待排查列表。
from collections import Counter edge_count = Counter() for _ in range(20): idx = rng.choice(len(data), len(data), replace=True) dag = k2_learn(data.iloc[idx], ['rain', 'sprinkler', 'mud'], max_parents=3) for node, pas in dag.items(): for pa in pas: edge_count[(pa, node)] += 1 print(edge_count)这段代码不需要跑太久，20 轮足够暴露哪些边是“换一批样本就消失”的噪声。如果你要更严谨，可以用 100 轮，但小数据上 20 轮已经能看出基本模式。记住 bootstrap 频率不等同于因果强度，它只说明网络结构对样本扰动的稳定性。
如果数据带时间戳，下一步再往动态贝叶斯网络走。做法是建时间片副本，把X_t和X_{t+1}都放进变量池，用同样的 K2 搜索逻辑学跨时间片边；同片内结构可以直接复用静态 K2 结果，跨时间片边要单独设 max_parents。我吃过不做事后校验的亏，曾经把一个只在 bootstrap 里出现 0.2 频率的边写进报告，后来被业务方当场质疑。从那以后，我的默认流程永远是“先跑 K2，再跑 bootstrap，最后才谈因果方向”。希望帮到你。
本文还有配套的精品资源，点击获取
