wolffd@0: function loglik = bethe_free_energy(engine, evidence)
wolffd@0: % BETHE_FREE_ENERGY Compute Bethe free energy approximation to the log likelihood
wolffd@0: % loglik = bethe_free_energy(engine, evidence)
wolffd@0: %
wolffd@0: % The Bethe free energy is given by an exact energy term and an approximate entropy term.
wolffd@0: % Energy
wolffd@0: %  E = -sum_f sum_i b(f,i) ln theta(f,i)
wolffd@0: % where b(f,i) = approximate Pr(family f = i) 
wolffd@0: % and theta(f,i) = Pr(f = i)
wolffd@0: % Entropy
wolffd@0: %  S = H1 - H2
wolffd@0: %  H1 = sum_f sum_p H(b(f))
wolffd@0: % where b(f) = belief on family f, H(.) = entropy
wolffd@0: %  H2 = sum_n (q(n)-1) H(b(n))
wolffd@0: % where q(n) = num. neighbors of n
wolffd@0: %
wolffd@0: % This function was written by Yair Weiss, 8/22/01.
wolffd@0: 
wolffd@0: hidden = find(isemptycell(evidence));
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: N = length(bnet.dag);
wolffd@0: 
wolffd@0: add_ev = 1;
wolffd@0: E=0;H1=0;H2=0;
wolffd@0: loglik=0;
wolffd@0: for n=1:N
wolffd@0:   ps=parents(bnet.dag,n);
wolffd@0:   if (length(ps)==0) % root node
wolffd@0:     qi=length(children(bnet.dag,n))-1;
wolffd@0:   else
wolffd@0:     qi=length(children(bnet.dag,n));
wolffd@0:   end
wolffd@0:   bf = marginal_family(engine, n, add_ev);
wolffd@0:   bf = bf.T(:);
wolffd@0:   e = bnet.equiv_class(n);
wolffd@0:   T = CPD_to_CPT(bnet.CPD{e});
wolffd@0:   T = T(:);
wolffd@0:   E = E-sum(log(T+(T==0)).*bf);
wolffd@0: 
wolffd@0:   if length(ps) > 0
wolffd@0:     % root nodes don't count as fmailies
wolffd@0:     H1 = H1+sum(log(bf+(bf==0)).*bf);
wolffd@0:   end
wolffd@0:   
wolffd@0:   bi = marginal_nodes(engine, n, add_ev);
wolffd@0:   bi = bi.T(:);
wolffd@0:   H2 = H2+qi*sum(log(bi+(bi==0)).*bi);
wolffd@0: end
wolffd@0: loglik=E+H1-H2;
wolffd@0: loglik=-loglik;