comparison toolboxes/FullBNT-1.0.7/bnt/learning/dirichlet_score_family.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function LL = dirichlet_score_family(counts, prior)
+% DIRICHLET_SCORE Compute the log marginal likelihood of a single family
+% LL = dirichlet_score(counts, prior)
+%
+% counts(a, b, ..., z) is the number of times parent 1 = a, parent 2 = b, ..., child = z
+% prior is an optional multidimensional array of the same shape as counts.
+% It defaults to a uniform prior.
+% 
+% We marginalize out the parameters:
+% LL = log \int \prod_m P(x(i,m) | x(Pa_i,m), theta_i) P(theta_i) d(theta_i)
+
+
+% LL = log[  prod_j gamma(alpha_ij)/gamma(alpha_ij + N_ij)  *
+%            prod_k gamma(alpha_ijk + N_ijk)/gamma(alpha_ijk)  ]
+% Call the prod_k term U and the prod_j term  V.
+% We reshape all quantities into (j,k) matrices
+% This formula was first derived by Cooper and Herskovits, 1992.
+% See also "Learning Bayesian Networks", Heckerman, Geiger and Chickering, MLJ 95.
+
+ns = mysize(counts);
+ns_ps = ns(1:end-1);
+ns_self = ns(end);
+
+if nargin < 2, prior = normalise(myones(ns)); end
+
+
+if 1
+  prior = reshape(prior(:), [prod(ns_ps) ns_self]);
+  counts = reshape(counts,  [prod(ns_ps) ns_self]);
+  %U = prod(gamma(prior + counts) ./ gamma(prior), 2); % mult over k
+  LU = sum(gammaln(prior + counts) - gammaln(prior), 2);
+  alpha_ij = sum(prior, 2); % sum over k
+  N_ij = sum(counts, 2);
+  %V = gamma(alpha_ij) ./ gamma(alpha_ij + N_ij);
+  LV = gammaln(alpha_ij) - gammaln(alpha_ij + N_ij);
+  %L = prod(U .* V);
+  LL = sum(LU + LV);
+else
+  CPT = mk_stochastic(prior + counts);
+  LL = sum(log(CPT(:) .* counts(:)));
+end
+