Bayesian Networks for Max-linear Models

TL;DR

This paper explores Bayesian networks based on max-linear models, analyzing their independence properties, estimation methods, and the potential for structure identification from observational data.

Contribution

It provides a comprehensive summary of independence properties, discusses generalized maximum likelihood estimation, and shows that minimal network structures can be identified asymptotically.

Findings

Distributions are generally not faithful to the independence model.

Generalized maximum likelihood estimation is applicable for coefficient estimation.

Minimal network structures can be asymptotically identified from observational data.

Abstract

We study Bayesian networks based on max-linear structural equations as introduced in Gissibl and Kl\"uppelberg [16] and provide a summary of their independence properties. In particular we emphasize that distributions for such networks are generally not faithful to the independence model determined by their associated directed acyclic graph. In addition, we consider some of the basic issues of estimation and discuss generalized maximum likelihood estimation of the coefficients, using the concept of a generalized likelihood ratio for non-dominated families as introduced by Kiefer and Wolfowitz [21]. Finally we argue that the structure of a minimal network asymptotically can be identified completely from observational data.

Equations

$$X_1=Z_1,X_2=\max (c_{21}{X}_1,Z_2),X_3=\max (c_{31}{X}_1,Z_3),\max (c_{42}{X}_2,c_{43}{X}_3,Z_4),$$

$$X_1$$

$$X_2$$

$$X_3$$

$$X_4$$

$$=c_{42}(c_{21}{Z}_1\vee Z_2)\vee c_{43}(c_{31}{Z}_1\vee Z_3)\vee Z_4$$

$$=(c_{42}{c}_{21}\vee c_{43}{c}_{31})Z_1\vee c_{42}{Z}_2\vee c_{43}{Z}_3\vee Z_4.$$

$$B=\left(\begin{array}{cccc}1& 0& 0& 0\\ c_{21}& 1& 0& 0\\ c_{31}& 0& 1& 0\\ c_{42}{c}_{21}\vee c_{43}{c}_{31}& c_{42}& c_{43}& 1\end{array}\right),$$

$$X_v=\underset{u\in \mathrm{pa}(v)}{\bigvee }{c}_{vu}{X}_k\vee c_{vv}{Z}_v,v=1,\dots ,d,$$

$$P(X_v\le x|X_{\mathrm{pa}(v)})=F(x|x_{\mathrm{pa}}(v)).$$

$$X_v=g_v(X_{\mathrm{pa}(v)},Z_v),v\in V,$$

$$X_v=\sum _{u\in \mathrm{pa}(v)}{c}_{vu}{X}_u+c_{vv}{Z}_v,v\in V,$$

$$\sum _{u\in \mathrm{pa}(v)}{c}_{vu}{X}_u+c_{vv}{Z}_v\approx \underset{u\in \mathrm{pa}(v)}{\bigvee }{c}_{vu}{X}_u\vee c_{vv}{Z}_v$$

$$A\perp \perp B|C⟺X_A{\perp \perp }_P{X}_B|X_C$$

$$v\perp \perp \mathrm{pr}(v)|\mathrm{pa}(v)$$

$$v\perp \perp (\mathrm{nd}(v)\setminus \mathrm{pa}(v))|\mathrm{pa}(v).$$

$$A{\perp }_{\mathcal{D}}B|C⟹A\perp \perp B|C.$$

$$\text{mbox}(\mathrm{O})⟺\text{mbox}(\mathrm{L})⟺\text{mbox}(\mathrm{G}).$$

$$A{\perp }_{\mathcal{D}}B|C⟺A{\perp \perp }_{P}B|C.$$

$$X_v=\underset{u\in \mathrm{pa}(v)}{\bigvee }{c}_{vu}{X}_u\vee c_{vv}{Z}_v,v=1,\dots ,d,$$

$$d_{vu}(\pi ):=\prod _{l=0}^{n-1}{c}_{k_{l+1}{k}_l}\text{mbox}and{b}_{vu}:=\underset{\pi \in {\Pi }_{uv}}{\bigvee }{d}_{vu}(\pi ),$$

$$\displaystyle b_{vu}=\bigvee\limits_{\pi\in\Pi_{uv}}d_{vu}(\pi)\mbox{ for $u\in{{\rm an}}(v)$};\,\,b_{vv}=c_{vv}=1;\,\,\mbox{$b_{vu}=0$ for $u\in{V\setminus{\rm An}(v)}$},$$

$$X_v=\underset{u\in \mathrm{An}(v)}{\bigvee }{b}_{vu}{Z}_u,v=1,\dots ,d.$$

$$\displaystyle(F\odot G)_{vu}=\Big{(}\bigvee\limits_{k=1}^{n}f_{vk}g_{ku}\Big{)}.$$

$$\displaystyle X=B\odot Z=\big{(}\bigvee_{u=1}^{d}b_{vu}Z_{j},i=1,\ldots,d\big{)}=\big{(}\bigvee_{u\in{\rm An}(v)}b_{vu}Z_{j},i=1,\ldots,d\big{)}.$$

$$B=\underset{k=0}{\overset{d-1}{\bigvee }}{C}^{\odot k}=C^{\odot (d-1)},$$

$$r_{vu}:=\{\begin{array}{ll}1,& \text{if there is a path from u to v, or if u=v},\\ 0,& \text{otherwise}.\end{array}$$

$$c_{vv}=b_{vv},c_{vs}=b_{vs}\text{ for s∈paB(v)},\text{ and }{c}_{vs}\in \left(0,b_{vs}\right)\text{ for s∈pa(v)∖paB(v)}.$$

$$c_{42}{c}_{21}=c_{43}{c}_{31},c_{42}{c}_{21}>c_{43}{c}_{31},\text{or}{c}_{42}{c}_{21}<c_{43}{c}_{31}.$$

$$\tilde{B}=\left(\begin{array}{cccc}1& 0& 0& 0\\ c_{21}& 1& 0& 0\\ 0& 0& 1& 0\\ c_{42}{c}_{21}& c_{42}& c_{43}& 1\end{array}\right)$$

$${\rho }_x(P,Q)=\frac{\mathrm{d}P}{\mathrm{d}(P+Q)}(x)$$

$$\rho_{x}(\hat{P},Q)\geq\rho_{x}(Q,\hat{P})\mbox{ for all $Q\in{\mathcal{P}}$,}$$

$${\rho }_x(c,c^{\ast })=\frac{\mathrm{d}{P}_c}{\mathrm{d}(P_c+P_{c^{\ast }})}(x_1,x_2)=\{\begin{array}{ll}1/2& \text{ for }{x}_2>c{x}_1\\ 1& \text{ for }{x}_2=c{x}_1\\ 0& \text{ for }{x}_2<c{x}_1\end{array}$$

Full text

1 Introduction

The type of model we are studying has been motivated by applications to risk analysis, where extreme risks play an essential role and may propagate through a network. For example, say, if an extreme rainfall happens on a specific location near a river network, it may effect water levels at other parts of the network in an essentially deterministic fashion. Similar phenomena occur in the analysis of risk for other complex systems.

Specifically, the model presented in (1) below arose in the context of technical risk analysis, more precisely, in an investigation of the “runway overrun” event of airplane landing. Numerous variables contribute to this event and extraordinary values of some variables lead invariably to a runway overrun (see Gissibl et al. (2017) for more details) naturally leading to questions about cause and effect of risky events. Other potential examples for risk-related cause and effect relations include chemical pollution of rivers (Ver Hoef et al. (2006)), flooding in river networks (Asadi et al. (2015)), financial risk (Einmahl et al. (2017)), and many others.

Statistical theory and applications of extreme value theory until the 1990s mainly focused on i.i.d. data as, for instance, yearly maximal water levels to predict future floodings or peaks over thresholds used to estimate the Value-at-Risk (e.g. Embrechts et al. (1997)). From this, both theory and applications moved on to multivariate data, modelling risks like joint wind and wave extremes as well as extreme risks in financial portfolios Beirlant et al. (2004). The investigation of extremes in time series models have proved useful in financial and environmental risk analysis, and also in telecommunication (see e.g. the book Finkenstädt and Rootzén (2004)). More recently, extreme space-time models have been suggested and applied to environmental risk data Buhl et al. (2016); Davis et al. (2013); Davison et al. (2012); Huser and Davison (2014).

The paper focuses on first steps reporting on the methodological development associated with a specific class of network models. We begin with introducing our leading example of a recursive max-linear model which is Example 2.1 of Gissibl and Klüppelberg (2015):

Example 1

Consider the network in the figure below:

1$$2$$3$$4

Each node $i$ in the network represents a random variable $X_{i}$ and the joint distribution of $X=(X_{1},X_{2},X_{3},X_{4})$ is determined by a system of max-linear structural equations

[TABLE]

where $Z_{1},Z_{2},Z_{3},Z_{4}$ are independent positive random variables and the coefficients $c_{ji}$ are all strictly positive.

The interpretation of a system like this is that each node in the network is subjected to a random shock $Z_{i}$ and the effect from shocks of other nodes pointing to it, the latter being attenuated or amplified by the coefficients $c_{ji}$. To simplify notation here and later we write $a\vee b$ for $\max(a,b)$. We can alternatively represent $X=(X_{1},X_{2},X_{3},X_{4})$ directly in terms of the noise variables as

[TABLE]

We may then summarize the above coefficients to the noise variables $Z_{1},\dots,Z_{4}$ in the matrix

[TABLE]

$\Box$

In greater generality we may write such a recursive max-linear model as

[TABLE]

where ${\rm pa}(v)$ denotes parents of $v$ in a directed acyclic graph (DAG) and $Z_{v}$ represent independent noise variables. The present article is concerned with such models and summarizes basic elements of Gissibl and Klüppelberg (2015) and Gissibl et al. .

In this setting, natural candidates for the noise distributions are extreme value distributions or distributions in their domains of attraction resulting in a corresponding multivariate distribution with dependence structure given by the DAG (for details and background on multivariate extreme value models see e.g. de Haan and Ferreira (2006); Resnick (1987, 2007)). The paper is structured as follows.

In Section 2 we establish the necessary terminology (Section 2.1), introduce Bayesian networks (Section 2.2), and basic properties of conditional independence (Section 2.3). In Section 2.4 we establish basic Markov properties of Bayesian networks. In Section 3 we study the specific Markov properties of Bayesian networks given by max-linear structural equations as in (1) and in Section 4 we study statistical properties of the models.

2 Preliminaries

2.1 Graph terminology

A graph as we use it here is determined by a finite vertex set $V$, an edge set $E$, and a map that to each edge $e$ in $E$ associates its endpoints $u,v\in V$. Our graphs are simple so that there are no self-loops (edges with identical endpoints) and no multiple edges. Therefore we can identify an edge $e$ with its endpoints $u,v$ so we can write $e=uv$. An edge $uv$ of a directed graph points from $u$ to $v$ and we write $u\to v$. Then $u$ is a parent of $v$ and $v$ is a child of $u$. The set of parents of $v$ is denoted ${\rm pa}(v)$ and the set of children of $u$ is ${\rm ch}(u)$. If $uv$ is an edge we also say that $u$ and $v$ are adjacent and write $u\sim v$ whether or not the edge is directed.

A walk $\omega$ from $u$ to $v$ of length $n$ is a sequence of vertices $\omega=[u=u_{0},u_{1},\ldots,u_{n}=v]$ so that $u_{i-1}\sim u_{i}$ for all $i=1,\ldots,n$. A walk is a cycle if $u=v$. A path is a walk with no repeated vertices. The walk is directed from $u$ to $v$ if $u_{i-1}\to u_{i}$ for all $i$. If all edges in a graph $\mathcal{D}=(V,E)$ are directed, $\mathcal{D}$ is a directed graph. A directed graph is acyclic if it has no directed cycles. A DAG is a directed acyclic graph. A DAG is a tree if every vertex has at most one parent and a polytree if there is at most one path between two vertices $u$ and $v$.

If there is a directed path from $u$ to $v$ in $\mathcal{D}$ we say that $u$ is an ancestor of $v$ and $v$ a descendant of $u$ and write $u\leadsto v$ or v\mathrel{\mathchoice{\reflectbox{\displaystyle\rightsquigarrow}}{\reflectbox{\textstyle\rightsquigarrow}}{\reflectbox{\scriptstyle\rightsquigarrow}}{\reflectbox{\scriptscriptstyle\rightsquigarrow}}}u. The set of ancestors of $v$ is denoted ${\rm an}(v)$. A set $A\subseteq V$ is said to be ancestral if ${\rm an}(v)\subset A$ for all $v\in A$, or, alternatively, if ${\rm pa}(v)\subset A$ for all $v\in A$. For a subset $A$ of $V$ we let ${\rm An}(A)$ denote the smallest ancestral set containing $A$.

We say that the vertex set $V$ of a DAG $\mathcal{D}$ is well-ordered if $V=\{1,\ldots,d\}$ and all edges in $\mathcal{D}$ point from low to high, i.e. if $ij\in E\implies i<j$. Then the set of predecessors of a vertex $i$ is ${\rm pr}(i)=\{1,\ldots,i-1\}$.

For a DAG $\mathcal{D}$ we define its moral graph $\mathcal{D}^{m}$ as the simple, undirected graph with the same vertex set but with $u$ and $v$ adjacent in $\mathcal{D}^{m}$ if and only if either $u\sim v$ in $\mathcal{D}$ or if $u$ and $v$ have a common child. For further general graph terminology we refer the reader to West (2001) but some of the concepts above are illustrated in Fig. 1.

2.2 Bayesian networks

A real-valued Bayesian network associated to a given DAG $\mathcal{D}=(V,E)$ is determined by specifying random variables $X=(X_{v},v\in V)$ and the conditional distribution of each of these, given values of their parent variables; for example as

[TABLE]

Because there are no directed cycles in $\mathcal{D}$ there is a unique joint distribution corresponding to this specification.

Alternatively, as in Example 1, we can specify these conditional distributions through structural equations which describe the conditional distribution of $X_{v}$ conditionally on $X_{{\rm pa}(v)}=x_{{\rm pa}(v)}$ in a functional form. More precisely a system of equations of the form

[TABLE]

where $(Z_{v})_{v\in V}$ are independent noise variables and $g_{v}$ suitable functions.

A system of structural equations as above is sometimes referred to as a data generating mechanism, interpreting each equation as a way of generating random variables with the desired conditional distribution.

An important instance of these models are linear structural equation models where the functions $g_{v}$ are linear and hence

[TABLE]

where $c_{vu},u\in{\rm pa}(v),c_{vv}$ are structural coefficients, see for example Bollen (1989). In general, a structural equation system need not be associated with a DAG, but if it is, the equation system is said to be recursive.

If the distributions of $Z_{v}$ have heavy tails and all structural coefficients are non-negative, the sum tends to be dominated by the largest term:

[TABLE]

and hence for such cases, the max-linear variant in (4) as described in more detail in Section 3 below.

2.3 Conditional independence

The notion of conditional independence is at the heart of graphical models, including Bayesian networks. For three random variables $(X,Y,Z)$ we say that $X$ is conditionally independent of $Y$ given $Z$ if the conditional distribution of $X$ given $(Y,Z)$ does not depend on $Y$ and we then write $X{\perp\!\!\!\!\perp}Y\,|\,Z$ or $X{\perp\!\!\!\!\perp}_{P}Y\,|\,Z$ if we wish to emphasize the dependence on the joint distribution $P$ of $(X,Y,Z)$.

The notion of conditional independence has a number of important properties, see e.g. Dawid (1980) or Lauritzen (1996).

Proposition 1

Let $(\Omega,\mathbb{F},P)$ be a probability space and $X$, $Y$, $Z$, $W$ random variables on $\Omega$. Then the following properties hold.

(C1)

If $X{\perp\!\!\!\!\perp}Y\,|\,Z$ then $Y{\perp\!\!\!\!\perp}X\,|\,Z$ (symmetry);

(C2)

If $X{\perp\!\!\!\!\perp}Y\,|\,Z$ and $W=\phi(Y)$ then $X{\perp\!\!\!\!\perp}W\,|\,Z$ (reduction);

(C3)

If $X{\perp\!\!\!\!\perp}(Y,Z)\,|\,W$ then $X{\perp\!\!\!\!\perp}Y\,|\,(Z,W)$ (weak union);

(C4)

If $X{\perp\!\!\!\!\perp}Z\,|\,Y$ and $X{\perp\!\!\!\!\perp}W\,|\,(Y,Z)$ then $X{\perp\!\!\!\!\perp}(Z,W)\,|\,Y$ (contraction);

It is occasionally important to abstract the notion of conditional independence away from necessarily being concerned with probability measures. An (abstract) independence model ${\perp_{\sigma}}$ over $V$ is a ternary relation over subsets of a finite set $V$. The independence model is a semi-graphoid if the following holds for mutually disjoint subsets $A$, $B$, $C$, $D$:

(S1)

If $A{\perp_{\sigma}}B\,|\,C$ then $B{\perp_{\sigma}}A\,|\,C$ (symmetry);

(S2)

If $A{\perp_{\sigma}}(B\cup D)\,|\,C$ then $A{\perp_{\sigma}}B\,|\,C$ and $A{\perp_{\sigma}}D\,|\,C$ (decomposition);

(S3)

If $A{\perp_{\sigma}}(B\cup D)\,|\,C$ then $A{\perp_{\sigma}}B\,|\,(C\cup D)$ (weak union);

(S4)

If $A{\perp_{\sigma}}B\,|\,C$ and $A{\perp_{\sigma}}D\,|\,(B\cup C)$, then $A{\perp_{\sigma}}(B\cup D)\,|\,C$ (contraction);

Further, the independence model is a graphoid if it also satisfies

(S5)

If $A{\perp_{\sigma}}B\,|\,(C\cup D)$ and $A{\perp_{\sigma}}C\,|\,(B\cup D)$ then $A{\perp_{\sigma}}(B\cup C)\,|\,D$ (intersection).

We shall in particular be interested in distributions on product spaces ${\mathcal{X}}=\times_{v\in V}{\mathcal{X}}_{v}$ where $V$ is a finite set. For $A\subseteq V$ we write $x_{A}=(x_{v},v\in A)$ to denote a generic element in${\mathcal{X}}_{A}=\times_{v\in A}{\mathcal{X}}_{v}$, and similarly $X_{A}=(X_{v})_{v\in A}$.

If $P$ is a probability distribution on ${\mathcal{X}}$, we can now define an independence model ${\perp\!\!\!\!\perp}$ by the relation

[TABLE]

and it follows from Proposition 1 that ${\perp\!\!\!\!\perp}$ is a semi-graphoid; in general ${\perp\!\!\!\!\perp}$ is not a graphoid without further assumptions on $P$.

Another important independence model is determined by separation in an undirected graph. More precisely, if ${\mathcal{G}}=(V,E)$ is an undirected graph we can define an independence model $\perp_{\mathcal{G}}$ by letting $A\perp_{\mathcal{G}}B\,|\,S$ mean that all paths in ${\mathcal{G}}$ from $A$ to $B$ intersect $S$. Then it is easy to see that $\perp_{\mathcal{G}}$ is always a graphoid; indeed the term graphoid refers to this fact.

For a directed graph, the relevant notion of separation is more subtle. A vertex $u$ is a collider on a path $\pi$ if two arrowheads meet on the walk at $u$, i.e. if the following situation occurs $\pi=[\cdots\rightarrow u\leftarrow\cdots]$.

We say that a path $\pi$ from $u$ to $v$ in a DAG $\mathcal{D}$ is connecting relative to $S$, if all colliders on $\pi$ are in the ancestral set ${\rm An}(S)$, and all non-colliders are outside $S$. A path that is not connecting relative to $S$ is said to be blocked by $S$. We then define an independence model $\perp_{\mathcal{D}}$ relative to a directed graph $\mathcal{D}$ as follows:

Definition 1

For three disjoint subsets $A$, $B$, and $S$ of the vertex set $V$ of a graph ${\mathcal{G}}=(V,E)$ we say that $A$ and $B$ are $\mathcal{D}$-separated by $S$ if all paths from $A$ to $B$ are blocked by $S$ and we then write $A\perp_{\mathcal{D}}B\,|\,S$.   $\Box$

Example 2

Consider the network in the figure below, only slightly more complex than in Example 1:

1$$2$$3$$4$$5

We have $2\perp_{\mathcal{D}}3\,|\,1$ since the path $2\leftarrow 1\to 3$ is blocked as the non-collider $1$ is in $S=\{1\}$ whereas the path $2\to 4\to 3$ is blocked because the collider $4$ is not an ancestor of $S=\{1\}$; on the other hand it holds that $\neg(2\perp_{\mathcal{D}}3\,|\,\{1,5\})$ since now the second path is rendered active as the collider $4$ is in ${\rm An}(\{1,5\})$.   $\Box$

Note that this definition in a natural way extends that of $\perp_{\mathcal{G}}$ for an undirected graph, as an undirected graph does not have colliders. The independence model $\perp_{\mathcal{D}}$ also satisfies the graphoid axioms, see e.g. Lauritzen and Sadeghi (2017).

There is an alternative method for checking $\mathcal{D}$-separation in terms of standard separation in a suitable undirected graph, associated with the query. More precisely we say that $A$ is $m$-separated from $B$ by $S$ and we write $A\perp_{m}B\,|\,S$ if $S$ separates $A$ from $B$ in the moral graph $(\mathcal{D}_{{\rm An}(A\cup B\cup S)})^{m}$. We then have:

Proposition 2

Let $A$, $B$ and $S$ be disjoint subsets of the nodes of a directed acyclic graph ${\mathcal{G}}$. Then $A\perp_{\mathcal{D}}B\,|\,S\iff A\perp_{m}B\,|\,S$.

For a proof, see Richardson (2003), amending an inaccuracy in Lauritzen et al. (1990).

Example 3

To illustrate the alternative procedure, we again consider the network in Example 2.

If we wish to check whether $2\perp_{\mathcal{D}}3\,|\,1$ we consider the subgraph induced by the ancestral set of $\{1,2,3\}$ and moralize to obtain the graph to the left in the figure below. Since $1$ is a separator in this graph, we conclude that $2\perp_{\mathcal{D}}3\,|\,1$.

1$$2$$3$$4$$5$$1$$2$$3

On the other hand, if the query is whether $2\perp_{\mathcal{D}}3\,|\,\{1,5\}$ we have ${\rm An}(\{1,5\})=V$ and thus the relevant moral graph is given to the right in the figure above; in this graph, $2$ and $3$ are not separated by $\{1,5\}$ so we conclude $\neg(2\perp_{\mathcal{D}}3\,|\,\{1,5\})$.   $\Box$

2.4 Markov properties of Bayesian networks

It follows directly from the construction of a Bayesian network, that the joint distribution $P$ satisfies the well-ordered Markov property (O) w.r.t. $\mathcal{D}$ if for some well-ordering of $V$, every variable is conditionally independent of its predecessors given its parents

[TABLE]

for all $v\in V=\{1,\ldots,d\}$.

We further say that $P$ obeys the local Markov property (L) w.r.t. $\mathcal{D}$ if every variable is conditionally independent of its non-descendants, given its parents:

[TABLE]

And, finally, $P$ satisfies the global Markov property (G) w.r.t. $\mathcal{D}$ if

[TABLE]

Example 4

Consider the network in the figure below:

1$$2$$3$$4$$6$$5

The numbering of the nodes here constitute a well-ordering so, for example, (O) implies $5{\perp\!\!\!\!\perp}\{1,3,4\}\,|\,2$, whereas the local Markov property (L) implies $5{\perp\!\!\!\!\perp}\{1,3,4,6\}\,|\,2$; the global Markov property implies, for example, $5{\perp\!\!\!\!\perp}\{1,6\}\,|\,4$.   $\Box$

In the case of undirected graphs, the local and global Markov properties are different (Lauritzen, 1996, Section 3.2), but here we have

Theorem 2.1

Let $\mathcal{D}$ be a directed acyclic graph with $V=\{1,\ldots,d\}$ well-ordered and $P$ a probability distribution on ${\mathcal{X}}=\times_{v\in V}{\mathcal{X}}_{v}$. Then we have

[TABLE]

In words, if $P$ satisfies any of these Markov properties, it satisfies all of them.

Proof

This fact is established in (Lauritzen et al., 1990, Corollary 2) for any semi-graphoid independence model ${\perp_{\sigma}}$.   $\Box$

Note that in particular it is true that if $P$ satisfies (O) w.r.t. one well-ordering, it satisfies (O) w.r.t. all well-orderings.

The global Markov property gives a sufficient condition for conditional independence in terms of $\mathcal{D}$-separation. Another central concept is that of faithfulness, formally defined below

Definition 2

A probability distribution $P$ on ${\mathcal{X}}=\times_{v\in V}{\mathcal{X}}_{v}$ is said to be faithful to a DAG $\mathcal{D}$ if

[TABLE]

In other words, if $\mathcal{D}$-separation is also necessary for conditional independence.   $\Box$

Generally, most probability distributions are faithful (Meek, 1995), but we shall later see that this is not the case for the special Bayesian networks we study here.

Finally, we need to emphasize that two different DAGs can define exactly the same independence model. Consider two graphs $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ as well as their associated independence models $\perp_{\mathcal{D}_{1}}$ and $\perp_{\mathcal{D}_{2}}$. It may well happen that even though the graphs are different, their independence models might be identical, see for example Figure 2 below.

Here all independence models are the same although the graphs are different. This also means that any probability distribution $P$ which satisfies the global Markov property for any of them, automatically satisfies the global Markov property for all of them. We formally define

Definition 3

Two DAGs $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ are Markov equivalent if and only if their independence models coincide, i.e. if $A\perp_{\mathcal{D}_{1}}B\,|\,C\iff A\perp_{\mathcal{D}_{2}}B\,|\,C$.   $\Box$

The following result was shown by Frydenberg (1990) and Verma and Pearl (1990) and gives a necessary and sufficient condition for two DAGs to be Markov equivalent.

Theorem 2.2

Two directed acyclic graphs $\mathcal{D}_{1}=(V,E_{1})$ and $\mathcal{D}_{2}=(V,E_{2})$ are Markov equivalent if and only if they have the same skeleton ${\rm ske}(\mathcal{D}_{1})={\rm ske}(\mathcal{D}_{2})$ and the same unshielded colliders.

Here the skeleton ${\rm ske}(\mathcal{D})$ of a DAG $\mathcal{D}$ is the undirected graph with $u\sim v$ in ${\rm ske}(\mathcal{D})$ if $u\sim v$ in $\mathcal{D}$, and an unshielded collider is a triple $u\to w\leftarrow v$ with $u\not\sim v$.

3 Recursive max-linear structural equation models

We shall be interested in Bayesian networks defined through structural equation systems (2) where the functions $g_{v}$ are max-linear, i.e. the additions in (3) are replaced with the operation of forming the maximum.

Henceforth we assume that the vertex set of our DAG $\mathcal{D}=(V,E)$ is well-ordered so $V=\{1,\ldots,d\}$ and assume a data generating mechanism specified via a recursive max-linear structural equation model, which has representation

[TABLE]

where $Z_{1},\dots,Z_{d}$ are independent and identically distributed with a continuous distribution having support ${\mathbb{R}}_{+}=(0,\infty)$, and $c_{vu}>0$, $u\in{\rm pa}(v)$, $c_{vv}$ are structural coefficients in the equations or edge weights for the associated DAG $\mathcal{D}$.

Following Gissibl and Klüppelberg (2015) we say this is a recursive max-linear model. Note that our use of indices for edge weights here is the opposite of that used in Gissibl and Klüppelberg (2015).

For simplicity we assume throughout the rest of the paper that $c_{vv}=1$ for all $v\in V$. For a path $\pi=[u=k_{0}\rightarrow k_{1}\rightarrow\dots\rightarrow k_{n}=v]$ of length $n$ from $u$ to $v$, we define the quantities

[TABLE]

where $\Pi_{uv}$ denotes all paths from $u$ to $v$. In summary, we define

[TABLE]

where ${\rm An}(v)={\rm an}(v)\cup\{v\}$ is the smallest ancestral set containing vertex $v$. We then arrange these coefficients in the max-linear coefficient matrix $B=(b_{vu})_{d\times d}$ and find

[TABLE]

This equation represents $X$ as a max-linear model as defined for instance in Wang and Stoev (2011).

For two non-negative matrices $F$ and $G$, where the number $n$ of columns in $F$ is equal to the number of rows in $G$ we introduce the product $\odot$ as

[TABLE]

If we collect the noise variables into the column vector $Z=(Z_{1},\ldots,Z_{d})^{\prime}$, the representation (7) of $X$ can then be written as

[TABLE]

Given the DAG  $\mathcal{D}$ and the edge weights $c_{ik}$ with $c_{ii}=1$ for all $i=1,\dots,d$, the max-linear coefficient matrix $B$ can be found by iterating the weighted adjacency matrix $C=(c_{vu}\mathbf{1}_{{\rm Pa}(v)}(u))_{d\times d}$ of $\mathcal{D}$ using this matrix multiplication; here $\mathbf{1}_{{\rm Pa}(v)}$ denotes the indicator function of ${\rm Pa}(v)={\rm pa}(v)\cup\{v\}$) :

[TABLE]

cf. Butkovič (2010), Lemma 1.4.1. For more details see Gissibl and Klüppelberg (2015), Theorem 2.4.

By (6) the max-linear coefficient $b_{vu}$ of $X$ is different from zero if and only if $u\in{\rm An}(v)$. This information is contained in the reachability matrix $R=(r_{vu})_{d\times d}$ of $\mathcal{D}$, which has entries

[TABLE]

If the $vu$-th entry of $R$ is equal to one, then $v$ is reachable from $u$. In the context of a DAG $\mathcal{D}$ with its reachability matrix $R$ and a recursive max-linear model $X$ on $\mathcal{D}$ with max-linear coefficient matrix $B$ it will be useful to keep the following in mind.

Remark 1

Let $\mathcal{D}$ be a DAG with reachability matrix $R$.

(i)   The max-linear coefficient matrix $B$ is a weighted reachability matrix of $\mathcal{D}$; i.e., $R={\rm sgn}(B)$.

(ii)   Since $V$ is assumed well-ordered, $B$ and $R$ are lower triangular matrices.   $\Box$

From (6) and (7) we conclude that a path $\pi$ from $u$ to $v$, whose weight $d_{vu}(\pi)$ is strictly less than $b_{vu}$ does not have any influence on $X_{i}$. For $v\in V$ and $u\in{\rm an}(v)$ we call a path $\pi$ from $u$ to $v$ max-weighted, if $b_{vu}=d_{vu}(\pi)$, and investigate its relevance for the recursive max-linear model in further detail.

Firstly we note that we can remove an edge from $\mathcal{D}$ which is not part of a max-weighted path without changing the distribution of $X$. The DAG obtained in this way is termed the minimum max-linear DAG $\mathcal{D}^{B}$. In the special case where $\mathcal{D}$ is a polytree, all paths are necessarily max-weighted and we clearly have

Proposition 3

If $\mathcal{D}$ is a polytree, it holds that $\mathcal{D}^{B}=\mathcal{D}$.

The following result describes exactly all DAGs and edge weights possible for a given max-linear coefficient matrix. Recall that we have set $c_{vv}=1$.

Theorem 3.1

*[Gissibl and Klüppelberg (2015), Theorem 5.4]

Let $X$ be given by a recursive max-linear structural equation system with coefficient matrix $B$. Let further $\mathcal{D}^{B}=(V,E^{B})$ be the minimum max-linear DAG of $X$ and ${\rm pa}^{B}(v)$ the parents of $v$ in $\mathcal{D}^{B}$.*

1. (a)

$\mathcal{D}^{B}$* is the DAG  with the minimum number of edges such that $X$ satisfies (4). The weights in (4) are uniquely given by $c_{vv}=b_{vv}$ and $c_{vs}=b_{vs}$ for $v\in V$ and $s\in{\rm pa}^{B}(v)$.* 2. (b)

Every DAG  with vertex set $V$ that has at least the edges of $\mathcal{D}^{B}$ and the same reachability matrix as $\mathcal{D}^{B}$ represents $X$ in the sense of (4) with weights satisfying

[TABLE]

There are no further DAGs and weights such that $X$ has representation (4).

In general, recursive max-linear models are not faithful to their DAG, not even if $\mathcal{D}=\mathcal{D}^{B}$, see Remark 3.9 (ii) in Gissibl and Klüppelberg (2015). This is illustrated in Example 5 below.

Example 5

[Example 3.8 of Gissibl and Klüppelberg (2015) and continuation of Example 1:] We note that the paths $[1\to 2]$, $[1\to 3]$, $[2\to 4]$, and $[3\to 4]$ are max-weighted as they are the only directed paths between their endpoints. It therefore holds that $\mathcal{D}^{B}=\mathcal{D}$ since they are the unique max-weighted paths. Still, the distribution determined by this recursive system is never faithful to $\mathcal{D}$, as we shall see below.

Concerning the paths from node $1$ to $4$ we have three situations:

[TABLE]

In the first situation, both paths from $1$ to $4$, $[1\to 2\to 4]$ and $[1\to 3\to 4]$, are max-weighted whereas in the other situations only one of them is.

If the path $[1\to 2\to 4]$ is max-weighted, we can consider the subdag $\tilde{\mathcal{D}}$ obtained from $\mathcal{D}$ by removing the edge $1\to 3$:

1$$2$$3$$4

In other words, we are changing the edge weights by letting $\tilde{c}_{31}=0$, keeping the other edge weights unchanged. The new max-linear coefficient matrix becomes

[TABLE]

where we have exploited that $c_{ii}=1$. The max-linear coefficient matrix for the marginal distribution of $(X_{1},X_{2},X_{4})$ is obtained by ignoring the third row and since only entries in the third row have changed, we see that $(X_{1},X_{2},X_{4})$ has the same joint distribution in the model determined by $\mathcal{D}$ as it has in the model determined by $\tilde{\mathcal{D}}$.

But as we clearly have $1\perp_{\tilde{\mathcal{D}}}4\,|\,2$, we conclude that $X_{1}{\perp\!\!\!\!\perp}X_{4}\,|\,X_{2}$ in the model determined by $\tilde{\mathcal{D}}$ and hence also by $\mathcal{D}$. But since $\neg(1\perp_{\mathcal{D}}4\,|\,2)$, the distribution is not faithful to $\mathcal{D}$.

If $[1\to 3\to 4]$ is also max-weighted, the similar argument yields $X_{1}{\perp\!\!\!\!\perp}X_{4}\,|\,X_{3}$, so the distribution is not faithful to $\mathcal{D}$ for any allocation of edge weights.   $\Box$

We note that Gissibl and Klüppelberg (2015) suggest in their Remark 3.9(i) that additional conditional independence relations that are valid for a given DAG can be revealed by considering a system of submodels determined by appropriate subgraphs, but here we refrain from giving a complete description of all valid conditional independence relations.

4 Statistical properties

The statistical theory of recursive max-linear models is challenging because standard assumptions for smooth statistical models are not satisfied. For example, if we for a given DAG $\mathcal{D}$ consider the family $\cal P$ of distributions with coefficients adapted to $\mathcal{D}$, this family is not dominated by any measure on the space of observations, so standard likelihood theory does not apply. On the other hand, as we shall see, estimation of coefficients and identification of the network structure for recursive max-linear models can be made in a simple fashion and procedures are more efficient than usual in that estimates of coefficients and structures converge at exponential rates to the true values. Here we shall give a summary of the most important findings in Gissibl et al. .

Throughout the following we consider a sample $\mathbf{x}=(X^{1}=x^{1},\ldots,X^{n}=x^{n})$ from a distribution $P$ given by the recursive max-linear model (4).

4.1 Estimation of coefficients

We first consider the situation where the DAG $\mathcal{D}=(V,E)$ and for the sake of simplicity we assume the distribution of noise variables $Z_{v},v\in V,$ is completely known, the coefficients $c_{vv}$ are all equal to one, whereas the edge weights $C=\{c_{vu},u\in{\rm pa}(v),v\in V\}$ are all strictly positive, but otherwise unknown. We let ${\mathcal{C}}$ denote the set of all possible coefficients and $P_{C}$ denote the distribution of $X$ determined by the corresponding recursive model (4).

The family ${\cal P}=P_{C},C\in{\mathcal{C}},$ is not dominated by any fixed $\sigma$-finite measure $\mu$ on ${\mathcal{X}}$, as the support of $P_{C}$ varies strongly with the coefficients; more precisely, the distributions have disjoint atomic components. This is a disadvantage in the sense that we cannot define a standard likelihood function; but, as we shall see, an advantage since these atomic components help identifying $P_{C}$ from a given sample. We illustrate this by a simple example.

Example 6

[Estimation from the atoms]

Consider the simple DAG $1\to 2$ with just two nodes and a single directed edge, and let $c=c_{21}$ be the corresponding coefficient.

Then $P_{c}$ has support on the cone given as $x_{2}\geq cx_{1}\geq 0$ and the line $A_{c}=\{x_{2}=cx_{1}\}$ is an atom for $P_{c}$ because $P_{c}(A_{c})=P(Z_{2}\leq cX_{1})=P(Z_{2}\leq cZ_{1})>0$ (cf. Remark LABEL:rem(ii)).

Still, since then $\{c\}$ is the only atom in $P_{c}$ for $Y=X_{2}/X_{1}$, the sample will for large $n$ with high probability have repeated values of $Y$ and $c$ will be the only value that is repeated. In other words, $\hat{c}=\min\{y^{\nu}=x_{2}^{\nu}/x^{\nu}_{1},\nu=1,\ldots,n\}$ will be exactly equal to the true parameter with high probability. A similar estimator has been considered by Davis and Resnick (1989) in a time-series framework.   $\Box$

Although most likelihood theory is concerned with dominated families, Kiefer and Wolfowitz (1956) considered the non-dominated case. Their formulation has been used rarely — an exception being Johansen (1978); see also Scholz (1980) and Gill et al. (1989), for example. This formulation turns out to be exactly what we need to discuss estimation of $C$ in a formal way.

For two probability measures $P$ and $Q$ on a measurable space $({\mathcal{X}},\mathbb{E})$, we define the generalized likelihood ratio $\rho_{x}(P,Q)$ at the observation $x$ as

[TABLE]

where ${\mathrm{d}P}/{\mathrm{d}(P+Q)}$ is the density of $P$ w.r.t. $P+Q$; the density always exists as, clearly, $P(A)+Q(A)=0\implies P(A)=0$ so $P$ is absolutely continuous w.r.t. $P+Q$.

The idea here is that if $\rho_{x}(P,Q)>\rho_{x}(Q,P)$, then $P$ is a more likely explanation of $x$ than $Q$. We note in particular that if $P$ and $Q$ have densities $f$ and $g$ w.r.t. a $\sigma$-finite measure $\mu$, we have $\rho_{x}(P,Q)=f(x)/\{f(x)+g(x)\}$ so then $\rho_{x}(P,Q)>\rho_{x}(Q,P)$ if and only if $f(x)>g(x)$. Hence $\rho_{x}$ extends the standard likelihood ratio in a natural way.

Clearly, the generalized likelihood ratio suffers from the same problem as the usual likelihood ratio: the densities are only defined almost surely, so can be changed on $P+Q$-null sets; therefore, a version of ${\mathrm{d}P}/{\mathrm{d}(P+Q)}$ must be chosen independently of the observation $x$.

Next we say that if ${\mathcal{P}}$ is a family of probability distributions, $\hat{P}$ is a generalized maximum likelihood estimate (GMLE) of $P$ based on $x\in{\rm supp}(\hat{P})$ if

[TABLE]

i.e. if $\hat{P}$ explains $x$ at least as well as any other member of ${\mathcal{P}}$.

Example 7

[Continuation of Example 6: GMLE]

We illustrate use of the generalized maximum likelihood ratio for the model described in Example 6. To identify the density, we consider two values $c>c^{*}$ where we have

[TABLE]

and

[TABLE]

If $c=c^{*}$ we may let

[TABLE]

Thus, if we consider a full sample, let $\hat{c}=\min\{y^{\nu}=x_{2}^{\nu}/x^{\nu}_{1},\nu=1,\ldots,n\}$ and $n_{+}(c,\mathbf{x})=\#\{\nu:y^{\nu}>c\}$, we get:

[TABLE]

whereas

[TABLE]

Clearly, $\rho_{\mathbf{x}}(\hat{c},c)\geq\rho_{\mathbf{x}}(c,\hat{c})$ showing that $\hat{c}$ is the unique GMLE of $c$.   $\Box$

Indeed, it holds in general for a recursive max-linear model that

[TABLE]

is a GMLE of the edge weights. We refer to Gissibl et al. for further details but should point out that in the general case, the GMLE is not unique. Since the distribution of $X$ only depends on the edge weights through the max-linear coefficient matrix $B$, only $B$ is uniquely estimable from a sample. We clearly have by (9) for the GMLE that

[TABLE]

An alternative estimate of the max-linear coefficient matrix is given as

[TABLE]

Although this estimate is also sensible and asymptotically consistent, it is less efficient than the GMLE as $X^{\nu}_{i}/X^{\nu}_{j}$ only attends its minimum value when all noise variables on the path from $j$ to $i$ are smaller than $b_{ij}X^{\nu}_{j}$ for the same $\nu$, whereas the minima for the $X^{\nu}_{v}/X^{\nu}_{u}$ on the path from $j$ to $i$ can be attained for different $\nu$s.

4.2 Identification of structure

General methods for identifying the structure of DAG $\mathcal{D}$ from a sample are often based on an assumption of faithfulness, so that observed conditional independence relations can be translated back to the structure of the DAG since then any observed conditional independence must correspond to a separation in $\mathcal{D}$, see for example Spirtes et al. (2000). Also, as noted in Theorem 2.2, two DAGs can be different but still Markov equivalent and thus any method based on observed direct conditional independence relations cannot distinguish between DAGs that are Markov equivalent.

As shown in Example 5, faithfulness is violated for max-linear Bayesian networks whenever $\mathcal{D}$ is not a polytree. However, as we shall see below, the minimal DAG $\mathcal{D}^{B}$ of a max-linear Bayesian network can still be completely recovered from observations.

This fact conforms with recent developments where the recursive linear structural equation systems have been shown to be completely identifiable if the errors follow a non-Gaussian distribution (Shimizu et al. (2006)) and it has been shown that the faithfulness assumption can be considerably weakened also in other situations (Spirtes and Zhang (2014); Peters and Bühlmann (2014)).

To explain why the structure $\mathcal{D}^{B}$ is identifiable, we consider the statistics

[TABLE]

and note that $Y_{ij}$ has support $[b_{ij},\infty)$ and an atom in $b_{ij}$ if and only if $j\in{\rm an}(i)$. Using this property one can show that the following estimate $\check{B}$ eventually identifies the max-linear coefficient matrix $B$.

[TABLE]

Then $\mathcal{D}^{B}$ is identifiable from $B$; we refer the reader to Gissibl et al. for further details.

5 Conclusion

We have reviewed basic elements of Bayesian networks based on recursive max-linear structural equations and some of their statistical properties. We conclude this article by pointing out some natural extensions of this work that we hope to address in the future.

Firstly, it would be of interest to have a simple and complete description of all independence properties which hold for a distribution determined by a recursive max-linear equation system, i.e. a global Markov property for max-linear Bayesian networks.

Secondly, it appears that a consequent use of algebraic theory; see e.g. Butkovič (2010), based on properties of the max-times semiring $\mathbb{S}=([0,\infty],\vee,\cdot)$ would be able to simplify the theory of these models.

Finally, we should emphasize that the models heuristically can be seen as limiting cases of standard linear recursive models where error distributions have heavy tails and therefore the maximal element of any sum will almost completely dominate the sum; a rigorous study of this limiting process will enhance the understanding of this class of models.

Acknowledgements

The authors have benefited from discussions with Nadine Gissibl and financial support from the Alexander von Humboldt Stiftung.