简介K2算法是一种基于贪心搜索的贝叶斯网络结构学习方法在给定节点顺序的条件下通过逐步为每个节点选择最优父节点集合来优化网络评分。面向有一定概率图模型基础的研究生、科研人员及算法工程师该压缩包提供一套从数据中学习贝叶斯网络结构的完整可运行实现。压缩包共7个文件包含4个M函数、1个C文件、1个TXT说明和1个MAT样例数据M函数覆盖主流程、控制中心与辅助子函数逻辑层次清晰C文件通过mex接口编译后可用于加速评分计算MAT文件提供现成数据便于快速验证算法效果。整个压缩包仅10KB代码精炼无冗余已有2168人学习下载是该领域常用的学习参考。借助这套实现读者可以深入理解K2算法中节点排序、父节点贪心选择及评分函数的实现细节也可将mex接口迁移到自己的项目中作为贝叶斯网络结构学习的工具模块或教学示例。1. K2算法不是黑匣子它解决的是贝叶斯网络结构从哪来的问题拿到一张满是变量的业务表最想做的第一件事通常是搞清楚谁影响谁。贝叶斯网络能给出答案但网络结构不会自己冒出来——靠专家画图成本高、周期长主观性还强。K2算法是贝叶斯网络结构学习里最常用的一类打分搜索方法给一组离散观测数据和变量顺序它自动算出一个有向无环图告诉你变量之间依赖关系长什么样。它不依赖专家手工连线也不需要你做干预实验只要数据够、顺序大致合理几分钟就能跑出第一版结构。这篇文章带着你把K2从打分公式到代码实现过一遍数据怎么准备、节点顺序怎么定、伪计数为何不能乱调、连续变量怎么喂进去以及哪些坑会让你学出一个“看起来合理、实际上没用”的网络。新手能照着代码跑通熟手也能在这里找到调参边界和分析盲区。2. 从数据到网络K2的贝叶斯评分与贪心搜索为什么能work2.1 贝叶斯评分先算一个节点在给定父节点下的拟合度K2的核心是把“哪个结构好”变成“哪个结构后验概率高”。给定数据D结构G的后验正比于P(G)P(D|G)其中P(D|G)是在该结构下观测数据出现的边际似然。对离散变量每个节点的局部条件概率用多项分布建模参数配Dirichlet先验那么P(D|G)可以写成每个节点家族得分乘积的形式。换句话说全局结构的好坏被分解成了局部变量间的依赖强度K2算法就靠这个分解把结构搜索变成一个个节点的独立操作。常见的做法是直接比较对数得分。某个节点vi在父节点集pa(vi)下的K2家族分可以表达为score(vi, pa(vi)) ∑_{j} [ log Γ(α_ij) − log Γ(α_ij N_ij) ∑_{k} ( log Γ(α_ijk N_ijk) − log Γ(α_ijk) ) ]这里的j遍历父状态组合k遍历当前变量的取值状态N_ij与N_ijk是统计频数α_ij与α_ijk是由等价样本大小α推导出的伪计数。伪计数不是可调优的“锦上添花”它决定了先验有多强α越大越倾向于少加边结构越稀疏α越小越容易被数据中的偶然共现带偏。我一般把α当作稀疏力度来理解而不是当作什么玄学参数。这个公式不必手算因为每个变量只依赖它的潜在父节点不同变量的得分计算完全独立。这个独立性正是K2贪心搜索能跑起来的理由你不需要枚举全网结构只需要对每个变量从空集开始逐步往父节点集合里塞入选的候选节点观察家族分有没有变高。家族分布在所有可能父组合上的真实验证后就能给出“谁的加入更合理”的排序不需要依赖任何全局图约束。这是K2和成对互信息这类启发式最大的区别它评估的不是“两个变量有没有相关”而是“在给定其他父节点的条件下加入这个父节点是否改善整体拟合”。2.2 贪心搜索节点顺序和父节点数上限是两个抓手K2假设用户提供一个节点顺序所有父节点必须出现在子节点之前。这个顺序是算法成立的前提不是可选配置。有了顺序变量vi的候选父节点只能从排在vi前面的节点里选。贪心搜索从空父集开始依次尝试把每一个候选父节点加入当前父集计算得分增量每次选增量最大的那个保留直到没有正向增量或父节点数达到上限max_parents。与一般爬山搜索不同的是K2这个策略一旦接受一个父节点就不会再移除它也不会尝试翻转边的方向。这是贪心的短视之处也是它高效的原因K2只做局部最优复杂度大致为O(n^2 · r_max^2 · N)n是节点数r_max是变量最大状态数N是样本量。在几百个节点、几千条样本的规模下通常几十秒就能出一版结构比穷举和结构EM都快一个数量级。节点顺序和max_parents是K2真正需要调的两个抓手。顺序错的后果比参数错更严重如果把真正的因排在果的后面算法只能给出反边或者干脆漏掉这条边。max_parents则控制结构的复杂度业务上通常设为3或4因为高维依赖在观测数据里很难被稳定估计即使在真实系统中存在某节点有6个父节点的情况用有限样本也学不出来。在样本量不够大时宁可接受漏边也好过被虚假依赖带偏。手写实现时我会把两个参数放在最显眼的位置order和max_parents。不要把它们藏在配置文件的角落里因为它们就是K2结果质量的生死线。下一章开始写代码时你会看到这两个参数如何直接决定搜索空间的大小。3. 手写 K2 结构学习数据预处理与核心代码3.1 离散化与节点顺序写代码前的两步准备K2只能用离散数据。连续变量必须先分箱每个变量的取值要映射成从0开始的整数编码。常见做法是先用等频分箱把连续列变成类别列再用因子编码转成整数。等频分箱的好处是每一箱样本量接近不会出现某箱只有几条样本导致后续伪计数失衡。箱数推荐控制在3到8之间箱数太多会让父状态组合数指数膨胀样本很快被稀释。节点顺序的处理我一般用单独一列order数组保存它决定候选父节点的范围。业务上有明确时间先后就用时间先后比如“订单创建支付发货”没有明显时序就用领域经验里最可能的因果方向。顺序写进代码时要注释来源不然两天后回看这段代码自己都不知道顺序当初是根据什么拍的。下面是数据预处理的代码输入是一张原始DataFrame输出是离散化与编码后的numpy矩阵和节点顺序。import pandas as pd import numpy as np def prepare_discrete_data(df, order_columns, bins_by_colNone): 把原始DataFrame转成K2可用的离散整数矩阵。 df: 原始数据一列一个变量 order_columns: 按父-子排列的列名列表 bins_by_col: 每个连续列的分箱数None表示该列已经是类别列 df df.copy() # 1. 连续变量等频分箱并映射为整数 if bins_by_col is not None: for col, n_bins in bins_by_col.items(): if col in df.columns and df[col].dtype.kind in fc: # qcut做等频分箱duplicatesdrop防止连续值集中在同一箱 df[col] pd.qcut(df[col], qn_bins, labelsFalse, duplicatesdrop) # 2. 所有列统一转成整数编码 for col in order_columns: df[col] df[col].astype(category).cat.codes # 3. 保证顺序合法order_columns里出现的列必须都在df中 missing [c for c in order_columns if c not in df.columns] if missing: raise ValueError(f缺少列: {missing}) # 4. 返回整数矩阵和顺序下标 X df[order_columns].to_numpy(dtypeint) order list(range(len(order_columns))) return X, order, order_columns代码里的qcut是关键等频比等距在K2里稳得多因为K2统计的是状态组合频数样本分布不均会让某些组合整段空缺。最后返回的order是整数下标数组对应order_columns的顺序。如果后续要解释学出来的边记得保留order_columns这个名字列表输出边时要把整数节点还原成列名。这一步完成后数据已经是K2需要的形态一个二维整数矩阵和一份节点顺序。接下来就可以写家族分函数了。3.2 K2 评分函数伪计数和 log 空间一起算评分函数是K2实现的核心不要直接按概率乘法算会下溢。正确做法是把所有项移到log空间用lgamma函数处理Γ。家族分函数接收当前变量编号、候选父节点列表、全局数据矩阵和伪计数强度α返回一个浮点数。这段代码可以完整复用到你自己项目里。from math import lgamma import numpy as np def k2_family_score(X, var, parents, alpha1.0): 计算变量var在给定parents下的K2家族分。 X: (N, M)整数矩阵每列取值0..r-1 var: 当前变量列号 parents: 候选父节点列号列表必须全部排在var前面 alpha: 等价样本大小K2原文默认1.0调大更稀疏 r int(X[:, var].max()) 1 # 当前变量的状态数 if not parents: parent_states np.zeros(X.shape[0], dtypeint) q 1 else: # 把多个父节点的状态组合编码成一个整数 r_p [int(X[:, p].max()) 1 for p in parents] q int(np.prod(r_p)) parent_states X[:, parents[0]].copy() for j, p in enumerate(parents[1:], start1): parent_states parent_states * r_p[j] X[:, p] score 0.0 for j in range(q): mask_j parent_states j n_j int(mask_j.sum()) alpha_j alpha / q # 每个父状态组合分到的伪计数 # 父组合j自身的边际似然项 score lgamma(alpha_j) - lgamma(alpha_j n_j) # 当前变量在该父组合下的条件分布项 counts np.bincount(X[mask_j, var], minlengthr) for k in range(r): n_ijk int(counts[k]) # 每个取值状态再分到 alpha_j / r 的伪计数 score lgamma(alpha_j / r n_ijk) - lgamma(alpha_j / r) return scoreparents为空时q1相当于只对变量的边际分布打分这个分数是所有候选结构的比较基准。不同的父节点组合之间可以直接相减得分差值就是“加入某个父节点带来的对数边际似然提升”这个差值正是下一章贪心搜索的决策依据。alpha在这里的作用可以从参数名直接看出alpha / q是每组父状态组合的伪计数缩小r倍后进入每个取值状态的lgamma。alpha越大n_ijk在lgamma里的主导作用越弱新边被接受的门槛越高。实际跑数据时发现用float64存score没问题但n_j为0的父组合仍然会进入循环这一项的lgamma相互抵消不影响结果。bincount的minlength参数保证了状态k即使没出现也返回0不会漏项。评分函数本身没有for循环逐样本遍历瓶颈只在q的规模上。所以max_parents一旦设为5q可能冲到几千评分计算时间会明显拉长。3.3 贪心搜索主循环三个参数怎么设有了家族分函数搜索就按K2的原始策略来按给定顺序遍历每一个变量候选父节点只从顺序更靠前的集合里取每次加入使增益最大的那一个。写完这个主循环K2结构学习就完成了90%。from itertools import combinations def k2_structure_learning(X, order, max_parents3, alpha1.0, min_gain1e-4): K2贪心结构搜索主流程。 X: 离散整数数据矩阵 (N, M) order: 节点顺序列表父节点必须排在子节点前面 max_parents: 单个节点最多允许的父节点数建议3~4 alpha: 等价样本大小越大越稀疏 min_gain: 得分增量阈值用于过滤数值噪声 n_vars X.shape[1] parents {v: [] for v in range(n_vars)} for var in order: # 候选父节点顺序在var之前且未被使用 candidates [v for v in order[:order.index(var)]] while len(parents[var]) max_parents: base_score k2_family_score(X, var, parents[var], alpha) best_gain 0.0 best_cand None for cand in candidates: if cand in parents[var]: continue new_score k2_family_score(X, var, parents[var] [cand], alpha) gain new_score - base_score if gain best_gain min_gain: best_gain gain best_cand cand if best_cand is None: break parents[var].append(best_cand) return parents主循环的逻辑和推导过程一一对应顺序靠前的节点可以为后面的节点当父节点反过来不行每次迭代以当前父集为基准尝试所有未选候选记录增益最大者如果该增益不能超过min_gain就停止为这个节点加父节点。min_gain是一个容易被忽略的参数很多人只调max_parents和alpha结果在数据量很大时加进一条边际贡献极小的边纯属数值噪声。我一般把min_gain设到1e-4到1e-3之间样本量上万时可以放宽到1e-3。三个参数之间的配合是这样的max_parents控制结构搜索的“上限空间”alpha控制先验对稀疏的偏好min_gain控制数值上的接受门槛。做参数扫描时先固定max_parents3然后看alpha从0.5到2之间的结构变化如果边数变化太剧烈说明数据支持度不足。跑完以后把parents字典转成有向边三元组输出顺手把节点编号映射回列名一份可读性高的结构图就出来了。4. 用现成库快速跑通K2 打分与爬山搜索落地4.1 现成库的本质一样需要顺序与离散数据手写实现适合理解原理和定制逻辑但实际项目里更多用现成库加速迭代。常见的Python贝叶斯网络库都内置了K2评分函数和结构搜索接口设计通常是把评分器和搜索器分开你可以自己组合。用库的好处是底层有大量优化比如边缘得分的缓存、多项式展开的加速几百个节点也不至于慢到不可用。但必须记住库实现的K2同样要求数据离散化、同样依赖传入节点顺序。顺序丢失是库使用中最常见的错误——很多人只传数据框不传顺序库默认按列顺序把先出现的列当作先验父节点结果学出来的结构完全不是业务预期。下面是直接用库跑结构学习的示例。数据就是上一章的离散矩阵这里只需要还原成DataFrame并显式传入顺序。import pandas as pd from pgmpy.estimators import K2Score, HillClimbSearch # 假设 prepare_discrete_data 已经处理好了 X 和 order_columns df_discrete pd.DataFrame(X, columnsorder_columns) # K2评分器内部就是家族分的实现alpha 参数对应上一章的等价样本大小 k2 K2Score(df_discrete, equivalent_sample_size1.0) # 爬山搜索默认会尝试加边、删边、翻转边三种操作 hc HillClimbSearch(df_discrete) model hc.estimate( scoring_methodk2, max_indegree3, white_listNone ) print(model.edges())这段代码的注意点在于HillClimbSearch默认不接收“变量顺序”参数它靠评分器内部的优先关系判断方向而不是靠用户给定顺序。因此如果你手头只有乱序的数据先手动限定white_list或black_list把“顺序靠前的节点指向顺序靠后的节点”以外的边全部禁止。否则爬山搜索会在等价类网络之间反复横跳结果不稳定。库代码看起来比手写实现少但参数控制反而更隐晦我用库跑实验时一般会先固定一个简单的三变量数据把库输出和手写实现对照一遍确认两边边集一致再上真实数据。4.2 结构学习之后参数学习与结构评估K2结构学习的产物只是DAG的骨架每条边上的条件概率表还没有数值。要想让这个网络具备推理能力还需要参数学习这一步。参数学习简单得多数据完整的情况下直接用极大似然估计统计频数再归一化得到条件概率表。这一步骤和结构学习完全解耦你先决定图再填参数。用现成库做参数估计的代码如下。from pgmpy.models import BayesianNetwork from pgmpy.estimators import MaximumLikelihoodEstimator # 从学到的边构造贝叶斯网络对象 model BayesianNetwork(list(model.edges())) # 用同一份离散数据做MLE参数估计 mle MaximumLikelihoodEstimator(model, df_discrete) cpd_map mle.get_parameters() for cpd in cpd_map: print(cpd)参数估计完成后模型就可以回答条件概率查询了。但结构评估不能只靠看一眼边集下结论。我每次拿到K2学出的网络至少会做两个检查一是对数似然或BIC在不同数据划分上是否稳定二是把边集和业务先验核对一遍。分数对比用库接口最快拿同一份数据分别用K2评分和BIC评分跑同一搜索器看两者给出的网络结构差异。若K2学出的边在BIC视角下明显过多基本可以判定是样本量不足导致的过拟合。除此之外还可以计算汉明距离或者精确率召回率前提是你有一个参考网络或者业务专家确认过的边集。别指望评分值本身告诉你“结构对不对”边际似然只衡量数据拟合度不度量因果正确性。一个容易忽略的事实是K2和BIC在有限样本下会收敛到不同的稀疏度K2默认α1偏向于加入更多边而BIC的惩罚项更重。所以在真实项目里我更习惯把K2的结果当作候选集再用BIC排序或专家删选收窄到最终结构。5. K2 算法实战避坑5 条踩过的记录5.1 节点顺序给错学出来的边全是反向现象数据真实关系是A → B但顺序里把B放在A前面。跑完K2后A → B的边完全消失B → A反而被加了进去有时还会连带多出一条虚假的中间节点边。整个网络结构看起来自洽但对业务解释完全没用。原因K2的搜索空间被顺序锁定候选父节点只能来自顺序靠前的一方。B排在A前A就不可能是B的子节点A的所有依赖只能从B那边漏过来。评分函数本身无法识别错误顺序它只会在受限的DAG空间里找拟合最好的那个。解决顺序必须来自业务因果而不是数据相关性。没有先验时先用无向依赖分析或互信息排序粗筛一遍再结合业务时间线定顺序。如果实在没有头绪保守做法是分别跑K2和BIC爬山把两组结果的公共边保留方向相反的边全部标记为“不可判定”不要强行采信某一条。5.2 连续变量直接塞进 K2评分震荡得很厉害现象把没离散化的数值型变量直接传给评分函数程序不报错但结果完全乱来。同一条数据跑两遍边集不同换一个变量排序方式结构天翻地覆。评分值有时为正有时为负看不出规律。原因K2内部把所有变量当作有限状态多项分布。连续值几乎每条样本都是一个独立状态父状态组合数量爆炸绝大多数组合只有1条甚至0条样本伪计数被稀释到失去平滑作用。边际似然在这种稀疏计数下变得极其敏感微小扰动都会改变得分排序。解决做分箱箱数控制在3~8之间。优先等频分箱而不是等距分箱等距分箱在长尾分布下会出现空箱等频则保证每个取值状态的统计量至少是几十条。分箱完成后用pd.cut cat.codes再检查每一列的取值个数如果某列只剩1个取值说明这个变量在数据里几乎没变异建议直接删掉否则只会给模型增添噪声。5.3 max_parents 设太大学到一张密集“关系网”现象max_parents6或更大时学出的网络每个节点都有5、6条入边几乎像一个完全图。这样的网络看起来很“丰富”但条件概率表大得没法看而且交叉验证对样本外的预测能力并不比稀疏网络好。原因K2打分有一个特性增加父节点几乎总会小幅提升训练数据的边际似然因为模型复杂度变高了。只要max_parents允许搜索就会持续往里塞父节点直到阈值拦不住为止。打分函数虽然带先验惩罚但惩罚力度在大量样本面前不足以抵消复杂度收益。解决先把max_parents固定为3跑一版作为基线再把max_parents增加到4比较两者在验证集上的BIC。如果BIC没有显著提升就退回3。另外把min_gain设高一点例如1e-3能过滤大量“只提升0.0001”的垃圾边。项目里稳定的做法是记录不同max_parents下的边数变化曲线挑曲线变平缓的位置作为取值点。5.4 零计数把小样本组合变成黑洞现象学出的结构里出现一条完全说不通的边比如“支付方式影响用户性别”。检查数据后发现这条边对应的父状态组合里样本量极少有的组合只有两三条记录但评分却给了很大的正增益。原因K2家族分里伪计数α_j / r是均匀分到每个状态组合的。组合本身样本少时观察频数N_ijk极低lgamma项的变化率反而最大。也就是说小样本组合对评分的边际贡献被放大了模型倾向于用一条边来“解释”这些零星的共现。解决调大alpha。alpha从1调到5到10伪计数增强小样本组合的过拟合收益会被抑制。同时检查每个父状态组合的样本覆盖度如果大量组合n_j 10考虑减少分箱数或者直接删掉该变量。我一般在跑完结构后写一段检查代码统计每条边对应的最小父组合样本量低于20的边全部打上“需复核”标签。5.5 把 K2 学出的有向边当成因果结论现象K2学出A → B业务人员直接拿去做决策认为干预A就能改变B。结果上线实验后B纹丝不动。这不是算法错了是使用方式错了。原因K2学习的是数据中的条件依赖关系不是因果机制。在同一个马尔可夫等价类中A → B和B → A可能有完全相同的边际似然搜索器只是按顺序随机选了一个方向。只有在存在v结构两个节点的共同子节点时方向才是可由数据识别的其他边的方向本质上是先验顺序或搜索偏好的产物。解决把K2的输出理解为“变量之间的依赖骨架”方向只表示概率分解方式不表示可干预的因果关系。要获得因果意义上的方向必须结合实验信息、时间顺序或者领域知识做人工约束。做因果推断时正确的打开方式是先用K2生成候选结构再标出哪些边方向不可识别最后用专家知识或A/B实验补全。这个边界想清楚K2就能避开“拿相关性当因果”的大坑变成高效建模工具。6. 把 K2 用得更稳先验顺序、连续数据与结构回放验证6.1 先用业务顺序当先验再用数据修正项目起步阶段没有确定顺序可以先把时间属性和流程顺序拼一份粗糙顺序创建时间早的变量排前面链路下游排后面。跑完K2后不要急着接受全部边把每条边和业务顺序核对凡是违反先验顺序的边先标出来看数据支撑度如何。常见做法是跑三组实验业务顺序、互信息排序顺序、倒序三组结果里方向一致且与业务顺序不冲突的边可靠性最高。这份“方向一致的公共边”列表比任何单一K2输出都有说服力。6.2 连续变量的离散化跟着分类任务走离散化直接决定结构学习的质量上限。多数人会陷入“箱数越多信息保留越多”的直觉误区。K2的状态组合数会随着箱数指数放大箱数超过8后样本覆盖度明显下降结构学习退化成噪声拟合。我用得比较稳的方式是跟着预测任务走如果这个连续变量后续要参与分类或回归就按目标变量做有监督分箱如果没有明确目标变量就做等频分箱然后人工检查分箱边界是否对应业务含义。不要用聚类后标签当离散取值聚类边界对结构学习基本没有帮助反而引入不必要的额外变量。6.3 结构回放用已知网络校验 K2 的调参把K2用在真实数据前花半小时做一个“结构回放”能省掉后面几天的返工从一个已知的小网络结构出发按该网络采样离散数据再用K2学习比较学到的边和真实边的差异。这是我的固定习惯每次拿到新数据先造一个小规模仿真数据把max_parents、alpha、箱数扫一遍观察结构恢复率变化选恢复率最高的参数组合再套到真实数据上。这样做的好处是避免了真实数据里“永远不知道正确答案”的尴尬。我曾经在一份样本量三千的数据上只调max_parents不调alpha结果结构在两条边之间反复横跳怎么都稳定不下来。后来做了回放实验才发现alpha3、max_parents3的组合能把恢复率稳定的边集找出来而alpha1时总是多带一条假边。从那以后每次换数据集第一件事就是拿已知结构做基准回放把参数确定下来再上正式模型希望帮到你。本文还有配套的精品资源点击获取 贝叶斯网络结构学习：K2算法原理、代码实现与避坑指南
简介K2算法是一种基于贪心搜索的贝叶斯网络结构学习方法在给定节点顺序的条件下通过逐步为每个节点选择最优父节点集合来优化网络评分。面向有一定概率图模型基础的研究生、科研人员及算法工程师该压缩包提供一套从数据中学习贝叶斯网络结构的完整可运行实现。压缩包共7个文件包含4个M函数、1个C文件、1个TXT说明和1个MAT样例数据M函数覆盖主流程、控制中心与辅助子函数逻辑层次清晰C文件通过mex接口编译后可用于加速评分计算MAT文件提供现成数据便于快速验证算法效果。整个压缩包仅10KB代码精炼无冗余已有2168人学习下载是该领域常用的学习参考。借助这套实现读者可以深入理解K2算法中节点排序、父节点贪心选择及评分函数的实现细节也可将mex接口迁移到自己的项目中作为贝叶斯网络结构学习的工具模块或教学示例。1. K2算法不是黑匣子它解决的是贝叶斯网络结构从哪来的问题拿到一张满是变量的业务表最想做的第一件事通常是搞清楚谁影响谁。贝叶斯网络能给出答案但网络结构不会自己冒出来——靠专家画图成本高、周期长主观性还强。K2算法是贝叶斯网络结构学习里最常用的一类打分搜索方法给一组离散观测数据和变量顺序它自动算出一个有向无环图告诉你变量之间依赖关系长什么样。它不依赖专家手工连线也不需要你做干预实验只要数据够、顺序大致合理几分钟就能跑出第一版结构。这篇文章带着你把K2从打分公式到代码实现过一遍数据怎么准备、节点顺序怎么定、伪计数为何不能乱调、连续变量怎么喂进去以及哪些坑会让你学出一个“看起来合理、实际上没用”的网络。新手能照着代码跑通熟手也能在这里找到调参边界和分析盲区。2. 从数据到网络K2的贝叶斯评分与贪心搜索为什么能work2.1 贝叶斯评分先算一个节点在给定父节点下的拟合度K2的核心是把“哪个结构好”变成“哪个结构后验概率高”。给定数据D结构G的后验正比于P(G)P(D|G)其中P(D|G)是在该结构下观测数据出现的边际似然。对离散变量每个节点的局部条件概率用多项分布建模参数配Dirichlet先验那么P(D|G)可以写成每个节点家族得分乘积的形式。换句话说全局结构的好坏被分解成了局部变量间的依赖强度K2算法就靠这个分解把结构搜索变成一个个节点的独立操作。常见的做法是直接比较对数得分。某个节点vi在父节点集pa(vi)下的K2家族分可以表达为score(vi, pa(vi)) ∑_{j} [ log Γ(α_ij) − log Γ(α_ij N_ij) ∑_{k} ( log Γ(α_ijk N_ijk) − log Γ(α_ijk) ) ]这里的j遍历父状态组合k遍历当前变量的取值状态N_ij与N_ijk是统计频数α_ij与α_ijk是由等价样本大小α推导出的伪计数。伪计数不是可调优的“锦上添花”它决定了先验有多强α越大越倾向于少加边结构越稀疏α越小越容易被数据中的偶然共现带偏。我一般把α当作稀疏力度来理解而不是当作什么玄学参数。这个公式不必手算因为每个变量只依赖它的潜在父节点不同变量的得分计算完全独立。这个独立性正是K2贪心搜索能跑起来的理由你不需要枚举全网结构只需要对每个变量从空集开始逐步往父节点集合里塞入选的候选节点观察家族分有没有变高。家族分布在所有可能父组合上的真实验证后就能给出“谁的加入更合理”的排序不需要依赖任何全局图约束。这是K2和成对互信息这类启发式最大的区别它评估的不是“两个变量有没有相关”而是“在给定其他父节点的条件下加入这个父节点是否改善整体拟合”。2.2 贪心搜索节点顺序和父节点数上限是两个抓手K2假设用户提供一个节点顺序所有父节点必须出现在子节点之前。这个顺序是算法成立的前提不是可选配置。有了顺序变量vi的候选父节点只能从排在vi前面的节点里选。贪心搜索从空父集开始依次尝试把每一个候选父节点加入当前父集计算得分增量每次选增量最大的那个保留直到没有正向增量或父节点数达到上限max_parents。与一般爬山搜索不同的是K2这个策略一旦接受一个父节点就不会再移除它也不会尝试翻转边的方向。这是贪心的短视之处也是它高效的原因K2只做局部最优复杂度大致为O(n^2 · r_max^2 · N)n是节点数r_max是变量最大状态数N是样本量。在几百个节点、几千条样本的规模下通常几十秒就能出一版结构比穷举和结构EM都快一个数量级。节点顺序和max_parents是K2真正需要调的两个抓手。顺序错的后果比参数错更严重如果把真正的因排在果的后面算法只能给出反边或者干脆漏掉这条边。max_parents则控制结构的复杂度业务上通常设为3或4因为高维依赖在观测数据里很难被稳定估计即使在真实系统中存在某节点有6个父节点的情况用有限样本也学不出来。在样本量不够大时宁可接受漏边也好过被虚假依赖带偏。手写实现时我会把两个参数放在最显眼的位置order和max_parents。不要把它们藏在配置文件的角落里因为它们就是K2结果质量的生死线。下一章开始写代码时你会看到这两个参数如何直接决定搜索空间的大小。3. 手写 K2 结构学习数据预处理与核心代码3.1 离散化与节点顺序写代码前的两步准备K2只能用离散数据。连续变量必须先分箱每个变量的取值要映射成从0开始的整数编码。常见做法是先用等频分箱把连续列变成类别列再用因子编码转成整数。等频分箱的好处是每一箱样本量接近不会出现某箱只有几条样本导致后续伪计数失衡。箱数推荐控制在3到8之间箱数太多会让父状态组合数指数膨胀样本很快被稀释。节点顺序的处理我一般用单独一列order数组保存它决定候选父节点的范围。业务上有明确时间先后就用时间先后比如“订单创建支付发货”没有明显时序就用领域经验里最可能的因果方向。顺序写进代码时要注释来源不然两天后回看这段代码自己都不知道顺序当初是根据什么拍的。下面是数据预处理的代码输入是一张原始DataFrame输出是离散化与编码后的numpy矩阵和节点顺序。import pandas as pd import numpy as np def prepare_discrete_data(df, order_columns, bins_by_colNone): 把原始DataFrame转成K2可用的离散整数矩阵。 df: 原始数据一列一个变量 order_columns: 按父-子排列的列名列表 bins_by_col: 每个连续列的分箱数None表示该列已经是类别列 df df.copy() # 1. 连续变量等频分箱并映射为整数 if bins_by_col is not None: for col, n_bins in bins_by_col.items(): if col in df.columns and df[col].dtype.kind in fc: # qcut做等频分箱duplicatesdrop防止连续值集中在同一箱 df[col] pd.qcut(df[col], qn_bins, labelsFalse, duplicatesdrop) # 2. 所有列统一转成整数编码 for col in order_columns: df[col] df[col].astype(category).cat.codes # 3. 保证顺序合法order_columns里出现的列必须都在df中 missing [c for c in order_columns if c not in df.columns] if missing: raise ValueError(f缺少列: {missing}) # 4. 返回整数矩阵和顺序下标 X df[order_columns].to_numpy(dtypeint) order list(range(len(order_columns))) return X, order, order_columns代码里的qcut是关键等频比等距在K2里稳得多因为K2统计的是状态组合频数样本分布不均会让某些组合整段空缺。最后返回的order是整数下标数组对应order_columns的顺序。如果后续要解释学出来的边记得保留order_columns这个名字列表输出边时要把整数节点还原成列名。这一步完成后数据已经是K2需要的形态一个二维整数矩阵和一份节点顺序。接下来就可以写家族分函数了。3.2 K2 评分函数伪计数和 log 空间一起算评分函数是K2实现的核心不要直接按概率乘法算会下溢。正确做法是把所有项移到log空间用lgamma函数处理Γ。家族分函数接收当前变量编号、候选父节点列表、全局数据矩阵和伪计数强度α返回一个浮点数。这段代码可以完整复用到你自己项目里。from math import lgamma import numpy as np def k2_family_score(X, var, parents, alpha1.0): 计算变量var在给定parents下的K2家族分。 X: (N, M)整数矩阵每列取值0..r-1 var: 当前变量列号 parents: 候选父节点列号列表必须全部排在var前面 alpha: 等价样本大小K2原文默认1.0调大更稀疏 r int(X[:, var].max()) 1 # 当前变量的状态数 if not parents: parent_states np.zeros(X.shape[0], dtypeint) q 1 else: # 把多个父节点的状态组合编码成一个整数 r_p [int(X[:, p].max()) 1 for p in parents] q int(np.prod(r_p)) parent_states X[:, parents[0]].copy() for j, p in enumerate(parents[1:], start1): parent_states parent_states * r_p[j] X[:, p] score 0.0 for j in range(q): mask_j parent_states j n_j int(mask_j.sum()) alpha_j alpha / q # 每个父状态组合分到的伪计数 # 父组合j自身的边际似然项 score lgamma(alpha_j) - lgamma(alpha_j n_j) # 当前变量在该父组合下的条件分布项 counts np.bincount(X[mask_j, var], minlengthr) for k in range(r): n_ijk int(counts[k]) # 每个取值状态再分到 alpha_j / r 的伪计数 score lgamma(alpha_j / r n_ijk) - lgamma(alpha_j / r) return scoreparents为空时q1相当于只对变量的边际分布打分这个分数是所有候选结构的比较基准。不同的父节点组合之间可以直接相减得分差值就是“加入某个父节点带来的对数边际似然提升”这个差值正是下一章贪心搜索的决策依据。alpha在这里的作用可以从参数名直接看出alpha / q是每组父状态组合的伪计数缩小r倍后进入每个取值状态的lgamma。alpha越大n_ijk在lgamma里的主导作用越弱新边被接受的门槛越高。实际跑数据时发现用float64存score没问题但n_j为0的父组合仍然会进入循环这一项的lgamma相互抵消不影响结果。bincount的minlength参数保证了状态k即使没出现也返回0不会漏项。评分函数本身没有for循环逐样本遍历瓶颈只在q的规模上。所以max_parents一旦设为5q可能冲到几千评分计算时间会明显拉长。3.3 贪心搜索主循环三个参数怎么设有了家族分函数搜索就按K2的原始策略来按给定顺序遍历每一个变量候选父节点只从顺序更靠前的集合里取每次加入使增益最大的那一个。写完这个主循环K2结构学习就完成了90%。from itertools import combinations def k2_structure_learning(X, order, max_parents3, alpha1.0, min_gain1e-4): K2贪心结构搜索主流程。 X: 离散整数数据矩阵 (N, M) order: 节点顺序列表父节点必须排在子节点前面 max_parents: 单个节点最多允许的父节点数建议3~4 alpha: 等价样本大小越大越稀疏 min_gain: 得分增量阈值用于过滤数值噪声 n_vars X.shape[1] parents {v: [] for v in range(n_vars)} for var in order: # 候选父节点顺序在var之前且未被使用 candidates [v for v in order[:order.index(var)]] while len(parents[var]) max_parents: base_score k2_family_score(X, var, parents[var], alpha) best_gain 0.0 best_cand None for cand in candidates: if cand in parents[var]: continue new_score k2_family_score(X, var, parents[var] [cand], alpha) gain new_score - base_score if gain best_gain min_gain: best_gain gain best_cand cand if best_cand is None: break parents[var].append(best_cand) return parents主循环的逻辑和推导过程一一对应顺序靠前的节点可以为后面的节点当父节点反过来不行每次迭代以当前父集为基准尝试所有未选候选记录增益最大者如果该增益不能超过min_gain就停止为这个节点加父节点。min_gain是一个容易被忽略的参数很多人只调max_parents和alpha结果在数据量很大时加进一条边际贡献极小的边纯属数值噪声。我一般把min_gain设到1e-4到1e-3之间样本量上万时可以放宽到1e-3。三个参数之间的配合是这样的max_parents控制结构搜索的“上限空间”alpha控制先验对稀疏的偏好min_gain控制数值上的接受门槛。做参数扫描时先固定max_parents3然后看alpha从0.5到2之间的结构变化如果边数变化太剧烈说明数据支持度不足。跑完以后把parents字典转成有向边三元组输出顺手把节点编号映射回列名一份可读性高的结构图就出来了。4. 用现成库快速跑通K2 打分与爬山搜索落地4.1 现成库的本质一样需要顺序与离散数据手写实现适合理解原理和定制逻辑但实际项目里更多用现成库加速迭代。常见的Python贝叶斯网络库都内置了K2评分函数和结构搜索接口设计通常是把评分器和搜索器分开你可以自己组合。用库的好处是底层有大量优化比如边缘得分的缓存、多项式展开的加速几百个节点也不至于慢到不可用。但必须记住库实现的K2同样要求数据离散化、同样依赖传入节点顺序。顺序丢失是库使用中最常见的错误——很多人只传数据框不传顺序库默认按列顺序把先出现的列当作先验父节点结果学出来的结构完全不是业务预期。下面是直接用库跑结构学习的示例。数据就是上一章的离散矩阵这里只需要还原成DataFrame并显式传入顺序。import pandas as pd from pgmpy.estimators import K2Score, HillClimbSearch # 假设 prepare_discrete_data 已经处理好了 X 和 order_columns df_discrete pd.DataFrame(X, columnsorder_columns) # K2评分器内部就是家族分的实现alpha 参数对应上一章的等价样本大小 k2 K2Score(df_discrete, equivalent_sample_size1.0) # 爬山搜索默认会尝试加边、删边、翻转边三种操作 hc HillClimbSearch(df_discrete) model hc.estimate( scoring_methodk2, max_indegree3, white_listNone ) print(model.edges())这段代码的注意点在于HillClimbSearch默认不接收“变量顺序”参数它靠评分器内部的优先关系判断方向而不是靠用户给定顺序。因此如果你手头只有乱序的数据先手动限定white_list或black_list把“顺序靠前的节点指向顺序靠后的节点”以外的边全部禁止。否则爬山搜索会在等价类网络之间反复横跳结果不稳定。库代码看起来比手写实现少但参数控制反而更隐晦我用库跑实验时一般会先固定一个简单的三变量数据把库输出和手写实现对照一遍确认两边边集一致再上真实数据。4.2 结构学习之后参数学习与结构评估K2结构学习的产物只是DAG的骨架每条边上的条件概率表还没有数值。要想让这个网络具备推理能力还需要参数学习这一步。参数学习简单得多数据完整的情况下直接用极大似然估计统计频数再归一化得到条件概率表。这一步骤和结构学习完全解耦你先决定图再填参数。用现成库做参数估计的代码如下。from pgmpy.models import BayesianNetwork from pgmpy.estimators import MaximumLikelihoodEstimator # 从学到的边构造贝叶斯网络对象 model BayesianNetwork(list(model.edges())) # 用同一份离散数据做MLE参数估计 mle MaximumLikelihoodEstimator(model, df_discrete) cpd_map mle.get_parameters() for cpd in cpd_map: print(cpd)参数估计完成后模型就可以回答条件概率查询了。但结构评估不能只靠看一眼边集下结论。我每次拿到K2学出的网络至少会做两个检查一是对数似然或BIC在不同数据划分上是否稳定二是把边集和业务先验核对一遍。分数对比用库接口最快拿同一份数据分别用K2评分和BIC评分跑同一搜索器看两者给出的网络结构差异。若K2学出的边在BIC视角下明显过多基本可以判定是样本量不足导致的过拟合。除此之外还可以计算汉明距离或者精确率召回率前提是你有一个参考网络或者业务专家确认过的边集。别指望评分值本身告诉你“结构对不对”边际似然只衡量数据拟合度不度量因果正确性。一个容易忽略的事实是K2和BIC在有限样本下会收敛到不同的稀疏度K2默认α1偏向于加入更多边而BIC的惩罚项更重。所以在真实项目里我更习惯把K2的结果当作候选集再用BIC排序或专家删选收窄到最终结构。5. K2 算法实战避坑5 条踩过的记录5.1 节点顺序给错学出来的边全是反向现象数据真实关系是A → B但顺序里把B放在A前面。跑完K2后A → B的边完全消失B → A反而被加了进去有时还会连带多出一条虚假的中间节点边。整个网络结构看起来自洽但对业务解释完全没用。原因K2的搜索空间被顺序锁定候选父节点只能来自顺序靠前的一方。B排在A前A就不可能是B的子节点A的所有依赖只能从B那边漏过来。评分函数本身无法识别错误顺序它只会在受限的DAG空间里找拟合最好的那个。解决顺序必须来自业务因果而不是数据相关性。没有先验时先用无向依赖分析或互信息排序粗筛一遍再结合业务时间线定顺序。如果实在没有头绪保守做法是分别跑K2和BIC爬山把两组结果的公共边保留方向相反的边全部标记为“不可判定”不要强行采信某一条。5.2 连续变量直接塞进 K2评分震荡得很厉害现象把没离散化的数值型变量直接传给评分函数程序不报错但结果完全乱来。同一条数据跑两遍边集不同换一个变量排序方式结构天翻地覆。评分值有时为正有时为负看不出规律。原因K2内部把所有变量当作有限状态多项分布。连续值几乎每条样本都是一个独立状态父状态组合数量爆炸绝大多数组合只有1条甚至0条样本伪计数被稀释到失去平滑作用。边际似然在这种稀疏计数下变得极其敏感微小扰动都会改变得分排序。解决做分箱箱数控制在3~8之间。优先等频分箱而不是等距分箱等距分箱在长尾分布下会出现空箱等频则保证每个取值状态的统计量至少是几十条。分箱完成后用pd.cut cat.codes再检查每一列的取值个数如果某列只剩1个取值说明这个变量在数据里几乎没变异建议直接删掉否则只会给模型增添噪声。5.3 max_parents 设太大学到一张密集“关系网”现象max_parents6或更大时学出的网络每个节点都有5、6条入边几乎像一个完全图。这样的网络看起来很“丰富”但条件概率表大得没法看而且交叉验证对样本外的预测能力并不比稀疏网络好。原因K2打分有一个特性增加父节点几乎总会小幅提升训练数据的边际似然因为模型复杂度变高了。只要max_parents允许搜索就会持续往里塞父节点直到阈值拦不住为止。打分函数虽然带先验惩罚但惩罚力度在大量样本面前不足以抵消复杂度收益。解决先把max_parents固定为3跑一版作为基线再把max_parents增加到4比较两者在验证集上的BIC。如果BIC没有显著提升就退回3。另外把min_gain设高一点例如1e-3能过滤大量“只提升0.0001”的垃圾边。项目里稳定的做法是记录不同max_parents下的边数变化曲线挑曲线变平缓的位置作为取值点。5.4 零计数把小样本组合变成黑洞现象学出的结构里出现一条完全说不通的边比如“支付方式影响用户性别”。检查数据后发现这条边对应的父状态组合里样本量极少有的组合只有两三条记录但评分却给了很大的正增益。原因K2家族分里伪计数α_j / r是均匀分到每个状态组合的。组合本身样本少时观察频数N_ijk极低lgamma项的变化率反而最大。也就是说小样本组合对评分的边际贡献被放大了模型倾向于用一条边来“解释”这些零星的共现。解决调大alpha。alpha从1调到5到10伪计数增强小样本组合的过拟合收益会被抑制。同时检查每个父状态组合的样本覆盖度如果大量组合n_j 10考虑减少分箱数或者直接删掉该变量。我一般在跑完结构后写一段检查代码统计每条边对应的最小父组合样本量低于20的边全部打上“需复核”标签。5.5 把 K2 学出的有向边当成因果结论现象K2学出A → B业务人员直接拿去做决策认为干预A就能改变B。结果上线实验后B纹丝不动。这不是算法错了是使用方式错了。原因K2学习的是数据中的条件依赖关系不是因果机制。在同一个马尔可夫等价类中A → B和B → A可能有完全相同的边际似然搜索器只是按顺序随机选了一个方向。只有在存在v结构两个节点的共同子节点时方向才是可由数据识别的其他边的方向本质上是先验顺序或搜索偏好的产物。解决把K2的输出理解为“变量之间的依赖骨架”方向只表示概率分解方式不表示可干预的因果关系。要获得因果意义上的方向必须结合实验信息、时间顺序或者领域知识做人工约束。做因果推断时正确的打开方式是先用K2生成候选结构再标出哪些边方向不可识别最后用专家知识或A/B实验补全。这个边界想清楚K2就能避开“拿相关性当因果”的大坑变成高效建模工具。6. 把 K2 用得更稳先验顺序、连续数据与结构回放验证6.1 先用业务顺序当先验再用数据修正项目起步阶段没有确定顺序可以先把时间属性和流程顺序拼一份粗糙顺序创建时间早的变量排前面链路下游排后面。跑完K2后不要急着接受全部边把每条边和业务顺序核对凡是违反先验顺序的边先标出来看数据支撑度如何。常见做法是跑三组实验业务顺序、互信息排序顺序、倒序三组结果里方向一致且与业务顺序不冲突的边可靠性最高。这份“方向一致的公共边”列表比任何单一K2输出都有说服力。6.2 连续变量的离散化跟着分类任务走离散化直接决定结构学习的质量上限。多数人会陷入“箱数越多信息保留越多”的直觉误区。K2的状态组合数会随着箱数指数放大箱数超过8后样本覆盖度明显下降结构学习退化成噪声拟合。我用得比较稳的方式是跟着预测任务走如果这个连续变量后续要参与分类或回归就按目标变量做有监督分箱如果没有明确目标变量就做等频分箱然后人工检查分箱边界是否对应业务含义。不要用聚类后标签当离散取值聚类边界对结构学习基本没有帮助反而引入不必要的额外变量。6.3 结构回放用已知网络校验 K2 的调参把K2用在真实数据前花半小时做一个“结构回放”能省掉后面几天的返工从一个已知的小网络结构出发按该网络采样离散数据再用K2学习比较学到的边和真实边的差异。这是我的固定习惯每次拿到新数据先造一个小规模仿真数据把max_parents、alpha、箱数扫一遍观察结构恢复率变化选恢复率最高的参数组合再套到真实数据上。这样做的好处是避免了真实数据里“永远不知道正确答案”的尴尬。我曾经在一份样本量三千的数据上只调max_parents不调alpha结果结构在两条边之间反复横跳怎么都稳定不下来。后来做了回放实验才发现alpha3、max_parents3的组合能把恢复率稳定的边集找出来而alpha1时总是多带一条假边。从那以后每次换数据集第一件事就是拿已知结构做基准回放把参数确定下来再上正式模型希望帮到你。本文还有配套的精品资源点击获取