{'fraction_non_alphanumeric': 0.053102161405691184, 'fraction_numerical': 0.014266832529933136, 'mean_word_length': 3.9433128362797847, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Consistency of Maximum Likelihood Estimators for Causal Network Identification

We consider the problem of identifying parameters of a particular class of Markov chains, called Bernoulli Autoregressive (BAR) processes. The structure of any BAR model is encoded by a directed graph. Incoming edges to a node in the graph indicate that the state of the node at a particular time instant is influenced by the states of the corresponding parental nodes in the previous time instant. The associated edge weights determine the corresponding level of influence from each parental node. In the simplest setup, the Bernoulli parameter of a particular node's state variable is a convex combination of the parental node states in the previous time instant and an additional Bernoulli noise random variable. This paper focuses on the problem of edge weight identification using Maximum Likelihood (ML) estimation and proves that the ML estimator is strongly consistent for two variants of the BAR model. We additionally derive closed-form estimators for the aforementioned two variants and prove their strong consistency.


I. INTRODUCTION
The spreading of ideas and information, the propagation of viruses and diseases, and the fluctuation of stock prices are examples of processes evolving over social, information or other types of networks [1]- [8]. Identifying the underlying network structure in these systems motivates the so-called network inference problem, which aims at recovering the underlying connectivity between entities or nodes in the system based on observed data. The dependencies, correlations or causal relationships between network entities can be modeled as undirected or directed edges in a graph. The associated dependency strengths can be described as edge weights. Many algorithms have been proposed to identify the network structure and edge weights from time series data for various processes. Clearly, efficient algorithms in terms of sample complexity are desired.
The network inference problem for various dynamic processes has been recently studied in both the machine learning literature and the system identification literature. Among the relevant studies in the machine learning community, the socalled continuous-time independent cascade model (CICE) considered in [5], [9], [10], presents a typical model for capturing the dynamics of virus or information spreading in networks. A discrete-time version of CICE is studied in [11]. In [7], the Generalized Linear Model is formulated, which is a class of diffusion models encompassing both the Xiaotian Xie, Carolyn L. Beck  System identification is a well-studied model estimation approach in the context of control theory; e.g., see [12]- [14] and references therein. The traditional task of system identification focuses on estimating system parameters from measured input and output data. Among recent emerging applications in the system identification area, two main extensions are relevant. The first involves cases where the state matrix of a state-space model represents a directed graph, in other words, the dynamic process can be factorized into subprocesses for each node [15]- [19]. In such cases, system identification methods are found effective in tackling dynamic network inference problems. The second extension is to consider state vectors evolving in discrete state spaces.
In this paper, we consider the Bernoulli Autoregressive (BAR) model, which is a parameterized discrete-time Markov chain initially introduced in [20]. In this model, the state of each node is a Bernoulli random variable with probability of success equal to a convex combination of the parental node states (or their flipped states) in the previous time step and an additional binary noise term ensuring persistence of excitation. The BAR model can be used to approximate opinion dynamics, biological and financial times series, and similar processes [20]. Another relevant discrete-time binary process is the ALARM model proposed in [21]. In contrast to the BAR model, the ALARM model defines the transition probabilities via a logistic function.
Relying on well-established statistical principles, we first formulate and study the consistency properties of the Maximum Likelihood (ML) parameter estimator for the BAR model in which every parental node causally influences each descendant node positively; the notion of positive correlations is formalized in [20]. The consistency of ML estimators in the case of independent and identically distributed (i.i.d.) random variables has been studied extensively; see e.g., [22], [23] and references therein. The consistency of ML estimators for Markov chains appears to be less well studied, see [24] for a reference.
To establish the (strong) consistency of the ML estimator for the BAR model, we prove that the vectorized transition probability matrix is an injective mapping of the model parameters. In the rest of the paper, we call the injectivity of this mapping identifiability of the BAR model. The strong consistency of the ML estimator is then shown by leveraging the injectivity and the continuity of the transition probabilities with respect to the parameters, as well as the compactness of the parameter space. By relying on the ML principle, a closed-form estimator is subsequently provided. Strong consistency is also shown to hold for this estimator. The identifiability proof is then extended to the generic BAR model with both positive and negative correlations, where the notion of negative correlations is also formalized in [20]. This identifiability extension establishes the strong consistency of the ML estimator for the general BAR model class. The closed-form estimator and its consistency are also extended to the generic BAR model. These analytical results provide a complement to the prior work [20]. Finally, numerical simulations are provided to demonstrate the sample complexity gain achieved by the derived estimators in this paper over other existing algorithms in the literature for the BAR model when focusing on the structure identification subproblem. We note here that structural inference is an identification subproblem that can be tackled by parameter estimation in processes with underlying network structures.
The rest of this paper is organized as follows. In Section II, the BAR model with positive correlations only is introduced. In Sections III and IV, the identifiability of the BAR model with positive correlations only and the strong consistency of the corresponding ML estimator, as well as a closedform estimator and its strong consistency are derived. The generic BAR model with positive and negative correlations and also, the identifiability and strong consistency of the corresponding ML estimator are provided in Section V. The extension of the closed-form estimator to the generic BAR model and its strong consistency are derived in Section VI. Finally, simulation results are provided in Section VII and Section VIII concludes the paper.
Notation: Matrices and vectors are denoted by bold upper and lowercase letters, respectively. Probability distributions in vector form may be either denoted by bold upper or lowercase letters. Random vectors are also denoted by uppercase bold letters, while their corresponding realizations are denoted by lowercase bold letters. Scalar random variables are denoted by uppercase letters. The i-th entry of a vector x is denoted by x i . For a matrix A, a i j corresponds to its (i, j)th entry. Depending on the context, vector and matrix entries may be indexed more generally, e.g., by state elements.

II. THE BAR MODEL WITH POSITIVE CORRELATIONS
The BAR model is a special form of a Markov chain defined on a directed graph G = (V , E ) with |V | = p nodes. Let X i (k) ∈ {0, 1} be the state of node i ∈ [p] at time instant k and let X(k) ∈ {0, 1} p be the associated BAR process state vector at the same time instant. The most natural BAR model, with positive correlations only, is described by where a i ∈ [0, 1] p , b i ∈ [0, 1], i = 1, . . . , p are parameters of the BAR model and Ber(ρ) represents the Bernoulli distribution with parameter ρ. Additionally, are independent noise random variables, also independent of X(t) for any t < k + 1, where ρ w i ∈ [ρ min , ρ max ] for all i ∈ [p] with 0 < ρ min < ρ max < 1. Moreover, the initial distribution is P X(0) , i.e., X(0) ∼ P X(0) . The interpretation here is that the entries of X(k + 1) are conditionally independent Bernoulli random variables given X(k).
To ensure that the Bernoulli random variables in (1) are well-defined, we require that For persistent excitation, we further assume , the BAR Markov chain will get absorbed in 0 p or 1 p upon visiting the state 0 p or 1 p , respectively. Furthermore, we assume that a i encodes a part of the graph structure through the equivalence where the ordered pair ( j, i) denotes a directed edge from node j to node i. The notion of positive correlations in (1) relies on the fact that a i j > 0 increases the probability of the event {X i (k + 1) = 1} when X j (k) = 1. A more general form of the BAR model with both positive and negative correlations is introduced in Section V. We now let A = [a 1 , a 2 , · · · , a p ] T , i.e., a T r corresponds to the r-th row of is an irreducible and aperiodic Markov chain with finite state space {0, 1} p . Moreover, for any vectors specifies the transition probability from state u to state v. We denote by π ∈ R 2 p the associated stationary distribution with component π u corresponding to the state u ∈ {0, 1} p and by P = (p uv ) ∈ R 2 p ×2 p the BAR transition probability matrix. The goal is to recover the model parameters from an observed sequence {X(k) = x(k)} T k=0 . Clearly, by inferring A, estimates of b and the underlying network structure are direct per (2) and (3), respectively.

III. MAXIMUM LIKELIHOOD ESTIMATION
In this section, we consider recovering the BAR model parameters via ML estimation and we establish the strong consistency of the ML estimator. Suppose that {x(k)} T k=0 is a sequence of observations generated by the BAR model (1). Let θ = (A, b, ρ w ) with the implicit relationship b = 1 p − A1 p . Clearly, b is a redundant parameter, but it is preserved here to facilitate the subsequent analysis. From (4) the rescaled log-likelihood function is given by In the rest of the paper, we assume that P X(0) is independent of the model parameters, which is well-aligned with the realistic scenario of arbitrarily initializing the Markov chain.
For any states u, v ∈ {0, 1} p , we denote by N uv the number of one-step transitions from state u to state v in the observed sequence and we let N u = ∑ v N uv be the amount of time spent in state u in a horizon of T time steps. Then (5) can be also written as Let θ 0 be the true parameter tuple (A, b, ρ w ). An application of the Ergodic Theorem [25] for Markov chains reveals that which is the negative of the entropy rate of the corresponding BAR chain with parameter set θ 0 and is always finite since the BAR model has a finite state space. Let θ 0 ∈ Θ, where Θ is a compact set appropriately defined on the basis of (2) and additional assumptions in Section II. The ML estimatorθ T of θ 0 satisfieŝ In the rest of this section, we will show the strong consistency ofθ T . The key idea of the proof is along the lines of the proof of Theorem 2.1 in [24] using techniques for general discrete-time Markov chains. To summarize, we first prove Here, vec(·) denotes the vectorization of a matrix. Finally, we complete the proof by leveraging the compactness of the parameter space Θ and the continuity of the components of p(θ ) = vec(P(θ )) or equivalently, of the transition probabilities with respect to the model parameters.
Remark: In the following, we will say that the BAR model is identifiable when (8) holds.
The main result of this section can be now stated. Theorem 1: The ML estimatorθ T of θ 0 , defined in (7), for the BAR model in (1) is strongly consistent.
Proof: We break up the proof into three parts. Proof of P(θ T ) a.s.
For each u ∈ {0, 1} p , we define the (row) vector Q u = (N uv /N u ) v∈{0,1} p ∈ R 2 p with the convention Q u = 2 −p 1 2 p for N u = 0 and we let P u = (P uv ) v∈{0,1} p ∈ R 2 p denote the transition distribution out of state u, which is also a row vector in the transition matrix P. In particular, it is well-known that the set {Q u } u∈{0,1} p is the ML estimator of the transition matrix P, assuming no further parameterization of the transition probabilities. Consider the fact that the Kullback-Leibler divergence from P u (θ T ) to Q u is nonnegative, i.e., Multiply both sides of the above inequality by N u T and sum over u to obtain where the last inequality is due to (6) and the definition of the ML estimator. From (9), we can further obtain By the Ergodic Theorem for Markov chains, By (11) and (12) we have that This together with (10) yields Employing Pinsker's inequality [26] and the fact that the total variation distance between two discrete measures q, r with vector forms q, r, respectively, is q − r TV = (1/2) q − r 1 , we have that for each u ∈ {0, 1} p Multiplying again with N u T and summing over u gives Now employing again the Ergodic Theorem for Markov chains, i.e., and combining (13) and (14) yields We then end up with Consider now the transition probability from u = 0 to v , where v 1 = 1 and v j = v j for j = 2, . . . , p. Since p uv (θ ) = p uv (θ ), Combining (16) and (17), it is easy to see that c 1 = c 1 , which is a contradiction. Thus, c = c =⇒ P(θ ) = P(θ ) or equivalently, p(θ ) = p(θ ). Now we consider the second case where c = c and A = A . Without loss of generality, let a 11 = a 11 . Consider u = e p,1 and the same v, v as before. By our assumption that p u v (θ ) = p u v (θ ) and p u v (θ ) = p u v (θ ), the contradiction a 11 = a 11 arises. Thus, c = c , A = A =⇒ P(θ ) = P(θ ) or equivalently, p(θ ) = p(θ ).

IV. A CLOSED-FORM ESTIMATOR
In this section, we provide a closed-form estimator for the BAR model parameters in (1), i.e., for (A, c), since by knowing (A, c) we can recover (b, ρ w ). Recall that c = diag(b)ρ w . Considering the log-likelihood function for θ = (A, c), we have Observe that P (v r = 1|u) and P (v r = 0|u) are independent of v. We can therefore define ϑ u,r,l = P ((·) r = l|u), for l ∈ {0, 1}. Furthermore, we define N u,r,l = ∑ T −1 k=0 I (x(k) = u, x r (k + 1) = l), which is the number of times the BAR chain transitions from state u to a state with r-th entry being equal to l. Moreover, N u,r,0 + N u,r,1 = N u = ∑ v N uv , ∀u ∈ {0, 1} p and r ∈ [p]. With these introductions we have the following theorem: .., u m such that p ≤ m ≤ 2 p . Let U m ∈ R m×p be a matrix with k-th row equal to u k for k ∈ [m] and y m,r = [N u 1 ,r,1 /N u 1 , · · · , N u m ,r,1 /N u m ] . Then, whenever U m is full-column rank,Â is an estimator of A, wherê Finally, to obtain a valid estimate of the parameter set for any T ≥ 1, we let where [·] + corresponds to a projection onto the parameter space Θ. Proof: First, let us rewrite the log-likelihood function in (18) as (N u,r,0 log ϑ u,r,0 + N u,r,1 log ϑ u,r,1 ) Instead of maximizing this function with respect to θ = (A, c), we maximize it with respect to the choice of the marginal conditional probabilities {ϑ u,r,0 , ϑ u,r,1 } u,r . Consider the constrained ML estimation problem Forming the Lagrangian and setting the gradient (with respect to {ϑ u,r,0 , ϑ u,r,1 } u,r ) to zero, we obtain Recall that ϑ u,r,1 is defined as the probability of transitioning from state u to some state with r-th component equal to 1. We can therefore require that or y m,r −ĉ r 1 m = U mâr , where the estimatesĉ r are provided in the statement of the theorem. Under the assumption that U m is full-column rank, U m U m is nonsingular and (19) follows. The proof is then concluded by projecting the obtained estimate onto Θ. Remark: (21) is reminiscent of the invariance property for ML estimation.
Theorem 3: The closed-form estimator in Theorem 2 is strongly consistent.
Proof: It is sufficient to show that (Â,ĉ) given by Theorem 2 is strongly consistent. Since the BAR chain is finitestate, the stationary probabilities satisfy π u > 0, ∀u ∈ {0, 1} p . Moreover, for any initial measure, the Ergodic Theorem for Markov chains implies that 1  1 The following convergences can be easily justified by straightforward embeddings of U m into R 2 p ×p in the first case and of y m,r into R 2 p in the second case via zero padding.
Combining these observations with (22) we obtain s. (23) and the strong consistency of the closed-form estimator in Theorem 2 follows if U 2 p is full-column rank or equivalently if U 2 p U 2 p is nonsingular. It is easy to see that U 2 p U 2 p has diagonal entries equal to 2 p−1 and off-diagonal entries equal to 2 p−2 . Thus, we can write U 2 p U 2 p = 2 p−2 1 p 1 p + 2 p−2 I p . This matrix is invertible for every p < ∞, since 1 + 2 p−2 1 p 2 p−2 I p −1 1 p = 1 + p > 0 as the condition in the Sherman-Morrison formula [27] dictates.

V. THE GENERIC BAR MODEL
Motivated by modeling positive and negative influences from parental nodes, an extension of the BAR model in (1) has been introduced in [20]. We first reformulate this generic BAR model.
. The nodes in S + i and S − i are said to have positive and negative influence on i, respectively. The generic BAR model, parameterized byθ = A,Ã, b, ρ w , is defined as for all i ∈ [p], where a i andã i are the i-th rows of A ∈ R p×p andÃ ∈ R p×p , respectively. Furthermore, we assume that S + i = supp(a i ) and S − i = supp(ã i ). Here, supp(·) denotes the support of a vector. As in the previous case, the constraints are also required in this case. Similarly, we assume that . Therefore, the parameter space is defined as , and a i j ,ã i j ≥ 0, Remark: The parameter spaceΘ is compact. To see this, first note that the associated constraints imposed onΘ are defined row-wise. This implies thatΘ can be viewed as the Cartesian product of row spaces, i.e.,Θ = ∏ p i=1Θ i , wherẽ ThenΘ is compact if and only ifΘ i is compact for every i ∈ [p]. It is straightforward to see thatΘ i is compact without the constraints a i jãi j = 0, ∀ j ∈ [p]. Also observe that adding a constraint a i jãi j = 0 leads to the coordinate projection of a compact set in R 2p+2 onto two (2p + 1)-dimensional subspaces of R 2p+2 with the corresponding images of the projected compact set being also compact sets. Denote the union of these images asΘ i, j . ThenΘ i = ∩ p j=1Θ i, j , i.e.,Θ i is the intersection of p compact sets and is therefore compact.
The ML estimator is a maximizer of the rescaled loglikelihood function, that is, The ML estimator for the generic BAR model can be shown to be strongly consistent via a direct extension of the analysis in Section III. More precisely, it is sufficient to establish identifiability. Theorem 4: For the generic BAR model in (24),θ = θ =⇒ p(θ ) = p(θ ), ∀(θ ,θ ) ∈Θ ×Θ withθ =θ . Proof: For two different sets of parameters (A,Ã, b, ρ w ), (A ,Ã , b , ρ w ) and by letting c i = b i ρ w i , i = 1, . . . , p, we can reparameterize the generic BAR model to obtainθ = (A,Ã, c),θ = (A ,Ã , c ) by recalling the oneto-one correspondence between the initial sets of parameters and the later ones as in the case of the BAR model with only positive correlations. The corresponding relations in this case are b = 1 p − (A +Ã)1 p and ρ w = (diag(b)) −1 c. We will examine different cases for whichθ =θ and by assuming that p(θ ) = p(θ ) we will arrive to contradictions.
• First Case: Suppose that there exists some i ∈ [p] such that S + i = S + i , where S + i and S + i correspond to the parental neighborhoods with positive influence on the ith node, as these neighborhoods are encoded inθ and θ , respectively. Without loss of generality, we assume that S + 1 = S + 1 . Translating this structural difference into the model parameters, suppose again without loss of generality that 1 ∈ S + 1 , i.e., a 11 = 0 and therefore, a 11 = 0 and 1 / ∈ S + 1 , i.e., a 11 = 0 andã 11 = 0 or a 11 > 0.
We first consider the transition probabilities from u = ∑ i∈S − 1 e p,i to states v and v , where v and v only differ in the first element with v 1 = 1 and v 1 = 0. Then by letting p uv (θ ) = p uv (θ ) and p uv (θ ) = p uv (θ ), we obtain Further by considering the probabilities of transitioning from u = e p,1 + ∑ i∈S − 1 e p,i to v and v , we obtain where the assumption a 11 = 0 has been used. By (26) and (27) we have that a 11 +ã 11 = 0, which implies that a 11 = 0 and contradicts our assumption. The case of , c = c and A,Ã = A ,Ã . Without loss of generality, we assume that a 11 = 0, a 11 = 0 and a 11 = a 11 . In this case, the contradiction a 11 = a 11 arises when selecting u = e p,1 + ∑ i∈S − 1 e p,i and v, v as before.

VI. A CLOSED-FORM ESTIMATOR FOR THE GENERIC BAR MODEL
Extending the closed-form estimator for the BAR model with only positive correlations, we now introduce a closedform estimator for the generic BAR model. Consider the same introductions as before and also the Bernoulli argument for X i (k + 1) in (24). We can rewrite and we note that due to the nonoverlapping supports of a i and a i for every i ∈ [p], the vector a i −ã i contains the entries of a i and the entries ofã i with flipped signs, each at a different location. We further note that . With these introductions, we can reparameterize the generic BAR model asθ = (Ā,c), whereĀ = [ā 1 , . . . ,ā p ] T andc = [c 1 , . . . ,c p ] T . Clearly, by knowingĀ, we can immediately separate A andÃ based on the underlying signs of the entries. Furthermore, by also knowingc, we can compute the products c i = b i ρ w i , i = 1, . . . , p and finally, we can compute ρ w i via the formula ρ w i = c i /(1 − ∑ p j=1 (a i j +ã i j )) for i ∈ [p]. The parameterizationθ = (Ā,c) is of the same form as the parameterization θ = (A, c) of the BAR model with only positive correlations and therefore, the same unprojected closed-form estimator as before is suitable.
The proofs of this theorem and of the strong consistency of the proposed closed-form estimator are straightforward based on the proofs of Theorems 2 and 3, respectively.

VII. SIMULATION RESULTS
To validate our analysis, experiments on synthetic networks of various sizes are performed. Focusing only on the structure identification of the underlying BAR graphs, we compare the ML estimators with the closed-form estimators and the BAR structure observer (BARobs) proposed in [20]. Our simulations show that the ML estimator either for the BAR model with positive correlations only or for the generic BAR model can perfectly recover the graph structure with sufficient data and outperforms the other estimators in terms of sample complexity.
We implement all experiments in MATLAB. In particular for the ML estimators, we use the "fmincon" function in the "optimization toolbox" to solve the maximization problems. The ground truth networks (A, b, ρ w ) or A,Ã, b, ρ w are randomly generated in such a way that all constraints specifying the parameter spaces Θ orΘ are satisfied. Moreover, to facilitate the simulation, we create networks with minimum edge weight a min ∈ (0, 1), i.e., a i j ≥ a min , ∀( j, i) ∈ E in the case of the BAR model (1) and either a i j ≥ a min orã i j ≥ a min , ∀( j, i) ∈ E in the case of the generic BAR model (24). The data {x(k)} T k=0 are generated according to the created ground truth models and then are used for network inference. We use the F 1 score, defined by  as the criterion to evaluate the performance of the algorithms. We note that recall is the fraction of correctly recovered edges among the ground truth edges and precision is the fraction of correctly recovered edges over all correct and incorrect edges identified by the algorithms. More specifically, an edge ( j, i) is viewed as being inferred ifâ i j ≥ ca min or a i j ≥ ca min for some empirically selected c ∈ (0, 1).
In Fig. 1, we compare the performance of ML and closedform estimators for the BAR model (1) with p = 10 nodes and maximum in-degree d max = 5. In Fig. 2, the performance of the ML, closed-form and BARobs estimators is demonstrated for a synthetic network with size p = 20 and maximum in-degree d max = 5 corresponding to the generic BAR model (24). The x-axis corresponds to the number of observations T and the y-axis to the F 1 score. In both cases, the ML estimator can achieve F 1 scores equal to 1 for a sufficient sample size, e.g., T = 1200 in Fig. 1 and T = 4000 in Fig. 2. In Fig. 2, the BARobs algorithm also achieves an F 1 score equal to 1 for a sufficiently large sample size. However, compared to the ML estimator, the BARobs requires a larger sample size. Finally, in Fig. 2 both the ML and BARobs estimators outperform the closed-form estimator in (30).
Finally, focusing on the ML estimator which outperforms the other two estimators in terms of sample complexity, we present a real network experiment in a biological application in which small sample sizes are critical [28]. The underlying biological network consists of 43 nodes. The state of each node is binary and is updated according to a boolean rule defined by the states of some nodes in the network. We approximate the network by the generic BAR model (24) and we generate pseudo-real data. Fig. 3 illustrates the performance of the ML estimator on this data set.

VIII. CONCLUSIONS
In this paper, we studied the problem of estimating the parameters of a class of Markov chains called BAR models. ML estimation for BAR chains was shown to be strongly consistent. Strong consistency was also established for closed-form estimators of these BAR models.