简介针对贝叶斯网络结构学习任务这份资源围绕K2算法提供从数据中自动学习网络结构的Matlab实现适用于概率图模型初学者、研究生及相关算法开发者可作为课程作业与算法研究的入门参考。包内共7个文件4个.m源码分别承担主控流程、子函数和图对象构造等任务1个.c文件为MEX接口经编译后可提升核心计算速度1个.mat样例数据可直接加载用于验证算法效果1个.txt为许可说明。整体压缩包仅10KB紧凑清晰便于按行阅读和二次改造文件划分简洁便于按功能定位代码。目前已有2168人学习过该资源。借助源码、MEX接口和样例数据读者可深入理解K2算法如何利用评分函数并通过贪心搜索逐步生成贝叶斯网络结构同时也能借鉴Matlab与C混合编程的实际用法为课程设计、科研实验或算法改进提供直接参考。1. 贝叶斯网络结构学习是什么K2算法为什么还有人在用一位做风险分析的同行曾跟我抱怨业务方指着二十多个指标问他哪些会相互影响他只能靠领域经验手动画依赖图画到一半自己都信不过。这其实就是贝叶斯网络结构学习要解决的事——从观测数据里自动推出变量之间的依赖关系图而K2算法是这条路上历史最久、也最容易落地的入口之一。它的名字来自上世纪90年代初一篇经典文献核心思路是给变量一个大致先后顺序用评分函数衡量「某个节点加上某几个父节点后模型和数据的匹配程度」再一步步贪心地补边。下面按原理、代码、参数、避坑的顺序一步步拆适合刚接触贝叶斯网络、手里有离散或可离散化数据、想快速得到一张可解释依赖图的从业者。2. 从评分到搜索K2算法的工作机制与选型理由2.1 K2解决的问题结构学习该走哪条路贝叶斯网络由两部分组成结构DAG即有向无环图和参数每个节点在给定父节点取值下的条件概率表。结构学习就是只从观测数据里推出这张图长什么样参数学习是在结构确定之后去估计每个节点的条件概率。K2算法管的是前半段而且只针对离散数据——它的统计基础是多项式分布和Dirichlet先验数据里出现小数或者NaN第一步就会卡住。结构学习业内分两条路线。一条是基于约束的典型代表是PC算法反复做条件独立性检验保留那些“在控制了某些变量之后仍然相关”的边最后再定向。这条路的问题在于高阶独立性检验非常吃样本量样本不够时检验结果互相矛盾算法只能靠某种启发式去调和稳定性比较让人头疼。另一条就是K2所在的基于评分的路线定义一个评分函数衡量“这张图数据”有多匹配然后在图空间里搜索高分结构。评分函数给出方向搜索策略决定效率两件事分开逻辑上比约束法好理解得多。K2在这条路线里的定位很特别。它是最早把评分和搜索组合成完整方案的算法之一至今仍是很多结构学习论文的baseline。它的关键设计——输入节点顺序、限制每个节点的父节点数、每步只贪心选增益最大的父节点——让它在几十个变量的规模下依然跑得动结果可复现还方便把领域经验塞进去。作为入门第一个结构学习算法K2的性价比很高。2.2 K2评分为什么局部打分能代表全局好坏K2的评分函数不是随便定义的经验值它来自贝叶斯后验。假设所有可能的结构先验相等那么结构G的后验概率正比于边际似然P(D|G)也就是对参数求积分之后的数据似然。离散数据配合Dirichlet先验时这个积分有解析解。K2原始论文用的先验是每个参数都取1的均匀先验跟后来常用的BDeu评分把等价样本量设为1完全等价——这一点在实现上非常讨巧因为主流库里的BDeu评分可以直接当K2评分用。把它取对数之后单个节点的局部评分长这样log g(i, π_i) Σ_j [ log Γ(r_i) - log Γ(N_ij r_i) Σ_k log Γ(N_ijk 1) ]其中符号的含义是符号含义r_i节点X_i的状态数也就是这个变量取几个离散值N_ijkX_i取第k个状态、父节点组合为第j种组合时样本里出现这种组合的个数N_ij父节点组合为第j种时所有样本数等于Σ_k N_ijkπ_i节点X_i的父节点集合总评分就是所有节点局部评分之和。这个性质决定了整条路线的走向既然总评分能按节点拆开那么搜图的时候就可以一个节点一个节点地处理不需要每次候选动作都重新计算整张图的分。K2算法的每一步都是在当前节点的候选父集合里找一个能让这个局部评分涨得最多的父节点——评分可分解贪心搜索才站得住脚。另一个值得注意的点是评分函数天然包含对复杂度的惩罚。父节点越多、状态组合越细数据被切得越碎Γ(N_ij r_i)这项就会把分数拉低。所以K2不会无限加父节点即便不设上限加过头了评分也会掉头向下。这个内在惩罚机制让K2对过拟合有一定自愈能力。2.3 贪心搜索K2怎样把指数级搜索变成线性近似DAG的数量随节点数超指数增长5个节点的DAG已经有几万种10个节点以上完全不可能穷举。K2没有去搜全空间而是用两个约束把空间砍到可控第一用户必须提供一个节点顺序算法只允许从“顺序靠前”的节点指向“顺序靠后”的节点第二每个节点的父节点数上限为u超过固定值就停。这两个约束听起来粗暴但正是它们让K2在几十个节点下依然有实用价值。搜索过程非常朴素。按给定顺序遍历每个节点对当前节点X把所有排在它前面的节点作为候选父集合初始父集合为空。每一轮逐个尝试把某个候选父节点加入父集合用局部评分算出增益保留增益最大的那个如果最大增益为正就真的加进去否则停止继续处理下一个节点。每节点最多u轮每个候选父至多试一次整体计算量大致在O(n²u)量级和穷举不可同日而语。用三个节点A、B、C做例子如果顺序是A在前、B居中、C最后那么合法的边只有A→B、A→C、B→C三种可能。算法永远不会生成C→A这样的边不管数据里这两个变量的相关性有多强。这就是K2把先验变成硬约束的方式——顺序一旦给定搜索空间就被圈死了代价是顺序错了所有结果都会跟着错。这个特性既是K2的高效来源也是它最大的坑后面避坑章会专门展开。3. 用pgmpy跑通K2结构学习离散化、代码与参数调优3.1 数据准备离散化与缺失值处理K2的评分模型要求数据里每列都是离散整数取值最好是从0开始的连续整数。如果拿到的是连续变量必须先在外部做离散化。我一般用pgmpy自带的KBinsDiscretizer做等频分箱分箱数先给3到5不要一上来就设10。等频分箱对长尾分布比较稳等距分箱碰上偏态分布很容易把多数样本塞进一个箱子里。缺失值也必须提前处理K2Score不会容忍NaN。常见做法是缺失率高的列直接删缺失率低的行先删行或按众数补数据集行数很紧张时再用插补import pandas as pd from pgmpy.preprocessing import KBinsDiscretizer data pd.read_csv(your_data.csv) data data.dropna().reset_index(dropTrue) discretizer KBinsDiscretizer(data, n_bins5, strategyquantile) discrete_data discretizer.transform(data).astype(int) for col in discrete_data.columns: print(col, sorted(discrete_data[col].unique()))这段代码先把带缺失的行丢掉再用等频分箱把连续列变成5档整数。astype(int)很关键不同版本的pgmpy分箱结果类型并不一致有的会返回浮点K2Score会直接报错。最后的循环是检查每个变量的实际取值确认它们是连续的整数集合比如0、1、2而不是0、2、4这种缺值的状态集。如果发现取值不连续说明数据里混入了异常编码先用factorize或者LabelEncoder重编码再往下走。3.2 K2结构学习的核心代码手写循环与库搜索两种姿势pgmpy的estimators模块里提供了K2Score也就是等价样本量为1的BDeu评分但并没有一个叫K2的现成搜索器。想跑严格意义上的K2常见做法是自己写贪心循环逻辑非常短。下面的函数就是标准实现from pgmpy.estimators import K2Score def k2_search(discrete_data, node_order, max_parents3): scorer K2Score(discrete_data) learned {} for i, node in enumerate(node_order): candidates [v for v in node_order[:i] if v ! node] parents [] while len(parents) max_parents and candidates: current scorer.local_score(node, parents) best_gain float(-inf) best_parent None for cand in candidates: if cand in parents: continue new_score scorer.local_score(node, parents [cand]) gain new_score - current if gain best_gain: best_gain, best_parent gain, cand if best_gain 1e-6: break parents.append(best_parent) candidates.remove(best_parent) learned[node] parents return learned核心逻辑只有两层循环外层按节点顺序逐个处理内层在当前节点的候选父集里找最大增益。local_score(node, parents)返回的是对数评分越大越好所以gain 0就说明这个父节点值得收编。1e-6的阈值是防止浮点噪声导致的多加一条无意义边实际调参时可以放宽到1e-4少一些边、结果更干净。候选只取node_order里排在前面的节点这就把K2的顺序先验落到了代码里。另一条路是用HillClimbSearch加上黑名单来模拟K2搜索空间。HillClimbSearch是全局贪心搜索允许加边、删边、反边比原版K2更激进通常结果也更好但已经不是严格意义上的K2from pgmpy.estimators import HillClimbSearch, K2Score node_order list(discrete_data.columns) black_list [ (later, earlier) for i, earlier in enumerate(node_order) for later in node_order[i 1:] ] hc HillClimbSearch(discrete_data) dag hc.estimate( scoring_methodK2Score(discrete_data), black_listblack_list, max_indegreemax_parents, epsilon1e-4, ) print(dag.edges())黑名单把所有“后序节点指向先序节点”的边全部禁掉搜索空间就和K2一致了。区别在于HillClimbSearch每走一步都会重估整张图还可能删掉之前加错的边而原版K2只加边不删边。对比文献复现用前面的手写版本实际建模图省事、想冲更高评分用HillClimbSearch更划算。我个人在项目里通常先跑手写K2拿一个稳妥的基线再用HillClimbSearch做增强两边结果一对比就能看出数据里有没有足够信号支撑更复杂的结构。提示node_order如果不是手动指定的直接用list(discrete_data.columns)会按列名排序这个顺序完全随机实战里必须用领域知识或者排序算法单独确定。3.3 参数说明max_parents、节点顺序、epsilon 怎么给K2需要设置的参数不多但每一个都直接决定结果质量。按重要程度排顺序是节点顺序、max_parents、epsilon、分箱数。参数含义和调整思路放在一起看参数作用常见设置调整思路node_order节点先后顺序K2的硬先验领域时间顺序/因果层次排错全盘皆错务必做敏感性检查max_parents每个节点最多父节点数3~5起步小样本从3开始别硬加epsilon评分增益阈值1e-4~1e-6越小搜得越狠越大图越稀疏n_bins连续变量分箱数3~5箱数少丢信息多则参数爆炸节点顺序是唯一没有后悔药的参数。如果没有领域知识常见做法是先算变量两两之间的互信息用最大生成树建一个粗略的拓扑顺序再把这个顺序喂给K2。这样不一定对但比瞎排强很多而且能跑通整个流程。顺序到底靠不靠谱有个很实用的诊断方法把顺序随机打乱10次每次跑一遍K2看学出来的边集波动有多大。如果10次结果SHD波动剧烈说明数据对顺序先验高度敏感任何单次K2输出都不能当结论必须回头去找更可靠顺序或者换PC这类不依赖顺序的算法交叉验证。max_parents的默认值我会控制在3到5。不是越多越好父节点组合数按乘积增长两个三值父节点就有9种组合再来一个就27种样本稍微小一点每个组合下的计数稀疏得像筛子评分全是噪声。epsilon在1e-4到1e-6之间调调大一点图更稀疏调小一点更容易塞进弱相关边。实际使用中epsilon1e-4在几千样本上效果通常已经足够干净。4. 模拟数据验证用已知结构反推K2的可靠性4.1 构造已知贝叶斯网络并生成样本结构学习是个没有标准答案的任务拿到真实数据你永远不知道学出来的图对不对。所以在正式跑业务数据之前先做一个合成数据验证是必须的步骤手工搭一个已知结构的贝叶斯网络从里面生成数据再把数据喂给K2看它能不能把原结构找回来。这个过程能在十几分钟内告诉你当前的数据质量、顺序和参数设定到底行不行。下面搭一个五个节点的疾病诊断模型X1和X2共同影响X3X3和X4共同影响X5。用pgmpy的TabularCPD逐个定义条件概率表from pgmpy.models import BayesianNetwork from pgmpy.factors.discrete import TabularCPD model BayesianNetwork([(X1, X3), (X2, X3), (X3, X5), (X4, X5)]) cpd_x1 TabularCPD(X1, 2, [[0.6], [0.4]]) cpd_x2 TabularCPD(X2, 2, [[0.5], [0.5]]) cpd_x3 TabularCPD( X3, 2, evidence[X1, X2], values[[0.9, 0.7, 0.6, 0.1], [0.1, 0.3, 0.4, 0.9]] ) cpd_x4 TabularCPD(X4, 2, [[0.8], [0.2]]) cpd_x5 TabularCPD( X5, 2, evidence[X3, X4], values[[0.9, 0.8, 0.5, 0.2], [0.1, 0.2, 0.5, 0.8]] ) model.add_cpds(cpd_x1, cpd_x2, cpd_x3, cpd_x4, cpd_x5) assert model.check_model() sim_data model.simulate(n_samples5000, random_state42) discrete_data sim_data.astype(int)这里有个新手常踩的坑TabularCPD的values列顺序是按evidence列表从后往前遍历的。以X3为例四列依次对应X10,X20X10,X21X11,X20X11,X21。写概率表之前先把这个顺序想清楚不然check_model()会报概率不归一化或者更隐蔽地直接构造出一个语义完全错掉的分布。simulate负责从网络里采样5000个样本对五节点二值网络来说足够学出稳定结构。4.2 学习恢复并计算结构距离沿用第3章手写的k2_search直接把合成数据丢进去学。节点顺序这里用真实生成顺序也就是X1、X2、X3、X4、X5——先有原因后有结果完全符合K2的顺序假设node_order [X1, X2, X3, X4, X5] learned k2_search(discrete_data, node_order, max_parents3) learned_edges [(p, c) for c, ps in learned.items() for p in ps] true_edges list(model.edges()) print(learned edges:, learned_edges) print(true edges:, true_edges) def compute_shd(edges_learned, edges_true): true_und {frozenset(e) for e in edges_true} learn_und {frozenset(e) for e in edges_learned} shd len(true_und.symmetric_difference(learn_und)) common true_und learn_und for pair in common: true_dir [e for e in edges_true if frozenset(e) pair][0] learn_dir [e for e in edges_learned if frozenset(e) pair][0] if true_dir ! learn_dir: shd - 1 return shd print(SHD:, compute_shd(learned_edges, true_edges))SHD结构汉明距离就是把学到的图变成真实图需要的最少操作次数加一条边、删一条边、反一条边各算1。反向边在这种实现里需要做一次删除加一次新增但我在公共边检测里把它扣了1次所以最终反向边按1次计。正确的顺序加上5000样本这个五节点例子通常能学到SHD0的结果偶尔有一两条边因为采样噪声丢失SHD落在1左右也正常。如果把这个顺序故意打乱比如把X5放到最前面SHD会明显上升而且错的边大多是那些“后被K2当作先验原因”的边。这一步对比能直观感受到顺序对结果的支配作用比任何原理说明都更入味。4.3 结果解读SHD、评分和方向一致性的判断线拿到SHD后怎么判断能不能收工我一般用三个维度交叉看。第一个是SHD本身五节点小网络SHD0到1属于完美恢复2到4说明数据噪声大或顺序有局部错位超过4就得回去查数据质量而不是怀疑K2——它只是在给定的错误约束下忠实执行搜索。第二个是总评分对比用K2Score(discrete_data).total_score(model)分别算真实图和学到的图两张图的分数应该非常接近学到图偶尔略高也不意外因为K2本就是在最大化这个分数。第三个维度是方向一致性这个最容易被忽略。K2学出的有向边不能直接理解成因果方向尤其在链式结构里X→Y和X←Y可能属于同一个马尔可夫等价类评分完全相同K2只是按黑名单约束选了一个方向。真正能从数据里识别方向的地方是V结构也就是某个节点有多个父节点且这些父节点之间没有边相连。解读结果时优先看V结构处的位置和方向链式部分只当作依赖关系用不要写进结论当成因果。模拟验证的意义就在于你先知道真图再看K2的哪些边是可靠的哪些只是等价类里的一个代表这个经验可以直接迁移到真实数据上。5. K2落地避坑节点顺序、离散化与等价类的五个翻车现场5.1 节点顺序排错整个DAG方向跟着乱现象跑出来的DAG确实没有环边也少但方向看起来就是不对劲。比如业务上明明是A导致B学出来却是B导致A而且评分还挺高。更迷惑的是如果把顺序彻底反过来再跑一遍得到的图评分几乎一样方向却整体反转。原因K2把节点顺序当成硬约束给它一个错的空间它就只能在错的空间里找最优解。A和B之间不管是A→B还是B→A只要两者没有形成V结构评分函数根本分不出哪个方向更对黑名单里没禁掉的那个方向就留下了。顺序排反等于先入为主把因果方向写错了数据再强也救不回来。解决没有后悔药只能重新跑。实战里有顺序直接排没有顺序先做互信息排序或者用PC算法先学无向骨架再定序。更务实的方式是做顺序敏感性分析随机生成10组顺序跑10次K2如果学出来的边结构大范围漂移就说明结果完全被顺序绑架任何单次输出都不可信必须回到数据探索阶段把顺序定稳再往下走。5.2 连续变量直接喂评分函数K2Score直接罢工现象DataFrame里明明有值代码报错说数值不是整数或者数据带小数时评分结果全都是NaN图学出来是空的。原因K2评分是建立在多项式计数上的底层逻辑是数样本落入哪个状态组合小数或者浮点类型根本没法计数。这是很多第一次用pgmpy的人会撞上的错误日志信息往往还比较隐晦只提示score计算失败不告诉你问题出在分箱。解决分箱是必经步骤。先用KBinsDiscretizer做等频分箱箱数控制在3到5。分箱数也不是越大越好5个箱子和2个父节点已经有25种组合样本量稍微不够就全是稀疏计数学出来的边带着大量虚高评分看起来结构丰富换个数据集立刻崩塌。分箱之后顺手astype(int)杜绝类型问题。5.3 缺失值没处理程序直接报错或静默丢行现象数据里有NaNK2Score直接抛异常换成某些旧版本库会静默把包含缺失的行整行丢弃样本量悄悄缩水学出的结构越来越稀。原因K2的边际似然计算要求每条样本在所有变量上都有完整观测缺失值在数学上直接没法计数。静默丢行比报错更危险因为没有人会发现样本量变了几千行的数据丢掉几百行结构学习结果可能完全漂移。解决在进K2之前就统一处理。缺失率超过20%的变量先删列个别缺失行按删除或众数插补处理处理完要看一眼每列的行数、取值集合。这条看似是数据预处理的基本功却是K2实践里翻车频率最高的一步因为报错时机往往不在读取时而在评分那一行。业界所谓「趁早踩坑」就是这个意思数据不进评分函数模型不会替你兜底。5.4 评分等价类学出来的方向不等于因果方向现象同一个数据集用不同的初始化顺序跑K2得到的图评分相差无几但某些边的方向完全相反看起来像是算法不稳定。原因这是评分函数的数学性质不是bug。结构学习只能识别到马尔可夫等价类这个粒度同一等价类里的所有DAG共用一个评分值。X→Y和X←Y在没有任何其他相邻结构时属于同一个等价类K2选哪个纯粹看黑名单先验。解决解读时区分两类结构。V结构如X→Z←Y且X与Y无直接边是可以从数据里识别方向的多路条件依赖会打破评分等价链式和分叉结构的方向则完全是先验说了算。交付结论时把K2学出的图称作「条件依赖图」而不是「因果图」只在V结构处谈方向倾向。如果业务上必须因果方向那就只能引入干预数据或者更严格的领域约束算法本身已经到极限了。5.5 max_parents设得太大评分表直接爆炸现象max_parents从3提到6之后运行时间暴涨几十倍内存也开始吃紧最后学出来的图反而多了好几条明显不靠谱的弱边。原因K2评分要对每一种父节点组合做计数组合数是各父节点状态数的乘积。两个三值父节点是9种组合四个三值父节点直接到81种叠加5000个样本每个组合平均只有60条数据评分对这60条数据的随机波动极其敏感弱相关边就趁虚而入。解决max_parents从3开始最多不要超过5这是绝大多数实际业务数据集的安全边界。如果确实怀疑某个节点受很多因素影响先用3跑通全流程、看一眼核心边再单独放开这个节点的上限做对比而不是全局盲调。小样本场景还有一个土办法先跑500行的子集做pipeline验证确定参数组合和代码没有性能问题再用全量数据正式跑能省下大量试错的等待时间。6. 进阶让K2结果更稳的BDeu参数与模型平均技巧前文用的K2Score本质上是等价样本量设为1的BDeu评分。这个等价样本量是可以调的把它调大相当于给参数先验加了更多虚拟样本评分对数据里的偶然波动更钝感学出来的结构更稀疏、更保守。数据量小、噪声大的场景我一般把equivalent_sample_size调到5到10能明显减少那些只在某一个数据集里成立、换个样本就消失的假边from pgmpy.estimators import BDeuScore scorer BDeuScore(discrete_data, equivalent_sample_size10)另一个比调参更能救场的习惯是模型平均。单次K2学出来的边不一定稳尤其是样本量几千这种常见规模随机采样波动就能让个别边从有到无。我一般从原始数据里做有放回抽样跑50到100次K2统计每条边出现的频率只保留出现率超过0.8的边作为最终结构。这样做会把那些评分略高但脆弱的边自动过滤掉剩下的基本都是数据里的稳定信号。这个技巧在真实业务数据上比任何提参技巧都管用代价只是多跑几十次搜索对K2这种量级的算法来说完全能接受。from sklearn.utils import resample edge_count {} n_bootstrap 50 for seed in range(n_bootstrap): boot resample(discrete_data, replaceTrue, n_sampleslen(discrete_data), random_stateseed) dag k2_search(boot, node_order, max_parents3) for child, parents in dag.items(): for parent in parents: edge_count[(parent, child)] edge_count.get((parent, child), 0) 1 stable_edges [e for e, c in edge_count.items() if c / n_bootstrap 0.8] print(stable_edges)和PC算法怎么分工也顺带说一句K2适合你已经掌握变量先后顺序的场景比如按时间发生的流水线数据、按因果层级组织的指标体系PC不依赖顺序但依赖独立性检验样本量大时才稳。两者不是替代关系真实项目中我会先用PC出一个无向骨架做交叉验证再用K2定向两边都对得上的边才敢写进结论。我有一次在本地模拟数据上把顺序排反了学出来的V结构方向整体反转要不是用模型平均看到边频率对半开差点就拿着这张图去交付了。从那之后我跑结构学习默认加两道工序先验证后解读。希望帮到你。本文还有配套的精品资源点击获取