# a_gibbs_sampler_for_learning_dags__d0680ece.pdf

Journal of Machine Learning Research 17 (2016) 1-39 Submitted 11/14; Revised 11/15; Published 4/16

A Gibbs Sampler for Learning DAGs

Robert J. B. Goudie robert.goudie@mrc-bsu.cam.ac.uk Medical Research Council Biostatistics Unit Cambridge CB2 0SR, UK

Sach Mukherjee sach.mukherjee@dzne.de German Centre for Neurodegenerative Diseases (DZNE) Bonn 53175, Germany

Editor: Peter Spirtes

We propose a Gibbs sampler for structure learning in directed acyclic graph (DAG) models. The standard Markov chain Monte Carlo algorithms used for learning DAGs are random-walk Metropolis-Hastings samplers. These samplers are guaranteed to converge asymptotically but often mix slowly when exploring the large graph spaces that arise in structure learning. In each step, the sampler we propose draws entire sets of parents for multiple nodes from the appropriate conditional distribution. This provides an eﬃcient way to make large moves in graph space, permitting faster mixing whilst retaining asymptotic guarantees of convergence. The conditional distribution is related to variable selection with candidate parents playing the role of covariates or inputs. We empirically examine the performance of the sampler using several simulated and real data examples. The proposed method gives robust results in diverse settings, outperforming several existing Bayesian and frequentist methods. In addition, our empirical results shed some light on the relative merits of Bayesian and constraint-based methods for structure learning. Keywords: structure learning, DAGs, Bayesian networks, Gibbs sampling, Markov chain Monte Carlo, variable selection

1. Introduction

We consider structure learning for graphical models based on directed acyclic graphs (DAGs). The basic structure learning task can be stated simply: given data X on p variables, assumed to be generated from a graphical model based on an unknown DAG G, the goal is to make inferences concerning the edge set of G (the vertex set is treated as ﬁxed and known). Structure learning appears in a variety of applications (for an overview see Korb and Nicholson, 2011) and is a key subtask in many analyses involving graphical models, including causal inference. In a Bayesian framework, structure learning is based on the posterior distribution P(G | X) (Koller and Friedman, 2009). The domain of the distribution is the space G of all DAGs with p vertices. The size of the space grows super-exponentially with p, precluding exhaustive enumeration for all but the smallest problems. Markov chain Monte Carlo (MCMC) methods are widely used to sample DAGs and thereby approximate the posterior P(G | X). Available methods include MC3 (Madigan and York, 1995) and related samplers (Giudici and Castelo, 2003; Grzegorczyk and Husmeier, 2008) and algorithms that

c 2016 Robert J. B. Goudie and Sach Mukherjee.

Goudie and Mukherjee

sample from the space of total orders (Friedman and Koller, 2003). Order-space sampling entails some restrictions on graph priors (Eaton and Murphy, 2007; Ellis and Wong, 2008), because the number of total orders with which a DAG is consistent is not constant. Orderspace approaches have more recently led to exact methods based on dynamic programming (Koivisto and Sood, 2004; Parviainen and Koivisto, 2009; Tian and He, 2009; Tamada et al., 2011; Parviainen and Koivisto, 2013; Nikolova et al., 2013). These are currently feasible only in small domain settings (typically p < 20 but p < 30 is feasible on large cluster computers) and usually share the same restrictions on graph priors as order-space samplers. Frequentist constraint-based methods such as the PC-algorithm (Spirtes et al., 2000; Colombo and Maathuis, 2014) and the approach of Xie and Geng (2008) are an alternative. These methods make ﬁrm decisions about the DAG structure via a series of tests of conditional independence. In this paper we propose a Gibbs sampler for structure learning of DAGs that ameliorates key deﬁciencies in existing samplers. Random-walk Metropolis-Hastings algorithms currently used for structure learning propose small changes at each iteration, which can leave the algorithms unable to escape local modes. The Gibbs sampler proposed here considers the parents of a set of nodes as a single component (a so-called block ), sampling an entire parent set for each node in the block in one step. These large moves are sampled from the appropriate conditional posterior distribution. This enables the sampler to eﬃciently locate and explore areas of signiﬁcant posterior mass. The method is based on the simple heuristic that the parents of a node are similar to the covariates or inputs in variable selection, with the node as the output variable, but accounts for acyclicity exactly so that the equilibrium distribution is indeed the correct posterior distribution over DAGs. The sampler does not impose restrictions on priors or graph space beyond those (maximum in-degree, modular prior) common to most samplers for structure learning (e.g. Friedman and Koller, 2003; Ellis and Wong, 2008; Grzegorczyk and Husmeier, 2008). The maximum in-degree restriction formally precludes large-sample consistency in the general case, but facilitates eﬀective and robust inference by reducing the size of the model space (see Discussion).

2. Notation and Model

Let G denote a DAG with vertex set V (G) = {1, . . . , p}, and directed edge set E(G) V (G) V (G); often we will refer to vertex and edge sets simply as V and E respectively, leaving G implicit. The binary adjacency matrix corresponding to G is denoted AG, with entries speciﬁed as AG uv = 1 (u, v) E(G) and diagonal entries equal to zero. For inference, G is treated as a latent graph in the space G of all possible DAGs with p vertices. When a DAG is used to deﬁne a probabilistic graphical model (a Bayesian network), each vertex (or node; we use both terms interchangeably) is associated with a component of a p-dimensional random vector X = (X1 . . . Xp)T. XZ denotes the random variables (RVs) corresponding to variable indices Z {1, . . . , p}. We use bold type for the corresponding data; n samples are collected in the n p matrix X, with Xu denoting the column corresponding to variable u and XZ the submatrix corresponding to variable index set Z {1, . . . , p} (for notational simplicity we assume columns are ordered by increasing variable index).

Gibbs Sampler for Learning DAGs

0 1 1 0 0 0 1 0 0 0 0 0 0 0 1 0

pa G 2 = {1}

pa G 3 = {1, 2, 4} pa G 4 =

Figure 1: Illustration of notation. A graph G, with vertex set V = {1, 2, 3, 4}, edge set E = {(1, 2), (1, 3), (2, 3), (4, 3)} and adjacency matrix AG as shown, can be speciﬁed using parent sets as G = pa1 = , pa2 = {1}, pa3 = {1, 2, 4}, pa4 = or abbreviated to G = , {1}, {1, 2, 4}, (see text for description of notation). For the vertex subset Z = {2, 4}, pa G Z = ({1}, ) and pa G Z = ( , {1, 2, 4}) and thus the graph can be speciﬁed as G = pa Z = pa G Z, pa Z = pa G Z . We abbreviate this as G = pa G Z, pa G Z .

The proposed Gibbs sampler changes entire parent sets en masse, not just for one node but for a set of nodes, and so it will often be convenient and natural to specify a DAG G in terms of the parents of each node. Let pa G v = {u : (u, v) E(G)} denote the set of nodes that are parents of node v in G. A tuple of p parent sets pa1 = pa G 1 , . . . , pap = pa G p fully speciﬁes the edge set E, since u pa G v (u, v) E(G). Thus, structure learning can be viewed as inference of parent sets pav. The parent set pav takes values in the power set P(V \ {v}) of nodes excluding node v, subject to acyclicity. We will usually suppress the labels, and simply write pa G 1 , . . . , pa G p when specifying parent sets. Since pa G v is a subset of the variable indices, we can use Xpa G v to denote the corresponding data submatrix.

We denote the tuple of parent sets of a set of nodes Z = {z1, . . . , zs} V by pa G Z = (pa G z )z Z. This is a tuple of s components (for notational convenience we take these to be ordered by increasing vertex index), each of which is the parent set for the corresponding node in G. We wish to make inferences about pa Z, which takes values in P(V \{z1}) . . . P(V \ {zs}), subject to acyclicity. The tuple of parent sets for the complement ZC = V \ Z is denoted pa G Z = (pa G z )z ZC. Clearly, a tuple pa G Z1, . . . , pa G Zs of tuples of parent sets speciﬁes the parents of every node in a graph whenever {Z1, . . . , Zs} forms a partition of V ; the entire edge set can thus be speciﬁed by such a tuple. In particular, note that any graph G can be speciﬁed as pa G Z, pa G Z for any Z V . Some of the notation we use is illustrated in Figure 1. Our statistical formulation for Bayesian networks is standard. We brieﬂy summarize the main points and refer the reader to the references below for details. Given a DAG G, the joint distribution of X is p(X | G, {θv}) = Q v V p(Xv | Xpa G v , θv), i.e. the joint distribution factors over nodes, and each node is conditioned on its parents in G, parameterized by θv. For structural inference, interest focuses on the posterior distribution P(G | X). This is proportional to the product of the marginal likelihood p(X | G) and a graph prior π(G). Our sampler is compatible with essentially any speciﬁc model, but inherits computational costs associated with evaluation of the relevant quantities. In all examples we assume conjugate priors for θv, as well as local parameter independence and modularity; this leads to a closed-form marginal likelihood in both multinomial (Heckerman et al., 1995) and Gaussian

Goudie and Mukherjee

(Geiger and Heckerman, 1994) cases. Computations are simpliﬁed if the graph prior is modular (Friedman and Koller, 2003), meaning it factorises as π(G) = Q v V πv(pa G v ); we assume modular priors in all examples. Note that the prior is not speciﬁed over the space of orders and so is not subject to the restrictions this entails (Ellis and Wong, 2008; Eaton and Murphy, 2007). Under these assumptions, the posterior factorises across nodes as

P(G | X) = P(pa1 = pa G 1 , . . . , pap = pa G p | X)

v V p(Xv | Xpa G v )πv(pa G v ),

where p(Xv | Xpa G v ) is the marginal likelihood for node v given the graph G = pa G 1 , . . . , pa G p . The distribution P(G | X) is the target distribution for our sampler.

3. A Gibbs Sampler for Structure Learning

In this section we provide a high-level description of the Gibbs sampler for structure learning. We ﬁrst recall the standard random-walk Metropolis-Hastings sampler (known as MC3) and then describe a na ıve Gibbs sampler, which oﬀers no gains over MC3, but prepares the ground for the introduction of ideas from the Gibbs sampling literature. Speciﬁcally we discuss a strategy known as blocking that can improve mixing in Gibbs samplers and show how to use blocking for structure learning of DAGs. Several points relating to computation are central to developing a practical Gibbs sampler in this setting, but for clarity of exposition we defer discussion of computational aspects to Section 4.

3.1 MC3 Sampler

The standard sampler for structure learning of DAGs is MC3 (Madigan and York, 1995). This is a classical Metropolis-Hastings sampler with proposal G drawn uniformly at random from the set neigh(G) of DAGs that diﬀer from the current DAG G by the addition or removal of a single edge. The acceptance probability is min(1, r(G , G)), where

r(G , G) = min

1, P(G | X) |neigh(G )| 1

P(G | X) |neigh(G)| 1

Variants of MC3 include single-edge direction reversal proposals (Giudici and Castelo, 2003).

3.2 A Na ıve (and Ineﬃcient) Gibbs Sampler

Constructing a Gibbs sampler that is analogous to MC3 is straightforward. The posterior distribution on DAGs is a joint distribution over the oﬀ-diagonal entries in the adjacency matrix AG, i.e. over the p(p 1) binary RVs AG uv, u = v. At each iteration, MC3 can be thought of as proposing to toggle the value AG uv for some u = v, subject to the restriction that the resulting graph must be acyclic. A simple Gibbs sampler works in a similar way. At each step, a move from graph G to a new graph G is chosen by sampling from the conditional distribution of AG uv (for some u = v) given the rest of the graph. If (u, v) E(G), deﬁne the graph G(u,v) to be identical to G, and G (u,v) to be the graph that diﬀers from G only in lacking an edge from u to

Gibbs Sampler for Learning DAGs

v; conversely, if (u, v) / E(G), deﬁne G (u,v) to be identical to G, and G(u,v) to be the graph that diﬀers from G only in including an edge from u to v. If G(u,v) is cyclic, G (u,v)

is sampled as the new graph G with probability 1. If G(u,v) is acyclic, the conditional distribution of AG uv is Bernoulli, i.e.

P(AG uv = a | AG uv) =

P(G (u,v) | X) P(G (u,v) | X) + P(G(u,v) | X) a = 0,

P(G(u,v) | X) P(G (u,v) | X) + P(G(u,v) | X) a = 1.

The choice of u and v can either be made sequentially (systematically) or randomly (Roberts and Sahu, 1997); in this paper, random-scan Gibbs samplers are used throughout. We prove convergence of the Gibbs sampler to the target distribution in Appendix A.

3.3 Ineﬃciency in Gibbs Sampling

The mixing of Metropolis-Hastings samplers depends upon the proposal distribution, which for convenience is often chosen as a random-walk. In contrast, Gibbs samplers make moves according to conditional distributions that reﬂect local structure of the target distribution. Nonetheless, Gibbs sampling is not always eﬃcient. In particular, correlation between the components being sampled can lead to ineﬃciency. To see this, consider a Gibbs sampler for a multivariate continuous distribution with highly correlated components. At each step, a single component is sampled according to its conditional distribution, but since it is strongly correlated with other component(s), the conditional is concentrated on only a small part of the support. This means the sampler is likely to make only small moves. Analogous issues arise with discrete distributions.

For graphical models based on DAGs, there may be strong dependence between the edge indicators AG uv, particularly for the collections of RVs corresponding to parent sets. For example, there may be RVs Xu and Xv that in combination score highly as parents of node w, but not individually. Then, AG uw and AG vw will be correlated. In addition, the acyclicity restriction may induce strong dependence. For example, suppose two RVs Xu and Xv are strongly correlated and both edges (u, v) and its reverse (v, u) have high posterior probability. If the edge (u, v) is present, the probability of it being removed is low, but its presence precludes the reversed edge (v, u) from being added. This possibility motivates the edge reversal move used in variants of MC3 (Giudici and Castelo, 2003). Figure 2 shows a three node scenario. Here, both graphs (a) and (b) have high probability. Since reversing the edge (w, x) forms a cycle in (a), moves that consider the parents of only the pair w and x at the same time will not move between graphs (a) and (b) easily. Samplers that alter only a single edge indicator, such as MC3, will also fail. In contrast, sampling the parents of all three nodes jointly would make it easy to move between graphs (a) and (b). Decorrelating transformations could alleviate such problems, but in general ﬁnding a suitable transformation is diﬃcult, and for DAGs must of course encapsulate the requirement for acyclicity.

Goudie and Mukherjee

Figure 2: Illustrative scenario in which small local moves limit transitions between two regions of high probability. If both graphs (a) and (b) have high probability, the near-cyclic nature of the graphs makes transitions between (a) and (b) diﬃcult.

3.4 A Gibbs Sampler with Blocking

We address the deﬁciencies of the na ıve Gibbs sampler by grouping a number of the components together and sampling from their joint conditional distribution. In Gibbs sampling, this is known as blocking. Sampling from such a joint conditional can ameliorate diﬃculties caused by correlations between components because the joint conditional naturally accounts for the correlation structure. In the multivariate normal case, Roberts and Sahu (1997) have shown that blocking improves convergence for random-scan Gibbs sampling.

In principle, any group of components can be taken as a block, but for the algorithm to be useful in practice, sampling from a block s joint conditional distribution must be feasible, and ideally simple. The blocks that we consider correspond to groups (speciﬁcally tuples) of parent sets, so that the parent sets of several nodes are considered simultaneously. It is natural that each block is a tuple of parent sets because, as described in Section 2, we can specify a graph G using a tuple pa G 1 , . . . , pa G p of parent sets. It is therefore convenient to describe the sampler using this alternative graph speciﬁcation, in which G = pa G 1 , . . . , pa G p .

We denote the set of q nodes whose parent sets will be sampled together as a block by W = {w1, . . . , wq} V . Suppose the current graph is G = pa G W , pa G W (recall that any graph can be written in this way with respect to any partition of the node set). A move to a new graph G = pa G W , pa G W is formed by changing the parents of the nodes in W from pa G W = (pa G w1, . . . , pa G wq) to pa G W = (pa G w1, . . . , pa G wq) and setting pa G W = pa G W (i.e. leaving the parents of nodes not in W unchanged). We sample the tuple pa G W of parent sets jointly, conditional on the tuple pa G W of parent sets of nodes not in W (that remain unchanged). In terms of the adjacency matrix AG, each block consists of the indicators specifying the parents of the nodes in W, i.e. AG vw for v V, w W, v = w.

To construct a Gibbs sampler using these blocks, we need to ﬁnd the conditional posterior distribution on the tuple pa W of parent sets, given that the remaining tuple pa W of parent sets is set to pa G W . The conditional distribution depends on whether the graph G formed using the proposed parent sets is acyclic. We therefore introduce Pa G W to denote the set of permissible tuples, i.e. tuples pa G W of parent sets such that G = pa G W , pa G W is acyclic. For tuples pa G W / Pa G W , the conditional probability is 0. For pa G W Pa G W , the

Gibbs Sampler for Learning DAGs

Algorithm 1 A Gibbs sampler for learning DAGs, with blocks

Initialise starting point G0 = pa G0 1 , . . . , pa G0 p for t in 1 to N do

Sample q nodes uniformly at random from V , and call this set of nodes W Sample pa G W from P(pa W = pa G W | pa Gt 1 W , X)

Set Gt G = pa G W , pa Gt 1 W end for

conditional probability is

P(pa W = pa G W | pa W = pa G W , X) = P( pa G W , pa G W | X) P(pa G W | X)

= P(G | X) P

pa G W Pa G W P(pa G W , pa G W | X). (1)

This is the conditional distribution needed to specify the blocked Gibbs sampler; Algorithm 1 outlines the procedure at a high level, with the set W of nodes chosen uniformly at random at each step. Asymptotic convergence follows from the argument in Appendix A (in fact the requirements on the graph prior will be weaker than in the na ıve case).

3.5 Tuning Parameters

To reduce the computational costs of structure learning it is common to set a maximum in-degree (e.g. Friedman and Koller, 2003; Grzegorczyk and Husmeier, 2008). We set a maximum in-degree κ = 3 in all empirical examples, except where stated otherwise (Section 5). This facilitates sampling from the conditional distribution in Equation 1 by reducing the computational cost of evaluating the normalising constant. We set the block size q = |W| = 3. While q can be chosen freely in principle, evaluating the conditional distribution in Equation 1 becomes unmanageable when q is large. Thus, both κ and q act as tuning parameters (and are in addition to any hyper-parameters in the Bayesian formulation).

3.6 Structure Learning and Variable Selection

It is interesting to note that when q = 1 (and thus W = w for some w V ) if no choice of parent set induces a cycle then Pa G W = P(V \ {w}), the power set of the remaining nodes. In this case, the conditional distribution in Equation 1 can be viewed as the posterior distribution of a variable selection problem with response or output variable w, and the other variables as covariates or inputs. If the addition of particular nodes introduces a cycle, we have a constrained problem. In particular, suppose adding any of the nodes Z V as parents of node w would create a cycle. Then Pa G W = P(V \ ({w} Z)) and the conditional distribution in Equation 1 is a constrained variable selection problem in which the nodes Z are excluded from the set of candidate inputs.

Goudie and Mukherjee

4. Computational Aspects

Designing a computationally eﬃcient sampler is not straightforward in this setting. To see why, note that to sample from the conditional in Equation 1 we need to be able to identify the set Pa G W of tuples of parent sets that is permissible (in the sense of maintaining acyclicity). The cardinality of this set is typically large, and the interdependence between the parent sets of each node in W makes decomposing the problem into subproblems nontrivial. A simple but na ıve approach would list all possible parent sets for each node in W and check each such combination for cyclicity, but this approach is slow and cumbersome. We propose a partitioning scheme that leads to a two-stage sampler. The key idea is to choose the partition of Pa G W so that, conditional on a component of the partition, the parents of each node are independent. This enables an eﬃcient two-stage sampling method that we describe in Section 4.3. Acyclicity constraints are met by an eﬃcient dynamic algorithm.

4.1 A Partition on Permissible Tuples of Parent Sets

We partition the set Pa G W of permissible tuples of parent sets as Pa G W = {Pa G,H1 W , . . . , Pa G,Hη W }. Each component Pa G,Hh W is a set of permissible tuples of parent sets for nodes in W. It is convenient to label the partition components using secondary (unrelated to G) DAGs Hh H, where H is the space of DAGs with q = |W| vertices; η is the cardinality of H. We describe the relationship between each Pa G,Hh W and DAG Hh using the following elements, illustrated in Figure 3 and Table 1:

The reduced graph G = pa G 1 , . . . , pa G p , which is a function of the graph G, and is identical to G except that edges directed into nodes in W are removed.

pa G w w / W

The (reﬂexive) transitive closure, which is the directed graph TG on nodes V with edges ET, where (u, v) ET if and only if u = v or a (directed) path from u to v exists in G (i.e. there exists a sequence of nodes z1, . . . , zs V (G), with z1 = u and zs = v, such that (z1, z2), (z2, z3), . . . , (zs 1, zs) E(G)).

The descendant nodes de G w = {v V : TG wv = 1} and the non-descendant nodes nd G w = {v V : TG wv = 0} of w W in the reduced graph G. Note that w de G w and w / nd G w by deﬁnition.

The nodes nd G = T w W nd G w that are not descendants in the reduced graph G of any node in W.

The nodes de G,H w = S x pa H w de G x that are descendants in the reduced graph G of any node x pa H w , for a given node w W. Each node x pa H w is a parent node of the node w in the graph H = pa H w1, . . . , pa H wq .

Gibbs Sampler for Learning DAGs

Original graph G

Reduced graph G

An example H

Figure 3: An illustrative example of the relevant graphs and sets, with W = {3, 4, 7} shown in red. The edges into W are removed from the original graph G, to form G. Here H is the set of all DAGs on the nodes {3, 4, 7} and thus η = |H| = 25. For H as shown, nd G = {1, 2}.

de G w nd G w de G,H w de G,H w QG,H w Condition (B)

w = 3 {3, 5, 6} {1, 2, 4, 7, 8} {3, 4, 5, 6, 7, 8} {1, 2} No restriction

w = 4 {4, 6} {1, 2, 3, 5, 7, 8} {3, 4, 5, 6, 7, 8} {1, 2} No restriction

w = 7 {7, 8} {1, 2, 3, 4, 5, 6} {4, 6} {3, 5, 6, 7, 8} {1, 2, 4} pa G 7 {4} =

Table 1: The relevant sets and requirements of conditions (A) and (B) for the illustrative graph G shown in Figure 3, with H as shown. Recall nd G = {1, 2}. Condition (B) depends on RG,H 7,4 = {4}, but it imposes no restriction when w {3, 4} because

pa H 3 = pa H 4 = . We ﬁnd that Pa G,H W contains all tuples of parent sets for which pa G 3 {1, 2}, pa G 4 {1, 2} and pa G 7 = {4} Z where Z {1, 2}.

The nodes de G,H w = S x W\pa H w de G x that are descendants in the reduced graph G of any node x W \ pa H w , for a given node w W. Each node x W \ pa H w is not a parent node of w in the graph H = pa H w1, . . . , pa H wq .

Using this notation, we deﬁne Pa G,H w , for given G G, W V , w W and H = pa H w1, . . . , pa H wq H, as the set of parent sets pa G w that satisfy the following conditions.

(A) pa G w QG,H w = nd G de G,H w \ de G,H w

(B) pa G w RG,H w,x = for all nodes x pa H w , with RG,H w,x = de G x \ de G,H w .

Note that (B) depends on de G x not de G,H x . We deﬁne Pa G,H W as the set of tuples formed by the Cartesian product of the sets Pa G,H w of parent sets for w W; in other words pa G W = (pa G w1, . . . , pa G wq) Pa G,H W if and only if pa G w1 Pa G,H w1 , . . . , pa G wq Pa G,H wq .

Lemma 1 {Pa G,H1 W , . . . , Pa G,Hη W } is a partition of Pa G W .

Proof See Appendix B.

Goudie and Mukherjee

Condition (A) ensures all graphs G = pa G W , pa G W , with parents pa G W Pa G,H W , are acyclic. It requires that all parents of w W are either not descendants in the reduced graph G of any node in W, or a node whose ancestors in the reduced graph G are all parents in the graph H of node w. In particular, no descendant of w is added as an ancestor of w. Condition (B) ensures each DAG G is a member of Pa G,H W for only a single H H for any given G G and thus that {Pa G,H1 W , . . . , Pa G,Hη W } is a partition of the set of permissible tuples of parent sets. The condition checks that each edge in the graph H is used i.e. that pa G W = (pa G w1, . . . , pa G wq) Pa G,H W would not be allowed under any H = H, H H. Speciﬁcally, it ensures that each parent set pa G w contains at least one descendant node v V of each parent pa H w of node w in the graph H, and that v is not a descendant in G of any node x W that is not a parent in H of w W. To see why this condition is required, consider a graph G in which all nodes in W have no parents. Then pa G w = for all nodes w W and so condition (A) holds for all H H. Thus if condition (B) is disregarded, pa G W Pa G,H W for all H H, implying that {Pa G,H1 W , . . . , Pa G,Hη W } is not a partition of Pa G W .

4.2 Fast Identiﬁcation of Partition Components

Parent sets in each set Pa G,H w can be identiﬁed easily using simple set operations, and consequently the partition components Pa G,H W can be easily identiﬁed via Cartesian products of these sets. Let Paall w = P(V \ {w}) denote the set of all possible parent sets of node w, subject to maximum in-degree κ, and Paall w (v) = Paall w \ P(V \ {w, v}) denote the subset of Paall w containing only parent sets that contain node v V . Parent sets in Pa G,H w must (A) not include any nodes not in QG,H w , and (B) include at least one node in RG,H w,x for each x pa H w , and thus

( Pa(B) w \ Pa(A) w if pa H w =

Paall w \ Pa(A) w if pa H w = ,

where Pa(A) w = S

v / QG,H w Paall w (v) and Pa(B) w = T x pa H w S

r RG,H w,x Paall w (r).

We ﬁx q = |W| as a small constant for all p and set a maximum in-degree κ. This means permissible tuples of parents sets can be identiﬁed in O(pκ+1) time by storing the lookup tables Paall w (v) as bit maps, assuming that O(1) querying of the descendants and non-descendants is available. This can be achieved using a dynamic algorithm that provides access to the transitive closure. Giudici and Castelo (2003) previously used a similar approach. For a graph with transitive closure TG with edges ET, the descendants and non-descendants of node u are {v : (u, v) ET} and {v : (u, v) / ET} respectively. Thus, the descendants and nondescendants of a node can be identiﬁed in O(1) time when the transitive closure is known. The transitive closure for an arbitrary directed graph with p vertices can be determined in O(pω) time (Munro, 1971), where ω is the best known exponent for matrix multiplication (Coppersmith and Winograd, 1990, show ω < 2.376). However, in the present setting only incremental changes are made to the current state G of the sampler, so it is not necessary to use an oﬄine algorithm at each iteration. Instead, we can use a fully dynamic transitive closure algorithm (Demetrescu et al., 2010) that provides procedures both for querying the

Gibbs Sampler for Learning DAGs

transitive closure, and for incrementally updating it when an edge is added or removed from the graph. We choose to implement the algorithm introduced by King and Sagert (2002), which performs queries in O(1) time, and updates in O(p2) worst-case time (see reference for details of required assumptions). A trade-oﬀexists between the performance of these two operations, but this bound for updates is thought to be the best possible whilst retaining O(1) queries (Demetrescu and Italiano, 2005) and the algorithm is simple to implement. The algorithm maintains a p p path count matrix CG whose elements CG uv are the number of distinct paths from node u to node v in the graph G for u = v, and CG uv = 1 for u = v. The transitive closure can be derived from the path count matrix by noting that TG uv = 1 if and only if CG uv > 0. Thus query operations are performed in O(1) by simply checking whether the relevant component of CG is positive. When an edge is added or removed CG is updated as follows. First consider adding an edge (u, v) to a graph G to form a new graph G . The increase in the number of distinct paths between any two nodes x and y is given by the (x, y) element of CG u CG v , where CG u denotes the uth column of CG, CG v denotes the vth row and denotes the outer product. Thus, the path count matrix is updated simply as CG = CG +CG u CG v . Similarly, when an edge (u, v) is removed from the graph the update is CG = CG CG u CG v .

4.3 Two-stage Sampling Method

We draw a new graph G = pa G W , pa G W starting from a graph G using a two-stage method:

we sample ﬁrst a component Pa G,H W of the partition of permissible tuples of parent sets and then a tuple pa G W of parent sets from the selected partition component. In the ﬁrst stage, Pa G,H W is drawn from the conditional distribution, given the tuple of parent sets of nodes not in W.

P(Pa G,H W | pa G W , X) = P(Pa G,H W , pa G W | X) P(pa G W | X)

pa G W Pa G,H W Q w W p(Xw | Xpa G w )πw(pa G w ) P H H P

pa G W Pa G,H W Q w W p(Xw | Xpa G w )πw(pa G w )

pa G w Pa G,H w p(Xw | Xpa G w )πw(pa G w ) P H H Q w W P

pa G w Pa G,H w p(Xw | Xpa G w )πw(pa G w )

The ﬁnal equality follows by an interchange of sum and product that is proved in Lemma 2. This makes evaluation more eﬃcient by allowing the sums to be evaluated separately for each node. Friedman and Koller (2003) used a similar interchange.

Lemma 2 The following identity holds for any H H, W V and G G. X

pa W Pa G,H W

w W p(Xw | Xpaw)πw(paw) = Y

paw Pa G,H w

p(Xw | Xpaw)πw(paw)

Proof See Appendix C.

Goudie and Mukherjee

Algorithm 2 A Gibbs sampler for learning DAGs, with general blocks

Initialise starting point G0 = pa G0 1 , . . . , pa G0 p Compute initial path count matrix CG0

for t in 1 to N do

Sample q nodes uniformly at random from V , and call this set of nodes W Let G = Gt 1 Form G as deﬁned in Section 4.1 Evaluate CG from CG as described in Section 4.2 for w W do

Evaluate de G w = {v : CG wv 1} Evaluate nd G w = {v : CG wv = 0} end for for H H do

Evaluate Pa G,H w as described in Section 4.2 Let KH w = P

pa G w Pa G,H w p(Xw | Xpa G w )πw(pa G w ) end for Let KH = Q w W KH w end for Sample Pa G,H W according to P(Pa G,H W | pa G W , X) = KH P H H KH for w W do

Sample pa G w according to P(pa G w | Pa G,H W , pa G W , X) end for Set Gt G = pa G W , pa Gt 1 W Update CGt

In the second stage, we sample new parents pa G W from the selected partition component, and form the new graph G = pa G W , pa G W . The parents of each node w W, conditional on H , are independent, and so can be sampled separately from the following conditional distribution:

P(pa G w | Pa G,H W , pa G W , X) = p(Xw | Xpa G w )πw(pa G w ) P

pa G w Pa G,H w p(Xw | Xpa G w )πw(pa G w ).

This step is straightforward because this distribution is simply the posterior distribution of a constrained variable selection with response w and Pa G,H w as the set of possible active sets (i.e. selected covariates).

Algorithm 2 outlines the complete algorithm. The methods in Sections 4.2 enable fast identiﬁcation of each partition component. Run-time is a function of p, maximum in-degree κ, and the number q of nodes in the block. We choose a small, ﬁxed q for all p, so the run-time is determined by the evaluation of Pa G,H w , which is O(pκ+1).

Gibbs Sampler for Learning DAGs

In this section, we empirically assess the performance of the proposed sampler, comparing it with existing samplers as well as frequentist methods for structure learning of DAGs.

5.1 Evaluation Setup

We compare the Gibbs sampler with the MC3 sampler (Madigan and York, 1995) and the REV sampler (Grzegorczyk and Husmeier, 2008), a variant of MC3 that uses a more extensive edge reversal move. We also compare with two frequentist constraint-based methods: the PC-algorithm (Spirtes et al., 2000), and the Xie and Geng (2008) method, that is shown by its authors to outperform the PC-algorithm in some settings. Tuning parameters for each method were set as follows. For the constraint-based methods, the signiﬁcance level was α = 0.05 by default, but we also show some results for α = 0.00005, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.1. The Gibbs sampler we use is a randomscan sampler, with q = 3 (i.e. the parent sets of three nodes are sampled jointly at each iteration). To permit a fair comparison, for MC3 we use the same fast online transitive closure algorithm (Section 4.2), and pre-computation and caching of local marginal likelihoods used in the Gibbs sampler (REV uses a similar pre-computation and caching scheme). We constrain all of the samplers to maximum in-degree κ = 3 and set the graph prior as π(G) 1. We use conjugate formulations throughout, speciﬁcally multinomial-Dirichlet (Heckerman et al., 1995) for discrete data and Normal with a g-prior (with g = n) (Geiger and Heckerman, 1994; Zellner, 1986) for continuous data. We consider six examples: a small domain example, where comparison with the exact posterior (Tian and He, 2009) is possible; data simulated from the ALARM network and from randomly-generated networks of varying sparsity; data sets from social science and biology; and a pathological 4-node example designed to highlight a failure case for our method. In practical use samplers can only be run for a ﬁnite number of iterations (depending on available time and computational resources). We set the maximum number of iterations as follows. In total, we drew 106 iterations of REV (retaining only every 10th iteration to reduce storage requirements). Following Grzegorczyk and Husmeier (2008), 85% of proposals within REV were MC3 proposals (without MC3 proposals, the REV sampler is not irreducible). In our implementation, the computational costs of the Gibbs sampler are an order of magnitude lower than REV s (accounting for MC3 moves), but we nevertheless treated the computational costs as the same for the purposes of comparison and drew 106 iterations of the Gibbs sampler (again retaining every 10th iteration). The computational costs of our implementation of MC3 are roughly 1/10 of a Gibbs iteration, and so we performed 107

iterations of MC3 (retaining every 100th). For each sampler, 10 independent runs starting from diﬀerent initial graphs were performed. We discarded the ﬁrst 1/4 of samples, but as an additional ﬁlter we considered trace plots (log marginal likelihood versus iteration) for each run of each sampler and discarded further samples if they appeared to be from a pre-convergent phase. Speciﬁcally, if there was a clear step change in the trace for a certain run we discarded all samples in that run drawn prior to the highest step being reached. We note that this simple ﬁlter cannot be regarded as detecting convergence because the highest step may correspond to a region

Goudie and Mukherjee

near a local rather than a global mode. We therefore also study convergence in detail via other metrics.

5.2 Evaluation Metrics

Each sampler gives an estimated posterior distribution over DAGs. From this we obtain a point estimate either as the maximum a posteriori (MAP) graph GMAP, or as a thesholded graph Gτ formed by including all directed edges whose marginal posterior edge probabilities are at least τ (note that Gτ need not be acyclic). When τ = 0.5 we get the median probability model Gmed (for normal linear models Barbieri and Berger, 2004, show that in some settings this is the optimal model for prediction). We also consider thresholding to match sparsity levels seen in frequentist point estimates, i.e. setting τ such that |E(Gτ)| = |E( ˆG)| for a point estimate ˆG. By GPC τ and GXie τ we denote thresholded graphs whose sparsity is matched to the PC and Xie-Geng estimates respectively. We use the structural Hamming distance (SHD) to quantify diﬀerences between DAGs. This is the minimum number of edge insertions and deletions needed to transform one graph into another. Comparisons are made using completed partially directed acyclic graphs (Chickering, 2002a). To show the behaviour of the samplers at shorter MCMC run lengths some metrics are shown against iteration t. These are based on the ﬁnal 3/4 of samples drawn up to t (except for log marginal likelihoods), and so may incorporate some samples drawn before convergence. We assess convergence and stability via the following metrics:

Trace plots of log marginal likelihood against iteration t.

Between-run agreement between posterior edge probabilities. These are visualized as scatter plots. When edge probabilities agree between a pair of runs, all points in the plot will lie close to the y = x line. We show two variants: a hexagonally-binned version, to avoid over-plotting (Carr et al., 1987); and a panelled plot, in which each panel shows one pair of runs. We also consider the number of edges with estimated posterior probability greater than 0.9 in one run, and less than 0.1 in another run. Such major discrepancies represent serious Monte Carlo artefacts.

Potential scale reduction factor (PSRF) is a multiple-chain convergence diagnostic, developed by Gelman and Rubin (1992), that compares withinand between-chain means and variances of a suitable summary statistic. PSRF < 1.1 is often taken to indicate that a sampler has run suﬃciently long. The natural summary statistics in this context are the posterior edge probabilities. However, we ﬁnd that the resulting PSRF is not a reliable indicator of convergence. In simulations (not shown) using a 20 node network for which we could calculate true posterior edge probabilities (Tian and He, 2009), we found little association between single-edge PSRF and the (absolute) error in posterior edge probability. In particular, the diagnostic appeared to be very conservative for edges with true posterior probability near 0 or 1. We therefore suggest that a reasonable use of PSRF in this context may be to assess relative convergence between diﬀerent samplers, rather than as an absolute indicator of convergence. Below, we compare samplers in terms of the proportion of edges with PSRF < 1.1 (a larger proportion suggests better mixing). To calculate PSRF we consider the 10 runs

Gibbs Sampler for Learning DAGs

of each sampler as a collection of 5 pairs of runs and calculate PSRF separately for each edge using the ﬁnal three quarters of the samples drawn up to that point.

For the real data, we assess stability under resampling ( shaking the data ) by comparing estimates across bootstrap samples.

For experiments in which the true data-generating graph G is known we use the following metrics to assess accuracy:

Structural Hamming distance (SHD) between G and the estimated graphs.

Receiver-operating characteristic (ROC) curves. These show agreement between G and an estimate Gτ by plotting true positive against false positive rates parameterized by threshold τ. We consider also the area under the ROC curve (AUROC), focusing in particular on the small false positive rate region of the curve that is often of interest in applications.

Finally, when the posterior edge probabilities can be computed exactly (Tian and He, 2009) (i.e. in small p settings) we also consider the maximum and average absolute error in posterior edge probability, calculated across the set of all possible edges. We note that while SHD and ROC scores are useful in assessing performance, they are not convergence measures, as they do not assess accuracy of the posterior distribution per se (to see why, consider a degenerate sampler that samples only G , giving perfect scores on SHD and ROC, but an incorrect posterior distribution).

5.3 Small Domain Comparison to Exact Posterior

We applied the methods to the Zoo data set (Newman et al., 1998) that records p = 17 (discrete) characteristics of n = 101 animals. The maximum log marginal likelihood found by the Gibbs and REV samplers was consistent across runs, but less so for the MC3 sampler, and the Gibbs sampler reached a plateau of high probability after the fewest iterations (Figure A13a). REV required about ten times as many iterations as Gibbs to achieve the same proportion of edges with PSRF < 1.1 while MC3 needed about 100 times as many (Figure A14a). The estimated edge probabilities given by the Gibbs sampler were stable between runs (Figure A15a). The results from the REV sampler were also stable, but MC3

less so, with major disparities between some runs. Figure 4 shows the error in posterior edge probabilities as a function of iterations. Convergence was quickest for the Gibbs sampler, followed by REV and then MC3. All runs of the Gibbs sampler reached a point at which the maximum absolute error was 0.05 (after 67,000 iterations on average). In contrast, only 5/10 runs of REV and 6/10 runs of MC3 reached the same level of accuracy in their complete run. Similarly, the Gibbs sampler achieved an average absolute error of less than 0.01 in the fewest iterations.

5.4 The ALARM Network

The ALARM network (Beinlich et al., 1989) consists of 37 discrete nodes and 46 edges and has been widely used in studying structure learning (e.g. Friedman and Koller, 2003; Grzegorczyk and Husmeier, 2008). We simulated data sets from ALARM with sample sizes

Goudie and Mukherjee

Figure 4: For the Zoo data set, maximum and average (across all possible 272 edges) error in posterior edge probabilities versus iteration number (on a log10 scale). For each MCMC algorithm, 10 independent runs, initialised at disparate initial values, are shown.

n = 100, 500, 1000, 2500, 5000. Figure 5 shows trace plots for the n = 1000 case as an illustrative example. The maximum log marginal likelihood found by the Gibbs sampler across runs was 10,608.20 0.8 (mean standard deviation), whereas for REV it was 10,627.2 40.4 and for MC3 11,176.9 273.7. The highest scoring graph discovered by any of the 10 runs of the MC3 sampler had log marginal likelihood 10,856.6, far below the Gibbs maximum, suggesting non-convergence in all 10 runs. REV appeared to reach convergence in all but two runs (runs 9 and 10). The Gibbs sampler consistently (and rapidly) reached a high scoring plateau and appeared to have converged in all 10 runs. PSRF results were in line with these ﬁndings (Figure A14b); more than 10 times as many iterations of REV were needed to give the same proportion of edges with PSRF < 1.1 as the Gibbs sampler. Figure 6 compares pairs of runs for n = 1000; this pair of runs is typical of all 10 runs (except for runs 9 and 10 of REV which disagree considerably; all runs shown in Figure A15b). There were no major between-run discrepancies (as deﬁned in Section 5.2) for the Gibbs sampler at any sample size. The mean number of major discrepancies (across pairs of runs) increased from 0 (n = 100) to 8 and 91 for MC3 and REV respectively (when n = 5000). Figure 7 shows ROC curves for false positive rates < 0.05. The Gibbs sampler performs better and with less variability than the other methods. Sample size is inﬂuential. For small n the Bayesian methods outperform the constraint-based methods. However, counter to the increase in statistical information, REV and MC3 perform less well with larger n: e.g. at n = 100 for a false positive rate (FPR) of 0, the Gibbs sampler found 21.9 1.9 (mean standard deviation) true edges; REV 20.1 1.4; and MC3 21.7 2.1. But at n = 5000, for the same FPR, Gibbs found 43.0 0.0 true edges; REV 13.2 21.3; and MC3 did not ﬁnd any true edges (the true graph has 46 edges). The constraint-based methods performed well for large sample sizes, as anticipated by the asymptotic consistency of the PC-algorithm (Kalisch and B uhlmann, 2007). The Xie-Geng method performed particularly well for n = 5000. Figure 8 shows AUROC as a function of iterations; the Gibbs sampler performed

Gibbs Sampler for Learning DAGs

Figure 5: Log marginal likelihood of the graphs visited by the three MCMC samplers in 10 independent runs, initialised at disparate initial conditions, on the ALARM data, with data sample size n = 1000. Iteration number is displayed on a log10 scale. The dashed lines indicate where the burn-in phase ended. In all cases the log marginal likelihood of the graphs reached by the MC3 sampler are below the displayed range.

Figure 6: Convergence diagnostic plots for each MCMC sampler for the ALARM data, with n = 1000. The posterior edge probabilities given by two independent MCMC runs are plotted against each other, binned into hexagonal areas to avoid over-plotting. When the edge probabilities of the two runs agree, all of the points in the plot will lie close to the y = x line; strongly oﬀ-diagonal points indicate extreme discrepancies between the runs.

Goudie and Mukherjee

Figure 7: ALARM data, receiver operating characteristic (ROC) curves for each of 10 independent MCMC runs for each MCMC sampler, and point estimates for the constraint-based methods. Note that the horizontal axis shows only false positives rates < 0.05, corresponding to the case in which interest focuses on highranked edges (the complete curves are shown in Figures A16). Point estimates from Xie-Geng s constraint-based method and the PC-algorithm are also shown for all 8 signiﬁcance levels considered.

best after all run lengths and sample sizes considered. These results are supported by SHDs between point estimates and the true graph (Table A2).

5.5 Networks of Varying Sparseness

Sparsity can be used to statistical and computational advantage but in practice it may be hard to know what level of sparsity is reasonable for a given application. We therefore sought to investigate the eﬀect of varying network sparsity, including scenarios where the data-generating graph can violate the in-degree constraint we impose. We simulated data following a procedure described in Kalisch and B uhlmann (2007) that we outline below. We ﬁrst generated a DAG G via its adjacency matrix AG, by drawing entries as AG uv Bernoulli(ρ), where ρ (0, 1) is a parameter controlling sparsity. Entries were drawn independently for each u < v, with AG uv = 0 otherwise. Larger ρ gives a denser graph and the expected number of neighbours (parents or immediate children) of a node is ρ(p 1). We set p = 25, n = 1000 and considered 5 values of the sparseness parameter ρ (corresponding to expected neighbourhood sizes 2, 3, 4, 5 and 6). For each ρ we drew 10 DAGs, simulating a data set from each DAG by ancestral sampling using a Normal linear model (see Kalisch and B uhlmann, 2007, for full details).

Gibbs Sampler for Learning DAGs

Figure 8: ALARM data, area under the receiver operating characteristic curve (AUROC) against iteration number (log10 scale) for each of the 10 independent runs for each MCMC algorithm. Each panel shows results for a particular sample size n = 100, . . . , 5000.

Figure 9 shows boxplots over SHDs. We see that accuracy decreases with increasing density, echoing results in Kalisch and B uhlmann (2007) for the PC algorithm. As throughout, all the samplers had an in-degree constraint (maximum in-degree κ = 3), while the frequentist methods did not. Nevertheless, the performance of the Bayesian methods did not appear to deteriorate any more rapidly than the frequentist methods. At all ρ s, the median probability graph Gmed outperformed the frequentist methods which in turn outperformed the MAP graph. In this context, we draw attention to a diﬀerence between the median probability and MAP graphs: the former, although obtained by averaging over DAGs satisfying the in-degree constraint, may itself have in-degree greater than κ, while the latter is necessarily subject to the constraint.

5.6 Survey Data

The publicly available Behavioral Risk Factor Surveillance System Survey (BRFSS) (Centers for Disease Control and Prevention, 2008) is a household-level random-digit telephone survey, collected by the U.S. National Center for Chronic Disease Prevention and Health, that has been conducted throughout the United States since 1984. We considered (discrete) responses to p = 24 questions (see Appendix D), spanning most of the topics covered in the survey. We considered the responses from New York in the 2008 survey, removing samples for respondents who refused or were unsure of their response, or whose response was

Goudie and Mukherjee

Figure 9: Synthetic data, accuracy of estimation for networks with diﬀerent levels of sparsity. The panels correspond to simulations in which the expected number of neighbours for each node in the data-generating graph is respectively 2, 3, 4, 5, and 6. Accuracy is quantiﬁed by structural Hamming distance (SHD) between estimates and data generating graphs (smaller SHDs correspond to more accurate estimates). Box plots are over the 10 independent networks/data sets simulated for each sparsity level. For MCMC methods, results from the MAP graph GMAP

and median probability model Gmed are shown.

missing, to any of the 24 questions. The resulting sample size was n = 4, 197. The median probability graph Gmed estimated by the Gibbs sampler is shown in Figure A17a. Figure 10 shows between-run agreement. The Gibbs runs showed better agreement than the REV runs and the MC3 runs disagreed considerably (these pairs of runs were typical; all runs shown in Figure A18a). Indeed there were no major between-run discrepancies (in the sense of Section 5.2) for the Gibbs sampler, whereas there were on average 2.4 major discrepancies for REV and 10.9 for MC3. The Gibbs sampler also had the highest proportion of edges with PSRF < 1.1 (Figure A14c). The maximum log marginal likelihoods found by each of the three samplers diﬀered considerably as did the number of iterations needed to reach a plateau (Figure A13b). The Gibbs sampler typically reached a plateau after around 500 samples (although in one run 103

samples were needed). REV took longer to reach a (usually lower) plateau. MC3 appeared to become stuck in a region with even lower log marginal likelihood. To investigate stability under resampling of the data we applied the methods to 10 bootstrap resamples of the data set. The Gibbs and REV samplers were more stable than the frequentist methods as well as than MC3. For example, using GPC τ as a point estimate

Gibbs Sampler for Learning DAGs

Figure 10: Convergence diagnostic plots for each MCMC sampler for the survey data. The posterior edge probabilities given by two independent MCMC runs are plotted against each other, binned into hexagonal areas to avoid over-plotting. When the two runs give the same estimates of the posterior edge probabilities all of the points appear on the line y = x; strongly oﬀ-diagonal points indicate extreme discrepancies between the runs. This pair of runs was typical of all pairs.

for the Bayesian methods (i.e. thresholding to give the same number of edges as PC), the mean SHD between results from pairs of bootstrap data sets was 21.7 (Gibbs), 22.3 (REV), 33.4 (MC3) and 30.8 (PC-algorithm). Results for GMAP and GXie τ are shown in Figure A19a.

5.7 Large-sample, Single-cell Molecular Data

We used single-cell molecular data from Bendall et al. (2011) with p = 34 continuous variables (see Appendix D) and n = 21, 691. The median probability graph Gmed estimated by the Gibbs sampler is shown in Figure A17b. Figure 11 shows between-run agreement for the samplers for a typical pair of runs; the Gibbs sampler shows better agreement than REV or MC3. All but one of the 10 runs of the Gibbs sampler were consistent with each other (Figure A18b) and the one inconsistent run nonetheless showed better agreement with the other Gibbs runs than any pair of runs of the other samplers. The Gibbs sampler had no major discrepancies (as deﬁned in Section 5.2) between any pairs of runs, while there were on average 26 and 87 major discrepancies for REV and MC3 respectively. Smaller discrepancies followed a similar pattern: the average number of edges diﬀering by 0.1 or more in posterior probability between runs were 6.4 (Gibbs sampler), 56.6 (REV), and 151.8 (MC3). The Gibbs sampler also had the highest proportion of edges with PSRF < 1.1 after any run length (Figure A14d). Trace plots are shown in Figure A13c. The Gibbs sampler found a region of higher log marginal likelihood than the other samplers and did so consistently. Indeed, the maximum log marginal likelihood achieved across all REV runs was exceeded in every Gibbs run and after only around 5000 iterations. Finally, we considered stability under resampling of the data (Figure A19b), including also the frequentist point estimators for comparison. The mean SHD between PC estimates from pairs of bootstrap samples was 214.6; the corresponding mean SHDs for the samplers (using GPC τ point estimates) were 128.8 (Gibbs), 140.1 (REV) and 225.2 (MC3). We note that the SHDs are particularly high here because the PC-algorithm estimates a

Goudie and Mukherjee

Figure 11: Convergence diagnostic plots for each MCMC sampler for the single-cell molecular data. The posterior edge probabilities given by two independent MCMC runs are plotted against each other, binned into hexagonal areas to avoid overplotting. When the two runs give the same estimates of the posterior edge probabilities all of the points appear on the line y = x; strongly oﬀ-diagonal points indicate extreme discrepancies between the runs. This pair of runs was typical of all pairs of runs.

network with a high density (recall that under GPC τ all methods choose the same number of edges as PC).

5.8 A Pathological 4-node Example

Our ﬁnal example demonstrates the potential sensitivity of the Gibbs sampler to choice of q = |W| (the number of nodes whose parent sets are sampled together in a single iteration). The example was constructed to highlight a situation in which the sampler can become stuck in a local mode unless q is large enough, potentially leading to extremely slow convergence. We used as data 105 simulated samples from a 4-node network in which the parents of node 2 were nodes 1 and 4; the parents of node 3 were nodes 1 and 2; and nodes 1 and 4 had no parents. Each node was Bernoulli distributed. The probability of success was 0.6 for nodes 1 and 4, while nodes 2 and 3 were noisy XORs with probability of success 0.9 when either parent (but not both) was true and 0.1 otherwise. For the purposes of demonstration, we set q = 2. As shown in Figure 12, after 106

iterations the maximum error in edge probability for the Gibbs sampler was 0.016 0.015 (mean standard deviation across 10 runs). Given that there are only 543 DAGs with p = 4 nodes, this magnitude of error is unexpectedly large. The slow convergence is due to the concentration of posterior mass on two graphs that the sampler cannot easily move between. The MAP graph (posterior probability 0.61) is the graph in which the parents of node 2 are nodes 1 and 3; the parents of node 4 are nodes 1 and 2; and nodes 1 and 3 have no parents. The data-generating graph has posterior probability 0.38. The graphs diﬀer in the parents of 3 nodes, and so with q = 2 the sampler cannot move between them in a single step. At the same time, the large sample size leads to all other graphs having low posterior probability (< 0.002), making multi-step transitions between the two graphs unlikely. Thus, the sampler becomes stuck on one of the two graphs. In this example,

Gibbs Sampler for Learning DAGs

Figure 12: For the pathological 4-node example, maximum and average (across all possible 12 edges) error in posterior edge probabilities by iteration number (on a log10 scale). For each MCMC algorithm, 10 independent runs, initialised at disparate initial values, are shown.

setting q = 3 improved convergence (maximum error 0.0011 0.0009 after 106 iterations), but in real-world examples it is diﬃcult to rule out the possibility that analogous issues may arise.

6. Discussion

We introduced a Gibbs sampler for structure learning of DAGs that can converge more easily than existing samplers due to its ability to eﬃciently make large moves in DAG space. We showed that it provides often substantial gains in accuracy and stability in comparison with existing (Bayesian and frequentist) methods in a range of settings. The formulation of the sampler develops and exploits the connection between variable selection and structure learning. This connection has been widely studied for undirected graphs (e.g. Meinshausen and B uhlmann, 2006), but for DAGs is complicated by the acyclicity requirement. Our approach accounts for acyclicity but further work in this area may ease the adaptation of results and methods from Bayesian variable selection to the case of DAGs. In the proposed sampler, existing variable selection methods could be of direct utility in sampling from the conditional distribution P(pa W | pa G W , X). When this sampling step is diﬃcult, a Metropolis-within-Gibbs approach (i.e. substituting a Metropolis step in place of the Gibbs step) could be considered. With q = |W| = 1, the conditional P(pa W | pa G W , X) is identical to the posterior of the corresponding variable selection problem. Then, a Metropolis-within-Gibbs move could directly exploit existing variable selection methods. In common with most of the structure learning literature, we used an in-degree constraint. This gives a smaller DAG space and in addition controls the number of parameters needed to specify conditional distributions. However, it is a strong constraint and a natural concern is whether it excludes higher in-degree models that could be appropriate for the data. This possibility cannot be ruled out, but the use of model averaging and marginal posterior summaries may ameliorate the concern to some extent, because even when the

Goudie and Mukherjee

true number of parents of a node exceeds the restriction, important candidate parents may still have high edge scores (and appear in a summary graph). For example, suppose the true in-degree of a node is 4, but at most 3 incoming edges are permitted. In this case, although no model including all 4 parents can be considered, provided the signal can be detected considering only 3 nodes at a time, the posterior probability of all 4 edges may nonetheless be high.

In our empirical examples the REV sampler of Grzegorczyk and Husmeier (2008) showed impressive gains over MC3. The Gibbs sampler generally seemed to outperform both. The use of conditional distributions in REV is a point of similarity with the Gibbs sampler proposed here. The performance gains of Gibbs could be explained by two key diﬀerences from REV. First, REV does not use the natural conditional distribution and requires an accept-reject step. Second, REV must include at least some MC3 proposals (otherwise the sampler is not irreducible), and these steps are not tailored to the shape of the posterior distribution (Grzegorczyk and Husmeier, 2008, make REV proposals with probability 1/15 and so most steps are based on MC3 proposals).

If there is strong correlation between the parent sets of more than q = |W| nodes, the Gibbs sampler may not mix well. In this situation, constraint-based methods may be useful. Alternatively, q could be chosen at each step according to some distribution, so that a mixture of diﬀerent block sizes is used. This would in particular allow larger blocks to be used without increasing the computational demands of the algorithm excessively. In this paper we ﬁxed q = 3, and found this simple choice gave a well-behaved and eﬀective sampler. But there is a trade-oﬀ: increasing q increases the time taken to evaluate P(pa W | pa G W , X), but also increases move size, with the potential to improve convergence.

Practical use of the Gibbs sampler in the form described here requires exact sampling from the conditional P(pa W | pa G W , X) and there are situations related to this requirement in which other methods may be more suitable. First, when an appropriate maximum indegree cannot be set, MC3 or a variant could be more appropriate (although convergence could be very slow). Alternatively, search procedures such as GES (Chickering, 2002b) or HCMC (Castelo and Koˇcka, 2003) could be used. Second, exact sampling from the conditional distribution is challenging in settings with thousands of nodes. In this case, eﬃcient constraint-based methods (such as the PC-algorithm) may be a better choice, particularly in the large sample setting. As noted by a referee, an interesting area for future research would be combining the Gibbs sampler with some aspects of other methods such as PC that are relatively well suited to the truly high-dimensional setting.

Acknowledgments

The authors are grateful to Karen Sachs for providing the single-cell molecular data; to Marco Grzegorczyk and Dirk Husmeier for their inspiring work in this area and for providing their implementation of the REV sampler; to the referees for numerous suggestions that improved the article; and to Frances Griﬃths, David Lunn and Paul Newcombe for helpful discussions. Part of this work was carried out while the authors were at the University of Warwick whose excellent scientiﬁc environment the authors warmly acknowledge and

Gibbs Sampler for Learning DAGs

supported by the Economic and Social Research Council (ESRC) and Engineering and Physical Sciences Research Council (EPSRC).

Appendix A. Convergence Conditions for Gibbs Samplers

Convergence of a Gibbs sampler for DAGs does not follow from the usual justiﬁcation of Gibbs sampling that relies on the Hammersley-Cliﬀord theorem (Besag, 1974). The theorem gives a positivity condition that is suﬃcient to prove that the univariate conditional distributions used by the Gibbs sampler uniquely deﬁne the joint distribution. The required condition is that the support of the joint distribution is given by the Cartesian product of the supports of the marginal distributions. An example of when this condition does not hold is the joint density p(x, y) for a pair of random variables X and Y with support only on [0, 1] [0, 1] and [2, 3] [2, 3]. Clearly p(x) and p(y) are both positive on [0, 1] and [2, 3] but neither [0, 1] [2, 3] or [2, 3] [0, 1] are in the support of the joint distribution (Hobert et al., 1997).

In the DAG setting, the acyclicity requirement means that this positivity condition is not satisﬁed. Consider a DAG consisting of two correlated random variables X1 and X2. The correlation means that both the graph with a single edge (1, 2) and the graph with a single edge (2, 1) have positive posterior probability. Thus P(AG 12 = 1 | X) > 0 and P(AG 21 = 1 | X) > 0 in the posterior marginal distributions. However, in the joint posterior distribution P(AG 12 = 1, AG 21 = 1 | X) = 0 because the corresponding graph (the complete graph) is cyclic. The complete graph is thus not in the support of the joint distribution but is in the Cartesian product of the supports of the marginal distributions.

An alternative suﬃcient condition for uniqueness of the joint distribution and convergence of the Gibbs sampler when positivity is not satisﬁed is given by Besag (1994), which was extended for continuous settings by Hobert et al. (1997). In the present context, the condition requires that for every pair G , G G of DAGs with p nodes there exists a ﬁnite sequence G1, . . . , GN, with G1 = G , GN = G and N N, and such that Gt and Gt 1 diﬀer in only a single component (in this context, a single edge), and that the joint posterior distribution P(Gt | X) > 0 for all t = 0, . . . N. When the graph prior π(G) > 0 for all G G, this condition is clearly satisﬁed: as an example, one such ﬁnite sequence removes every edge of G , one at a time, and then adds every edge of G , one at a time. Each graph in the sequence is clearly acyclic, since the sequence is composed of subgraphs of the acyclic G and G , and so has positive probability in the joint distribution when the graph prior is positive everywhere in G. A similar proof follows if the graph prior has support on all subgraphs of graphs with support in the graph prior, as is true for most widely used priors.

Appendix B. Proof of Lemma 1

Lemma 1 {Pa G,H1 W , . . . Pa G,Hη W } is a partition of Pa G W .

Proof. We show: (i) S h=1,...η Pa G,Hh W Pa G W ; (ii) S h=1,...η Pa G,Hh W Pa G W ; and (iii)

Pa G,Hh1 W Pa G,Hh2 W = for Hh1 = Hh2, Hh1, Hh2 H.

Goudie and Mukherjee

(i) We prove the relationship by showing that Pa G,Hh W Pa G W for each h = 1, . . . , η. By the deﬁnition of Pa G W in Section 3.4 we need to show that G = pa G W , pa G W is acyclic for all tuples pa G W Pa G,Hh W of parent sets, for all h = 1, . . . , η.

We proceed by contradiction. Suppose some graph G = pa G W , pa G W with pa G W Pa G,Hh W is cyclic. First note that G is acyclic because G is a subgraph of the acyclic G. Since the graph G diﬀers from the acyclic graph G only in the parents of nodes in W, each cycle in G must include at least one node in W. Denote by x1 x2 the existence of a path (that obeys the edge directions) in the graph G from node x1 W to node x2 W that does not include any nodes in W (except x1 and x2). Let w1, . . . , wr W be the (minimal) complete set of nodes in W included in some cycle in G . Without loss of generality suppose that w1 w2 . . . wr in G . Note that since w1, . . . , wr is the complete set of nodes in W in the cycle, no node between wi and wi+1 in the path can be in W, i {1, . . . , r 1}.

We now show that for x1, x2 W, x1 x2 only if an edge (x1, x2) links node x1 to x2 in the graph Hh. Since x1 x2, there must exist a node v V that is a parent of x2 in G such that v is a descendant in G of x1. Note that if the edge (x1, x2) is in the graph G then v = x1. Since v is a parent of node x2 in the graph G and x2 W, then v nd G de G,Hh x2 \ de G,Hh x2 by condition (A). Also since x1 x2 does

not include any nodes in W, v is also a descendant in G of x1 because G and G diﬀer only in which nodes are parents of nodes in W. We proceed by contradiction. Suppose no edge (x1, x2) exists in Hh. Then v is a descendant of x1 in the graph G, but x1 is not a parent of x2 in the graph Hh. So v de G,Hh x2 , which by condition (A) is a contradiction. Thus x1 x2 only if (x1, x2) is present in Hh, for x1, x2 W.

Now, recall that w1 w2 . . . wr. Since a cycle is formed we must in addition have a path in G from node wr to node w1. Since w1, . . . , wr is the complete set of nodes in W involved in the cycle, no node on the path from wr to w1 can be in W. Thus wr w1. However, this implies that the edges (w1, w2), (w2, w3), . . . , (wr 1, wr), (wr, w1) are all in Hh, which implies Hh is cyclic. But Hh is acyclic by assumption, and so we have a contradiction.

(ii) Suppose we start from a graph G = pa G W , pa G W . We want to show that for each

DAG G = pa G W , pa G W there is some H H such that pa G W Pa G,H W . We will show

that pa G W Pa G,H W , where H = pa H 1 , . . . , pa H q H is a DAG on nodes in W and

pa H w = {x W : there exists some v pa G w such that v de G x }. As usual, pa G w is the parent set in the graph G of the node w; and de G x is the descendants in G of the node x. Note that H is a subgraph (on the nodes in W) of the transitive closure TG . By deﬁnition, G is a DAG, so TG is also a DAG, and thus H is a DAG.

We show that pa G W Pa G,H W by showing that for each node w W, both conditions (A) and (B) that specify membership of Pa G,H w are satisﬁed.

Gibbs Sampler for Learning DAGs

For (A) we need pa G w nd G de G,H w \ de G,H w . Let v pa G w meaning that v is a

parent of the node w in the graph G . First we show that v / de G,H w , and then show that v nd G de G,H w .

To show that v / de G,H w , note that if v de G,H w then v must be a descendant in G of some node y W that is not in pa H w . However, every such y is in pa H w by the deﬁnition of pa H w , thus v / de G,H w .

To show that v nd G de G,H w , we suppose v / nd G and show this implies that v de G,H w . This follows because if v / nd G then it must be the descendant in G of some node y W. Then y pa H w by deﬁnition of H . Therefore v de G,H w , as required.

For (B), we need pa G w de G x \ de G,H w = for all x pa H w . Consider x pa H w . By

the deﬁnition of pa H w , this means that there exists some node v pa G w such that v de G x . Additionally, any y W for which v de G y is such that y pa H w , by deﬁnition of pa H w . Thus v is not a descendant in G of any node in W that is not in pa H w . Then v pa G w de G x \ de G,H w , and thus the condition is satisﬁed.

(iii) Since Hh1 = Hh2, at least one node has diﬀerent parents in Hh1 and Hh2. Suppose that the node w W is such a node, and that x W is a parent of w in Hh1 but not in Hh2. Consider a parent set pa G1 w Pa G,Hh1 w . We will prove that pa G1 w / Pa G,Hh2 w , and the result follows.

By condition (B), there must exist some v de G x \ de

G,Hh1 w such that v pa G1 w . We

will show that v / pa G2 w for every parent set pa G2 w Pa G,Hh2 w . This follows because v de G x and so is a descendant in G of x, which is not a parent of w in Hh2. Thus

v de Hh2 w , and so v / nd G de

G,Hh2 w . Therefore v cannot be a parent of

Appendix C. Proof of Lemma 2

Lemma 2 The following identity holds for any H H, W V and G G. X

pa W Pa G,H W

w W p(Xw | Xpaw)πw(paw) = Y

paw Pa G,H w

p(Xw | Xpaw)πw(paw)

Proof. Deﬁne λ(i) w = p(Xw | Xpa(i) w )πw(pa(i) w ), i {1, . . . P}, where P is the cardinality of

Pa G,H W , and where pa(i) w is the parent set of node w for the ith member of Pa G,H W i.e. we

have that Pa G,H W = n (pa(1) w1, . . . , pa(1) wq ), . . . , (pa(P) w1 , . . . , pa(P) wq ) o . Similarly, deﬁne Pa G,H w =

Goudie and Mukherjee

n pa(1) w , . . . , pa(Pw) w o , and thus Pw is the number of tuples of parent sets in Pa G,H w . Recall

that Pa G,H W is the Cartesian product of the sets Pa G,H w of parent sets for w W, thus

pa(i) w Pa G,H w

λ(1) w + + λ(Pw) w

i1 {1,...,P1}, ..., iq {1,...,Pq} λ(i1) w1 . . . λ(iq) wq

(pa(i) w1,...,pa(i) wq ) Pa G,H W

λ(i) w1 . . . λ(i) wq

(pa(i) w1,...,pa(i) wq ) Pa G,H W

w W λ(i) w .

Appendix D. Details of Data Analysed

We included in our analysis the following variables from the survey data (Centers for Disease Control and Prevention, 2008): SEX, AGE G, RACEGR2, MARITAL, CHLDCNT, INCOMG, USEEQUIP, HCVU65, MEDCOST, SMOKER3, ASTHMST, RFDRHV3, RFBING4, QLREST2, RFSEAT3,

TOTINDA, BMI4CAT, DIABETE2, EMTSUPRT, LSATISFY, EXTETH2, AIDTST2, DENVST1, IMONTH. In our analysis of the single-cell molecular data (Bendall et al., 2011) we included the following quantities, including the binding of antibodies, viability and DNA content: 191-DNA, 193-DNA, 103-Viability, 115-CD45, 139-CD45RA, 141-p PLCgamma2, 142-CD19, 144-CD11b, 145-CD4, 146-CD8, 148-CD34, 150-p STAT5, 147-CD20, 152-Ki67, 154-p SHP2, 151-p ERK1/2, 153-p MAPKAPK2, 156-p MAPKAPK2AP70/Syk, 158-CD33, 160-CD123, 159-p STAT3, 164-p SLP-76, 165-p NFk B, 166-Ik Balpha, 167-CD38, 168-p H3, 170-CD90, 169-p P38, 171-p Btk/Itk, 172-p S6, 174-p Src FK, 176-p CREB, 175-p Crk L, 110 114-CD3.

Gibbs Sampler for Learning DAGs

Appendix E. Additional Table

Method n = 100 250 500 1000 2500 5000 Gibbs 69.2 6.5 32.0 3.7 25.8 2.8 22.1 3.2 14.5 3.3 8.7 2.4 23.8 0.4 8.3 0.5 8.9 0.3 15.1 0.6 6.0 0.0 10.1 0.3 MC3 67.4 4.9 35.8 6.6 44.7 7.8 47.8 13.3 63.7 11.2 78.5 16.5 23.9 0.9 12.5 3.1 27.3 8.8 43.8 14.2 65.4 11.6 82.3 17.0 REV 68.2 5.5 34.7 5.8 27.4 4.6 28.1 4.0 24.4 9.5 20.4 6.7 25.9 2.1 8.3 1.1 4.0 0.0 11.5 2.0 11.6 9.9 13.5 8.0 PC 48.0 39.0 38.0 28.0 22.0 14.0 Xie-Geng 54.0 40.0 34.0 43.0 32.0 17.0

Table A2: ALARM data, structural Hamming distances (SHD) between estimated graphs and the true data-generating graph. For Bayesian methods we compare to both the maximum a posteriori (top line) and the median probability graph (bottom line), and report mean SHD over 10 independent Monte Carlo runs, along with the corresponding standard deviation. For constraint-based methods, we report results for α = 0.05.

Goudie and Mukherjee

Appendix F. Additional Figures

Figure A13: Log marginal likelihood of the graphs visited by the three MCMC samplers in 10 independent runs, initialised at disparate initial conditions. In (A) Zoo data; (B) survey data; (C) single-cell molecular data. Iteration number is displayed on a log10 scale. The dashed lines indicate where the burn-in phase ended.

Gibbs Sampler for Learning DAGs

Figure A14: Proportion of edges with PSRF < 1.1 against iteration number. (A) Zoo data; (B) ALARM data, n = 1000; (C) survey data; (D) single-cell molecular data.

Goudie and Mukherjee

Figure A15: Convergence diagnostic plots for all pairs of 10 independent MCMC runs for each sampler for (A) Zoo data and (B) ALARM data, n = 1000. In each cell, the posterior edge probabilities given by two independent runs are plotted against each other. Each point represents a single edge. The lower half of both panels compares runs of the Gibbs sampler; the upper half compares runs of the MC3 and the REV sampler respectively. When the two runs give the same estimates of the posterior edge probabilities, all of the points appear on the line y = x. The blue to orange colour scale represents the distance from this line, with orange points the furthest away.

Gibbs Sampler for Learning DAGs

Figure A16: ALARM data, receiver operating characteristic (ROC) curves for each of 10 independent MCMC runs for each MCMC sampler, and point estimates for the constraint-based methods. Point estimates from Xie-Geng s constraint-based method and the PC-algorithm are also shown for all 8 signiﬁcance levels.

Goudie and Mukherjee

Figure A17: The median probability model Gmed estimated by the Gibbs sampler, for (A) survey data and (B) single-cell molecular data. Darker shading indicates higher posterior edge probability. Note that no hard constraints were speciﬁed to ensure, for example, an in-degree of 0 for Age Group in (A); such constraints were omitted to keep the implementations of the various methods simple.

Gibbs Sampler for Learning DAGs

Figure A18: Convergence diagnostic plots for all pairs of 10 independent MCMC runs for each sampler for (A) survey data and (B) single-cell molecular data. In each cell, the posterior edge probabilities given by two independent runs are plotted against each other. Each point represents a single edge. The lower half of both panels compares runs of the Gibbs sampler; the upper half compares runs of the MC3 and the REV sampler respectively. When the two runs give the same estimates of the posterior edge probabilities, all of the points appear on the line y = x. The blue to orange colour scale represents the distance from this line, with orange points the furthest away.

Goudie and Mukherjee

Figure A19: Stability of estimators under bootstrapping, as quantiﬁed by structural Hamming distance (SHD) between estimates obtained from pairs of bootstrap iterations for (A) survey data and (B) single-cell molecular data. In each case 10 bootstrap data sets were drawn and each estimator was run on each bootstrap data set. Smaller SHDs indicate stable estimation in the sense of graph estimates that are robust to resampling. Boxplots are shown over all pairs of bootstrap iterations. The edge probabilities were thresholded so that the resulting graphs had the same number of edges as the graph indicated in the panel title: GPC τ (the graph estimated by the PC-algorithm), GXie τ (Xie-Geng method), or GMAP (MAP graph).

Gibbs Sampler for Learning DAGs

M. M. Barbieri and J. O. Berger. Optimal predictive model selection. Annals of Statistics, 32(3):870 897, 2004.

I. A. Beinlich, H. J. Suermondt, R. M. Chavez, and G. F. Cooper. The ALARM monitoring system: A case study with two probabilistic inference techniques for belief networks. In Second European Conference on Artiﬁcial Intelligence in Medicine, pages 247 256. Springer-Verlag, Berlin, 1989.

S. C. Bendall, E. F. Simonds, P. Qiu, E. D. Amir, P. O. Krutzik, R. Finck, R. V. Bruggner, R. Melamed, A. Trejo, O. I. Ornatsky, R. S. Balderas, S. K. Plevritis, K. Sachs, D. Pe er, S. D. Tanner, and G. P. Nolan. Single-cell mass cytometry of diﬀerential immune and drug responses across a human hematopoietic continuum. Science, 332(6030):687 696, 2011.

J. E. Besag. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B (Methodological), 36(2):192 236, 1974.

J. E. Besag. Discussion of Markov chains for exploring posterior distributions . Annals of Statistics, 22(4):1734 1741, 1994.

D. B. Carr, R. J. Littleﬁeld, W. L. Nicholson, and J. S. Littleﬁeld. Scatterplot matrix techniques for large N. Journal of the American Statistical Association, 82(398):424 436, 1987.

R. Castelo and T. Koˇcka. On inclusion-driven learning of Bayesian networks. Journal of Machine Learning Research, 4:527 574, 2003.

Centers for Disease Control and Prevention. Behavioral Risk Factor Surveillance System Survey Data. U.S. Department of Health and Human Services, Atlanta, Georgia, 2008.

D. M. Chickering. Learning equivalence classes of Bayesian-network structures. Journal of Machine Learning Research, 2:445 498, 2002a.

D. M. Chickering. Optimal structure identiﬁcation with greedy search. Journal of Machine Learning Research, 3:507 554, 2002b.

D. Colombo and M. H. Maathuis. Order-independent constraint-based causal structure learning. Journal of Machine Learning Research, 15:3921 3962, 2014.

D. Coppersmith and S. Winograd. Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation, 9(3):251 280, 1990.

C. Demetrescu and G. F. Italiano. Trade-oﬀs for fully dynamic transitive closure on DAGs: breaking through the O(n2) barrier. Journal of the ACM, 52(2):147 156, 2005.

C. Demetrescu, D. Eppstein, Z. Galil, and G. F. Italiano. Dynamic graph algorithms. In M. J. Atallah and M. Blanton, editors, Algorithms and Theory of Computation Handbook: General Concepts and Techniques, pages 9.1 9.28. Chapman and Hall/CRC, Boca Raton, FL, 2010.

Goudie and Mukherjee

D. Eaton and K. Murphy. Bayesian structure learning using dynamic programming and MCMC. In Proceedings of the Twenty-Third Conference on Uncertainty in Artiﬁcial Intelligence, pages 101 108, Corvallis, Oregon, 2007. AUAI Press.

B. Ellis and W. H. Wong. Learning causal Bayesian network structures from experimental data. Journal of the American Statistical Association, 103(482):778 789, 2008.

N. Friedman and D. Koller. Being Bayesian about network structure. A Bayesian approach to structure discovery in Bayesian networks. Machine Learning, 50(1-2):95 125, 2003.

D. Geiger and D. Heckerman. Learning Gaussian networks. In Proceedings of the Tenth Conference on Uncertainty in Artiﬁcial Intelligence, pages 235 240, San Francisco, CA, 1994. Morgan Kaufmann.

A. Gelman and D. B. Rubin. Inference from iterative simulation using multiple sequences. Statistical Science, 7(4):457 472, 1992.

P. Giudici and R. Castelo. Improving Markov chain Monte Carlo model search for data mining. Machine Learning, 50(1-2):127 158, 2003.

M. Grzegorczyk and D. Husmeier. Improving the structure MCMC sampler for Bayesian networks by introducing a new edge reversal move. Machine Learning, 71(2-3):265 305, 2008.

D. Heckerman, D. Geiger, and D. M. Chickering. Learning Bayesian networks: the combination of knowledge and statistical data. Machine Learning, 20(3):197 243, 1995.

J. P. Hobert, C. P. Robert, and C. Goutis. Connectedness conditions for the convergence of the Gibbs sampler. Statistics & Probability Letters, 33(3):235 240, 1997.

M. Kalisch and P. B uhlmann. Estimating high-dimensional directed acyclic graphs with the PC-Algorithm. Journal of Machine Learning Research, 8:613 636, 2007.

V. King and G. Sagert. A fully dynamic algorithm for maintaining the transitive closure. Journal of Computer and System Sciences, 65(1):150 167, 2002.

M. Koivisto and K. Sood. Exact Bayesian structure discovery in Bayesian networks. Journal of Machine Learning Research, 5:549 573, 2004.

D. Koller and N. Friedman. Probabilistic Graphical Models: Principles and Techniques. MIT Press, Cambridge, MA, 2009.

K. B. Korb and A. E. Nicholson. Bayesian Artiﬁcial Intelligence. Chapman & Hall/CRC Press, Boca Raton, FL, 2011.

D. Madigan and J. C. York. Bayesian graphical models for discrete data. International Statistical Review / Revue Internationale de Statistique, 63(2):215 232, 1995.

N. Meinshausen and P. B uhlmann. High-dimensional graphs and variable selection with the Lasso. Annals of Statistics, 34(3):1436 1462, 2006.

Gibbs Sampler for Learning DAGs

I. Munro. Eﬃcient determination of the transitive closure of a directed graph. Information Processing Letters, 1(2):56 58, 1971.

D. Newman, S. Hettich, C. Blake, and C. Merz. UCI Repository of Machine Learning Databases. University of California, Department of Information and Computer Science, Irvine, CA, 1998. URL http://www.ics.uci.edu/~mlearn/MLRepository.html.

O. Nikolova, J. Zola, and S. Aluru. Parallel globally optimal structure learning of Bayesian networks. Journal of Parallel and Distributed Computing, 73(8):1039 1048, 2013.

P. Parviainen and M. Koivisto. Exact structure discovery in Bayesian networks with less space. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artiﬁcial Intelligence, pages 436 443, Corvallis, Oregon, 2009. AUAI Press.

P. Parviainen and M. Koivisto. Finding optimal Bayesian networks using precedence constraints. Journal of Machine Learning Research, 14:1387 1415, 2013.

G. O. Roberts and S. K. Sahu. Updating schemes, correlation structure, blocking and parameterization for the Gibbs sampler. Journal of the Royal Statistical Society: Series B (Methodological), 59(2):291 317, 1997.

P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. The MIT Press, Cambridge, MA, 2000.

Y. Tamada, S. Imoto, and S. Miyano. Parallel algorithm for learning optimal Bayesian network structure. Journal of Machine Learning Research, 12:2437 2459, 2011.

J. Tian and R. He. Computing posterior probabilities of structural features in Bayesian networks. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artiﬁcial Intelligence, pages 538 547, Corvallis, Oregon, 2009. AUAI Press.

X. Xie and Z. Geng. A recursive method for structural learning of directed acyclic graphs. Journal of Machine Learning Research, 9:459 483, 2008.

A. Zellner. On assessing prior distributions and Bayesian regression analysis with g-prior distributions. In P. K. Goel and A. Zellner, editors, Bayesian Inference and Decision Techniques: Essays in Honour of Bruno de Finetti, pages 233 243. North-Holland, Amsterdam, 1986.