|
wolffd@0
|
1 function loglik = bethe_free_energy(engine, evidence)
|
|
wolffd@0
|
2 % BETHE_FREE_ENERGY Compute Bethe free energy approximation to the log likelihood
|
|
wolffd@0
|
3 % loglik = bethe_free_energy(engine, evidence)
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % The Bethe free energy is given by an exact energy term and an approximate entropy term.
|
|
wolffd@0
|
6 % Energy
|
|
wolffd@0
|
7 % E = -sum_f sum_i b(f,i) ln theta(f,i)
|
|
wolffd@0
|
8 % where b(f,i) = approximate Pr(family f = i)
|
|
wolffd@0
|
9 % and theta(f,i) = Pr(f = i)
|
|
wolffd@0
|
10 % Entropy
|
|
wolffd@0
|
11 % S = H1 - H2
|
|
wolffd@0
|
12 % H1 = sum_f sum_p H(b(f))
|
|
wolffd@0
|
13 % where b(f) = belief on family f, H(.) = entropy
|
|
wolffd@0
|
14 % H2 = sum_n (q(n)-1) H(b(n))
|
|
wolffd@0
|
15 % where q(n) = num. neighbors of n
|
|
wolffd@0
|
16 %
|
|
wolffd@0
|
17 % This function was written by Yair Weiss, 8/22/01.
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 hidden = find(isemptycell(evidence));
|
|
wolffd@0
|
20 bnet = bnet_from_engine(engine);
|
|
wolffd@0
|
21 N = length(bnet.dag);
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 add_ev = 1;
|
|
wolffd@0
|
24 E=0;H1=0;H2=0;
|
|
wolffd@0
|
25 loglik=0;
|
|
wolffd@0
|
26 for n=1:N
|
|
wolffd@0
|
27 ps=parents(bnet.dag,n);
|
|
wolffd@0
|
28 if (length(ps)==0) % root node
|
|
wolffd@0
|
29 qi=length(children(bnet.dag,n))-1;
|
|
wolffd@0
|
30 else
|
|
wolffd@0
|
31 qi=length(children(bnet.dag,n));
|
|
wolffd@0
|
32 end
|
|
wolffd@0
|
33 bf = marginal_family(engine, n, add_ev);
|
|
wolffd@0
|
34 bf = bf.T(:);
|
|
wolffd@0
|
35 e = bnet.equiv_class(n);
|
|
wolffd@0
|
36 T = CPD_to_CPT(bnet.CPD{e});
|
|
wolffd@0
|
37 T = T(:);
|
|
wolffd@0
|
38 E = E-sum(log(T+(T==0)).*bf);
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 if length(ps) > 0
|
|
wolffd@0
|
41 % root nodes don't count as fmailies
|
|
wolffd@0
|
42 H1 = H1+sum(log(bf+(bf==0)).*bf);
|
|
wolffd@0
|
43 end
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 bi = marginal_nodes(engine, n, add_ev);
|
|
wolffd@0
|
46 bi = bi.T(:);
|
|
wolffd@0
|
47 H2 = H2+qi*sum(log(bi+(bi==0)).*bi);
|
|
wolffd@0
|
48 end
|
|
wolffd@0
|
49 loglik=E+H1-H2;
|
|
wolffd@0
|
50 loglik=-loglik;
|
