K2算法详解：贝叶斯网络结构学习与节点顺序优化

简介：一套基于K2算法从数据中学习贝叶斯网络结构的MATLAB/C实现资源，面向机器学习、生物信息及概率图模型方向的学生和工程师。它解决在给定节点顺序下，利用贪心搜索构建有向无环图（DAG）并计算K2评分的问题，完整演示了K2算法以节点顺序为先验、逐步为每个节点挑选最优父节点的结构学习流程。资源共7个文件，以4个MATLAB脚本为主，含主程序、子函数和控制流程模块；另有C源码用于编译MEX加速，mat文件提供可直接运行的样例数据，txt为许可证说明，压缩包仅10KB，轻量易部署。已有2167人学习下载，适合希望理解K2算法细节并在自身数据上开展结构学习的读者。通过这套代码可掌握LG对象构造、闭合函数调用、C-MEX接口封装等关键环节，参照样例数据即可运行出DAG与K2Score，为进一步扩展参数学习或因果推断奠定基础。

1. K2算法是什么：一个靠“节点顺序”吃饭的结构学习算法

如果你接过贝叶斯网络的项目，大概率会卡在第一步：结构从哪来？人工画图主观又费劲，纯数据驱动又怕学到一堆解释不通的边。K2算法正是贝叶斯网络结构学习里最经典的一条落地路径：它先把所有变量排成一个先后顺序，再沿着这个顺序用贪心搜索为每个变量挑选父节点，把组合爆炸的图搜索压缩成一个个局部“挑父节点”的小决策。真正上手后你会发现，K2算法实现并不难，但决定结果上限的往往不是打分函数，而是那个看起来不起眼的节点顺序——顺序给错，学出来的网络方向可能全反。这篇文章把K2的原理、代码、参数调整和常见翻车场景一起讲透，适合正在做结构学习、又不想把模型当黑匣子的从业者。

2. 为什么是K2：打分搜索路线与BDeu的数学直觉

2.1 结构学习的两条路线，K2站在哪一边

贝叶斯网络结构学习大体分两条路线。第一条是约束测试路线，代表是PC算法：通过一轮又一轮条件独立性检验来判断两个变量之间是否有边、方向如何。这条路在小样本场景下经常翻车，因为每一次独立性检验都在消耗样本，几十个变量十几条边测下来，p值早就不稳了。第二条是打分搜索路线，代表就是K2：定义一个能同时反映“拟合得好不好”和“结构复杂不复杂”的得分函数，然后到所有可能的图里去搜索得分最高的那个。

K2站在打分搜索这边，并且做了一个关键简化：它不搜索所有图，而是要求使用者先给定一个节点顺序。父节点只能从顺序靠前的节点里选出来，这等于把搜索空间从2^(n(n-1)/2)直接压到每个节点的父集组合数量级，工程上立刻可解。这也是K2在结构学习里长盛不衰的原因——它把NP-hard问题拉回了一个普通笔记本能跑完的规模。

和全局枚举比，K2唯一的妥协是贪心。全局枚举在节点数超过15以后基本跑不动；K2在100个节点、每个节点最多3个父节点的设定下，十几分钟能出一版结构，这个量级正好是工业场景里最常见的大小。

2.2 BDeu打分函数：为什么它能同时奖励拟合和惩罚复杂度

打分函数最通用的是BDeu（Bayesian Dirichlet equivalent uniform）。它计算的是“给定这个父节点组合，子节点的观测数据在Dirichlet先验下出现的边际似然”。单个节点Xi的BDeu对数得分如下：

score(Xi|πi) = Σ_j [ lnΓ(αij) − lnΓ(αij + Nij) + Σ_k ( lnΓ(αijk + Nijk) − lnΓ(αijk) ) ]

其中πi是Xi的父节点集合；qi是父节点状态组合数；ri是Xi的状态数；Nijk是“父组合为j且Xi取状态k”的样本数；α是等效样本量，决定先验的强度，默认取1.0；αij = α/qi，αijk = α/(qi·ri)。Γ是伽马函数，代码里用math.lgamma取对数再算，避免大数溢出。

直觉上，第一项随样本量增大而减小，给“额外多出来的父节点组合”施压；后一项随拟合精度提高而增大，奖励父节点把子节点状态“解释”得更集中。两者对抗，就避免了纯似然打分那种“边加越多越好”的毛病。参数α越大，先验越强，越不爱加边；α越小，越倾向于相信数据、加边更积极。这个参数后面调参会用到。

顺带说一句，原始K2论文里用的K2评分是BD家族在超参数取1时的特例，社区实现一般直接拿BDeu替代，本文代码也用BDeu，公式更通用，效果等价。

2.3 贪心搜索：K2怎么把搜索空间切碎

K2的搜索策略可以一句话说清：从空父集开始，逐个尝试加入“还没加入、且排在当前节点之前”的候选父节点，选得分增量最大的那个加进去，重复到没有正增益或父节点数到上限为止。每个节点的父节点搜索是独立的，互不干扰，所以整个算法天然可以并行，把每个节点的父节点搜索丢到多进程里就行。

停止条件有两个：增益不再为正，或者父节点数量到了上限max_parents。前者常见于弱依赖，后者常见于样本量大、信号强的场景。贪心的代价是可能错过“加一条边得分变差、再加一条反而更好”的组合，这个问题第五部分专门讲。

2.4 K2和PC、MMHC怎么选：一张选型表

注意这表的“要不要节点顺序”不是优点和缺点，而是使用条件。有可靠的领域顺序就用K2，没有就先用Tabu或PC出一个粗略方向，再回过头来定序，这也是第四部分要展开的做法。

3. 用可运行的Python实现K2：离散化、BDeu打分与贪心搜索

3.1 数据准备：连续变量离散化与样本量底线

K2吃的是离散数据。连续变量要么等宽切段，要么分位数切段，我一般优先分位数，因为能保证每档样本量均匀，打分时不至于出现空组合。下面这段把pandas的DataFrame切成整型状态表：

import pandas as pd import numpy as np df = pd.read_csv("your_data.csv") # 常量列对打分没有贡献，先剔除 df = df.loc[:, df.nunique() > 1] # 分位数离散化：每列切成3档，产出0/1/2整型 for col in df.columns: df[col] = pd.qcut(df[col], q=3, labels=[0, 1, 2], duplicates="drop") X = df.values.astype(int) n_states = [len(np.unique(X[:, i])) for i in range(X.shape[1])] print("状态数:", n_states)

pd.qcut按分位数切，比等宽切对离群值更稳；duplicates="drop"处理某列重复值过多导致分位点重合的情况。样本量底线我一般按经验卡两条：总样本数不少于节点数的100倍；学完后每个父节点组合的平均样本数不少于5，否则打分里的计数矩阵会大量出现0，分数值失真。

3.2 BDeu打分函数先落地

import math import itertools import numpy as np def bdeu_score(data, node, parents, n_states, alpha=1.0): """计算单个节点的BDeu对数得分，parents是父节点下标列表""" N = data.shape[0] r_i = n_states[node] if len(parents) == 0: # 无父节点时 q=1, 直接统计节点自身状态分布 counts = np.bincount(data[:, node], minlength=r_i) aij = alpha aijk = alpha / r_i score = math.lgamma(aij) - math.lgamma(aij + N) for k in range(r_i): score += math.lgamma(aijk + counts[k]) - math.lgamma(aijk) return score # 把父节点状态组合编码为整数key，枚举所有组合 strides = [] stride = 1 for p in parents: strides.append(stride) stride *= n_states[p] q = stride keys = np.zeros(N, dtype=int) for p, s in zip(parents, strides): keys += data[:, p] * s combo_counts = np.zeros((q, r_i), dtype=int) for row, key in enumerate(keys): combo_counts[key, data[row, node]] += 1 combo_total = combo_counts.sum(axis=1) aij = alpha / q aijk = alpha / (q * r_i) score = 0.0 for j in range(q): if combo_total[j] == 0: continue score += math.lgamma(aij) - math.lgamma(aij + combo_total[j]) for k in range(r_i): if combo_counts[j, k] > 0: score += math.lgamma(aijk + combo_counts[j, k]) - math.lgamma(aijk) return score

这段代码的核心是用累乘stride把父节点状态组合压成一个整数下标，避免嵌套循环枚举组合。alpha是等效样本量，默认1.0；调大它会让算法更不愿加边，适合样本量偏小的场景。combo_total[j] == 0时整项为0，跳过不计算，省时间也避免lgamma(0)问题。打分是分解式的，每个节点独立评分，整网得分等于所有节点得分之和，所以后续贪心搜索可以按节点逐个做。

3.3 贪心搜索主循环：节点顺序、候选父节点与停止条件

def k2_fit(data, order, n_states, alpha=1.0, max_parents=3): """ K2主算法 order: 节点顺序列表，例如 [0,2,1] 表示节点0最早，节点1最晚 只有排在node前面的节点才允许成为node的父节点 """ n = data.shape[1] dag = {i: [] for i in range(n)} pos = {node: idx for idx, node in enumerate(order)} for node in range(n): allowed = [p for p in order if pos[p] < pos[node]] cur_score = bdeu_score(data, node, [], n_states, alpha) improved = True while improved and len(dag[node]) < max_parents: improved = False best_p = None best_gain = 0.0 for p in allowed: if p in dag[node]: continue cand_score = bdeu_score(data, node, dag[node] + [p], n_states, alpha) gain = cand_score - cur_score if gain > best_gain: best_gain = gain best_p = p if best_p is not None and best_gain > 1e-8: dag[node].append(best_p) cur_score += best_gain improved = True return dag order = list(range(X.shape[1])) # 先用自然顺序跑一版，第4章专门讲定序 dag = k2_fit(X, order, n_states, alpha=1.0, max_parents=3) for node, parents in dag.items(): print(f"节点{node}的父节点: {parents}")

主循环里每一轮只保留增益最大的那个候选父节点，然后重新评估剩余候选。best_gain > 1e-8这层容差避免浮点噪声导致的无效加边。max_parents=3是工业场景比较稳妥的上限，超过3个父节点后参数数量指数上涨，打分开始惩罚复杂度，而且业务上解释起来也费劲。注意这里allowed列表只依赖order，K2强顺序假设的全部含义就在这一行。

3.4 结果校验：BIC到底够不够好

学完结构别急着交差。我一般同时算两个指标：整网对数似然LL，以及BIC = LL − (d/2)·lnN，其中d是参数个数，公式为d = Σ_i (r_i − 1)·Π_{p∈Pa_i} r_p。BIC和BDeu的排序在样本量足够时基本一致，所以BIC可以作为K2结果的快速体检表——你不需要额外重学，只需把K2学出的图拿去算BIC，再和Tabu搜索算出的BIC对比。如果K2的BIC比Tabu低太多，说明贪心大概率停在局部最优，回到第2.3节那个已知代价。

4. 节点顺序怎么定：三种方案与我的效果排序

4.1 为什么说节点顺序是K2的命门

K2算法里，节点顺序不是锦上添花，而是把搜索空间砍掉大半的硬约束：一条边只能从顺序靠前的节点指向顺序靠后的节点。这个假设如果和真实因果方向一致，算法很快收敛；如果反了，哪怕数据信号再强，K2也会强行给出一个与事实相反的方向。我见过不少项目在打分函数、离散化上调来调去，最后发现问题出在最初顺手写的那个order上。所以定序这件事值得花整块时间做，下面三个方案按可靠程度从高到低排。

4.2 方案一：领域专家直接给顺序，最可靠但最容易被低估

故障诊断、基因调控、用户转化这类场景里，专家往往能直接说出“先有因、后有果”的链条：比如设备温度升高在前，报警触发在后；某个基因表达变化在前，蛋白浓度变化在后。把这些链条汇总成一个全序或偏序，就是K2最理想的食物。做法是把变量按“原因优先”排成列表，原因排前面，结果排后面。专家给序时经常只给偏序，这时候我一般保留偏序，把无先后关系的变量随机放在中间位置，多跑几轮取总分最高的一组顺序，而不是硬凑一个全序。

4.3 方案二：先学一个粗略DAG再拓扑排序，数据定序的兜底

没有专家知识时，常见做法是先用Tabu搜索或PC算法跑一个低配结构，然后对这个DAG做拓扑排序，把拓扑序作为K2的order。这个方案的好处是顺序完全从数据来，坏处是Tabu或PC本身的错误会带进order，形成误差传递。我的建议是把得到的DAG先人工看一眼，删掉明显不符合常识的边，再做拓扑排序；如果DAG里有环，拓扑排序会失败，这时需要用贪心破环或直接换方案三。

4.4 方案三：互信息启发式排序，快速但粗糙

还有一种完全不用学结构的办法：计算每对变量的互信息，把每个变量与其他所有变量的互信息之和当成“依赖强度”，依赖强度大的排前面。这个假设是“越像因的变量，跟全系统的关联越强”，在部分场景下成立，但不保证方向正确。代码实现很短：

def discrete_mi(x, y, n_states_x, n_states_y): pxy = np.zeros((n_states_x, n_states_y)) for i in range(len(x)): pxy[x[i], y[i]] += 1 pxy /= len(x) px = pxy.sum(axis=1, keepdims=True) py = pxy.sum(axis=0, keepdims=True) return float((pxy * np.log(pxy / (px * py) + 1e-12)).sum()) n = X.shape[1] mi_mat = np.zeros((n, n)) for i in range(n): for j in range(i + 1, n): v = discrete_mi(X[:, i], X[:, j], n_states[i], n_states[j]) mi_mat[i, j] = mi_mat[j, i] = v strength = mi_mat.sum(axis=1) heuristic_order = np.argsort(-strength).tolist() print("互信息启发式顺序:", heuristic_order)

这个顺序只作初始值，我一般会拿它跑一次K2，再把结果和专家序的结果对比。对比时主要看BIC和边的方向是否符合业务预期，两者都过不了再换方案二。三种方案效果排序上，我的经验是专家序明显强于数据序，数据序里方案二强于方案三；但方案三胜在快，适合快速试错阶段。

5. K2算法避坑指南：5个让结构学习翻车的真实场景

5.1 样本量不足时BDeu也会“瞎连”

现象：学出来的网络几乎全连接，每个节点都挂了三四个父节点，但看数据根本没那么强的关系。原因很简单：样本量小的时候，打分里的拟合项增长比复杂度惩罚快，BDeu被“数据太少、计数稀疏”带偏，相当于每多一条边都能靠少数几个样本白捡似然增益。解决：把alpha从1.0往大调，比如alpha=10甚至alpha=N/10，先验强度大了就不爱加边；同时把max_parents压到2；最后用Bootstrap看边频率，频率低于0.5的边直接删。调整alpha后重新学的结构，边数通常会明显下降。

5.2 离散化边界一变，结构天翻地覆

现象：同样的数据，用等宽切3段和分位数切3段，学出来的边有近一半不一样；切4段又变一版。这不是K2的锅，是所有离散结构学习方法的通病：离散化把连续变量的信息截断了，边界移动改变了互信息和条件独立关系。解决：固定分位数策略，别换来换去；每列状态数尽量一致，不要有的列切3段有的列切8段，状态数差异过大会让打分比较失衡；对最终结构做敏感性分析——换个bins数重学一遍，只保留在多次离散化下都出现的边，这部分边才是真正稳的。

5.3 节点顺序给反，有向边整体掉头

现象：已知业务上是A→B，但K2学出B→A，而且BIC还提示学出来的方向更优。原因不是打分错，而是order里B排在A前面，K2只允许B做A的父节点，数据也支持B→A这个方向，于是结构就翻了。这类问题最难排查，因为BIC正常、边也合理，只有方向反了。解决：定序之前把已知的因果链条写出来，一条条比对新学结构；有争议方向的边单独拿出来做干预实验或A/B验证，不要指望纯观测数据能定方向。

5.4 独立节点变成孤儿：是好事还是数据问题

现象：某个节点从头到尾没有父节点也没有子节点，整张图它是孤立点。原因有两种：它确实独立于其他变量，这是正确结果；但也可能是离散化后该列信息量太低，比如90%样本落在同一档，互信息趋近于0，算法自然找不到边。解决：先查该列的状态分布，如果某一档占比超过95%，这个变量在离散化后基本废了，考虑换离散化参数或扔掉；如果分布正常仍是孤立点，保留它作为独立节点，这也是结构学习给出的有效结论。

5.5 贪心走不到“加两条边更优”的局部解

现象：K2学完的结构，把某条缺失的边手工加上去之后BIC反而涨了，但K2当时没选它。原因是贪心顺序问题：单独加那条边时增益为负，先加另一条边再一起加才是正的，贪心每次只走一步，永远到不了这个局面。解决：把K2学出的结构作为初始解，再用Tabu搜索或随机重启的Hill Climbing继续优化，让禁忌表能接受“暂时变差”的移动；或者多跑几组随机order，取BIC最高的结果。K2在这个流程里充当的角色是给后续优化提供一个高质量起点。

6. 把K2学到的结构用进动态贝叶斯网络：先验网络与Bootstrap验证

6.1 DBN的两段式建模：先验网络与转移网络

动态贝叶斯网络（DBN）把静态结构扩展到时序上：t时刻的变量可以依赖t−1时刻的变量，同时保留片内的静态依赖。工程上常见的做法是先做静态结构学习，把K2学到的网络直接作为先验网络固定下来，再去学片间转移边。这样既避免DBN全量搜索的组合爆炸，也让静态结构学习的结果有了实际去向。转移网络的处理可以这样：

# X_prev: 上一时刻整型状态矩阵，X_now: 当前时刻整型状态矩阵 # 两者行数相同，列数都是n X_pair = np.hstack([X_prev, X_now]) # 左边是t-1片，右边是t片 n_states_pair = n_states + n_states pair_order = list(range(n)) + [n + i for i in range(n)] # t-1片全部排前面 pair_parents = k2_fit(X_pair, pair_order, n_states_pair, alpha=1.0, max_parents=3) # 只保留跨片边：父节点在左半，子节点在右半 transfer_edges = [] for node in range(n, 2 * n): for p in pair_parents[node]: if p < n: transfer_edges.append((p, node - n)) print("转移边:", transfer_edges)

这里的pair_order保证片间边只能从右往左，也就是从t−1时刻指向t时刻，满足一阶马尔可夫假设。

6.2 用Bootstrap检验边的可信度：阈值取多少

K2跑一次的结构不能直接当结论用，我习惯用Bootstrap重采样给每条边打可信度分：

from sklearn.utils import resample n = X.shape[1] edge_freq = np.zeros((n, n)) B = 200 for _ in range(B): idx = resample(range(X.shape[0]), replace=True) boot_X = X[idx] boot_dag = k2_fit(boot_X, order, n_states, alpha=1.0, max_parents=3) for node, parents in boot_dag.items(): for p in parents: edge_freq[p, node] += 1 edge_prob = edge_freq / B stable_edges = [(p, node) for node in range(n) for p in range(n) if node != p and edge_prob[p, node] > 0.7] print("可信边:", stable_edges)

阈值取值经验：样本量在几千以下时用0.5保边数量，样本过万且信号强时用0.7以上。只把Bootstrap后仍然稳定的边用于下游模型，剩下的边留给业务方去确认。

6.3 我现在的落地习惯

我现在的落地流程是：先用专家知识或粗糙Tabu定序，K2跑200轮Bootstrap，频率低于0.5的边直接砍掉，剩下的结构交给评审看方向；涉及时序数据时，再用这个静态结构做DBN的先验网络去学片间边。K2最大的价值不是一次给出最终答案，而是给一个稳定、可解释、能迭代的起点，后续所有结构优化都从它出发。这个流程救过我很多次，希望帮到你。

本文还有配套的精品资源，点击获取