简介：这份资源围绕K2算法从数据中学习贝叶斯网络结构这一经典问题展开，面向具备一定概率图模型基础、希望动手实现结构学习的研究生与算法工程师。包内共7个文件，以4个m脚本为主体，配合1个c源码、1个mat数据文件和1个txt说明，压缩包约10KB，体量轻便，便于快速阅读与二次开发。其中m脚本承担主流程控制、评分函数与对象构造等职责，c源码给出K2核心计算的底层实现，mat文件提供可直接加载的样本数据，txt则补充许可与使用说明。已有2168人学习下载，说明该实现具备一定的参考价值。读者可借此理解K2算法在给定节点顺序与最大父节点数约束下逐步搜索最优父集的过程，掌握评分函数与候选父集评估的配合方式，并对照源码完成结构学习实验、结果复现与参数调整，为后续研究或工程落地提供可运行的起点。
1. 从一张“跑不通”的评分表说起：K2 算法到底在算什么
很多人第一次接触贝叶斯网络结构学习，是从一份评分函数表开始的。表里每个节点、每个候选父节点集合都对应一个分数，看起来只要挑分数最高的组合就行。但真把数据丢进去，跑出来的网络要么边多到像毛线团，要么关键依赖全被吞掉。问题往往不在数据，而在搜索策略和评分口径没对齐。K2 算法就是在这个环节里被反复提起的名字：它不靠条件独立性检验逐对“审问”变量，而是用贝叶斯评分给候选结构打分，再配合一个固定的节点顺序做贪心搜索。适合谁？手里有离散变量、样本量不算夸张、又需要一个能复现、能解释的结构学习基线的人。它不保证全局最优，但胜在流程透明、参数少、每一步都能打印出来看。下面按“评分怎么算、顺序怎么定、代码怎么落、坑怎么排”四段推进，中间给可直接抄的 Python 实现。
2. K2 评分与贪心搜索：把结构学习拆成可打印的步骤
2.1 评分函数在算什么：从边缘似然到可加性
K2 评分本质上是贝叶斯评分的一种简化形式。对每个节点 $X_i$ 和它的父节点集合 $\pi_i$，评分假设参数先验是均匀的，然后对参数做积分，得到该节点在给定父节点下的边缘似然。把所有节点的局部评分加起来，就是整个结构的 K2 分数。写成式子大致是：
$$ \text{K2}(G) = \sum_{i=1}^{n} \sum_{j=1}^{q_i} \left[ \log \frac{(r_i - 1)!}{(N_{ij} + r_i - 1)!} + \sum_{k=1}^{r_i} \log N_{ijk}! \right] $$
其中 $r_i$ 是节点 $X_i$ 的取值个数，$q_i$ 是父节点组合数，$N_{ijk}$ 是父节点取第 $j$ 种组合时 $X_i$ 取第 $k$ 个值的样本数，$N_{ij}$ 是这些样本数的和。这个式子看起来吓人，但落到代码里就是数频次、算阶乘对数。关键点在于：评分是可分解的，每个节点独立算自己的局部分数，改一条边只影响两个节点的局部分数。这也是贪心搜索能高效跑起来的前提。
常见做法是先把所有变量离散化，再统计频次表。如果变量是连续的，要么先分箱，要么换别的评分函数，K2 本身不处理连续变量。分箱的边界怎么定，后面避坑章节会专门说。
2.2 节点顺序为什么是命门：先验顺序与搜索空间
K2 算法要求你提供一个节点的拓扑顺序。这个顺序决定了每个节点只能从它前面的节点里选父节点，不能反向选。换句话说，你给的不是“建议”，而是硬约束。顺序对了，搜索空间被砍掉一大半，贪心几步就能收敛；顺序错了，正确的边永远进不了候选集，评分再准也白搭。
那顺序从哪来？常见有三条路：一是领域知识，比如时间上先发生的变量放前面；二是用其他方法先跑一个粗略排序，比如基于互信息或条件独立性做一版；三是直接按变量名或列顺序硬上，然后看结果是否合理。我一般会先用互信息排一版，再人工调整个别明显反直觉的位置。注意，K2 对顺序敏感，但也不是完全不能错——如果两个变量真实存在依赖，但顺序反了，算法会漏掉这条边，而不是报错。所以顺序验证要单独做一步。
2.3 贪心加边：每一步只做一件事
有了评分函数和节点顺序，搜索过程很直白：从空图开始，按顺序遍历每个节点，尝试把当前节点之前的某个节点加为父节点，算局部评分增量，如果增量为正就保留，否则跳过。每个节点只往前看，不回头改已经定下的父节点集合。伪代码大致是：
for i in range(n): parents = [] score_old = local_score(i, parents) while True: best_gain = 0 best_parent = None for j in range(i): # 只能选顺序在前的节点 if j in parents: continue score_new = local_score(i, parents + [j]) gain = score_new - score_old if gain > best_gain: best_gain = gain best_parent = j if best_parent is None: break parents.append(best_parent) score_old = local_score(i, parents) final_parents[i] = parents这段逻辑里有两个容易忽略的点。第一，best_gain初始为 0，意味着只有严格提升才加边，评分持平不加。第二，每加一条边就重新算局部评分，而不是一次性算所有候选。这样写虽然多算几次，但逻辑清晰，也方便打印每一步的增量。实际工程里可以缓存频次表，避免重复统计。
2.4 一个可复现的最小实现
下面给一个不依赖第三方贝叶斯网络库的最小实现，只用 numpy 和 pandas。数据假设是离散的，每列取值从 0 开始连续编号。
import numpy as np import pandas as pd from math import lgamma def count_table(data, node, parents): """统计 N_ijk：父节点组合 × 当前节点取值""" if not parents: counts = np.bincount(data[:, node]) return counts.reshape(1, -1) parent_vals = data[:, parents] # 把父节点组合编码成单一整数 comb = np.zeros(len(data), dtype=int) for p in parent_vals.T: comb = comb * (p.max() + 1) + p n_comb = comb.max() + 1 r = data[:, node].max() + 1 table = np.zeros((n_comb, r)) for c, v in zip(comb, data[:, node]): table[c, v] += 1 return table def k2_local_score(data, node, parents): """计算单个节点的 K2 局部评分""" table = count_table(data, node, parents) r = table.shape[1] score = 0.0 for row in table: n_ij = row.sum() score += lgamma(r) - lgamma(n_ij + r) # 先验均匀时的常数项 for n_ijk in row: score += lgamma(n_ijk + 1) return score def k2_structure(data, order): """K2 贪心搜索主流程""" n = data.shape[1] final_parents = [[] for _ in range(n)] for idx, node in enumerate(order): parents = [] score_old = k2_local_score(data, node, parents) while True: best_gain = 0.0 best_parent = None for j in order[:idx]: if j in parents: continue score_new = k2_local_score(data, node, parents + [j]) gain = score_new - score_old if gain > best_gain: best_gain = gain best_parent = j if best_parent is None: break parents.append(best_parent) score_old = k2_local_score(data, node, parents) print(f"节点 {node} 加入父节点 {best_parent}，增量 {best_gain:.4f}") final_parents[node] = parents return final_parents参数说明：data是二维整数数组，每列一个变量，取值从 0 开始；order是节点索引的列表，表示拓扑顺序。count_table里用乘加编码父节点组合，要求每个父节点取值是紧凑整数，如果原始数据有缺失或跳号，先做映射。k2_local_score里用lgamma代替阶乘对数，避免大数溢出。k2_structure打印每一步的增量，方便你判断是不是加边加过头。
跑完之后，final_parents就是每个节点的父节点列表，转成边集就能画图。如果发现某个节点父节点特别多，先别急着高兴，大概率是顺序或分箱出了问题，下一章会展开。
3. 数据准备与顺序验证：让 K2 跑出可解释的结构
3.1 离散化不是随便分箱：等频、等宽与领域切点
K2 评分建立在频次统计上，连续变量必须先离散化。常见做法有三种：等宽分箱、等频分箱、按领域切点分箱。等宽简单，但对偏态分布不友好，极端值会把大部分样本挤进一两个箱；等频每箱样本数接近，但切点会随数据变动，跨数据集复现时要注意；领域切点最稳，但需要你先验知识。
我一般会先看每个变量的分布直方图，如果明显双峰，就按谷底切；如果近似均匀，等频三到五箱就够。箱数不是越多越好，箱数多会让频次表稀疏，K2 评分里的阶乘项对稀疏计数很敏感，容易出现“假高分”。一个经验规则：每个父节点组合下，平均样本数最好不低于 5。如果变量多、样本少，宁可减少箱数，也别硬上高分辨率。
def discretize_equal_freq(series, bins=3): """等频离散化，返回 0 起始的整数编码""" labels = pd.qcut(series, q=bins, labels=False, duplicates='drop') return labels.fillna(0).astype(int).values这段代码用qcut做等频切分，duplicates='drop'处理重复边界。注意返回的编码可能不连续，如果后续count_table要求紧凑整数，需要再做一次映射。
3.2 顺序从哪来：互信息排序与人工校正
没有领域顺序时，我通常先用互信息给变量排一版。互信息高的一对变量，谁先谁后不好定，但可以把互信息总和高的变量往前放。具体做法是算两两互信息矩阵，按行和降序排。这只是一个启发式，不保证正确，但比随机顺序强。
def mutual_info_order(data): """按互信息和降序给出一个初始顺序""" n = data.shape[1] mi = np.zeros((n, n)) for i in range(n): for j in range(i + 1, n): table = count_table(data, i, [j]) # 简化的互信息估计，实际可用 sklearn 的 mutual_info_score p = table / table.sum() pi = p.sum(axis=1, keepdims=True) pj = p.sum(axis=0, keepdims=True) with np.errstate(divide='ignore', invalid='ignore'): term = p * np.log(p / (pi * pj)) mi[i, j] = mi[j, i] = np.nansum(term) return list(np.argsort(-mi.sum(axis=1)))排完之后，人工扫一眼：如果某个变量明显是结果变量却排在前面，手动往后挪。顺序调整后重新跑 K2，对比边集变化。如果挪一个变量导致大量边消失，说明原来的顺序可能靠巧合撑住了评分，需要更谨慎地验证。
3.3 评分增量打印与边集过滤
K2 跑完会给出每个节点的父节点列表，但并不是所有边都值得保留。常见过滤手段有两种：一是设一个最小增量阈值，低于阈值的边即使评分为正也砍掉；二是用 bootstrap 重采样，看每条边出现的频率，频率低的边视为不稳定。
我一般先看增量打印。如果某条边的增量只有 0.01 量级，而其他边都在 1 以上，这条边大概率是噪声。可以设一个相对阈值，比如最大增量的 5%。另外，如果两个变量在顺序上相邻且互信息很高，K2 可能给出双向边中的一条，这时候要结合领域判断方向是否合理。
def filter_edges(parents, min_gain=0.05, gains=None): """按增量阈值过滤边，gains 是 (child, parent, gain) 列表""" if gains is None: return parents kept = [[] for _ in parents] max_gain = max(g for _, _, g in gains) if gains else 1.0 for child, parent, gain in gains: if gain >= min_gain * max_gain: kept[child].append(parent) return kept这段过滤逻辑需要你在搜索过程中记录每条边的增量。如果不想改主流程，也可以在跑完后用局部评分重新算一遍每条边的增量，再过滤。
3.4 用模拟数据验证：已知真值下的对照
验证 K2 实现是否正确，最直接的办法是生成一批已知结构的模拟数据，看能不能恢复出来。下面生成一个简单的链式结构 A→B→C，每个变量三取值，样本量 2000。
np.random.seed(42) n_samples = 2000 A = np.random.randint(0, 3, n_samples) B = (A + np.random.randint(0, 3, n_samples)) % 3 C = (B + np.random.randint(0, 3, n_samples)) % 3 data = np.column_stack([A, B, C]) order = [0, 1, 2] parents = k2_structure(data, order) print(parents)如果实现正确，输出应该是[[], [0], [1]]或接近的形式。如果 B 的父节点里没有 A，或者 C 的父节点里多了 A，说明评分或顺序有问题。这个对照实验建议每改一次代码就跑一遍，比看理论公式快得多。
4. 避坑与排查：K2 跑不通时先看这五条
4.1 现象：所有节点都没有父节点，边集为空
原因通常是评分增量始终为负或零。K2 评分在样本量小、箱数多的时候，先验项会压过数据项，导致加边不划算。另外，如果数据里变量之间确实独立，空图就是正确结果，不一定是 bug。
解决：先检查样本量和箱数的比例，每个父节点组合下平均样本数低于 5 就减少箱数。然后手动算一对明显相关变量的局部评分增量，如果增量仍为负，检查频次表编码是否正确，特别是父节点组合编码有没有溢出或错位。
4.2 现象：某个节点父节点特别多，边集密集
原因多半是节点顺序把结果变量放在了前面，导致后面的变量拼命往前连。另外，如果离散化箱数过多，评分对稀疏频次的惩罚不够，也会鼓励多加边。
解决：调整顺序，把疑似结果变量往后挪。同时降低箱数，或者引入最小增量阈值。还可以用 bootstrap 看边频率，把频率低于 0.6 的边砍掉。我一般会先跑一版不加阈值的结果，再跑一版加阈值的结果，对比哪些边是稳定的。
4.3 现象：正确边方向反了
K2 只能按给定顺序选父节点，如果真实方向是 B→A，但顺序里 A 在 B 前面，算法只会考虑 A 作为 B 的父节点，不会反向。结果就是方向反了或者边消失。
解决：这不是算法 bug，是顺序约束的必然结果。要么调整顺序，要么换用不依赖顺序的搜索策略。如果领域上确实无法确定方向，可以在结果里标注“方向未验证”，不要强行解释。
4.4 现象：评分增量打印正常，但最终边集和预期差很多
原因可能是频次表统计时父节点组合编码有误，比如两个父节点取值组合被映射到同一个编码，导致计数合并。另外，如果数据里有缺失值被当成一个合法取值，也会污染频次表。
解决：在count_table里加断言，检查编码后的组合数是否等于各父节点取值数的乘积。缺失值在离散化之前就要处理掉，不要留到频次统计阶段。可以用data.isna().sum()先扫一遍。
4.5 现象：换一批数据后结果完全变样
原因可能是离散化切点依赖当前数据分布，换数据后切点漂移，导致频次表结构变化。等频分箱尤其明显。
解决：如果要做跨数据集对比，尽量用固定切点，比如按领域阈值或第一版数据的切点固化下来。等频分箱只适合探索性分析，不适合需要复现的场景。另外，样本量差异大时，评分尺度也会变，增量阈值要重新校准。
5. 进阶技巧：用 bootstrap 边频率给结构加一层稳定度
K2 跑一次的结果受顺序和分箱影响很大，单次边集只能当参考。我后来养成的习惯是：不管多信任当前顺序，都跑 200 次 bootstrap，统计每条边出现的频率，只保留频率高于 0.7 的边。这一步多花几分钟，但能过滤掉大量靠巧合撑起来的边。
具体做法是对原始数据有放回重采样，每次重采样后跑一遍 K2，记录边集。然后按边聚合频率。下面是一个可复用的 bootstrap 封装：
def bootstrap_k2(data, order, n_boot=200, seed=0): """bootstrap 重采样跑 K2，返回边频率字典""" rng = np.random.default_rng(seed) edge_count = {} n = data.shape[1] for b in range(n_boot): idx = rng.integers(0, len(data), len(data)) sample = data[idx] parents = k2_structure(sample, order) for child in range(n): for parent in parents[child]: edge = (parent, child) edge_count[edge] = edge_count.get(edge, 0) + 1 freq = {e: c / n_boot for e, c in edge_count.items()} return freq def stable_edges(freq, threshold=0.7): """按频率阈值筛选稳定边""" return {e: f for e, f in freq.items() if f >= threshold}参数说明：n_boot是重采样次数，200 次在中小规模数据上通常够用；threshold是频率阈值，0.7 偏保守，0.5 偏宽松，按你对假阳性的容忍度调。bootstrap_k2里每次重采样都重新跑完整搜索，计算量是单次的n_boot倍，如果数据大，可以先把频次表缓存起来，或者减少n_boot到 50 次先看趋势。
拿到稳定边集后，和单次结果对比。如果某条边单次出现但 bootstrap 频率只有 0.3，基本可以判定是噪声。反过来，如果某条边单次没出现但 bootstrap 频率很高，说明顺序或分箱可能压制了它，值得回头检查。
还有一个细节：bootstrap 重采样会改变样本分布，等频分箱的切点如果每次重新算，会引入额外波动。我一般会在 bootstrap 之前先把离散化切点固定下来，重采样只对已经离散化的数据做。这样边频率反映的是结构稳定性，而不是分箱稳定性。
从那以后我每次跑 K2，都强制走一遍“固定切点 → 互信息排序 → 单次搜索 → bootstrap 过滤”的流程，不再直接拿单次结果下结论。希望帮到你。
本文还有配套的精品资源，点击获取
