1、Matlab 贝叶斯网络建模1 FullBNT 简介基于 Matlab 的贝叶斯网络工具箱 BNT 是 kevin p.murphy 基于matlab 语言开发的关于贝叶斯网络学习的开源软件包，提供了许多贝叶斯网络学习的底层基础函数库，支持多种类型的节点（概率分布） 、精确推理和近似推理、参数学习及结构学习、静态模型和动态模型。1.1 贝叶斯网络表示BNT 中使用矩阵方式表示贝叶斯网络，即若节点 i 到 j 有一条弧，则对应矩阵中 值为 1，否则为 0。(， )1.2 结构学习算法函数BNT 中提供了较为丰富的结构学习函数，都有：1. 学习树扩展贝叶斯网络结构的 算法 ._()2. 数据完整条
2、件下学习一般贝叶斯网络结构学习算法表 1-1 数据完整条件下贝叶斯结构算法算法名称 调用函数K2 算法 learn_struct_k2()贪婪搜索 GS（greedy search ）算法 earn_struct_gs()爬山 HC（hill climbing）算法 learn_struct_hc() 3. 缺失数据条件下学习一般贝叶斯网络结构学习算法表 1-2 缺失数据条件下贝叶斯结构算法算法名称 调用函数最大期望 EM（expectation maximization）算法 learn_struct_EM（）MCMC（Markov Chain Monte Carlo） learn_stru
3、ct_mcmc() 1.3 参数学习算法函数1. BNT 中也提供了丰富的参数学习函数，都有：2. 完整数据时，学习参数的方法主要有两种：最大似然估计learn_params()和贝叶斯方法 bayes_update_params();3. 数据缺失时，如果已知网络拓扑结构，用 EM 算法来计算参数， learn_params_em ()。1.4 推理机制及推理引擎为了提高运算速度，使各种推理算法能够有效应用，BNT 工具箱采用了引擎机制，不同的引擎根据不同的算法来完成模型转换、细化和求解。这个推理过程如下：BNT 中提供了多种推理引擎，都有：表 1-3 BNT 推理引擎算法名称 调用函数联合
4、树推理引擎 jtree_inf_engine()全局联合树推理引擎 global_joint_inf_engine()信念传播推理引擎 belprop_inf_engine()变量消元推理引擎 var_elim_inf_engine()采样传播引擎 gibbs_sampling_inf_engine2 参数学习在 BNT 中，参数评估程序可以分为 4 类。分类依据是否目标是通过参数或仅仅一个点的估计来计算贝叶斯全部的后验概率，是否全部的节点是可观察的或者存在数据/隐含变量（局部可观察） 。Full obs Partial obsPoint learn_params learn_params_e
5、mBayes Bayes_update_params not yet supported2.1 节点类型Noisy-or 节点一个 Noisy-or 节点就像通常的“或”门一样，但有时父节点的效果将被抑制。受抑制的父节点 i 的概率用 来表示。一个节P（ i）点 C，有两个父节点 A 和 B，有如下 CPD，使用 F 和 T 来表达关和开，（在 BNT 中是 1 和 2） 。A B P(C=off) P(C=on)F F 1.0 0.0T F qA) 1-q(A)F T q(B) 1-q(B)T T q(A)q(B) 1-q(A)q(B)Softmax 节点神经网络节点 使用一个多层感知器实现
6、了从连续父节点向离散子节点的映射。高斯节点 将连续值的节点处理成一个离散的情况广义线性模型节点分类/回归树节点2.2 最大似然参数估计bnet3 = learn_params(bnet2, samples);2.3 先验参数分布tabular_CPD(bnet, i, prior_type, dirichlet, dirichlet_type, unif);B=1 B=2 B=3A=1 1 1 1A=2 1 1 1tabular_CPD(bnet, i, prior_type, dirichlet, dirichlet_type, BDeu);B=1 B=2 B=3A=1 1/6 1/6 1/
7、6A=2 1/6 1/6 1/6我们将 N/(q*r) 放入每个格； N 是等效的样本大小，r=|A|, q = |B|. 这可以按如上面方式创建：tabular_CPD(bnet, i, prior_type, dirichlet, dirichlet_type, .BDeu, dirichlet_weight, 10);这里 1 是等效样本大小，也是先验概率的强度。你可以使用上面面方式更改它，3 结构学习问题：以下两模型结构评分是否相等？表 3-1 算法概要贝叶斯模型选择算法1.建立模型 A-B，生成样本数据2.建立所有可能的结构: (1)A B, (2)BB 并计算先验概率3.模型 2
8、和模型 3 为 Markov equivalent4.B 节点使用 noisy Not gate5.正确的模型在 12 次后收敛代码示例% 贝叶斯选择模型示例.% 建立模型A-B，生成样本数据% 建立所有可能的结构: (1)A B, (2)BB 并计算先验概率% 模型2和模型3为Markov equivalent% B节点使用noisy Not gate% 正确的模型在12 次后收敛% ground truthN = 2;dag = zeros(N);A = 1; B = 2; dag(A,B) = 1;ntrials = 25;ns = 2*ones(1,N);true_bnet = mk_
9、bnet(dag, ns);true_bnet.CPD1 = tabular_CPD(true_bnet, 1, 0.5 0.5);pfail = 0.1;psucc = 1-pfail;true_bnet.CPD2 = tabular_CPD(true_bnet, 2, pfail psucc; psucc pfail); % NOT gateG = mk_all_dags(N);nhyp = length(G);A BA Bhyp_bnet = cell(1, nhyp);for h=1:nhyphyp_bneth = mk_bnet(Gh, ns);for i=1:Nfam = fami
10、ly(Gh, i);hyp_bneth.CPDi= tabular_CPD(hyp_bneth, i, prior_type, dirichlet, .CPT, unif);endend% 结构先验，假设为均匀分布 prior = normalise(ones(1, nhyp);% 保存结果并初始化训练结构init_hyp_bnet = hyp_bnet; init_prior = prior; rand(state, 0);hyp_w = zeros(ntrials+1, nhyp);hyp_w(1,:) = prior(:); data = zeros(N, ntrials);% 逐步计算
11、后验概率 LL = zeros(1, nhyp);ll = zeros(1, nhyp);for t=1:ntrialsev = cell2num(sample_bnet(true_bnet);data(:,t) = ev;for i=1:nhypll(i) = log_marg_lik_complete(hyp_bneti, ev);hyp_bneti = bayes_update_params(hyp_bneti, ev);endprior = normalise(prior .* exp(ll);LL = LL + ll;hyp_w(t+1,:) = prior;end% Plot，模型
12、后验概率% Red = model 1 (no arcs), blue/green = models 2/3 (1 arc)% Blue = model 2 (2-1)% Green = model 3 (1-2, “ground truth“)if 1figure;m = size(hyp_w, 1);h=plot(1:m, hyp_w(:,1), r-, 1:m, hyp_w(:,2), b-., 1:m, hyp_w(:,3), g:);axis(0 m 0 1)title(model posterior vs. time)drawnowend% 检验结果hyp_bnet2 = init
13、_hyp_bnet;prior2 = init_prior;图 3-1 贝叶斯模型选择后验概率对比BNT 中的结构学习程序可以按类似参数学习的情况分成四类：Full obs Partial obsPoint learn_struct_K2 not yet supportedBayes learn_struct_mcmc not yet supported3.1 Markov 等效如果两个 DAGs 编码同样的条件独立，它们被叫做 Markov 等效。所有 DAGs 的集合可以被分割成 Markov 等效类。同一类内的线图可以有方向，它们弧的颠倒不会改变任何 CI 关系。每一类都可以用一个 PD
14、AG（partially directed acyclic graph，局部有向非循环图）这种图被称为本质图或方向图。这个详细说明哪个边必须在某一个方位上被定向，哪个可能被颠倒。3.2 穷举搜索结构学习的强有力手段是列举 DAGs 的所有可能性，并对它们一一打分。这为其它算法的比较提供了一个“黄金标准” 。我们按如下做：dags = mk_all_dags(N);score = score_dags(data, ns, dags);默认的情况下，我们使用贝叶斯打分规则，并假定 CPDs 是用带有 BDeu 的先验表表示的。如果想是用一致的先验值，我们可以通过如下方法覆盖这些默认值。params
15、 = cell(1,N); for i=1:N paramsi = prior, unif; end score = score_dags(data, ns, dags, params, params, scoring_fn, bic); 实际上不能列举 N5 的所有可能的 DAGs。3.3 K2 算法K2 算法(Cooper and Herskovits, 1992)是一种按如下方式工作的贪婪搜索算法。每一个起始点没有父节点。然后增加结果结构打分最高时的父节点。当单独添加父节点再不能提高分数时，停止添加父节点。当我们使用固定的顺序时，我们不需要做循环检查，也不需要为每个节点单独选择父节点。B
16、NT 推广了这点允许使用任何种类的 CPD，无论贝叶斯打分规则还是 BIC，另外，你可以对每一个节点指定一个任意的父节点数量的上限。order = C S R W;max_fan_in = 2;dag2 = learn_struct_K2(data, ns, order, max_fan_in, max_fan_in);3.3 爬山法爬山算法从状态空间中的一个指定点开始，考虑所有最接近的邻节点，然后移向得分最高的相邻节点。当相邻节点得分没有高于当前节点时（例如到达了局部最大值。 ） ，算法停止。然后从空间的其它部分重新开始。“相邻”通常定义为所有的图可以通过从当前的图添加、删除或翻转一个单独的
17、弧得出，并服从无环的约束。其它相邻的可能详。learn_struct_hc()3.4 MCMC使用 Metropolis-Hastings (MH)的马尔可夫链蒙特卡尔算法来搜索所有的图空间。标准的分配提案是考虑移动所有最近的按上面定义的邻节点。这个函数可以按如下方法调用：sampled_graphs, accept_ratio = learn_struct_mcmc(data, ns, nsamples, 500, burnin, 10);图 3-2 精确后验概率和 MCMC 后验概率对比图 3-3 MCMC 接受率3.5 SEM 算法计算贝叶斯打分时，有部分是计算具有挑战性的观测，因为参数
18、学习的后验概率变成了多峰的状态（这是由于隐含节点导致了一个混合的分布） 。因此需要使用逼近算法，如 BIC。不幸的是搜索算法仍然是代价高昂的，因为我们需要在每一步运行 EM 算法来计算 MLE值，它需要对每一个模型进行计算打分。一个变换的方法是在每步进行局域搜索来替代第 M 步的 EM，当数据是“添满”状态时这种方法非常有效。以上被称为结构上的 EM 算法（Friedman 1997）它可以通过 BIC 打分收敛的局部最大值来证明。4 推断引擎创立好一个贝叶斯网络，我们现在可以用它来进行推断。贝叶斯网络中有许多不同的算法来作为推断的的工具，在速度、复杂性、普遍性和精确性上有不同的表现。BNT
19、因此提供了多种多样的不同的推断引擎。推理引擎是一个包含了 bnet (Bayesian net )支持 enter_evidence 和marginal_nodes 方法的对象。引擎设计者把 bnet 作为一个自变量，并且可以执行一些特殊处理的模型。当调用 enter_evidence，引擎可以处理一些经过特殊处理的证据。最后，当调用， marginal_nodes 引擎可以执行一些特殊处理的查询。4.1 联合树推断引擎（所有精确推断引擎的根本）engine = jtree_inf_engine(bnet);4.2 全局推理算法最简单的推理方法是直接构建所有结点的联合分布，然后得到边缘概率。这
20、已在 global_joint_inf_engine 中实现，但它仅适用于教学和调试。4.3 近似传播引擎likelihood_weighting_inf_engine，可以完成重要采样，并能够处理任意种类的结点。 gibbs_sampling_inf_engine，由 Bhaskara Marthi 写的. 目前它仅能处理表格条件概率分布（tabular CPDs） 。5 上市公司贝叶斯网络模型5.1 数据研究上市公司财务指标间的关系，本文共选取了 17 个财务指标的 240 条样本数据。17 个财务指标如下表所示。表 5-1 财务指标对应编号编号 名称 编号 名称X1 每股收益 X10 总资产增长率X2 每股净资产 X11 股东权益增长率X3 调整后净资产 X12 每股收益增长率X4 净资产收益率 X13 流动比率X5 股东权益比率 X14 速动比率X6 每股未分利润 X15 应收帐款周转率
