calculations do not rely on it.

Typically, a locally optimal DAG is found by starting from a random graph and randomly ascending the likelihood by adding, modifying, or removing one directed edge at a time [@Friedman:1999; @Friedman:2000; @koller2009]. Alternatively, the posterior distribution $p({\mathcal{G}}\mid{\mathbf{X}})$ can be estimated with Bayesian inference using Markov Chain Monte Carlo (MCMC) simulation, allowing us to estimate the significance levels at an extra computational cost. The parameter values of $\alpha$, $\beta$, and $\sigma$, as part of ${\mathcal{G}}$, can be estimated with maximum likelihood.

When Bayesian network is modified by a single edge, only the vertices that receive a change would require a recalculation, whilst all others remain intact. This significantly reduces the amount of computation needed for each random step. A further speedup is achievable if we constrain the maximum number of parents each vertex can have, either by using the same fixed number for all nodes, or by pre-selecting a variable number of potential parents for each node using, for instance, a preliminary $L_1$-regularisation step [@schmidt2007learning].

Two DAGs are called Markov equivalent if they result in the same PDF [@koller2009]. Clearly, using gene expression data alone, Bayesian networks can only be resolved up to Markov equivalence. 