How is the log likelihood calculated for bayesian networks? In structure learning, there are score-based methods which rely on information criteria such as BIC or AIC.
BIC, specifically, is defined as:
$$
BIC = k \ln(n) - 2\ln\left(\hat{L}\right)
$$
Where $k$ is the number of parameters in the model, $n$ is the number of training examples and $\hat{L}$ is the likelihood function associating the model itself with observed data $x.$
In the case of structure learning, we are trying to determine how well a bayesian network approximates the underlying causal structure in the observed data. The bnlearn R package implements such calculations in its methods and, as far as I can tell, the log-likelihood is usually the preferred likelihood function, as it is supposed to be easier to compute.
So my main question here is: how is $\hat{L}$ calculated in the context of bayesian networks? I understand the "concept" (I think), but I can't seem to visualize a practical way of implementing that. What would be the most straight-forward way of representing such calculation?
 A: Conceptually, one could use the "standard" way of calculating likelihood, see How to calculate the likelihood function in conjunction with the rule of product decomposition i.e., chain rule for the network. So the joint-probability for the entire network reads $$P(x_{1}, x_{2}, ..., x_{n}) = {\Large \Pi} P(x_{i}|x_{parent})$$
$x_{parent}$ is the parent of $x_{i}$. One would need to incorporate  this in formulating "standard" likelihood $L(\theta|x_{1}, x_{2}, ..., x_{n})$, $\theta$ being the parameters.
A: These slides introduce many popular scoring methods for BNs: http://www.lx.it.pt/~asmc/pub/talks/09-TA/ta_pres.pdf
As I understand, for categorical variables and dataset $D$,
$\text{Likelihood}(D) = \prod_{i=1}^n \prod_{j=1}^{q_i} \prod_{k=1}^{r_i} \left( \frac{N_{ijk}}{N_{ij}} \right)^{N_{ijk}}$,
where $n$ is number of variables, $q_i$ is the number of configurations of the parents of the $i$th variable, $r_i$ is number of categories for $i$th variable. $N_{ijk}$ is the number of instances in the dataset $D$, where the $i$th variable takes on the $k$th value while the parents are set in their $j$th configuration. And $N_{ij}$ marginalise out over $k$, so it's the total number of total instances with the $j$th configuration.
My interpretation of this is, denoting $x_i$ the $i$th variable.
$ = \prod_{i=1}^n \prod_{j=1}^{q_i} \prod_{k=1}^{r_i} Pr(x_i = k | \text{parents of } x_i \text{ have } j\text{th configuration})^{\text{number of times occurs}}$ .
