methodology for learning statistical dependency structures from possibly noisy data [@Friedman:1999; @Friedman:2000; @koller2009].

We adopt our previous convention in Section \[sec:genet-gene-expr\], where we have the gene expression data ${\mathbf{X}}$ and genetic markers ${\mathbf{G}}$. The model contains a total of $k$ vertices (i.e. random variables), $X_i$ with $i=1,\dots,k$, corresponding to the expression level of gene $i$. Given a DAG ${\mathcal{G}}$, and denoting the parental vertex set of $X_i$ by ${\mathbf{Pa}}^{({\mathcal{G}})}(X_i)$, the acyclic property of ${\mathcal{G}}$ allows to define the joint probability distribution function as $$\label{eq:2}
  p(X_1,\dots,X_k\mid{\mathcal{G}})
  =\prod_{i=1}^kp(X_i \mid {\mathbf{Pa}}^{({\mathcal{G}})}(X_i)).$$ In its simplest form, we model the conditional distributions as $$p\bigl(X_i\mid {\mathbf{Pa}}^{({\mathcal{G}})}(X_i)\bigl) = N\biggl(\alpha_i +
  \sum_{X_j\in{\mathbf{Pa}}^{({\mathcal{G}})}(X_i)} \beta_{ji}(X_j-\alpha_j),\sigma_i^2\biggr),$$ where $(\alpha_i,\sigma_i)$ and $\beta_{ji}$ are parameters for vertex $X_i$ and edge $X_j\rightarrow X_i$ respectively, as part of the DAG structure ${\mathcal{G}}$. Under such modelling, the Bayesian network is called a linear Gaussian network.

The likelihood of data ${\mathbf{X}}$ given the graph ${\mathcal{G}}$ is $$\begin{aligned}
  p({\mathbf{X}}\mid{\mathcal{G}}) =  \prod_{i=1}^k\prod_{l=1}^n 
  p(X_{il}\mid \{X_{jl}, X_j\in{\mathbf{Pa}}^{({\mathcal{G}})}(X_i)\}).\end{aligned}$$ Using Bayes’ rule, the log-likelihood of the DAG ${\mathcal{G}}$ based on the gene expression data ${\mathbf{X}}$ becomes [$$\log p({\mathcal{G}}\mid{\mathbf{X}})=\log
  p({\mathbf{X}}\mid {\mathcal{G}})+\log p({\mathcal{G}})-\log p({\mathbf{X}}),$$]{} where $p({\mathcal{G}})$ is the prior probability for ${\mathcal{G}}$, and $p({\mathbf{X}})$ is a constant when the expression data is provided, so the follow-up 