贝叶斯网络推理之最大可能解释问题（基于FullBNT-1.0.4的MATLAB实现）

题目：贝叶斯网络推理之最大可能解释问题（基于FullBNT-1.0.4的MATLAB实现）
        看本篇前，如果对最大可能解释的概念不汪楚，可以阅读《贝叶斯网络与最大可能解释(MPE)问题》，如果想简单了解贝叶斯网络推理算法，可以阅读《贝叶斯网络推理算法简单罗列》，尤其是常用的联结树算法（junction tree），在文中附录摘有详细解释。
        在BNT工具箱中，最大可能解释问题的求解函数是find_mpe，该函数在不同的推理引擎下进行重载：
      hmm_inf_engine/find_mpe
      smoother_engine/find_mpe
      var_elim_inf_engine/find_mpe
      jtree_mnet_inf_engine/find_mpe
      jtree_inf_engine/find_mpe
      global_joint_inf_engine/find_mpe
      belprop_mrf2_inf_engine/find_mpe
      belprop_inf_engine/find_mpe
      belprop_fg_inf_engine/find_mpe
        有关find_mpe的使用例子在目录\FullBNT-1.0.4\BNT\examples\static中有两个m文件：mpe1.m和mpe2.m。两个例子大同小异，以下就以mpe1.m为例来说明。
        例子mpe1.m的源码如下：seed = 1;
rand('state', seed);
randn('state', seed);

N = 4;
dag = zeros(N,N);
C = 1; S = 2; R = 3; W = 4;
dag(C,[R S]) = 1;
dag(R,W) = 1;
dag(S,W)=1;

false = 1; true = 2;
ns = 2*ones(1,N); % binary nodes

bnet = mk_bnet(dag, ns);
if 0
bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
else
for i=1:N, bnet.CPD{i} = tabular_CPD(bnet, i); end
end

evidence = cell(1,N);
onodes = [1 3];
data = sample_bnet(bnet);
evidence(onodes) = data(onodes);

clear engine;
engine{1} = belprop_inf_engine(bnet);
engine{2} = jtree_inf_engine(bnet);
engine{3} = global_joint_inf_engine(bnet);
engine{4} = var_elim_inf_engine(bnet);
E = length(engine);

clear mpe;
for e=1:E
mpe{e} = find_mpe(engine{e}, evidence);
end
for e=2:E
assert(isequal(mpe{1}, mpe{e}))
end
        整个程序不足50行。前面一半就是人为地生成了一个贝叶斯网络（第5~10行规定了贝叶斯网络结构，第15~23行规定了贝叶斯网络参数）。这个贝叶斯网络是“草地潮湿原因”模型（详见本文结尾链接），初始化的贝叶斯网络结构如下（C—Cloudy, S—Sprinkler, R—Rain, W—WetGrass）：


        这里再多说几句贝叶斯网络参数的设定，第16行到23行用了一个if…else…结构，第16行默认是if 0，即执行else部分内容（第22行），第22行实际上是随机初始化贝叶斯网络参数，当然，因为第1~3行对对伪随机数初始状态进行了限定，因此每次执行的结果都是一样的（具体可以参见：Matlab中rand('state',s)和rand('state',0)表示什么意思？，链接：http://www.ilovematlab.cn/thread-57952-1-1.html）；若将第16行改为if 1，则执行第17~20行参数初始化代码，这个初始化参数什么意思呢？
        


        第17行很简单，天气是不是Cloudy的概率均为0.5；
        第18行初始化Rain概率，如下：
[align=center]

	Rain=False	Rain=True
Cloudy=False	0.8	0.2
Cloudy=True	0.2	0.8

[/align]即天气是Cloudy时Rain的概率是0.8（不Rain的概率是0.2），天气不是Cloudy时Rain的概率是0.2（不Rain的概率是0.8）；
        第19行初始化打开Sprinkler的概率，如下：
[align=center]

	Sprinkler=False	Sprinkler=True
Cloudy=False	0.5	0.5
Cloudy=True	0.9	0.1

[/align]即天气是Cloudy时打开Sprinkler的概率是0.1（不打开Sprinkler的概率是0.9），天气不是Cloudy时打开Sprinkler的概率是0.5（不打开Sprinkler的概率也是0.5）；
        第20行初始化WetGrass的概率，如下：
[align=center]

	WetGrass=False	WetGrass=True
Rain=False, Sprinkler=False	1	0
Rain=False, Sprinkler=True	0.1	0.9
Rain=True, Sprinkler=False	0.1	0.9
Rain=True, Sprinkler=True	0.01	0.99

[/align](注：Rain和Sprinkler的先后顺序不太确定，虽然此处并没有影响)
即天气不下雨(Rain=False)也不洒水(Sprinkler=False)时草地肯定是干的(WetGrass=False)，天气下雨(Rain=True)同时又洒水(Sprinkler=True)时草地湿(WetGrass=True)的概率是0.99（草地不湿的概率是0.01），天气下雨和洒水器打开二者仅有一个为真时草地湿(WetGrass=True)的概率是0.9（草地不湿的概率是0.1）；
        证据变量由第27~30行给定，这个写法是很固定的，证据变量均初始化为cell类型，贝叶斯网络有多少个节点，cell变量就初始化为多大（evidence = cell(1,N);），然后对证据变量赋值即可（evidence(onodes)= data(onodes);），这里N=4，证据变量为第1个和第3个节点（onodes = [1 3];个人猜测 onodes应该是observed nodes的简写吧），证据变量的值是根据贝叶斯网生成了一组数据（data = sample_bnet(bnet);），然后将这组数据中的第1个和第3个节点的值赋给了证据变量（evidence(onodes) = data(onodes);）。
        


        推理引擎由第32~37行生成，共生成了四个推理引擎（第33~36行），当然具体使用时根据需要生成一个即可，这里是为了比较各推理引擎的推理结果是否相同。
        最大可能解释的求解是由39~42行完成的，分别使用四个推理引擎得到MPE结果。 
        第43~45行对比第2~4个推理引擎的推理结果与第1个引擎的推理结果是否相同，若不同则会报错。
        运行mpe1.m，不会出现任何结果，工作区变量如下：


        若想查看各引擎的推理结果，双击工作区的mpe变量：


        共有四个引擎，双击上图中的每个“1×4 cell”，可以发现都是一样的：


        这个结果什么意思呢？注意程序前一半对贝叶斯网络的设定，这里证据变量是节点1和节点3（值均为1），预测结果是节点2和节点4的值也均为1。前面已经提到，这个贝叶斯网络是“草地潮湿原因”模型，代码第7行的四个节点均为英文单词首字母C—Cloudy,S—Sprinkler, R—Rain, W—WetGrass，贝叶斯网络结构如下：


        现在已知天气并非Cloudy（多云or阴天？）且并没有Rain（下雨），预测结果是Sprinkler（洒水器）没有开且Grass不是Wet，注意程序设定的是值为1即false，值为2时即true。
        不明白为什么天气非Cloudy且没有下雨的前提下不开洒水器呢？
        OK，到此讲完了，总结一下find_mpe函数的使用方法：
        该函数需要输出两个参数：推理引擎和证据变量，因此首先要准备好推理引擎和证据变量；推理引擎只需要调用对应的函数即可，推理引擎函数输入是贝叶斯网络bnet（在第15行命名的，已知网络的结构和参数）；证据变量给定比较简单，初始化一个cell类型变量，把证据变量位置赋值即可。

        其实BNT工具箱里的例子有两个在网上流传较广的，一个就是“草地潮湿原因”，另一个是“是窃贼还是地震”：
        WetGrass相关：
                http://blog.csdn.net/a_302/article/details/17916569
                https://zhidao.baidu.com/question/1769951750149622540.html
                http://blog.sina.com.cn/s/blog_7600796a0100p5lx.html
                http://blog.csdn.net/yefengnidie/article/details/4247254
                https://wenku.baidu.com/view/379405dcd15abe23482f4d1c.html
        Burglary相关：
                https://wenku.baidu.com/view/4a775b70336c1eb91a375d33.html
                https://my.oschina.net/SnifferApache/blog/343756