{'fraction_non_alphanumeric': 0.050956422599442615, 'fraction_numerical': 0.01260450975424373, 'mean_word_length': 4.152063958231359, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Large-Scale Generalized Hypergeometric Distribution (GHD) DAG models

We introduce a new class of identifiable DAG models where the conditional distribution of each node given its parents belongs to a family of generalized hypergeometric distributions (GHD). A family of generalized hypergeometric distributions includes a lot of discrete distributions such as the binomial, Beta-binomial, negative binomial, Poisson, hyper-Poisson, and many more. We prove that if the data drawn from the new class of DAG models, one can fully identify the graph structure. We further present a reliable and polynomial-time algorithm that recovers the graph from finitely many data. We show through theoretical results and numerical experiments that our algorithm is statistically consistent in high-dimensional settings (p>n) if the indegree of the graph is bounded, and out-performs state-of-the-art DAG learning algorithms.


Introduction
Probabilistic directed acyclic graphical (DAG) models or Bayesian networks provide a widely used framework for representing causal or directional dependence relationships among many variables. One of the fundamental problems associated with DAG models is learning a large-scale causal structure given samples from the joint distribution P (G) over a set of nodes of a graph G.
Prior works have addressed the question of identifiability for different classes of joint distribution P (G). [5,6] show the Markov equivalence class (MEC) where graphs that belong to the same MEC have the same conditional independence relations. [2,18,19,23] show that the underlying graph of a DAG model is recoverable up to MEC under the faithfulness or some related conditions. However since many MECs contain more than one graph, a true graph cannot be determined.
Recently, many works show fully identifiable DAG models under stronger assumptions on joint distribution P (G). [14] proves that Gaussian structural equation models with equal or known error variances are identifiable. In addition, [17] shows that linear non-Gaussian models where each variable is determined by a linear function of its parents plus a non-Gaussian error term are identifiable. [7,11,16] relax the assumption of linearity and prove that nonlinear additive noise models where each variable is determined by a non-linear function of its parents plus an error term are identifiable under suitable regularity conditions. Instead of considering linear or additive noise models, [12,13] introduce discrete DAG models where the conditional distribution of each node given its parents belongs to the exponential family of discrete distributions such as Poisson, binomial, and negative binomial. They prove that the discrete DAG models are identifiable as long as the variance is a quadratic function of the mean. Learning DAG or causal discovery from count data is an important research problem because such count data are increasingly ubiquitous in big-data settings, including high-throughput genomic sequencing data, spatial incidence data, sports science data, and disease incidence data. However as we discussed, most existing methods focus on the continuous or limited discrete DAG models. Hence it is important to model complex multivariate count data using a broader family of discrete distributions.
In this paper, we generalize the main idea in [12,13] to a family of generalized hypergeometric distributions (GHD) that includes Poisson, hyper-Poisson, binomial, negative binomial, beta-binomial, hypergeometric, inverse hypergeometric and many more (see more examples in Table 1 and [3,8,9]). We introduce a new class of identifiable DAG models where the conditional distribution of each node given its parents belongs to a family of GHDs. In addition, we prove that the class of GHD DAG models is identifiable from the joint distribution P (G) using convex relationship between the mean and the r-th factorial moment for some positive integer r under the causal sufficiency assumption that all relevant variables have been observed.
However we do not assume the faithfulness assumption that can be very restrictive [21].
We also develop the reliable and scalable Moments Ratio Scoring (MRS) algorithm which learns any large-scale GHD DAG model. We provide computational complexity and statistical guarantees of our MRS algorithm to show that it has polynomial run-time and is consistent for learning GHD DAG models, even in the high-dimensional p > n setting when the indegree of the graph d is bounded. We demonstrate through simulations and a real NBA data that our MRS algorithm performs better than state-of-the-art GES [2], MMHC [20], and ODS [12] algorithms in terms of both run-time and recovering a graph structure.
The remainder of this paper is structured as follows: Section 2.1 summarizes the necessary notation, Section 2.2 defines GHD DAG models and Section 2.3 proves that GHD DAG models are identifiable.
In Section 3, we develop a polynomial-time algorithm for learning GHD DAG models and provide its theoretical guarantees and computational complexity in terms of the triple (n, p, d). Section 4 empirically evaluates our methods compared to GES, MMHC, and ODS algorithms on synthetic and real basketball data.

GHD DAG Models and Identifiability
In this section, we first introduce some necessary notations and definitions for directed acyclic graph (DAG) models. Then we propose novel generalized hypergeometric distribution (GHD) DAG models for multivariate count data. Lastly, we discuss their identifiability using a convex relation between the mean and r-th factorial moments.

Problem set-up and notation
A DAG G = (V, E) consists of a set of nodes V = {1, 2, · · · , p} and a set of directed edges E ∈ V × V with no directed cycles. A directed edge from node j to k is denoted by (j, k) or j → k. The set of parents of node k denoted by Pa(k) consists of all nodes j such that (j, k) ∈ E. If there is a directed path j → · · · → k, then k is called a descendant of j and j is an ancestor of k. The set De(k) denotes the set of all descendants of node k. The non-descendants of node k are Nd(k) := V \ ({k} ∪ De(k)). An important property of DAGs is that there exists a (possibly non-unique) ordering π = (π 1 , ...., π p ) of a directed graph that represents directions of edges such that for every directed edge (j, k) ∈ E, j comes before k in the ordering. Hence learning a graph is equivalent to learning an ordering and skeleton that is a set of edges without their directions.
We consider a set of random variables X := (X j ) j∈V with a probability distribution taking values in probability space X v over the nodes in G. Suppose that a random vector X has a joint probability density function P (G) = P (X 1 , X 2 , ..., X p ). For any subset S of V , let X S := {X j : j ∈ S ⊂ V } and X (S) := × j∈S X j . For any node j ∈ V , P (X j | X S ) denotes the conditional distribution of a variable X j given a random vector X S . Then, a DAG model has the following factorization [10]: where P (X j | X Pa(j) ) is the conditional distribution of a variable X j given its parents X Pa(j) .
We suppose that there are n i.i.d samples X 1:n := (X (i) ) n i=1 drawn from a given DAG models where is a p-variate random vector. We use the notation · to denote an estimate based on samples X 1:n . In addition, we assume the causal sufficiency that all variables have been observed.

Generalized Hypergeometric Distribution (GHD) DAG Models
In this section, we begin by introducing a family of generalized hypergeometric distributions (GHDs) defined by [8]. A family of GHDs includes a large number of discrete distributions and has a special form of probability generating functions expressed in terms of the generalized hypergeometric series. We borrow the notations and terminologies in [9] to explain detailed properties of a family of GHDs. Let a j = a(a + 1) · · · (a + j − 1) be the rising factorial, (a) j = a(a − 1) · · · (a − j + 1) be the falling factorial, and a 0 = (a) 0 = 1. In addition, generalized hypergeometric function is: [ 8,9] show that GHDs have probability generating functions of the following form: This class of distributions includes a lot of discrete distributions such as binomial, beta-binomial, Poisson, Poisson type, displaced Poisson, hyper-Poisson, logarithmic, and generalized log-series. We provide more examples with their probability generating functions in Table 1 (see also in [3,8,9]). Now we define the generalized hypergeometric distribution (GHD) DAG models: Hyper-Poisson (Bardwell and Crow) (GHD) DAG models if the conditional distribution of each node given its parents belongs to a family of generalized hypergeometric distributions and the parameter depend only on its parents: For each j ∈ V , X j | X Pa(j) has the following probability generating function where a(j) = (a j1 , ..., a jp j ), b(j) = (b j1 , ..., b jq j ), and θ : X Pa(j) → R.
A popular example of GHD DAG models is a Poisson DAG model in [12] where a conditional distribution of each node j ∈ V given its parents is Poisson and the rate parameter is an arbitrary positive function θ j (X Pa(j) ). Unlike Poisson DAG models, GHD DAG models are hybrid models where the conditional distributions have various distributions which incorporate different data types. In addition, the exponential family of discrete distributions discussed in [13] is also included in a family of GHDs. Hence, our class of DAG models is strictly broader than the previously studied identifiable DAG models for multivariate count data.
GHD DAG models have a lot of useful properties for identifying a graph structure. One of the useful properties is the recurrence relation involving factorial moments: Consider a GHD DAG model. Then for any j ∈ V and any integer r = 2, 3, ..., there exists a r-th factorial constant moments ratio (CMR) function as long as max X j ≥ r.
The detail of the proof is provided in Appendix. Prop. 2.2 claims that the GHD DAG models always satisfy the r-th constant moments ratio (CMR) property that the r-th factorial moment is a function of the Figure 1: Bivariate directed acyclic graphs of G 1 , G 2 and G 3 mean. The condition max X j ≥ r for r ≥ 2 rules out the DAG models with Bernoulli and multinomial distributions which are known to be non-identifiable [6]. We will exploit the CMR property for model identifiability in the next section.

Identifiability
In this section we prove that GHD DAG models are identifiable. To provide intuition, we show identifiability for the bivariate Poisson DAG model discussed in [12]. Consider all possible graphical models illustrated in Fig. 1: , where X 1 and X 2 are independent; G 2 : X 1 ∼ Poisson(λ 1 ) and X 2 | X 1 ∼ Poisson(θ 2 (X 1 )); and G 3 : X 2 ∼ Poisson(λ 2 ) and ) for arbitrary positive functions θ 1 , θ 2 : N ∪ {0} → R + . Our goal is to determine whether the underlying graph is G 1 , G 2 or G 3 from the probability distribution P (G).
as long as E(X 2 | X 1 ) is not a constant. The inequality follows from the Jensen's inequality.
Similarly for G 3 , E((X 2 ) r ) = E(X 2 ) r and E((X 1 ) r ) > E(X 1 ) r as long as E(X 1 | X 2 ) is not a constant.
Hence we can distinguish graphs G 1 , G 2 , and G 3 by testing whether a moments ratio E((X j ) r )/E(X j ) r is greater than or equal to 1.
Now we state the identifiability condition for the general case of p-variate GHD DAG models: . For a given GHD DAG model, the conditional distribution of each node given its parents is known. In other words, the r-th factorial CMR functions (f Prop. 2.2 and Assumption 2.3 enable us to use the following property: for any node j (E(X j | X Pa(j) ); a(j), b(j))), while for any non-empty Pa 0 (j) ⊂ Pa(j) and S j ⊂ Nd(j) \ Pa 0 (j), because the CMR function is strictly convex.
We state the first main result that general p-variate GHD DAG models are identifiable:

Theorem 2.4 (Identifiability). Under Assumption 2.3, the class of GHD DAG models is identifiable.
We defer the proof in Appendix. The key idea of the identifiability is to search a smallest conditioning set S j for each node j such that the moments ratio E(( j (E(X j | X S j ))) = 1. Thm. 2.4 claims that the assumption on nodes distributions is sufficient to uniquely identify GHD DAG models. In other words, the well-known assumptions such as faithfulness, non-linear causal relation, non-Gaussian additive noise assumptions are not necessary [7,11,14,16,17].

Algorithm
In this section, we present our Moments Ratio Scoring (MRS) algorithm for learning GHD DAG models.
Our MRS algorithm has two main steps: 1) identifying the skeleton (i.e., edges without their directions) using existing skeleton learning algorithms; and 2) estimating the ordering of the DAG using moments ratio scores, and assign the directions to the estimated skeleton based on the estimated ordering.
Although GHD DAG models can be recovered only using the r-th CMR property according to Thm. 2.4, our algorithm exploits the skeleton to reduce the search space of DAGs. From the idea of constraining the search, our algorithm achieves computational and statistical improvements. More precisely, Step 1) provides candidate parents set for each node. The concept of candidate parents set exploits two properties; (i) the neighborhood of a node j in the graph denoted by N (j) := {k ∈ V | (j, k) or (k, j) ∈ E} is a superset of its parents, and (ii) a node j should appear later than its parents in the ordering. Hence, the candidate parents set for a node j is the intersection of its neighborhood and elements of the ordering which appear before that node j, and is denoted by C mj := N (j) ∩ {π 1 , π 2 , ..., π m−1 } where m th element of the ordering is j (i.e., π m = j). The estimated candidate parents set is C mj := N (j) ∩ { π 1 , π 2 , ..., π m−1 } that is specified in This candidate parents set is used as a conditioning set for a moments ratio score in Step 2). If the candidate parents set is not applied, the size of the conditioning set for a moments ratio score could be p − 1.

Since
Step 2) computes the r-th factorial moments, the sample complexities depends significantly on the number of variables we condition on as illustrated in Sections 3.2. Therefore by making the conditioning set for a moments ratio score of each node as small as possible, we gain huge statistical improvements.
The idea of reducing the search space of DAGs has been studied in many sparse candidate algorithms.
Hence for Step 1) of our algorithm, any off-the-shelf candidate parents set learning algorithms can be applied such as the SC [4] and MMPC [19] algorithms. Moreover, any standard MEC learning algorithms such as GES ( [2]), PC [18] and MMHC [20] can be exploited because MEC provides the skeleton of a graph [22].
In Section 4, we provide the simulation results of the MRS algorithm where GES and MMHC algorithms are applied in Step 1).
Step 2) of the MRS algorithm involves learning the ordering by comparing moments ratio scores of nodes using Eqn. (1). The basic idea is to test which nodes have moment ratio score 1. The ordering is determined one node at a time by selecting the node with the smallest moment ratio score because the correct element of the ordering has the score 1, otherwise strictly greater than 1 in population.
Regarding the moments ratio scores, the score can be exploited for recovering the ordering only if the CMR property holds, which implies that the score should not be zero. Even if the zero value score is impossible in population, zero value scores often arise for a low count data such that all samples are less than r. Hence in order to avoid zero value scores due to a sample r-th factorial moment (i.e., E((X) r ) = 0), is Stirling numbers of the first kind. This alternative ratio score comes from (

Hence the moments ratio scores in
Step 2) of Alg.1 involve the following equations: , and (1) where C mj is the estimated candidate parents set of node j for the m th element of the ordering. In addition, , that refers to the truncated conditional sample size for x S , and n S := x S n(x S ) refers to the total truncated conditional sample size for variables X S . We discuss the choice of N min later in Sec 3.2. Lastly, we use the method of moments estimators Since there are many conditional distributions, our moments ratio score is the weighted average of the levels of how well each distribution satisfies the r-th CMR property. The score only contains the conditional expectations with n(x S ) ≥ N min for better accuracy because the accuracy of the estimation of a conditional expectation E(X j | x S ) relies on the sample size.
Finally, a directed graph is estimated combining the estimated skeleton from Step 1) and the estimated

Computational Complexity
The MRS algorithm uses any skeleton learning algorithms with known computational complexity for Step 1). Hence we first focus on our novel Step 2) of the MRS algorithm. In Step 2), there are (p − 1) iterations and each iteration has a number of moments ratio scores to be computed which is bounded by O(p). Hence Algorithm 1: Moments Ratio Scoring Input : n i.i.d. samples, X 1:n Output: Estimated ordering π and an edge structure, E ∈ V × V Step 1: Estimate the skeleton of the graph; Step 2: Estimate an ordering of the graph using r-th moments ratio scores; Calculate r-th moments ratio scores S r (m, j) using Eqn. (1); end The m th element of the ordering π m = arg min j S(m, j);

end
The last element of the ordering π p = {1, 2, · · · , p} \ { π 1 , π 2 , · · · , π p−1 }; Return: Estimate the edge sets: the total number of scores to be calculated is O(p 2 ). The computation time of each score is proportional to the sample size n, the complexity is O(np 2 ).
The total computational complexity of the MRS algorithm depends on the choice of the algorithm in Step 1). Since learning a DAG model is NP-hard [1], many state-of-the-art DAG learning algorithms such as GES [2], GDS [14,15], PC [18], and MMHC [20] are inherently heuristic algorithms. Although these algorithms take greedy search strategies, the computational complexities of greedy search based GES and MMHC algorithms are empirically O(n 2 p 2 ). In addition, PC algorithm runs in the worst case in exponential time. Hence, Step 2) may not the main computational bottleneck of the MRS algorithm. In Section 4, we compare the MRS to GES algorithm in terms of log run-time, and show that the addition of estimation of ordering does not significantly add to the computational bottleneck.

Statistical Guarantees
The MRS algorithm exploits well-studied existing algorithms for Step 1). Hence, we focus on theoretical guarantees for Step 2) of the MRS algorithm given that the skeleton is correctly estimated in Step 1). The main result is expressed in terms of the triple (n, p, d) where n is a sample size, p is a graph node size, and d is the indegree of a graph. Lastly, we discuss the sufficient conditions for recovering the graph via the MRS algorithm according to the chosen skeleton learning algorithm for Step 1).
We begin by discussing three required conditions that the MRS algorithm recovers the ordering of a graph.
is the predefined minimum sample size in the MRS algorithm.
The first condition is a stronger version of Assumption 2.3 since we move from the population to the finite sample setting. The second assumption is to control the tail behavior of the conditional distribution of each variable given its parents. It enables to control the accuracy of moments ratio scores in Eqn. (1) in high dimensional settings (p > n). The last assumption ensures that the score can be calculated.
We now state the second main result under Assumption 3.1. Since the true ordering π is possibly not unique, we use E(π) to denote the set of all the orderings that are consistent with the DAG G. for any > 0 such that if sample size is sufficiently large n > C log 2r+d (max (n, p))(log(p) + log(r)), the MRS algorithm with the r-th moments ratio scores recovers the ordering with high probability: The detail of the proof is provided in Appendix. Intuitively, it makes sense because the method of moment estimator converges to the true moment as sample size n increases by weak law of large number.
This allows the algorithm to recover a true ordering for the DAG consistently.
Thm. 3.2 claims that if the sample size n = Ω(log 2r+d (max (n, p)) log(p)), our MRS algorithm accurately estimates a true ordering with high probability. Hence our MRS algorithm works in high-dimensional settings (p > n) provided that the indegree of the graph d is bounded. This theoretical result is also consistent with learning Poisson DAG models shown in [12] where if n = Ω(log 4+d (max (n, p)) log(p)) their ODS algorithm recovers the ordering well. Since [12] uses the variance (the second order moments r = 2), both algorithms are expected to have the same performance of recovering graphs.
However the MRS algorithm performs significantly better than the ODS algorithm in general because the moments difference the ODS algorithm exploits is proportional to magnitude of the conditional mean while the moments ratio is not. For a simple Poisson DAG X 1 → X 2 , E((X 2 ) 2 ) − E(X 2 ) 2 = Var(E(X 2 | X 1 )). Hence if E(X 2 | X 1 ) ≈ 0, the score in ODS is inevitably close to 0, while the score in MRS, E((X 2 ) 2 )/E(X 2 ) 2 = 1 + Var(E(X 2 | X 1 ))/E(X 2 ) 2 is not necessarily close to 1. Hence, Assumption 3.1(A1) is much milder than the related assumption for the ODS algorithm. Now we discuss the sufficient conditions for recovering the true graph via the MRS algorithm according to the choice of the algorithm in Step 1). The GES [2], PC [18], and MMHC [20] algorithms require the Markov, faithfulness, and causal sufficiency or related assumptions to recover the skeleton of a graph.
Moreover GES and MMHC algorithms are greedy search based algorithms that are not guaranteed to recover the true skeleton of a graph. Therefore, the MRS algorithm may require strong assumptions or large sample size to recover the true graph based on the choice of the algorithm in Step 1). Although these assumptions can be very restrictive, we show through the simulations that the MRS algorithm can recover the directed edges well even in high dimensional settings (p > n).

Numerical Experiments
In this section, we provide various simulation results by comparing the MRS algorithm to state-of-the-art ODS, GES and MMHC algorithms in terms of recovering MECs and DAGs. Also, we support our theoretical results in Thm. 3.2 and computational complexity in Section 3.1 with synthetic and real basketball data. In addition, we show that our algorithm performs favorably compared to the ODS, GES, and MMHC algorithms in terms of recovering the directed graphs.

Synthetic Data
We conduct two sets of simulation study using 150 realizations of p-node random GHD DAG models with the indegree constraints d = 2: (1) Poisson DAG models where the conditional distribution of each node given its parents is Poisson; and (2) Hybrid DAG models where the conditional distributions are sequentially Poisson, Binomial with N = 3, hyper-Poisson with b = 2, and Binomial with N = 3.
We set the (hyper) Poisson rate parameter θ j (Pa(j)) = exp(θ j + k∈Pa(j) θ jk X k ) and the binomial probability p j (Pa(j)) = logit −1 (θ j + k∈Pa(j) θ jk X k ). For Poisson DAG models, the set of non-zero parameters θ jk ∈ R were generated uniformly at random in the range θ jk ∈ [−1.  provides that the MRS algorithm is more accurate than state-of-the-art ODS, GES and MMHC algorithms in both precision and recall. It makes sense because the moments ratio scores the MRS algorithm exploits are less sensitive to the magnitude of the moments than the score the ODS algorithm uses as discussed in Section 3.2, and because the GES and MMHC algorithms recover up to the MEC by leaving some arrows undirected. However it must be pointed out that our MRS algorithm apply to GHD DAG models while GES and MMHC apply to general classes of DAG models.
Not surprisingly, the MRS algorithm performs better than ODS, GES and MMHC algorithms in terms of recovering MECs since the ordering is well estimated by the MRS algorithm.

Computational Complexity
To validate the computational complexity discussed in Section 3.1, we show the log run-time of Step 1) and Step 2) of the MRS algorithm in Fig. 6 where the GES is applied for Step 1). We measured the run-time for

Deviations to the Distribution Assumptions
When the data are generated by a GHD DAG model where the conditional distribution of each node given its parents is unknown, our algorithm is not guaranteed to estimate the true graph and its ordering. Therefore, an important question is how well the MRS algorithm recovers graphs when incorrect distributions are used.
In this section, we heuristically investigate this question.
We use the same setting of the data generation for Poisson and Hybrid GHD DAG models with the node size p = 20. We consider (i) the true (conditional) distributions, and assume all nodes (conditional) distri-

Real Multi-variate Count Data: 2009/2010 NBA Player Statistics
We demonstrate the advantages of our graphical models for count-valued data by learning a 441 NBA player statistics from season 2009/2010 (see R package SportsAnalytics for detailed information) because we believe that some basketball statistics have causal or directional relationships. The original data set contains 24 covariates: player name, team name, players position, total minutes played, total number of field goals made, field goals attempted, threes made, threes attempted, free throws made, free throws attempted, offensive rebounds, rebounds, assists, steals, turnovers, blocks, personal fouls, disqualifications, technicals fouls, ejections, flagrant fouls, games started and total points. We eliminated player name, team name, number of games played, and players position, because our focus is to find the directional or causal relationships between statistics. We also eliminated ejections and flagrant fouls because both did not occur in our data set. Therefore the data set we consider contains 18 discrete variables.  played is likely to be positively correlated with other statistics, and some statistics have causal or directional relationships (e.g., the more shooting attempt implies the more shooting made).
The MRS and ODS algorithms are applied where GES algorithm is used in Step 1). We assume that the conditional distribution of each node given its parents is hyper-Poisson because most of NBA statistics we consider are the number of successes or attempts counted in the season. We emphasize that our method requires a known conditional distribution assumption to recover the true graph. However since we do not have prior node distribution information, we set b j = Var(X j )/ E(X j ) as we used in simulations that enables the MRS algorithm successfully recovers the directed edges. Fig. 9 show the estimated directed graph using the MRS and ODS algorithm. We provide the differences between the estimated DAGs in Proof. For any positive integer r ≥ 1, [8] shows that Then, the expectation can be obtained when r = 1.
By plugging this into Eqn.
As we discussed in the main body of the paper, for any node j ∈ V \{π 1 }, Prop. 2.2 and Assumption 2.3 Hence we can determine π 1 as the first element of the ordering.
For the (m − 1) th element of the ordering, assume that the first m − 1 elements of the ordering and their parents are correctly estimated. Now, we consider the m th element of the ordering and its parents. Again Prop. 2.2 and Assumption 2.3 yield that forj ∈ {π m+1 , π m+2 , · · · , π p }, j (E(X j | X Pa(j) )) | X 1:(m−1) )) = 1.
Hence we can estimate a valid m th component of the ordering π m and its parents by testing whether the r-th moments ratio is whether greater than or equal to 1. By the mathematical induction this completes the proof.

Proof for Theorem 3.2
Proof. We first reintroduce some necessary notations and definitions to make the proof concise. Without loss of generality, assume that the true ordering is unique and π = (π 1 , ..., π p ). In addition, we assume the true skeleton is provided. For ease of notation, we drop the r in the both r-th moments ratio scores and CMR function. Then, the element of the score can be written as: , and .
where C jk is the candidate parents set and s(r, k) is Stirling numbers of the first kind. Hence the r-th moments ratio score is We define the following important events: For each node j ∈ V and set S ⊂ V \ (De(j) ∪ {j}) and for any > 0, let ζ 1 := min j=1,...,p−1 min k=j,...,p S(j, π k ) − S(j, π k ) > M min 2 Here we use the method of moments estimators 1 n n i=1 (X (i) j ) k as unbiased estimators for E(X k j ) for all 1 ≤ k ≤ r.
By Assumption 3.1 (A1) S(j, π k ) > 1 + M min , the above lower bound of the probability is reduced to Next we introduce the following three propositions to show that the above lower bound converges to 1.
The first proposition proves that estimated score is accurate under some regularity conditions. For ease of notation, let (2E(X r j |X S )+(2−M min )g j (E(X j |X S ))) . The first condition (i) is satisfied if M min in Assumption 3.1 (A1) is sufficiently large and the second condition (ii) is satisfied if is sufficiently small. This means that if the estimated r-th factorial moment is close to the true r-th factorial moment, then ζ 1 is not satisfied with probability 1. Hence we discuss the error bound for the r-th factorial moment estimator in the next.
The following propositions show the error bound for the higher order moment X k j for 1 ≤ k ≤ r given the set ζ 5 , and therefore the error bound for the r-th factorial moment estimator: Proposition 5.2. For any node j ∈ V and any set S ⊂ V \ (De(j) ∪ {j}) and for any > 0, (ii) For ζ 3 , and m ∈ (E(X j | X S ) − /2, E(X j | X S ) + /2), P (ζ c 3 | ζ 5 ) ≤ 2 · p · exp − N min 2 8(max(f j (m))) 2 log 2 η .
Lastly, we show the relationship between the sample size n and N min to satisfy Assumption 3.1 (A3).
Suppose that d is the maximum number of parents of a node. Then the maximum size of the candidate parents set is d. The scenario is that a conditioning set has two possible cases. If there is only one element for the conditioning set, there is no difference between the conditional and marginal distributions.
In the best case, n = 2N min when the there are two conditional distributions |X C | = 2. Hence if n = O((log 2r (η)(log(p) + log(r))), our algorithm works in the high dimensional settings. In the worst case given ζ 5 , the sample size is n = (4 log(η) d − 2)(N min − 1) + 2N min = 4 log(η) d (N min − 1) + 2 where the number of all elements of {x ∈ X C | n i 1(X (i) C = x) ≥ N min } is two and all other elements of X C has N min − 1 repetitions. In this worst case, if n = O(log(η) (2r+d) (log(p) + log(r))) our algorithm correctly recovers the ordering with high probability.

Proof for Proposition 5.1
Proof. For ease of notation, let η = max{n, p} and the r-th moments ratio score: S(j, k) := x∈X C jk n(x) n C jk S r (j, k)(x).
In addition, let s(r, m) E(X m k | X C jk ).

Proof for Proposition 5.2
Since the proof for Prop. 5.2 (i) -(iii) are analogous, we provide the proof for (iii) and then we provide the proof for (ii).
Now we provide the proof for (ii).
Proof. By Mean value theorem, we obtain f j ( E(X j | X S )) − f j (E(X j | X S )) = f j (m) E(X j | X S ) − E(X j | X S ) .
where f j is the first derivative of f j andm is some point between E(X j | X S ) and E(X j | X S ).