Talk is cheap, I show you the code
第一章的作业主要是关于PGM的因子操作。实际上，因子是整个概率图的核心。对于有向图而言，因子对应的是CPD（条件分布）；对无向图而言，因子对应的是势函数。总而言之，因子是一个映射，将随机变量空间映射到实数空间。因子表现的是对变量之间关系的一种设计。每个因子都编码了一定的信息。
因子的数据结构：
phi = struct('var', [3 1 2], 'card', [2 2 2], 'val', ones(1, 8));
在matlab中，因子被定义为一个结构体。结构体中有三个变量，分别是 var : variable, card : cardinate, val : value.这三个变量构成了因子表。其中，var 后面是变量名 X_3，X_2,X_1. card 是每个变量所对应的取值范围。[2 2 2]表示这三个变量都是二值变量。如果是骰子，则应该写成[6 6 6]。val 后面接的是一个列向量。该向量的长度应该为 prod(card).
1、因子相乘
两个Scope相同或不同或部分不同的因子可以相乘，称为因子联合。对于有向图而言，因子相乘就是贝耶斯的链式法则。对于无向图而言，因子相乘就是合并两个相似的信息。因子相乘的原则是能够冲突，也就是只有对应项相乘（因子都是表）。
% FactorProduct Computes the product of two factors.
% C = FactorProduct(A,B) computes the product between two factors, A and B,
% where each factor is defined over a set of variables with given dimension.
% The factor data structure has the following fields:
% .var Vector of variables in the factor, e.g. [1 2 3]
% .card Vector of cardinalities corresponding to .var, e.g. [2 2 2]
% .val Value table of size prod(.card)
%
% See also FactorMarginalization.m, IndexToAssignment.m, and
% AssignmentToIndex.m function C = FactorProduct(A, B)
%A = struct('var', [1], 'card', [2], 'val', [0.11, 0.89]);
%B = struct('var', [2, 1], 'card', [2, 2], 'val', [0.59, 0.41, 0.22, 0.78]);
% Check for empty factors
if (isempty(A.var)), C = B; return; end;
if (isempty(B.var)), C = A; return; end; % Check that variables in both A and B have the same cardinality
[dummy iA iB] = intersect(A.var, B.var);
if ~isempty(dummy)
% A and B have at least 1 variable in common
assert(all(A.card(iA) == B.card(iB)), 'Dimensionality mismatch in factors');
end % Set the variables of C
C.var = union(A.var, B.var); % Construct the mapping between variables in A and B and variables in C.
% In the code below, we have that
%
% mapA(i) = j, if and only if, A.var(i) == C.var(j)
%
% and similarly
%
% mapB(i) = j, if and only if, B.var(i) == C.var(j)
%
% For example, if A.var = [3 1 4], B.var = [4 5], and C.var = [1 3 4 5],
% then, mapA = [2 1 3] and mapB = [3 4]; mapA(1) = 2 because A.var(1) = 3
% and C.var(2) = 3, so A.var(1) == C.var(2). [dummy, mapA] = ismember(A.var, C.var);
[dummy, mapB] = ismember(B.var, C.var); % Set the cardinality of variables in C
C.card = zeros(1, length(C.var));
C.card(mapA) = A.card;
C.card(mapB) = B.card; % Initialize the factor values of C:
% prod(C.card) is the number of entries in C
C.val = zeros(1, prod(C.card)); % Compute some helper indices
% These will be very useful for calculating C.val
% so make sure you understand what these lines are doing.
assignments = IndexToAssignment(1:prod(C.card), C.card);
indxA = AssignmentToIndex(assignments(:, mapA), A.card);
indxB = AssignmentToIndex(assignments(:, mapB), B.card); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE:
% Correctly populate the factor values of C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:size(indxA)
C.val(i) = A.val(indxA(i))*B.val(indxB(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end
2、变量边际化
变量边际化是概率图模型的核心操作，其目的是消除其他变量，获得单个变量的边缘分布。例如有CPD: struct('var', [3 1 2], 'card', [2 2 2], 'val', ones(1, 8)), 如果我们需要知道x_1 最本质的分布，则需要对x_2,x_3进行边际化操作。变量边际化最大的难点在于算法。比如要边际 x_2，那么应该把x_1,x_3取值相同的组合求和。此函数用了一个巧妙的方法来解决：先由var 求assignment,之后抽去assignment里需要消除的变量所对应列，在将assignment转成index.此时index里数字相同的编码就是需要求和的。
% FactorMarginalization Sums given variables out of a factor.
% B = FactorMarginalization(A,V) computes the factor with the variables
% in V summed out. The factor data structure has the following fields:
% .var Vector of variables in the factor, e.g. [1 2 3]
% .card Vector of cardinalities corresponding to .var, e.g. [2 2 2]
% .val Value table of size prod(.card)
%
% The resultant factor should have at least one variable remaining or this
% function will throw an error.
%
% See also FactorProduct.m, IndexToAssignment.m, and AssignmentToIndex.m function B = FactorMarginalization(A, V)
% A = Joint;
% V = [2];
% Check for empty factor or variable list
if (isempty(A.var) || isempty(V)), B = A; return; end; % Construct the output factor over A.var \ V (the variables in A.var that are not in V)
% and mapping between variables in A and B
[B.var, mapB] = setdiff(A.var, V); % Check for empty resultant factor
if isempty(B.var)
error('Error: Resultant factor has empty scope');
end; % Initialize B.card and B.val
B.card = A.card(mapB);
B.val = zeros(1, prod(B.card)); % Compute some helper indices
% These will be very useful for calculating B.val
% so make sure you understand what these lines are doing
assignments = IndexToAssignment(1:length(A.val), A.card);
indxB = AssignmentToIndex(assignments(:, mapB), B.card); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
% Correctly populate the factor values of B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% shunxu = 1:size(indxB);
hunhe = [indxB,shunxu'];
hunhe_paixu = sortrows(hunhe,1);
tmp_1 = hunhe_paixu(:,1);
tmp_2 = hunhe_paixu(:,2);
k = 1;
for i = 1:length(tmp_1)-1
if tmp_1(i) == tmp_1(i+1)
B.val(k) = B.val(k) + A.val(tmp_2(i));
else
B.val(k) = B.val(k) + A.val(tmp_2(i));
k = k+1;
end
end
B.val(k) = B.val(k) + A.val(i+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
3、变量观测
变量观测是将未知的变量置为已知。那么该变量其他不符合观测条件的取值都应该置为0. 该算法的核心思想是首先获得变量的assignment，之后选出被观测变量所占assignment的列。最后依次读取该列的值，将不符合观测结果的value置0.
% ObserveEvidence Modify a vector of factors given some evidence.
% F = ObserveEvidence(F, E) sets all entries in the vector of factors, F,
% that are not consistent with the evidence, E, to zero. F is a vector of
% factors, each a data structure with the following fields:
% .var Vector of variables in the factor, e.g. [1 2 3]
% .card Vector of cardinalities corresponding to .var, e.g. [2 2 2]
% .val Value table of size prod(.card)
% E is an N-by-2 matrix, where each row consists of a variable/value pair.
% Variables are in the first column and values are in the second column. function F = ObserveEvidence(F, E) % Iterate through all evidence for i = 1:size(E, 1)
v = E(i, 1); % variable
x = E(i, 2); % value % Check validity of evidence
if (x == 0),
warning(['Evidence not set for variable ', int2str(v)]);
continue;
end; for j = 1:length(F),
% Does factor contain variable?
indx = find(F(j).var == v); if (~isempty(indx)), % Check validity of evidence
if (x > F(j).card(indx) || x < 0 ),
error(['Invalid evidence, X_', int2str(v), ' = ', int2str(x)]);
end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
% Adjust the factor F(j) to account for observed evidence
% Hint: You might find it helpful to use IndexToAssignment
% and SetValueOfAssignment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:prod(F(j).card)
Assignment_ = IndexToAssignment(k,F(j).card);
if Assignment_(indx) ~= x
indx_ = AssignmentToIndex(Assignment_,F(j).card);
F(j).val(indx_) = 0;
end
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Check validity of evidence / resulting factor
if (all(F(j).val == 0)),
warning(['Factor ', int2str(j), ' makes variable assignment impossible']);
end; end;
end;
end; end
4、计算联合分布
计算联合分布就是将若干个变量进行相乘。联合分布的输入是因子集。将因子集中的变量依次相乘，所得最后结果则为联合分布。
%ComputeJointDistribution Computes the joint distribution defined by a set
% of given factors
%
% Joint = ComputeJointDistribution(F) computes the joint distribution
% defined by a set of given factors
%
% Joint is a factor that encapsulates the joint distribution given by F
% F is a vector of factors (struct array) containing the factors
% defining the distribution
% % FACTORS.INPUT(1) = struct('var', [1], 'card', [2], 'val', [0.11, 0.89]);
%
% % FACTORS.INPUT(2) contains P(X_2 | X_1)
% FACTORS.INPUT(2) = struct('var', [2, 1], 'card', [2, 2], 'val', [0.59, 0.41, 0.22, 0.78]);
%
% % FACTORS.INPUT(3) contains P(X_3 | X_2)
% FACTORS.INPUT(3) = struct('var', [3, 2], 'card', [2, 2], 'val', [0.39, 0.61, 0.06, 0.94]);
%
% F = FACTORS.INPUT; function Joint = ComputeJointDistribution(F) % Check for empty factor list
if (numel(F) == 0)
warning('Error: empty factor list');
Joint = struct('var', [], 'card', [], 'val', []);
return;
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE:
% Compute the joint distribution defined by F
% You may assume that you are given legal CPDs so no input checking is required.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_ = F;
%Joint = struct('var', [], 'card', [], 'val', []); % Returns empty factor. Change this.
for i = 2 : numel(F_)
F_(i) = FactorProduct(F_(i),F_(i-1));
end
Joint = F_(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
5、计算边缘分布
边缘分布是在联合分布的基础上更进一步的处理，对于给定联合分布，某些变量可能被观测到。某些变量可能需要边际掉，最后的结果就是边缘分布。边缘分布得名于其在因子表中处以边缘的位置。其关键操作在于得到联合分布后必须归一化。因为边缘分布的和总是1.
%ComputeMarginal Computes the marginal over a set of given variables
% M = ComputeMarginal(V, F, E) computes the marginal over variables V
% in the distribution induced by the set of factors F, given evidence E
%
% M is a factor containing the marginal over variables V
% V is a vector containing the variables in the marginal e.g. [1 2 3] for
% X_1, X_2 and X_3.
% F is a vector of factors (struct array) containing the factors
% defining the distribution
% E is an N-by-2 matrix, each row being a variable/value pair.
% Variables are in the first column and values are in the second column.
% If there is no evidence, pass in the empty matrix [] for E. % % FACTORS.INPUT(1) contains P(X_1)
% FACTORS.INPUT(1) = struct('var', [1], 'card', [2], 'val', [0.11, 0.89]);
%
% % FACTORS.INPUT(2) contains P(X_2 | X_1)
% FACTORS.INPUT(2) = struct('var', [2, 1], 'card', [2, 2], 'val', [0.59, 0.41, 0.22, 0.78]);
%
% % FACTORS.INPUT(3) contains P(X_3 | X_2)
% FACTORS.INPUT(3) = struct('var', [3, 2], 'card', [2, 2], 'val', [0.39, 0.61, 0.06, 0.94]);
%
% V = [3];
% F = FACTORS.INPUT;
% E = [];
function M = ComputeMarginal(V, F, E) % Check for empty factor list
if (numel(F) == 0)
warning('Warning: empty factor list');
M = struct('var', [], 'card', [], 'val', []);
return;
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE:
% M should be a factor
% Remember to renormalize the entries of M!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Joint = ComputeJointDistribution(F);
Obser = ObserveEvidence(Joint,E);
To_be_Mglzed = setdiff(Joint.var,V);
if ~isempty(To_be_Mglzed)
M = FactorMarginalization(Obser,To_be_Mglzed);
else
M = Obser;
end
M.val = M.val/sum(M.val);
% M = struct('var', [], 'card', [], 'val', []); % Returns empty factor. Change this. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M = StandardizeFactors(M);
end
最后，所有代码请点这里