title: Learning All Credible Bayesian Network Structures for Model Averaging
authors: Zhenyu A. Liao Charupriya Sharma James Cussens Peter van Beek
text: 
  
  
  
  
  
  
  
  
  Editor
  
  Learning All Credible Bayesian Network Structures for Model Averaging
      Zhenyu A. Liao z6liao@uwaterloo.ca
   
  Charupriya Sharma c9sharma@uwaterloo.ca
   
  David R. Cheriton School of Computer Science
   
  University of Waterloo
   
  Waterloo, ON N2L 3G1, Canada
  James Cussens james.cussens@bristol.ac.uk
  
  Department of Computer Science
  
  University of Bristol
  
  Bristol, BS8 1QU, United Kingdom
  Peter van Beek vanbeek@uwaterloo.ca
  
  David R. Cheriton School of Computer Science
   
  University of Waterloo
   
  Waterloo, ON N2L 3G1, Canada
      December 30, 2023
  ====================================================================================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  A Bayesian network is a widely used probabilistic graphical model with
  applications in knowledge discovery and prediction.  Learning a
  Bayesian network (BN) from data can be cast as an optimization
  problem using the well-known score-and-search approach.
  However, selecting a single model (i.e., the best scoring BN)
  can be misleading or may not achieve the best possible accuracy.
  An alternative to committing to a single model is to perform
  some form of Bayesian or frequentist model averaging, where
  the space of possible BNs is sampled or enumerated in some
  fashion. Unfortunately, existing approaches for model averaging
  either severely restrict the structure of the Bayesian network
  or have only been shown to scale to networks with fewer than 30 random
  variables. In this paper, we propose a novel approach to model
  averaging inspired by performance guarantees in approximation
  algorithms. Our approach has two primary advantages. First,
  our approach only considers credible models in that they
  are optimal or near-optimal in score. Second, our approach
  is more efficient and scales to significantly larger
  Bayesian networks than existing approaches.
  
  
  
  
  
  
  
    Bayesian Networks, Structure Learning, Bayes Factor, Unsupervised Learning
  
  
  
  
  
  § INTRODUCTION
  
  
  
  A Bayesian network is a widely used probabilistic graphical
  model with applications in knowledge discovery, explanation,
  and prediction <cit.>.
  A Bayesian network (BN) can be learned from data
  using the well-known score-and-search approach, where a
  scoring function is used to evaluate the fit of a proposed BN
  to the data, and the space of directed acyclic graphs (DAGs)
  is searched for the best-scoring BN.
  However, selecting a single model (i.e., the best-scoring BN)
  may not always be the best choice.
  When one is using BNs for knowledge discovery and explanation
  with limited data, selecting a single model may be
  misleading as there may be many other BNs that have scores that
  are very close to optimal and the posterior probability of
  even the best-scoring BN is often close to zero. As well, when one is using
  BNs for prediction, selecting a single model may not achieve
  the best possible accuracy. 
  
  An alternative to committing to a single model is to
  perform some form of Bayesian or frequentist model averaging  <cit.>.
  In the context of knowledge discovery, Bayesian model averaging
  allows one to estimate, for example, the posterior probability
  that an edge is present, rather than just knowing whether the
  edge is present in the best-scoring network.
  Previous work has proposed
  Bayesian and frequentist model averaging approaches to network structure
  learning that enumerate the space of all possible DAGs <cit.>, sample from the space of all
  possible DAGs <cit.>, consider
  the space of all DAGs consistent with a given ordering of
  the random variables <cit.>,
  consider the space of tree-structured or other restricted
  DAGs <cit.>, and consider
  only the k-best scoring DAGs for some given value of k <cit.>.
  Unfortunately, these existing approaches either severely
  restrict the structure of the Bayesian network, such as
  only allowing tree-structured networks or only considering a
  single ordering, or have only been shown to scale to small
  Bayesian networks with fewer than 30 random variables.
  
  In this paper, we propose a novel approach to
  model averaging for BN structure learning that is inspired
  by performance guarantees in approximation algorithms. Let
   be the score of the optimal BN and assume 
  without loss of generality that the optimization problem
  is to find the minimum-score BN.
  Instead of
  finding the k-best networks for some fixed value of k,
  we propose to find all Bayesian networks 𝒢 that
  are within a factor ρ of optimal; i.e.,
  
      ≤(𝒢) ≤ρ·,
  
  for some given value of ρ≥ 1, or equivalently,
  
      ≤(𝒢) ≤+ϵ,
  
  for ϵ = (ρ - 1) ·. Instead of choosing arbitrary 
  values for ϵ, ϵ≥ 0, we show that for the 
  two scoring functions BIC/MDL and BDeu, a good choice for the value of ϵ is 
  closely related to the Bayes factor, a model selection criterion 
  summarized in <cit.>.
  
  Our approach has two primary advantages. First, our approach
  only considers credible  models in that they are
  optimal or near-optimal in score. Approaches that enumerate
  or sample from the space of all possible models consider DAGs
  with scores that can be far from optimal; for example,
  for the BIC/MDL scoring function the ratio of worst-scoring
  to best-scoring network can be four or five orders of
  magnitude[<cit.>
  deem such models discredited when they make a similar
  argument for not considering models whose probability is
  greater than a factor from the most probable.]. A
  similar but more restricted case can be made against the
  approach which finds the k-best networks since there is no 
  a priori way to know how to set the parameter k such that only
  credible  networks are considered. Second, and perhaps most
  importantly, our approach is significantly more efficient and
  scales to Bayesian networks with almost 60 random variables. Existing
  methods for finding the optimal Bayesian network structure <cit.>
  rely heavily for their success on a significant body of
  pruning rules that remove from consideration many candidate
  parent sets both before and during the search. We show that many 
  of these pruning rules can be
  naturally generalized to preserve the Bayesian networks that are 
  within a factor of optimal. 
  We modify GOBNILP <cit.>, a state-of-the-art method for finding
  an optimal Bayesian network, to implement our generalized pruning rules
  and to find all near-optimal networks. We show in an experimental
  evaluation that the modified GOBNILP scales to significantly
  larger networks without resorting to restricting the structure of the Bayesian networks that are learned.
  
  
  
  
  § BACKGROUND
  
  
  
  In this section, we briefly review the necessary background in
  Bayesian networks and scoring functions, and define the Bayesian network structure
  learning problem <cit.>.
  
  
  
  
   §.§ Bayesian Networks
  
  
  
  A Bayesian network (BN) is a probabilistic graphical model that
  consists of a labeled directed acyclic graph (DAG), = (, ) in which the vertices = {V_1, …, V_n}
  correspond to n random variables, the edges E represent direct
  influence of one random variable on another, and each vertex
  i is labeled with a conditional probability
  distribution P(V_i|Π_i) that
  specifies the dependence of the variable V_i on its
  set of parents Π_i in the DAG. A BN can
  alternatively be viewed as a factorized representation of the
  joint probability distribution over the random variables and
  as an encoding of the Markov condition on the nodes; i.e., given its parents, every variable is conditionally independent of its non-descendants.
  
  Each random variable V_i has state space Ω_i = {v_i1, …, v_ir_i}[Our method works with continuous and mixed BNs, although the discussion focuses on the discrete case.], where r_i is the cardinality of Ω_i and typically r_i≥ 2. Each _i has state space Ω__i = {π_i1,…,π_ir__i}.  We use r__i to refer to the number of possible instantiations of the parent set _i of V_i (see Figure <ref>). 
  The set θ = {θ_ijk} for all i = {1,…, n }, j = {1,…,r__i} and  k = {1,…,r_i} represents parameter estimates in G obtained either from expert knowledge or from a dataset, where each θ_ijk estimates the conditional probability P(V_i = v_ik|Π_i = π_ij).
  
  
  
  The predominant method for Bayesian network structure learning (BNSL) from data is
  the score-and-search method.
  Let 
  
  I =
  {I_1, …, I_N} be a dataset where each instance I_i is an n-tuple that is a
  complete instantiation of the variables in . A scoring
  function σ( | I ) assigns a real value measuring 
  the quality of =(,) given the data I.
  Without loss of generality, we assume that a lower score
  represents a better quality network structure and omit I when the data is clear from context. 
  
  
  Given a non-negative constant ϵ and a dataset I = {I_1, …, I_N}, a credible  network G is a network that has a score  such that ≤≤ + ϵ, where  is the score of the optimal Bayesian network.  
  
  
  In this paper, we focus on solving a problem we call the 
  ϵ-Bayesian Network Structure Learning (). Note 
  that the BNSL for the optimal network(s) is a special case of 
   where ϵ=0.
  
  
  Given a non-negative constant ϵ, a dataset I = {I_1, …, I_N} 
  over random variables  = {V_1, …, V_n} and a scoring 
  function σ, the ϵ-Bayesian Network Structure Learning () 
  problem is to find all credible networks.
  
  
  
  
  
  
   §.§ Scoring Functions
  
  
  
  Scoring functions usually balance goodness of fit to
  the data with a penalty term for model complexity
  to avoid overfitting. Common scoring functions
  include BIC/MDL <cit.> and BDeu <cit.>. An important property of these
  (and most) scoring functions is decomposability,
  where the score of the entire network
  σ(  )
  can be rewritten as the sum of local scores associated to each vertex
  ∑_i = 1^nσ(i,_i )
  that only depends on i and its parent set _i in
  . The local score is abbreviated below as _i  when the local node i is clear from context.
  
  
  
  
  
  
  
  Pruning techniques can be used to reduce the
  number of candidate parent sets that need to be considered,
  but in the worst-case the number of candidate parent sets
  for each variable i is exponential in n, where n
  is the number of vertices in the DAG.
  
  In this work, we focus on the Bayesian Information Criterion (BIC) and the Bayesian Dirichlet, specifically BDeu, scoring functions. The BIC scoring function[We adopt the MDL notation that calculates a positive score.] in this paper is defined as,
  
      BIC :  = -max_θ L_G,I (θ) + t(G)· w = - ∑_i=1^n ∑_j=1^r_Π_i∑_k=1^r_i n_ijklogn_ijk/n_ij +  ∑_i=1^n r_Π_i(r_i -1)log N/2.
  
  Here, w = log N/2, t(G) is a penalty term and L_G,I (θ) is the log likelihood, given by,
  
      L_G,I (θ) = ∑_i=1^n ∑_j=1^r_Π_i∑_k=1^r_ilogθ_ijk^n_ijk ,
  
  where n_ijk is the number of instances in I where v_ik and π_ij co-occur.
  As the BIC function is decomposable, we can associate a score to Π_i, a candidate parent set of i as follows,
  
      BIC : Π_i = -max_θ_i L(θ_i) + t(Π_i)· w = - ∑_j=1^r_Π_i∑_k=1^r_i n_ijklogn_ijk/n_ij +  r_Π_i(r_i -1)log N/2.
  
  Here, L (θ_i) =∑_j=1^r_Π_i∑_k=1^r_ilogθ_ijk^n_ijk and t(Π_i) = r_Π_i(r_i -1).
  The BDeu scoring function[Our BDeu notation calculates a positive score that is consistent with the minimization setting.] in this paper is defined as,
  
      BDeu : G   =   - ∑_i=1^n ( ∑_j=1^r_Π_i( logΓ( α/r__i)/Γ( α/r__i + n_ij) + ∑_k=1^r_ilogΓ(α/r_i r__i+ n_ijk)/Γ(α/r_i r__i)) ),
  
  where α is the equivalent sample size and n_ij = ∑_k n_ijk.
  As the BDeu function is decomposable, we can associate a score to Π_i, a candidate parent set of i as follows,
  
      BDeu : Π_i   =   - ∑_j=1^r_Π_i( logΓ( α/r__i)/Γ( α/r__i + n_ij) + ∑_k=1^r_ilogΓ(α/r_i r__i+ n_ijk)/Γ(α/r_i r__i)).
  
  
  
  
  
  § THE BAYES FACTOR
  
  
  
  In this section, we show that a good choice for the value of ϵ for the
   problem is closely related to the 
  Bayes factor (BF), a model selection criterion summarized in <cit.>.
  
  The BF was proposed by Jeffreys as an alternative to significance tests <cit.>. It was thoroughly examined as a practical model selection tool in <cit.>. Let _0 and _1 be two DAGs (BNs) in the set of all DAGs  defined over a set of random variables . The BF in the context of BNs is defined as,
  
      BF(_0,_1)=P(|_0)/P(|_1),
  
  namely the odds of the probability of the data predicted by network _0 and _1. The actual calculation of the BF often relies on  Bayes' Theorem as follows,
  
      P(_0 |)/P(_1 |)=P(|_0)/P(|_1)·P(_0)/P(_1)=P(,_0)/P(,_1).
  
  Since it is typical to assume the prior over models is uniform, the BF can then be obtained using either P(|) or P(,), ∀∈. We use those two representations to show how BIC and BDeu scores relate to the BF.
  
  Using the Laplace approximation and other simplifications, <cit.> derived the following approximation to the logarithm of the marginal likelihood for network  <cit.>,
  
      logP(|) =   max_θ L_G,I (θ) - t(G)·log N/2+t(G) ·log2π/2
         -1/2log|J_G,I (θ)| + logP(θ|),
  
  where J_G,I (θ) is the Hessian matrix evaluated at the maximum likelihood estimate. It follows that,
  
      logP(|)=-BIC(,)+O(1).
  
  The above equation shows that the BIC score was designed to approximate the log marginal likelihood. 
  If we drop the lower-order term, we can then obtain the following equation,
  
      BIC(,_1)-BIC(,_0)   =logP(|_0)/P(|_1) = logBF(_0,_1).
  
  
  It has been indicated in <cit.> that as N→∞, the difference of the two BIC scores, 
  dubbed the Schwarz criterion,  approaches the true value of logBF such that,
  
      BIC(,_1)-BIC(,_0)-logBF(_0,_1)/logBF(_0,_1)→ 0.
  
  Therefore, the difference of two BIC scores can be used as a rough approximation to logBF. Note that some papers define BIC to be twice as large as the BIC defined in this paper, but the above relationship still holds albeit with twice the logarithm of the BF.
  
  Similarly, the difference of the BDeu scores can be expressed in terms of the BF. In fact, the BDeu score is the log marginal likelihood where there are Dirichlet distributions over the parameters <cit.>; i.e.,
  
      logP(,)=-BDeu(,),
  
  and thus,
  
      BDeu(,_1)-BDeu(,_0)   =logP(,_0)/P(,_1) = logBF(_0,_1).
  
  
  The above results are consistent with the observation in <cit.> that the logBF can be interpreted as a measure for the relative success of two models at predicting data, sometimes referred to as the “weight of evidence”, without assuming either model is true. The maximal acceptable distance from the optimal model, however, is often specific to a study and determined with domain knowledge; e.g., a BF of 1000 is more appropriate in forensic science. <cit.> proposed the following interpreting scale for the BF: a BF of 1 to 3 bears only anecdotal evidence, a BF of 3 to 20 suggests some positive evidence that _0 is better, a BF of 20 to 150 suggests strong evidence in favor of _0, and a BF greater than 150 indicates very strong evidence. If we deem 20 to be the desired BF in , i.e., _0=^* and ϵ=log(20), then any network with a score less than log(20) away from the optimal score would be credible, otherwise it would be discredited. Note that the ratio of posterior probabilities was defined as λ in <cit.> and was used as a metric to assess arbitrary values of k in finding the k-best networks.
  
  Finally, the  problem using the BIC or BDeu scoring function given a desired BF can be written as,
  
      ≤() ≤+logBF.
  
  
  
  
  
  § PRUNING RULES FOR CANDIDATE PARENT SETS
  
  
  
  To find all near-optimal BNs given a BF, the local score σ( _i ) for each candidate parent set _i⊆ 2^ - {i} and each random variable i must be considered. As this is very cost prohibitive, it is important that the search space of candidate parent sets be pruned, provided that global optimality constraints are not violated. In this section, we generalize existing pruning rules such that the generalized rules hold when solving the  problem.
  
  
  
  
  
   
  
  
  
  A candidate parent set Π_i can be safely pruned given a non-negative constant ϵ∈ℝ^+ if Π_i cannot be the parent set of V_i in any network in the set of credible networks. Note that for ϵ=0, the set of credible  networks just contains the optimal network(s). We discuss the original rules and their generalization below and proofs for each can be found in the appendix.
  
  
  
  
  <cit.> give a pruning rule for all decomposable scoring functions. This rule compares the score of a candidate parent set to those of its subsets. We give a relaxed version of the rule.
  
  
  
  
  
  	
  
  
  	
  	Given a vertex variable j, candidate parent sets
  	Π_j and  Π_j^', and some ϵ∈ℝ^+,  if Π_j ⊂Π_j^' and Π_j + ϵ≥Π_j',
  	Π_j^' can be safely pruned if σ is a decomposable scoring function. 
  	
  	
  
  
  
   §.§ Pruning with BIC/MDL Score
  
  
  
  
  
  
  
  
  
  
  
  A pruning rule comparing the BIC score and penalty associated to a candidate parent set to those of its subsets was introduced in <cit.>. The following theorem gives a relaxed version of that rule.
  
  Given a vertex variable j,  candidate parent sets
  	Π_j and Π_j', and some ϵ∈ℝ^+,
  	if Π_j ⊂Π_j' and Π_j -  Π_j' + ϵ < 0,
  	Π_j' and all supersets of Π_j^' can be safely pruned  if σ is the BIC scoring function. 
  	
  
  
  Another pruning rule for BIC appears in <cit.>. This provides a bound on the number of possible instantiations of subsets of a candidate parent set.
  
  
  
  
  
  
  
  
  
  
  	
  
  
  
  	Given a vertex variable V_i, and a candidate parent set Π_i such that r_Π_i> N/wlog r_i/r_i -1 + ϵ for some ϵ∈ℝ^+,  if _i ⊊_i' , then _i' can be safely pruned if σ is the BIC scoring function. 
  
  
  The following corollary of Theorem <ref> gives a useful upper bound on the size of a candidate parent set.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	Given a vertex variable i and candidate parent set 
  	Π_i, if Π_i has more than log_2 (N + ϵ) elements, for some ϵ∈ℝ^+, Π_i can be safely pruned if σ is the BIC scoring function. 
  	
  	
  Corollary <ref> provides an upper-bound on the size of parent sets based solely on the dataset size N. The following table summarizes such an upper-bound given different amounts of data N and a BF of 20.
  
  
  The entropy of a candidate parent set is also a useful measure for pruning. A pruning rule, given by <cit.>, provides an upper bound on the conditional entropy of candidate parent sets and their subsets. We give a relaxed version of their rule. First, we note that the sample estimate of entropy for a variable i is given by,
  
      H(i)    =  -∑^r_i_k=1n_ik/Nlogn_ik/N ,
  
  where n_ik represents how many instances in the dataset contain v_ik, where v_ik is an element in the state space Ω_i of i. Similarly, the sample estimate of entropy for a candidate parent set Π_i is given by,
  
      H(Π_i)    =  -∑^r_Π_i_j=1n_ij/Nlogn_ij/N .
  
  Conditional entropy is given by,
  
      H(X | Y ) = H(X ∪ Y) - H(Y) .
  
  
  
  
  
  
  
  
  
  
  
  
  
  Given a vertex variable V_i, and candidate parent set Π_i, let V_j ∉Π_i  such that N ·min{H(V_i |Π_i), H(V_j |Π_i)}≥ (1 - r_j) · t(Π_i) +ϵ for some ϵ∈ℝ^+. Then the candidate parent set Π_i' = Π_i ∪{V_j } and all its supersets can be safely pruned if σ is the BIC scoring function.
  
  
  
  
  
  
   §.§ Pruning with BDeu Score
  
  
  
  A pruning rule for the BDeu scoring function appears in <cit.> and a more general version is included in <cit.>. Here, we present a relaxed version of the rule in <cit.>.  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Given a vertex variable V_i and candidate parent sets Π_i and Π_i' such that Π_i ⊂Π_i' and Π_i ≠Π_i', let
  r_i^+(Π_i') := |{j : n_ij > 0 , j ∈Ω_Π_i'}| be the total number of  instantiations of Π_i' that appear in the dataset. If Π_i + ϵ < r_i^+(Π_i') log r_i, for some ϵ∈ℝ^+ then Π_i' and the supersets of Π_i' can be safely pruned if σ is the BDeu scoring function.
  
  
   
  
  
  
  
  § EXPERIMENTAL EVALUATION
  
  
  
  In this section, we evaluate our proposed BF-based method and compare its performance with published k-best solvers. 
  
  Our proposed method is more memory efficient comparing to the k-best based solvers in BDeu scoring and often collects more networks in a shorter period of time. With the pruning rules generalized above, our method can scale up to datasets with 57 variables in BIC scoring, whereas the previous best results are reported on a network of 29 variables using the k-best approach with score pruning <cit.>.
  
  The datasets are obtained from the UCI Machine Learning 
  Repository <cit.> and the Bayesian Network 
  Repository[<http://www.bnlearn.com/bnrepository/>].
  Both BIC/MDL <cit.> and 
  BDeu <cit.> scoring functions are used 
  where applicable. All experiments are conducted on computers with 2.2 GHz Intel E7-4850V3 processors. Each experiment is limited to 64 GB of memory and 24 hours of CPU time.
  
  We demonstrate the effect of applying pruning rules during scoring in Section <ref>. 
  Our generalized rules are able to eliminate the majority of the search space and therefore allow us to apply the ϵ-BNSL algorithm to medium sized networks. 
  We discuss the implementation details of the BF approach in Section <ref> and present experimental results with BIC scores on a wide range of datasets commonly used in BNSL. 
  We show the effect of varying sample sizes on our approach using data generated from synthetic BNs in Section <ref>. 
  Finally, We compare our approach with the k-best method in Section <ref>.
  
  
  
  
   §.§ The Effect of Pruning
  
  
  
  
  
  We modified the development version (9c9f3e6) of GOBNILP to apply the generalized pruning rules in Section <ref>. In particular, Lemma <ref> is applied to both BIC and BDeu; Theorem <ref> and Corollary <ref> are applied to BIC; Theorem <ref> is applied to BDeu. The combination of those rules effectively pruned more than 95% of the parent sets for almost all datasets. The worst pruning rate of 88.9% is observed on the letter dataset with a BF of 150 using BIC. We report the number of remaining candidate parent sets in Table <ref>. The generalized pruning rules allow us to scale up to medium sized networks, unlike previous approaches where the lack of effective pruning rules restricts them to small networks.
  
  
  
  
   §.§ The Bayes Factor Approach
  
  
  
  
  
  We modified the development version (9c9f3e6) of GOBNILP, denoted hereafter as _dev, and supplied appropriate parameter settings for collecting near-optimal networks[The modified code is available at: <https://www.cs.york.ac.uk/aig/sw/gobnilp/>]. The code is compiled with SCIP 6.0.0 and CPLEX 12.8.0.  extends the SCIP Optimization Suite <cit.> by adding a constraint handler for
  handling the acyclicity constraint for DAGs. If multiple BNs are
  required _dev just calls SCIP to ask it to collect feasible
  solutions.  In this mode, when SCIP finds a solution, the solution is stored, a
  constraint is added to render that solution infeasible and the search
  continues. This differs from (and is much more efficient than)
  the method used in the current stable version of  for finding k-best BNs where an entirely
  new search is started each time a new BN is found. A recent version of
  SCIP has a separate “reoptimization” method which might allow better
  k-best performance for  but we do not explore that here.
  By default when SCIP is asked to collect solutions it turns off all
  cutting plane algorithms. This led to very poor  performance
  since  relies on cutting plane generation. Therefore, this default setting is overridden in _dev to allow cutting planes when collecting solutions.
  To find only solutions with objective no worse than (+ ϵ), SCIP's
   function is used. Note that, for efficiency
  reasons, this is not effected by adding a linear constraint.
  
  
  
  We first use GOBNILP_dev to find the optimal score since GOBNILP_dev takes objective limit (+ ϵ) for enumerating feasible networks. For BIC, We set the limit on the size of the parent set based on Corollary <ref> that guarantees optimality, whereas for BDeu we set the number to -1 to allow all possible sizes. Then all networks falling into the limit are collected with a counting limit of 150,000. Finally the collected networks are categorized into Markov equivalence classes (MECs)[Our code can also collect only one DAG from each MEC.], where two networks belong to the same MEC iff they have the same skeleton and v-structures <cit.>. The proposed approach is tested on datasets with up to 57 variables. The search time T, the number of collected networks || and the number of MECs ℳ in the collected networks at BF = 3, 20 and 150 using BIC are reported in Table <ref>, where n is the number of random variables in the dataset and N is the number of instances in the dataset. The three thresholds are chosen according to the interpreting scale suggested by <cit.> where 3 marks the transition between anecdotal and positive evidence, 20 marks the transition between positive and strong evidence and 150 marks the transition between strong and very strong evidence. The search time mostly depends on a combined effect of the size of the network, the sample size and the number of MECs at a given BF. Some fairly large networks such as alarm, sponge and barley are solved much faster than smaller networks with a large sample size; e.g., msnbc and letter.
  
  The results also indicate that the number of collected networks and the number of MECs at three BF levels varies substantially across different datasets. In general, datasets with smaller sample sizes tend to have more networks collected at a given BF since near-optimal networks have similar posterior probabilities to the best network. Although the desired level of BF for a study, like the p-value, is often determined with domain knowledge, the proposed approach, given sufficient samples, will produce meaningful results that can be used for further analysis.
  
  
  
  
   §.§ Synthetic Data
  
  
  
  We use BNs up to 76 nodes from the Bayesian Network Repository to generate synthetic data in the form of 10 random samples for various sample sizes up to 1,000. Near optimal networks are collected following the same procedure outlined in Section <ref> using BIC. For the three largest datasets, hailfinder, hepar2, and win95pts, the scoring process failed to complete within 24 hours for the sample size of 1,000.
  
  The average number of networks in the credible sets is reported in Table <ref>. The results are consistent with Table <ref> in demonstrating that the number of collected networks are specific to the dataset. Collecting a large number of networks is not always ideal for large datasets, e.g., water, mildew, and barley all have a very small number of networks in the credible sets comparing to some other datasets with fewer nodes. We also note that increasing sample size does not always lead to smaller credible sets. Instead, the number of networks in the credible sets tend to peak around certain sample sizes. For example, the largest credible sets for water are collected with a sample size of 500, although the numbers are quite different sometimes across 10 trials as indicated by the large standard deviation.
  
  
  
  The average number of equivalence classes is reported in Table <ref>. The large number of networks can indeed be represented by a handful of equivalence classes. Increasing the amount of training data can both lead to a decrease in the number of networks and an increase in the complexity of the collected networks. Although the former is more evident for most datasets, water, mildew, and barley are examples of the latter. There are still peaks of sample sizes where large numbers of equivalence classes are collected, but the scenario occurs much less frequently than in Table <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
        
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
        
  
  
  
  
  
  
  
   §.§ Bayes Factor vs. K-Best
  
  
  
  In this section, we compare our approach with published solvers that are able to find a subset of top-scoring networks with the given parameter k. The solvers under consideration are KBest_12b[<http://web.cs.iastate.edu/ jtian/Software/UAI-10/KBest.htm>] from <cit.>, KBestEC[<http://web.cs.iastate.edu/ jtian/Software/AAAI-14-yetian/KBestEC.htm>] from <cit.>, and GOBNILP 1.6.3 <cit.>,
  
  referred to as KBest, KBestEC and GOBNILP below. The first two solvers are based on the dynamic programming approach introduced in <cit.>. 
  
  Due to the lack of support for BIC in KBest and KBestEC, only BDeu with a equivalent sample size of one is used in corresponding experiments.
  
  The most recent stable version of  is 1.6.3 that works with SCIP 3.2.1. The default configuration is used and experiments are conducted for both BIC and BDeu scoring functions. However, the k-best results are omitted here due to its poor performance. Despite that  can iteratively find the k-best networks in descending order by adding linear constraints, the pruning rules designed to find the best network are turned off to preserve sub-optimal networks. In fact, the memory usage often exceeded 64 GB during the initial ILP formulation, indicating that the lack of pruning rules posed serious challenge for GOBNILP. _dev, on the other hand, can take advantage of the pruning rules presented above in the proposed BF approach and its results compare favorably to KBest and KBestEC.
  
  The experimental results of KBest, KBestEC and _dev are reported in Table <ref>, where n is the number of random variables in the dataset, N is the number of instances in the dataset, and k is the number of top scoring networks. The search time T is reported for KBest, KBestEC and GOBNILP_dev (BF = 20). The number of DAGs covered by the k MECs |_k| is reported for KBestEC. In comparison, the last two columns are the number of found networks |_20| and the number of MECs |ℳ_20| using the BF approach with a given BF of 20 and BDeu scoring function.
  
  As the number of requested networks k increases, the search time for both KBest and KBestEC grows exponentially. The KBest and KBestEC are designed to solve problems of size fewer than 20, and so they have some difficulty with larger datasets.[Obtained through correspondence with the author.] They also fail to generate correct scoring files for msnbc. KBestEC seems to successfully expand the coverage of DAGs with some overhead for checking equivalence classes. However, KBestEC took much longer than KBest for some instances, e.g., nltcs and letter, and the number of DAGs covered by the found MECs is inconsistent for nltcs, letter and zoo. The search time for the BF approach is improved over the k-best approach except for datasets with very large sample sizes. The generalized pruning rules are very effective in reducing the search space, which then allows GOBNILP_dev to solve the ILP problem subsequently. Comparing to the improved results in <cit.>, our approach can scale to larger networks if the scoring file can be generated.[We are unable to generate BDeu score files for datasets with 30 or more variables.]
  
  
  
  Now we show that different datasets have distinct score patterns in the top scoring networks. The scores of the 1,000-best networks for some datasets in the KBest experiment are plotted in Figure <ref>. A specific line for a dataset indicates the deviation ϵ from the optimal BDeu score by the kth-best network. For reference, the red dash lines represent different levels of BFs calculated by ϵ=logBF (See Equation <ref>). The figure shows that it is difficult to pick a value for k a priori to capture the appropriate set of top scoring networks. For a few datasets such as adult and letter, it only takes fewer than 50 networks to reach a BF of 20, whereas zoo needs more than 10,000 networks. The sample size has a significant effect on the number of networks at a given BF since the lack of data leads to many BNs with similar probabilities. It would be reasonable to choose a large value for k in model averaging when data is scarce and vice versa, but only the BF approach is able to automatically find the appropriate and credible set of networks for further analysis.
  
  
  
  
  § CONCLUSION
  
  
  
  Existing approaches for model averaging for Bayesian network structure learning
  either severely restrict the structure of the Bayesian network
  or have only been shown to scale to networks with fewer than 30 random
  variables. In this paper, we proposed a novel approach to model
  averaging in Bayesian network structure learning that 
  finds all networks within a factor of optimal. 
  Our approach has two primary advantages. First,
  our approach only considers credible models in that they
  are optimal or near-optimal in score. Second, our approach
  is significantly more efficient and scales to much larger
  Bayesian networks than existing approaches.
  We modified GOBNILP <cit.>, a state-of-the-art method for finding
  an optimal Bayesian network, to implement our generalized pruning rules
  and to find all near-optimal networks. Our experimental results demonstrate
  that the modified GOBNILP scales to significantly
  larger networks without resorting to restricting the structure 
  of the Bayesian networks that are learned.
  
  
  
  
  § PROOFS OF PRUNING RULES
  
  
  
  We discuss the original pruning rules and prove their generalization below. A candidate parent set Π_i can be safely pruned given a non-negative constant ϵ∈ℝ^+ if Π_i cannot be the parent set of V_i in any network in the set of credible networks. Note that proofs of the original rules can be obtained by setting ϵ=0.
  
  
  
  
   §.§ Proof of Lemma <ref>
  
  
  
  <cit.> give a pruning rule that is applicable for all decomposable scoring functions.
  <cit.>
  	Given a vertex variable j, and candidate parent sets
  	Π_j and  Π_j^', if Π_j ⊂Π_j^' and Π_j < Π_j^',
  	Π_j^' can be safely pruned. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  We relax this pruning rule and prove Lemma <ref> below.
  
  	(Lemma <ref>)
  Consider networks  and ' that are the same except for the parent set of j, where  has the parent set Π_j for j and ' has the parent set Π_j^' for j. 
  
      =  Π_j + ∑_i≠ jΠ_j   [ is decomposable]
         ≤Π_j^' + ϵ+ ∑_i≠ jΠ_j   [given]
         ='.
  
  Thus, ' cannot be in the set of credible  networks.
  	
  	
  
  
  
   §.§ Proof of Theorem <ref>
  
  
  
  An additional pruning rule can be derived from Theorem <ref> that is applicable to the BIC/MDL scoring function. 
  
  <cit.>
  Given a vertex variable i, and candidate parent sets Π_i and Π_i^',
  	if Π_i ⊂Π_i^' and Π_i -  Π_i^' < 0,
  	Π_i^' and all supersets of Π_i^' can be safely pruned  if σ is the BIC/MDL scoring function.    
  
  
  
  
  
  
  
  
  
  
  
  
  
  Here, Π_i^' is the penalty term in the BIC scoring function. This pruning rule is relaxed as Theorem <ref> and we prove it below. 
  
  (Theorem <ref>)
      
      Π_i -  Π_i^'  + ϵ < 0     [given]
         ⇒ -Π_i + Π_i^' -ϵ > 0
         ⇒ -Π_i + Π_i^' - L(Π_i^')  - ϵ > 0    [L(Π_i^') < 0]
         ⇒Π_i^' > Π_i + ϵ.
  
      By Lemma <ref>, Π_i^' cannot be an optimal parent set. Using the fact that penalties increase with increase in parent set size, supersets of Π_i^' cannot be in the set of credible  networks. The result follows.
  
  
  
  
  
  
  
  
  
   §.§ Proof of Theorem <ref>
  
  
  
  
  
  <cit.>
  	 Given a vertex variable V_i and candidate parent set Π_i such that r_Π_i> N/wlog r_i/r_i -1,  if _i ⊊_i' , then _i' can be safely pruned if σ is the BIC scoring function. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	
  We relax the pruning rule given in Theorem <ref> as Theorem <ref> and prove it below.
  
  (Theorem <ref>)
  
      σ(Π_i') - σ(Π_i)
         []0= -max_θ_i L(Π_i') + t(Π_i')· w + max_θ_i L(Π_i) - t(Π_i)· w 
         []1≥ -max_θ_i L(Π_i) + t(Π_i')· w  - t(Π_i)· w 
         []2= -∑_j=1^r_Π_in_ij(-∑^r_i_i=1n_ijk/n_ijlogn_ijk/n_ij )+ t(Π_i')· w - t(Π_i)· w 
         []3≥ -∑_j=1^r_Π_in_ijH(θ_ij)- t(Π_i')· w  + t(Π_i)· w 
         []4≥ -∑_j=1^r_Π_in_ijlog r_i + r_Π_i·(r_e -1)· (r_i -1)· w 
         []5≥ -∑_j=1^r_Π_in_ijlog r_i + r_Π_i· (r_i -1)· w 
         []6=  -Nlog r_i + r_Π_i· (r_i -1)· w
         []7>ϵ.
  
      Step 0  uses the definition of BIC. Step 1 uses max_θ_i L(Π_i') is negative. Step 2 uses the fact that the maximum likelihood estimate, θ^∗_ijk = n_ijk/n_ij and n_ij = ∑_i=1^r_in_ijk. Step 3 uses the definition of entropy. Step 4 uses the definition of the penalty function t. Step 5 uses r_e ≥ 2.  Finally, the RHS in Step 6 follows because of the definition of n_ij. Step 7 uses the assumption of the theorem. 
      
      Using Lemma <ref>, we get the result as desired.  
  	
  	
  
  
  
   §.§ Proof of Corollary <ref>
  
  
  <cit.>
  	Given a vertex variable i and candidate parent set 
  	Π_i, if Π_i has more than log_2 N elements,
  	Π_i can be safely pruned if σ is the BIC scoring function. 
  
  
  Using Theorem <ref>, we generalize Corollary <ref> to Corollary <ref> and prove it below.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	(Corollary <ref>)
  	Assuming N>4, take a variable V_i and a parent set Π_i
  with |Π_i| = log_2 (N + ϵ) elements.  Because every variable has at least two states, we know that r_Π_i≥ 2^|Π_i|≥ N +ϵ > N/wlog r_i/r_i -1+ ϵ, because w =logN/2 gives us log r_i/w(r_i -1) < 1 , and by Theorem <ref> we know that no proper superset of Π_i can be an optimal parent set for i as desired.
  	 
  
  
  
  
   §.§ Proof of Theorem <ref>
  
  
  
  <cit.>
  Given a vertex variable V_i, and candidate parent sets Π_i, Π_i' such that Π_i = Π_i' ∪{V_j} for some variable V_j ∉Π_i', we have  L(Π_i) - L(Π_i') ≤  N ·min{H(V_i |Π_i'), H(V_j  |Π_i')}. 
  
  
  
  First, consider the definition of L_i(Π_i),
  
      L(Π_i) = max_θ∑_j=1^r_Π_i∑_k=1^r_i n_ijklogθ_ijk ,
  
  
  where the maximum likelihood estimate of θ_ijk is n_ijk/n_ij.
  This gives us N· H(V_i |Π_i) = -L(Π_i). Thus, we get,
  
      L(Π_i) - L(Π_i')    = N· (H(V_i |Π_i') - H(V_i |Π_i)[]1≤ N· H(V_i |Π_i') .
  
  We use the fact that entropy is positive. Now, consider the definition of mutual information,
  
      I(X,Y | Z) = H(X | Z) = H(X | Y ∪ Z).
  
  This gives us,
  
      L(Π_i) - L(Π_i')    = N · I(V_i, V_j  |Π_i')
         []2= N· (H(V_j |Π_i') - H(V_j |Π_i' ∪{V_i}))
      ⇒ L(Π_i) - L(Π_i')    []3≤ N ·min{H(V_i |Π_i'), H(V_j  |Π_i')}.
  
  Step 3 combines Steps 1 and 2. The result follows as desired.
  
  
  <cit.>
  Given a vertex variable V_i, and candidate parent set Π_i, let V_j ∉Π_i  such that N ·min{H(V_i |Π_i), H(V_j  |Π_i)}≥ (1 - r_j) · t(Π_i). Then the candidate parent set Π_i' = Π_i ∪{V_j } and all its supersets can be safely pruned if σ is the BIC scoring function. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We relax Theorem <ref> and prove its generalization below.
  
  (Theorem <ref>)
  
      Π_i'   []0= -L(Π_i')  +  t(Π_i')
         []1≥ -L(Π_i) - N ·min{H(V_i |Π_i); H(V_j |Π_i)} + t(Π_i')
         []2≥ -L(Π_i) + (1 - r_j)· t(Π_i) +ϵ  +t(Π_i') 
         []3= -L(Π_i) + t(Π_i) - r_j · t(Π_i) + ϵ + t(Π_i') 
         []4= -L(Π_i) + t(Π_i) - r_j · r_Π_i· (r_i -1) + ϵ + t(Π_i') 
         []5= -L(Π_i) + t(Π_i) - t(Π_i') + ϵ + t(Π_i') 
         []6=Π_i + ϵ.
  
  Step 1 uses Lemma <ref>. Step 2 uses the assumptions of the question. Step 4 uses the definition of t. Step 5 uses Π_i' = Π_i ∪{V_j }. Using Lemma <ref>, the result follows as desired.
  
  
  
  
  
  
  
   §.§ Proof of Theorem <ref>
  
  
  
  
  
  Let n_ij be a positive integer and α' be a positive real number. Then
  
      logΓ(n_ij + α')/Γ(α') = ∑_i=0^n_ij -1log (i + α')
   
  
  
  
  We start with the property that Γ(x+1) = xΓ(x) for any positive real number x. As α' > 0, this gives us,
  
      Γ(1 + α')/Γ(α')    []0=α' 
      Γ(2 + α')/Γ(1+ α')   []1= (1+ α') 
      ⇒Γ(1 + α')·Γ(2 + α')/Γ(1+ α')Γ(α')    []2=α' (1+ α')
      ⇒Γ(1 + α')⋯Γ(n_ij
          + α')/Γ(n_ij -1 + α')⋯Γ(α')   []3=α'⋯ (n_ij -1 + α') 
      ⇒Γ(n_ij + α') /Γ(α')   []4=α'⋯ (n_ij -1 + α')
      ⇒logΓ(n_ij+ α')/Γ(α')   []5=∑_i=0^n_ij -1log (i + α').
  
  Step 1 uses 1 + α'. Step 2 follows by multiplication of the equations in Step 1 and Step 0. Step 3 follows by repeated application of the identity. Step 4 cancels identical terms in the LHS. The result follows as desired.
  
  
  6B
  Let {n_ijk}_k=1,...r_i be non-negative integers with a positive sum, n_ij = ∑_k=1^r_i n_ijk and α” be a positive real number. Then
  
      ∑_k=1^r_ilogΓ(n_ijk + α”)/Γ(α”)≤logΓ(n_ij + α”)/Γ(α”)
   
  
  
  
  Consider allocation of {n_ijk}_k=1,...,r_i  items over the r_i bins. There are two cases.
  
      
    * Let there be some index k^∗ such that n_ijk^∗ = n_ij. This means that n_ijk=0 for all k≠ k^∗. It follows that ∑_k=1^r_ilogΓ(n_ijk + α”)/Γ(α”) = logΓ(n_ij + α”)/Γ(α”).
      
    * Let there be two indices k_1 and k_2 such that n_ijk_1 > 0 and n_ijk_2 > 0. Without loss of generality, we can assume that n_ijk_1≥ n_ijk_2. We move one item from bin k_1 to bin k_2. The sum n_ij remains constant. By Lemma <ref>, an increase in the RHS by log (n_ijk_1  + α”) - log (n_ijk_2  -1+ α”), results in a corresponding increase in the LHS. Note that the assumption n_ijk_1≥ n_ijk_2 means that this increase is positive. 
       By increasing counts at the expense of small counts in  this way a sequence of distributions of the fixed sum n_ij over the r_i bins can be constructed for which the LHS of Lemma <ref> is increasing. The sequence terminates when n_ijk^∗ = n_ij for some k^∗. The result follows.
   
  
  
  
  
  <cit.>
  
      ∑_j=1^r_Π_i( Γ( α')/Γ(n_ij+α') +  ∑_k=1^r_ilogΓ(n_ijk + α'/r_i)/Γ( α'/r_i)) ≤∑_i=0, j:n_ij>0^n_ijlog( i+ a'/r_i/i+α).
   
  
  
  
  
   
      ∑_j=1^r_Π_i( logΓ( α')/Γ(n_ij+α') +  ∑_k=1^r_ilogΓ(n_ijk + α'/r_i)/Γ( α'/r_i))
         []1≤∑_j=1^r_Π_i(logΓ( α')/Γ(n_ij+ α') + logΓ(n_ij + α'/r_i)/Γ( α'/r_i)) 
         []2≤∑_j=1^r_Π_i(logΓ( α')/Γ(n_ij+ α')Γ(n_ij + α'/r_i)/Γ( α'/r_i)) 
         []3≤∑_i=0, j:n_ij>0^n_ij-1( logi+ a'/r_i/i+α')  
         []4≤∑_i=0, j:n_ij>0^n_ijlog( i+ a'/r_i/i+α').
  
  
  Step 1 uses Lemma <ref>. Step 2 assumes n_ij > 0, and uses properties of the logarithm function. Step 3 uses Lemma <ref>. The result follows as desired.
  
  <cit.> 
  Given that r_i^+ := |{j: n_ij > 0}|, then 
  
      ∑_j=1^r_Π_ilogΓ( α')/Γ(n_ij + α') +  ∑_k=1^r_ilogΓ(n_ijk + α'/r_i)/Γ( α'/r_i)≤ -r_i^+log r_i.
  
  
  
  
  
   If n_ij > 0, then 
   
      ∑_i=0^n_ijlog( i+ a'/r_i/i+α')    = - log r_i ∑_i=1^n_ijlog( i+ a'/r_i/i+α') ≤ -log r_i.
  
  Note that as r_i ≥ 2, and α' > 0, it is clear that i + α'/r_i < i+ α'. This means that each term in ∑_i=1^n_ijlog( i+ a'/r_i/i+α') is negative.  This gives us the second inequality. The result then follows from Theorem <ref> as desired.
  
  <cit.>
  Given a vertex variable V_i and candidate parent sets Π_i and Π_i'
  such that Π_i ⊂Π_i' and Π_i ≠Π_i', let
  r_i^+(Π_i') be the number of positive counts in the contingency table for Π_i'. If Π_i < r_i^+(Π_i') log r_i then Π_i' and the supersets of Π_i' can be safely pruned.
  
  
  
  
  
   
  
  
  
  
  
   We generalize Corollary <ref> to Theorem <ref> and prove it below.
  
  (Theorem <ref>)
   Let G' be a Bayesian network where Π_i' or one of its supersets is a parent set for V_i.  Let G be another Bayesian network where Π_i is the parent set for V_i. 
   
  Consider the LHS of Corollary <ref>. It is the local BDeu score for a parent set Π_i' which has r_Π_i counts n_ij in its contingency table and counts n_ijk in the contingency table for Π_i' ∪{V_i}, where α' = α/r_Π_i for some ESS α. If r_i^+(Π_i') log r_i > Π_i + ϵ then Π_i + ϵ is
  lower than the local BDeu score for Π_i' due to Corollary <ref>. Take a candidate parent set Π_i”. If Π_i' ⊂Π_i” then r_i^+(Π_i”)  ≤ r_i^+(Π_i') and so r_i^+(Π_i”) log r_i ≤ r_i^+(Π_i') log r_i, as r_i ≥ 2. From this
  it follows that the local score for Π_i” must also be more than Π_i + ϵ. Using Lemma <ref>, the result follows as desired.
  
  
  0.2in
  
  
  