机器学习 —— 概率图模型（Homework: Exact Inference）

　　在前三周的作业中，我构造了概率图模型并调用第三方的求解器对器进行了求解，最终获得了每个随机变量的分布（有向图），最大后验分布（双向图）。本周作业的主要内容就是自行编写概率图模型的求解器。实际上，从根本上来说求解器并不是必要的。其作用只是求取边缘分布或者MAP，在得到联合CPD后，寻找联合CPD的最大值即可获得MAP，对每个变量进行边缘分布求取即可获得边缘分布。但是，这种简单粗暴的方法效率极其低下，对于MAP求取而言，每次得到新的evidance时都要重新搜索CPD，对于单个变量分布而言，更是对每个变量都要反复进行整体边缘化。以一个长度为6字母的单词为例，联合CPD有着多达26^6个数据，反复操作会浪费大量的计算资源。

　　团树算法背后的思路是分而治之。对于一组随机变量ABCDEFG,如果A和其他变量之间是独立的，那么无论是求边缘分布还是MAP都可以将A单独考虑。如果ABC联系比较紧密，CDE联系比较紧密，那么如果两个团关于C的边缘分布是相同的，则我们没有必要将ABCDE全部乘在一起再来分别求各个变量的边缘分布。因为反过来想，乘的时候也只是把对应的C乘起来，如果C的边缘分布相同，在相乘的时候其实两个团之间并没有引入其他信息，此时乘法不会对ABDE的边缘分布产生影响。团树算法的数学过程和Variable Elimination是相同的。

　　PGM在计算机中的表达是factorLists，factor的var(i)，var表示节点连接关系。val描述了factor中var的关系。cliqueTree其实是一种特殊的factorLists，它的var是clique，表示一堆聚类的var。它的val表示的还是var之间的关系。只不过此时var之间的连接不复存在了。所以clique由两个变量组成：1、cliqueTree 2、edges.

　　团树算法的初始化可以分为两个过程：1、将变量抱团；2、获取团的初始势；

　　变量抱团是一个玄学过程，因为有很多不同的抱法，而且还都是对的。比较常见的是最小边，最小割等...其实如果是人来判断很容易就能得到结果，但是使用计算机算法则要费一些功夫了。不过这不涉及我们对团树算法的理解，所以Koller教授代劳了。

　　团的初始势表示团里变量之间的关系。其算法如下，需要注意的是不能重复使用factor.因为一个factor表达了一种关系，如果两个团里都有同一个factor，那么就是...这个事情。。。你帮他重复一遍。。。等于你也有责任的，晓得吧？

 1 %COMPUTEINITIALPOTENTIALS Sets up the cliques in the clique tree that is
 2 %passed in as a parameter.
 3 %
 4 %   P = COMPUTEINITIALPOTENTIALS(C) Takes the clique tree skeleton C which is a
 5 %   struct with three fields:
 6 %   - nodes: cell array representing the cliques in the tree.
 7 %   - edges: represents the adjacency matrix of the tree.
 8 %   - factorList: represents the list of factors that were used to build
 9 %   the tree. 
10 %   
11 %   It returns the standard form of a clique tree P that we will use through 
12 %   the rest of the assigment. P is struct with two fields:
13 %   - cliqueList: represents an array of cliques with appropriate factors 
14 %   from factorList assigned to each clique. Where the .val of each clique
15 %   is initialized to the initial potential of that clique.
16 %   - edges: represents the adjacency matrix of the tree. 
17 %
18 % Copyright (C) Daphne Koller, Stanford University, 2012
19 
20 
21 
22 function P = ComputeInitialPotentials(C)
23 Input = C;
24 % number of cliques
25 N = length(Input.nodes);
26 
27 % initialize cluster potentials 
28 P.cliqueList = repmat(struct('var', [], 'card', [], 'val', []), N, 1);
29 P.edges = zeros(N);
30 
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 % YOUR CODE HERE
33 %
34 % First, compute an assignment of factors from factorList to cliques. 
35 % Then use that assignment to initialize the cliques in cliqueList to 
36 % their initial potentials. 
37 
38 % C.nodes is a list of cliques.
39 % So in your code, you should start with: P.cliqueList(i).var = C.nodes{i};
40 % Print out C to get a better understanding of its structure.
41 %
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 % N_factors = length(C.factorList);
44 for i = 1:N
45     k = 1;
46     clear clique_
47      N_factors = length(Input.factorList);
48     for j = 1:N_factors
49         if min(ismember(Input.factorList(j).var,Input.nodes{i}))
50             clique_(k) = Input.factorList(j);
51             k = k+1;
52             Input.factorList(j) =struct('var', [], 'card', [], 'val', []);
53         end
54     end
55     Joint_Dis_cliq = ComputeJointDistribution(clique_);
56     Joint_Dis_cliq_std = StandardizeFactors(Joint_Dis_cliq);
57     P.cliqueList(i) = Joint_Dis_cliq_std;
58 end
59 P.edges = Input.edges;
60 end