|
wolffd@0
|
1 function LL = dirichlet_score_family(counts, prior)
|
|
wolffd@0
|
2 % DIRICHLET_SCORE Compute the log marginal likelihood of a single family
|
|
wolffd@0
|
3 % LL = dirichlet_score(counts, prior)
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % counts(a, b, ..., z) is the number of times parent 1 = a, parent 2 = b, ..., child = z
|
|
wolffd@0
|
6 % prior is an optional multidimensional array of the same shape as counts.
|
|
wolffd@0
|
7 % It defaults to a uniform prior.
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % We marginalize out the parameters:
|
|
wolffd@0
|
10 % LL = log \int \prod_m P(x(i,m) | x(Pa_i,m), theta_i) P(theta_i) d(theta_i)
|
|
wolffd@0
|
11
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 % LL = log[ prod_j gamma(alpha_ij)/gamma(alpha_ij + N_ij) *
|
|
wolffd@0
|
14 % prod_k gamma(alpha_ijk + N_ijk)/gamma(alpha_ijk) ]
|
|
wolffd@0
|
15 % Call the prod_k term U and the prod_j term V.
|
|
wolffd@0
|
16 % We reshape all quantities into (j,k) matrices
|
|
wolffd@0
|
17 % This formula was first derived by Cooper and Herskovits, 1992.
|
|
wolffd@0
|
18 % See also "Learning Bayesian Networks", Heckerman, Geiger and Chickering, MLJ 95.
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 ns = mysize(counts);
|
|
wolffd@0
|
21 ns_ps = ns(1:end-1);
|
|
wolffd@0
|
22 ns_self = ns(end);
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 if nargin < 2, prior = normalise(myones(ns)); end
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26
|
|
wolffd@0
|
27 if 1
|
|
wolffd@0
|
28 prior = reshape(prior(:), [prod(ns_ps) ns_self]);
|
|
wolffd@0
|
29 counts = reshape(counts, [prod(ns_ps) ns_self]);
|
|
wolffd@0
|
30 %U = prod(gamma(prior + counts) ./ gamma(prior), 2); % mult over k
|
|
wolffd@0
|
31 LU = sum(gammaln(prior + counts) - gammaln(prior), 2);
|
|
wolffd@0
|
32 alpha_ij = sum(prior, 2); % sum over k
|
|
wolffd@0
|
33 N_ij = sum(counts, 2);
|
|
wolffd@0
|
34 %V = gamma(alpha_ij) ./ gamma(alpha_ij + N_ij);
|
|
wolffd@0
|
35 LV = gammaln(alpha_ij) - gammaln(alpha_ij + N_ij);
|
|
wolffd@0
|
36 %L = prod(U .* V);
|
|
wolffd@0
|
37 LL = sum(LU + LV);
|
|
wolffd@0
|
38 else
|
|
wolffd@0
|
39 CPT = mk_stochastic(prior + counts);
|
|
wolffd@0
|
40 LL = sum(log(CPT(:) .* counts(:)));
|
|
wolffd@0
|
41 end
|
|
wolffd@0
|
42
|
