1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
function mpe = find_mpe(engine, new_evidence, max_over)
% FIND_MPE Find the most probable explanation of the data (assignment to the hidden nodes)
% function mpe = find_mpe(engine, evidence, order)
%
% PURPOSE:
% CALC_MPE Computes the most probable explanation to the network nodes
% given the evidence.
%
% [mpe, ll] = calc_mpe(engine, new_evidence, max_over)
%
% INPUT:
% bnet - the bayesian network
% new_evidence - optional, if specified - evidence to be incorporated [cell(1,n)]
% max_over - optional, if specified determines the variable elimination order [1:n]
%
% OUTPUT:
% mpe - the MPE assignmet for the net variables (or [] if no satisfying assignment)
% ll - log assignment probability.
%
% Notes:
% 1. Adapted from '@var_elim_inf_engine\marginal_nodes' for MPE by Ron Zohar, 8/7/01
% 2. Only discrete potentials are supported at this time.
% 3. Complexity: O(nw*) where n is the number of nodes and w* is the induced tree width.
% 4. Implementation based on:
% - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference",
% UA1 96, pp. 211-219.
bnet = bnet_from_engine(engine);
ns = bnet.node_sizes;
n = length(bnet.dag);
evidence = cell(1,n);
if (nargin<2)
new_evidence = evidence;
end
onodes = find(~isemptycell(new_evidence)); % observed nodes
hnodes = find(isemptycell(new_evidence)); % hidden nodes
pot_type = determine_pot_type(bnet, onodes);
if pot_type ~= 'd'
error('only disrete potentials supported at this time')
end
for i=1:n
fam = family(bnet.dag, i);
CPT{i} = convert_to_pot(bnet.CPD{bnet.equiv_class(i)}, pot_type, fam(:), evidence);
end
% handle observed nodes: set impossible cases' probability to zero
% rather than prun matrix (this makes backtracking easier)
for ii=onodes
lIdx = 1:ns(ii);
lIdx = setdiff(lIdx, new_evidence{ii});
sCPT=struct(CPT{ii}); % violate object privacy
sargs = '';
for jj=1:(length(sCPT.domain)-1)
sargs = [sargs, ':,'];
end
for jj=lIdx
eval(['sCPT.T(', sargs, num2str(jj), ')=0;']);
end
CPT{ii}=dpot(sCPT.domain, sCPT.sizes, sCPT.T);
end
B = cell(1,n);
for b=1:n
B{b} = mk_initial_pot(pot_type, [], [], [], []);
end
if (nargin<3)
max_over = (1:n);
end
order = max_over; % no attempt to optimize this
% Initialize the buckets with the CPDs assigned to them
for i=1:n
b = bucket_num(domain_pot(CPT{i}), order);
B{b} = multiply_pots(B{b}, CPT{i});
end
% Do backward phase
max_over = max_over(length(max_over):-1:1); % reverse
maximize = 1;
for i=max_over(1:end-1)
% max-ing over variable i which occurs in bucket j
j = bucket_num(i, order);
rest = mysetdiff(domain_pot(B{j}), i);
%temp = marginalize_pot_max(B{j}, rest);
temp = marginalize_pot(B{j}, rest, maximize);
b = bucket_num(domain_pot(temp), order);
% fprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b);
sB=struct(B{b}); % violate object privacy
if ~isempty(sB.domain)
B{b} = multiply_pots(B{b}, temp);
else
B{b} = temp;
end
end
result = B{1};
marginal = pot_to_marginal(result);
[prob, mpe] = max(marginal.T);
% handle impossible cases
if ~(prob>0)
mpe = [];
ll = -inf;
%warning('evidence has zero probability')
return
end
ll = log(prob);
% Do forward phase
for ii=2:n
marginal = pot_to_marginal(B{ii});
mpeidx = [];
for jj=order(1:length(mpe))
%assert(ismember(jj, marginal.domain)) %%% bug
temp = find_equiv_posns(jj, marginal.domain);
mpeidx = [mpeidx, temp] ;
if isempty(temp)
mpeidx = [mpeidx, Inf] ;
end
end
[mpeidxsorted sortedtompe] = sort(mpeidx) ;
% maximize the matrix obtained from assigning values from previous buckets.
% this is done by building a string and using eval.
kk=1;
sargs = '(';
for jj=1:length(marginal.domain)
if (jj~=1)
sargs = [sargs, ','];
end
if (mpeidxsorted(kk)==jj)
sargs = [sargs, num2str(mpe(sortedtompe(kk)))];
if (kk<length(mpe))
kk = kk+1 ;
end
else
sargs = [sargs, ':'];
end
end
sargs = [sargs, ')'] ;
eval(['[val, loc] = max(marginal.T', sargs, ');'])
mpe = [mpe loc];
end
[I,J] = sort(order);
mpe = mpe(J);
mpe = num2cell(mpe);
%%%%%%%%%
function b = bucket_num(domain, order)
b = max(find_equiv_posns(domain, order));
|
