Daniel@0: function [S, LL] = bic_score(counts, CPT, ncases)
Daniel@0: % BIC_SCORE Bayesian Information Criterion score for a single family
Daniel@0: % [S, LL] = bic_score(counts, CPT, ncases)
Daniel@0: %
Daniel@0: % S is a large sample approximation to the log marginal likelihood,
Daniel@0: % which can be computed using dirichlet_score.
Daniel@0: %
Daniel@0: % S  = \log [ prod_j _prod_k theta_ijk ^ N_ijk ]  - 0.5*d*log(ncases) 
Daniel@0: % where counts encode N_ijk, theta_ijk is the MLE comptued from counts,
Daniel@0: % and d is the num of free parameters.
Daniel@0: 
Daniel@0: %CPT = mk_stochastic(counts);
Daniel@0: tiny = exp(-700);
Daniel@0: LL = sum(log(CPT(:)  + tiny) .* counts(:));
Daniel@0: % CPT(i) = 0 iff counts(i) = 0 so it is okay to add tiny
Daniel@0: 
Daniel@0: ns = mysize(counts);
Daniel@0: ns_ps = ns(1:end-1);
Daniel@0: ns_self = ns(end);
Daniel@0: nparams = prod([ns_ps (ns_self-1)]);
Daniel@0: % sum-to-1 constraint reduces the effective num. vals of the node by 1
Daniel@0: 
Daniel@0: S = LL - 0.5*nparams*log(ncases);