Non-Parametric Bayesian Networks: Improving

Non-Parametric Bayesian Networks: Improving Theory and Reviewing Applications Anca Hanea

a ∗

[email protected] a b

Oswaldo Morales Napolesb

Dan Ababeic

[email protected]

[email protected]

Centre of Excellence for Biosecurity Risk Analysis, University of Melbourne

Hydraulic Engineering, Faculty of Civil Engineering & Geosciences, Delft University of Technology c

Lighttwist Software

Abstract Applications in various domains often lead to high dimensional dependence modelling. A Bayesian network (BN) is a probabilistic graphical model that provides an elegant way of expressing the joint distribution of a large number of interrelated variables. BNs have been successfully used to represent uncertain knowledge in a variety of fields. The majority of applications use discrete BNs, i.e. BNs whose nodes represent discrete variables. Integrating continuous variables in BNs is an area fraught with difficulty. Several methods that handle discrete-continuous BNs have been proposed in the literature. This paper concentrates only on one method called non-parametric BNs (NPBNs). NPBNs were introduced in 2004 and they have been or are currently being used in at least twelve professional applications. This paper provides a short introduction to NPBNs, a couple of theoretical advances, and an overview of applications. The aim of the paper is twofold: one is to present the latest improvements of the theory underlying NPBNs, and the other is to complement the existing overviews of BNs applications with the NPNBs applications. The latter opens the opportunity to discuss some difficulties that applications pose to the theoretical framework and in this way offers some NPBN modelling guidance to practitioners. Keywords: Bayesian networks, model selection, quantification, inference, structured expert judgement, risk, safety, reliability, parameter estimation

1

Introduction

Understanding and representing multivariate distributions along with their dependence structure is a highly active area of research. A large body of scientific work treating multivariate models is available. This paper advocates graphical models in general and Bayesian networks (BNs) in particular to represent high dimensional distributions with complex dependence structures. A BN consists of a directed acyclic graph (DAG) and a set of (conditional) distributions. Each node in the graph corresponds to a random variable and the arcs represent direct qualitative dependence relationships. The absence of arcs guarantees a set of (conditional) independence facts. The direct predecessors (successors) of a node are called parents (children). A marginal distribution is specified for each node with no parents, and a conditional distribution is associated with each child node. The (conditional) distributions serve as the quantitative information about the strength of the dependencies between the variables involved. The graph with the conditional independence statements encoded by it, together with the (conditional) distributions represent the joint distribution over the random variables denoted by the nodes of the graph. The quantitative information needed can be either retrieved from data, when available, or from experts. The relatively simple visualization of the complicated relationships among the random variables is one of the most appealing features of a BN model. The main use of BNs is to update distributions given observations. This is referred to as inference in BNs. ∗ Correspondence

to: A.M.Hanea, CEBRA, University of Melbourne, Parkville, VIC 3010, Australia

1

BNs have been successfully used to represent uncertain knowledge, in a consistent probabilistic manner, in a variety of fields (Weber et al. 2010). Nevertheless, most of the applications use discrete BNs, i.e BNs whose nodes represent discrete random variables. These models specify marginal distributions for the nodes with no parents, and conditional probability tables for child nodes. Inference in discrete BNs can be performed using either exact algorithms (e.g. Zhang and Poole 1994, Pearl 1988), or approximate algorithms (e.g. Jordan et al. 1999; Murphy 2002b). Despite their popularity, discrete BNs suffer from severe limitations. Applications involving high complexity in data-sparse environments are severely limited by the excessive assessment burden which leads to rapid, informal and indefensible quantification. Furthermore, this type of representation, using only discrete variables, is inadequate for many problems of practical importance. Discrete BNs have been extensively studied, hence we omit their presentation here, but refer to (Weber et al. 2010) and references therein. Many domains require reasoning about the joint behavior of a mixture of discrete and continuous variables. These domains are often called hybrid domains. Hence, BNs involving both discrete and continuous variables will be called hybrid BNs (HBNs). Working with HBNs proves considerably more challenging than working with their discrete counterparts. Several methods for HBNs are available and discussed in the literature (e.g. Lerner 2002, Langseth et al. 2009). We will only briefly mention some of them, insist on the most recent ones, and present in more detail the method called non-parametric BNs (NPBNs) that is overlooked in many HBNs review papers. The emphasis here is on the application of the NPBN methods in different fields of science and engineering in the course of approximately 10 years since NPBNs were first presented. A study which is perhaps the most comprehensive overview study of BNs applications (Weber et al. 2010) incorporated information about most methods available for working with HBNs up to the date of its publication. Given the amount of research available on the theoretical and practical aspects of BNs, it is not surprising that even such an impressively large study overlooked a couple of methods and their applications. For instance only one application of NPBNs was mentioned, without any specification of the methodology applied. We would like to complement their presentation of NPBN applications. We aim to provide a short introduction to the NPBNs, present a couple of theoretical refinements and an overview of their use in practice. Moreover a number of difficulties and challenges occur when working with HBNs in general, and with NPBNs in particular. We address these challenges from the practitioner’s perspective, and in this way we provide guidance in using the NPBN framework. The remainder of the paper is organised as follows. Section 2 provides some background on methods available for HBNs. Section 3 describes the NPBN model and gives details about quantification and inference using NPBNs. The main theorem of NPBNs is reformulated and proved for more general settings. Two NPBN structure learning algorithms are detailed. The first one represents the initially proposed algorithm used in most of the applications presented later in the paper, whereas the second one is a new proposal that shows promise. NPBNs have been or are currently being used in at least twelve professional applications. Section 4 reviews these applications. Two of the applications are described in more detail with the intention of presenting some more theoretical details on real examples. The other applications are only very briefly presented, with more details given in the Appendix. Conclusions and directions for future research are given in Section 5.

2

Hybrid Bayesian Networks

One of the first ideas in dealing with HBNs is the use of the conditional Gaussian model (Pearl 1988; Shachter and Kenley 1989; Lauritzen 1992). The price of the discrete-normal HBNs is the restriction (of the continuous part) to the joint normal distribution. Some of the exact algorithms for inference designed for discrete BNs were extended to discrete-normal HBNs. Approximation algorithms are also available(e.g. Lerner 2002). Nevertheless uncertainty distributions need not conform to a particular parametric form. When the joint normality assumption is not appropriate, the most common method to deal with continuous variables in HBNs is discretization. When discretizing a continuous variable, a large number of partitions should be used in order to obtain a reasonable approximation. Unfortunately, in large, complex structures this leads to extremely large conditional probability tables that have to be quantified in a defendable way. Enough data for such a quantification task are very rarely available, hence the compromise is to use a small number of partitions (very often two per variable) to approximate continuous variables. More2

over, in many discretization techniques, the discretization is done locally, without taking into account the entire dependence structure of the variables represented by the nodes of the graph. This may lead to a poor control of the global error in the model (Langseth et al. 2009; Kozlov and Koller 1997). Some of these problems are resolved by using dynamic discretization (DD) (Neil et al. 2007). Nevertheless the quantification of a hybrid structure in the context of DD is not often discussed in the existing literature. Other methods to deal with HBNs include Markov Chain Monte Carlo simulations, variational methods, enhanced BN, and mixtures of truncated basis functions (MoTBFs). A review of some of these methods with their advantages and disadvantages is presented in Langseth et al. (2009) and Langseth et al. (2012). It can be argued that the MoTBFs method equates the need for good approximations in the tail of distributions with a relatively low computational complexity. MoTBFs are generalisations of the discretization technique, where an arbitrary density is approximated using generalized Fourier series approximations. Moreover, the MoTBFs framework also generalises the mixture of truncated exponentials (Langseth et al. 2010) and the mixtures of polynomials (Shenoy and West 2010) frameworks. In spite of its theoretical appeal, the MoTBFs can be extremely challenging in practice, and at the moment its applications are mainly hypothetical examples, rather than real life problems. The particular case of MTEs can handle small structures or simple ones (e.g. Aguilera et al. 2010; Cobb 2009; Cobb 2010). To our knowledge, all applications of the MTE method include networks containing less than eleven nodes. In practice, the quantification of MoTBFs requires large amounts of data (Fernandez et al. 2013). When data are not available, incorporating expert opinion could be somewhat cumbersome since the parameters of the MoTBF densities do not have a meaning. Instead the experts could provide full distributions and this information can be used in the quantification using the procedure proposed in Fernandez et al. (2013). Nevertheless, this procedure only applies to marginal densities, and it has to yet be extended to conditional densities, which in most BN models outnumber by far the marginal densities required. Another approach in dealing with HBNs combines BNs and structural reliability methods to create a new computational framework, called enhanced BN (eBN) (Straub and Der Kiureghian 2010). An eBN is transformed into a reduced BN with discrete nodes only, for which existing exact methods of inference can be used. Since inference can only be made on variables that are present in the reduced BN, the outcome space of continuous variable on which inference is desired must be discretized. Therefore, we are dealing with yet another instance of the discretization method. The quantification of an eBN is not discussed, neither from the perspective of data, nor from that of experts’ opinion. Even though the eBN methodology can handle inference in large structures - e.g. of a 115 variables (Straub and Der Kiureghian 2008; Straub and Der Kiureghian 2009) - at the moment its applications are again only hypothetical examples. All methods dealing with HBNs briefly mentioned above build upon a ”classical” construction of a joint distribution of a set of variables represented by a DAG that uses marginal distributions for the nodes with no parents, and conditional distributions for child nodes. The methods could very roughly be split into fully parametric methods (where the joint distribution of the continuous part belongs to a parametric family, e.g. the joint normal distribution) or discretization methods. It soon becomes apparent that the Achile’s heal of these methods is, in most cases, a transparent, reliable and defensible quantification process. Only after the quantification problem is solved, should one worry about performing inference using exact or approximate methods. A method that handles HBNs, but departs from the ”classical” approach is called non parametric Bayesian networks (NPBNs). It was introduced in Kurowicka and Cooke (2004) and extended in Hanea et al. (2006) and Hanea et al. (2010). NPBNs construct the joint distribution of a set of variables represented as a DAG by coupling marginal distributions of all variables with the dependence structure constructed from bivariate pieces of dependence. In NPBNs nodes are associated with arbitrary distributions and arcs with (conditional) one-parameter copulae (Joe 1997; Nelsen 1999). The quantification of NPBNs reduces to the quantification of a number of marginal distributions equal to the number of variables and a number of (conditional) dependence parameters equal to the number of arcs of the NPBN. These can be estimated from data or elicited from experts. Elicitation protocols are available, and have been successfully applied with domain experts in a number of applications presented in Section 4. A lengthier, yet very general exposition of NPBNs follows in the next section. For the mathematical details we refer to the original references. When it comes to applying the NPBN framework (or any other model obviously) difficulties are nevertheless unavoidable. Even though very often the ”classical” approaches encounter greater difficulties when applied to real

3

life problems, the mere-exposure effect makes practitioners prefer them anyway to completely different methodologies. The lack of personal experience with an approach could (for good reasons) be a great barrier in the decision to use it. Nevertheless documented experiences of other practitioners can provide a guide through the advantages and disadvantages of a fairly different approach. The following sections can be regarded as such a guide for NPBNs.

3

Non-Parametric Bayesian Networks

The NPBN methodology was initially developed for purely continuous BNs. NPBNs associate nodes with random variables for which no marginal distribution assumption is made, and arcs with one-parameter conditional copulae (Joe 1997), parameterised by Spearman’s rank correlations. The (conditional) copulae are assigned to the arcs of the NPBN according to a protocol that depends on a (non-unique) ordering of the parent nodes. The main result of NPBNs states that a particular choice of conditional copulae, together with the one-dimensional marginal distributions and the conditional independence statements implied by the graph uniquely determine the joint distribution, and every such specification is consistent (Hanea et al. 2006). The particularity of the chosen copula refers to having the so-called zero-independence property. A copula is said to have this property if zero correlation entails the independent copula. When using a copula with the zero independence property additional (conditional) independence statements are obtained by assigning rank correlation zero to an arc. In other words, there are (conditional) independence statements implied by the graphical structure and (conditional) independence statements implied by the specification. In many applications (e.g. in the financial sector) having to represent dependence only through copulae with the zero independence property places unrealistic assumptions on other features of the dependence structure (for example the t-copula that is often used in financial applications does not have the zero independence property, but it does represent tail dependence). For this reason, the necessity of this condition was further investigated. A number of definitions and equations, and the previous formulation of the NPBNs’ main result are needed in order to introduce the new formulation of the main theorem. Let us consider a BN on n variables. Then the factorization of the joint density in the standard way (following the sampling order 1, . . . , n) is: f1,...,n (x1 , ..., xn ) = f1 (x1 )

n Y

fi|P a(i) (xi |xP a(i) ),

(3.1)

i=2

where f1,...,n denotes the joint density of the n variables, fi denotes their marginal densities, and fi|j denotes conditional densities. Each variable Xi is represented by the node i. The parent nodes of i form the set P a(i). Assume P a(i) = {i1 ...ip(i) }. We associate the arcs ip(i)−k −→ i with the conditional rank correlations:  ri,ip(i) , k=0 (3.2) ri,ip(i)−k |ip(i) ,...,ip(i)−k+1 ,

1 ≤ k ≤ p(i) − 1.

The assignment is vacuous if {i1 ...ip(i) } = ∅. Assigning (conditional) rank correlations for i = 1, ..., n, as above results in associating every arc of the NPBN with a (conditional) rank correlation between parent and child. The following theorem shows that these assignments are algebraically independent and they uniquely determine the joint distribution for a particular choice of copulae. Theorem 3.1. Given: 1. a directed acyclic graph (DAG) with n nodes specifying conditional independence relationships in an NPBN; 2. n variables X1 , ..., Xn , assigned to the nodes, with invertible distribution functions F1 , ..., Fn , 3. the specification (3.2), i = 1, ..., n of conditional rank correlations on the arcs of the NPBN; 4

4. a copula realizing all correlations [−1, 1] for which correlation 0 entails independence; the joint distribution of the n variables is uniquely determined. This joint distribution satisfies the characteristic factorization (3.1) and the conditional rank correlations in (3.2) are algebraically independent. We will further propose a modification of Theorem 3.1 as to allow for various types of copulae (see Appendix 1 for the proof). If we did not specify the conditional independence statements as zero rank correlations, but as (conditional) independent copulae, then the (conditional) rank correlations associated with the arcs could be realized by any copula that realizes all correlations [−1, 1]. For instance, an NPBN could be quantified with a mixture of (conditional) independent copulae and t-copulae. Theorem 3.2. Given: 1. a directed acyclic graph with n nodes specifying conditional independence relationships in an NPBN, 2. n variables X1 , ..., Xn , assigned to the nodes, with invertible distribution functions F1 , ..., Fn , 3. the specification (3.2), i = 1, .., n, of (conditional) rank correlations on the arcs of the NPBN, 4. any copula realizing all correlations [−1, 1], 5. the (conditional) independent copula realizing all (conditional) independence relationships encoded by the graph of the NPBN; the joint distribution of the n variables is uniquely determined. This joint distribution satisfies the characteristic factorization (3.1) and the conditional rank correlations in (3.2) are algebraically independent. Of course this reformulated theorem still accommodates the use of copulae with the zero independence property, since such copulae reduce to the independent copula for correlation zero. What might appear as a slight modification of the original theorem may indeed bring a lot more flexibility for modelling multivariate distribution that allow more complex dependence structures. As mentioned earlier the possibility of using the multivariate t-copula, with different tail dependence for each pair of (conditional) variables allows the modeller to capture the phenomenon of dependent extreme values. It follows that in order to quantify NPBNs one needs to specify marginal distributions and arbitrary (conditional) copulae. The marginal distributions can be obtained from data or experts (Cooke 1991). Even though the empirical marginal distributions are used in most cases, parametric forms can be also fitted. The (conditional) copulae used in this method are parametrised by constant (conditional) rank correlations that can be calculated from data or elicited from experts (Morales et al. 2008). The conditional rank correlations need not (in general) be constant. Nevertheless they are usually assumed constant for convenience. Each parent-child influence is associated with a (conditional) rank correlation. For each child with more than one parent, the set of respective parents can be ordered in a non-unique way. When experts participate in building the structure of the NPBN they can agree upon an ordering of the parents that corresponds to their decreasing strength of influence on the child. Some other times, the order of the parents can be driven by the availability of data or another criterion. In general there is no unique procedure or criterion to build the qualitative part of the model which corresponds to the DAG of the NPBN. The rank correlation has a number of attractive properties, among which: it measures monotone dependence, rather than just linear, and it is independent of the marginal distributions. The (conditional) bivariate copulae parametrised by these correlations are then used as building blocks of the joint (multivariate) distribution. In principle, different copulae could be used for different arcs, but only the joint normal copula affords the advantages of rapid calculations/inference in large and complex problems. The theory of continuous NPBNs was extended to include ordinal discrete random variables, that is, variables that can be written as monotone transforms of uniform variables (Hanea et al. 2007). The practical limitations of this extension are discussed for a particular application. The name NPBN is somewhat inappropriate but it is used to stress the fact that the joint distribution is specified via marginal distributions, upon which no restrictions are placed, and the dependence structure is given in terms of a non parametric measure of dependence. Mixed continuous-discrete NPBNs are implemented in the UniNet software, initially developed by The Technical University of Delft (TU Delft). Uninet is freely downloadable for academic purposes from http://www.lighttwist.net/wp/uninet. 5

3.1

Quantification/Specification

Assume the DAG of an NPBN is known. Then the marginal distributions of the variables and the (conditional) rank correlations associated with the arcs need to be quantified. A copula choice needs to be made. When data is available all of the above could be estimated/fitted from data. In some cases, physical models are available, and based on them, synthetic data can be generated and used for the quantification of the NPBN. Otherwise structured expert judgement needs to be employed. 3.1.1

Data

There is not much to be said about estimating marginal distributions from data. One could either use the empirical marginals, or apply standard techniques to fit a parametric distribution. Fitting parametric distributions might improve models where estimating extremes is necessary. Conditional rank correlations are not always estimated from data directly. Rather, given a copula, these can be obtained from conditional exceedance probabilities (Morales et al. 2008). Let us now assume that the DAG is not known, and it is to be learned from data. The authors of (Hanea et al. 2010) propose a structure learning algorithm from an ordinal multivariate data set that may contain a large number of variables. The distinctive feature of this algorithm is that the onedimensional marginal distributions are taken directly from data, and the model assumes only that the joint distribution has a normal copula. Intuitively, when variables are joined by the normal copula, one could transform each marginal distribution to a standard normal to obtain a joint normal distribution with the same rank correlation structure as the original one. The concepts of learning and validation are closely connected, as indeed the goal is to learn an NPBN that is valid. The learning method performs two validation steps; it validates that the joint normal copula adequately represents the multivariate data, and that the NPBN is an adequate model of the saturated graph. Maybe it would be more appropriate for the second step to be called a verification (rather than validation) step, since the model is not tested on a different data set than the one the model was calibrated on. This should indeed be the aim of any validation technique. Nevertheless, in our experience there is rarely enough data available for performing both an adequate quantification and a full validation study. Most often the validation step is relaxed in favour of a better quantification. For both validation (verification) steps an overall measure of multivariate dependence (on which statistical tests can be based) is needed. A suitable measure in this case is the determinant of the rank correlation matrix (Hanea et al. 2010). The determinant is 1 if all variables are independent, and 0 if there is linear dependence between the variables transformed to standard normals. A measure that generalises the determinant in order to capture general dependence, and might sound more attractive to use for the validation steps, is the mutual information (Shannon and Weaver 1949; Cover and Thomas 1991). Nevertheless calculating the mutual information from data is cumbersome; moreover, for the normal copula there is a known relationship between the determinant and the mutual information, hence calculating one is equivalent with calculating the other. There are a number of essential concepts used in the calculations/simulations needed in the validation steps. Many are characteristic to the normal copula. Maybe the most important is the relationship between the rank correlation and the product moment correlation given by Pearson’s transformation (Pearson 1907) (see the Appendix 1 for the exact formulae). For the first validation step we distinguish 2 determinants. DER is the determinant of the empirical rank correlation matrix. DNR is the determinant of the empirical normal rank correlation matrix, which is the rank correlation matrix obtained by transforming the marginals to standard normals, and then transforming the product moment correlations to rank correlations. DNR will generally differ from DER because DNR assumes the normal copula, which may differ from the empirical copula. A statistical test for the suitability of DNR for representing DER is to obtain the sampling distribution of DNR and check whether DER is within the 90% central confidence band of DNR. One can do that as follows: 1. Compute DER by first transforming the marginals to uniforms and then calculating the product moment correlation of the transformed variables 2. Construct the normal version of each variable using the marginal distributions and the standard normal distribution (Zi = Φ−1 (FXi (Xi ))) and calculate the product moment correlation matrix

6

3. Compute DNR using the rank correlation and the product moment correlation 4. Re-sample the ”normal” data to obtain the empirical distribution of DNR and extract the 5th and the 95th quantiles of this distribution 5. Compare DER with the quantiles from step 4. If DER is within these bounds, do not reject the joint normal copula, otherwise reject If the normal copula assumption is not rejected on the basis of this test, one may attempt to build an NPBN which represents the DNR parsimoniously. For small data sets, the normal copula assumption is rarely rejected on the basis of this test. Nevertheless, for small dataset, this is true for any other copula assumption. On the other hand, for large datasets, the determinant validation test is rather severe. In our experience there are advantages in using the normal copula even when it is not a perfect fit, and would be rejected on large datasets; but this is only true for applications where tail dependence is not expected. Note that the saturated DAG will induce a joint distribution whose rank determinant is equal to DNR, since the NPBN uses the normal copula. However, many of the influences only reflect sample jitter and they will be eliminated from the model. The main reason for choosing the determinant of the rank correlation matrix as a multivariate dependence measure is that an approximation of it factorizes on the arcs of the NPBN. Small conditional rank correlations do not significantly contribute to the value of this approximated determinant, rather its value is driven by the largest conditional rank correlations. Moreover, the zero conditional rank correlations will correspond to the conditional independence statements encoded in the DAG structure. Let us return to the learning procedure. The NPBN will be built by adding arcs between variables only if their rank correlation is larger than a threshold value. Unfortunately, this threshold value is application dependent, and there is no universal value that we can advise. Rather, this value will be determined by the algorithm itself. The result of introducing arcs to capture causal or temporal relations is called a skeletal NPBN. Let DBN denote the determinant of the rank correlation matrix of an NPBN using the normal copula. The second validation (verification) step is similar to the first. 1. Construct a skeletal NPBN 2. Re-sample to obtain the distribution of DBN 3. If DNR is within the 90% central confidence band of DBN, then stop, else continue with the following steps 4. Find the pair of variables such that the arc between them is not in the DAG and their rank correlation is greater (in absolute value) than the rank correlation of any other pair not in the DAG. Add an arc between them and recompute DBN together with its 90% central confidence band 5. If DNR is within the 90% central confidence band of DBN, then stop, else repeat step 4 The resultant NPBN may contain nodes that have more than one parent. If the correlations between the parents of a node are neglected in the NPBN (i.e. if the parents are considered independent), then DBN will be different for different orderings of the parents. These differences will be small if the neglected correlations are also small. In general, there is no ”best” model; the choice of directionality may be made on the basis of non-statistical reasoning. Some influences may be included because the user wants to see these influences, even though they are small. The procedure for building an NPBN to represent a given data set is not fully automated, as it is impossible to infer all directionality of influence from multivariate data. Because of this fact, there are different NPBN structures that are wholly equivalent, and many non-equivalent NPBNs may provide statistically acceptable models of a given multivariate ordinal data set. More details about this algorithm are given in Section 4.1.2, where a particular application is presented. Nevertheless an improved, semi automated version of the learning algorithm is proposed in Blauw (2014). A sketch of this algorithm is as follows: 7

1. Start with an empty graph (nodes only) 2. Compute an independence matrix recording pairs of independent nodes. Use a kernel-based independence test (Zhang et al. 2012) to determine pairs of independent variables 3. Add N edges between variables with the largest absolute correlations. An initial value for N could be half of the number of variables. Make sure no edge is added between independent variables 4. Direct edges using four different rules (similar to those in the IC and PC algorithms) based on avoiding collides, avoiding directed cycles and representing (conditional) independences found with the kernel-based conditional independence test proposed by Zhang et al. (2012) 5. Direct all remaining edges randomly. A different approach for this step involves expert judgement and is proposed in (Sneller 2013) 6. Compute DBN together with its 90% central confidence band 7. If DNR is within the 90% central confidence band of DBN, then stop, else repeat step 3 for N/2 The above algorithm was evaluated on three datasets, one artificially created and two publicly available datasets. Comparisons with the well-known PC algorithm showed promise, hence we are confident to consider the above algorithm as an improved version of the original one for at least three reasons: • it uses a structured way to represent (conditional) independences along with the largest dependences • it is semi-automated • it explicitly incorporates the possibility of enhancing the structure learning process using expert opinion Both algorithms presented above are based on the assumption of the normal copula. Even though this assumption is extremely convenient from a computational point of view, and it is a very good option for smooth, symmetrical dependence structures, for applications where dependence between extremes is present and important, the normal copula is a poor choice. Nevertheless Theorem 3.2 permits the use of the t-copula as an alternative. The properties of the t-copula allow very similar calculations and sampling procedures to the ones used under the normal copula assumption (Jager 2013; Embrechts et al. 2003). Hence, the extension of the NPBN methodology is imminent. 3.1.2

Structured Expert Judgement

When field data is not available for the quantification of an NPBN, the alternative is to use data provided by experts. The subject of eliciting univariate margins from experts has been extensively studied. One of the most widely used techniques for this kind of elicitation is the classical model (Cooke 1991). The elicitation of dependence and the combination of experts’ individual estimates for dependence measures are based on the classical model. In practice, two techniques for eliciting (conditional) rank correlations have been used. The first technique consists in asking experts about conditional probabilities of exceedance and the second one about ratios of rank correlations. Besides questions about variables of interest and their dependence, experts are also asked calibration questions. These are questions for which the answer is known to the analyst, but not to the expert at the moment of the elicitation. As in the classical method, the answers to calibration questions are used to assign weights to individual expert opinions. The combined opinion consists of a weighted average of individual ones. This however does not necessarily preserve the conditional independence relationships encoded by the graph of the NPBN. In Morales-N´ apoles et al. (2014) the two methods for eliciting conditional rank correlations via related quantities are compared in order to obtain insight about which of the two renders more accurate estimates of conditional rank correlations. The results shows that good performance in uncertainty assessments (as measured by the calibration score in the classical model) does not automatically translate into good performance in dependence estimates. This observation relates to the d-calibration (dependencecalibration) score introduced in Morales-N´apoles et al. (2014). This score has an objective similar to the

8

calibration score from the classical model except that instead of evaluating expert opinions in terms of uncertainty it evaluates expert opinions in terms of statistical dependence estimation. It is observed there that analogously to uncertainty estimates, combining experts’ estimates of dependence according to their performance results in better estimates of the dependence structure. For details the reader is referred to Morales et al. (2008), Morales-N´ apoles (2009) and Morales-N´apoles et al. (2014). More application specific details are presented in Section 4.1.1.

3.2

Inference

After the model parameters are specified, the joint distribution becomes attainable. Since no analytical/parametric form of the joint distributions is available, the only way to stipulate it is by sampling it. A general sampling procedure for NPBNs is presented in Hanea et al. (2006). Any one-parameter copula could, in principle, be used. In large structures, that contain large undirected cycles, when arbitrary copulae are used, the sampling procedure involves numerical evaluations of multiple integrals. The bigger the undirected cycle is, the larger the number of multiple integrals that have to be numerically evaluated. This may constitute a big disadvantage for the use of NPBNs in (near) realtime decision support problems. Nevertheless, the disadvantage mentioned above vanishes when the normal copula is used. When the Normal copula is a fair assumption for a particular application one may operate on a joint normal distribution with the same rank correlation structure as the original one. Hence we can use the properties of the joint normal distribution and transform the (conditional) rank correlations, associated to the arcs of the NPBN, to conditional product moment correlations by inverting Pearson’s relationship (see appendix). For the joint normal distribution, the conditional and partial correlations are equal, and the correlation matrix can be calculated from the partial correlations (Yule and Kendall 1965). We sample from a joint normal distribution with given correlation matrix and transform back to the original margins. In this way we obtain a sample from the joint distribution of the initial variables, together with the dependence structure realized by the normal copula. This results in a dramatic decrease in the computational time. For examples and comparisons we refer to Hanea et al. (2006). An alternative sampling procedure based on the t-copula could follow similar steps (Jager 2013). There are several ways of performing inference in NPBNs. When arbitrary copulae are used, similar methods as for sampling could be employed after evidence becomes available. If the DAG contains undirected cycles, multiple integrals need to be evaluated for each sample, and for any new conditionalization. Hence, this approximate method of inference is suitable for structures without large undirected cycles. When this is not the case, the advantages of fast (exact or approximate) updating algorithms for discrete BNs are decisive. As mentioned in Section 2, the main criticism of the discretization methods is the extensive assessment burden that leads to indefeasible quantifications. In contrast, the quantification of NPBNs requires no more than the marginal distributions and a number of (conditional) rank correlations equal to the number of arcs in the DAG. Once the NPBN is quantified, the joint distribution can be extensively sampled (e.g. to produce hundreds of thousands of samples) and only after such a rich synthetic data set is produced the variables may be discretized using a large number of partitions. In this way, the reduced assessment burden and modelling flexibility of the NPBNs can be then combined with the fast updating algorithms of discrete BNs in a hybrid method presented in Hanea et al. (2006). The last and fastest way of performing inference in an NPBN is in the particular case in which the normal copula is used to realise the rank correlations. When the normal copula is assumed, by transforming the marginals to standard normals we are essentially arriving at a joint normal distribution. Any conditional distribution will then be normally distributed with known mean and variance (Whittaker 1990) and inference can be performed analytically. Finding the conditional distribution of the corresponding original variables will just be a matter of transforming back using the inverse distribution function of the variables and the standard normal distribution function (Hanea et al. 2006). This can be seen as an exact propagation method. Examples of inference on a specific problem are presented in Section 4.1.2.

9

4

Applications of NPBNs

This section introduces twelve of the largest projects that use the NPBN methodology. Nine of these applications are divided according to the different subjects they cover: risk analysis, reliability of structures, properties of materials. Two applications of temporal/dynamic NPBNs: one that tackles a parameter estimation problem, and another related to traffic prediction, are coupled in Section 4.4. The last part of this section presents three ongoing research projects. In all NPBNs, the (conditional) rank correlations are realised by the normal copula. This assumption was not restrictive in any way for the present applications. The respective dependence structures did not exhibit particularities that would make the use of the normal copula inappropriate. Hence the inference is exact and very fast. Nevertheless situations where the normal copula is not an appropriate assumption are easy to imagine. All it takes is an asymmetry in the dependence structure. A discussion on the consequences of erroneously choosing a Gaussian copula (not in NPBNs though) is presented in Huibregtse et al. (2015).

4.1

Risk analysis applications

The following five applications cover a variety of problems associated with the risks of the failure of dams, pollution, food consumption, air transport and building fires. Only the first two applications are discussed in detail, since they cover most of the techniques necessary in working with NPBNs, which were briefly introduced in the previous section. The other three applications are only sketched. Nevertheless we aim to point out the difficulties and the lessons learned from each of the application. A few more details about each of the applications can be found in the Appendix. 4.1.1

Earth Dams Safety in the State of Mexico

Motivated by the large flooding observed in the South East of Mexico in November 2007 (Aparicio et al. 2009), the Autonomous University of the State of Mexico, in cooperation with TU Delft, developed a model for earth dam safety in the State of Mexico (Delgado-Hern´andez et al. 2014; Morales-N´apoles et al. 2014). One of the central motivations for carrying out this investigation was the scarcity of systematic research to date, to comprehensively examine the variables that influence dam failures, describe their dependence and combine them into a model that allows to assess risks and to evaluate risk reduction alternatives for the case of Earth Dams in the State of Mexico. The application of the model is extensively discussed in Delgado-Hern´andez et al. (2014). In this section we highlight some of the main features of the model focussing on the quantification through structured expert judgement. Seven dams were considered in the analysis. The NPBN model contains eleven continuous variables and it is shown in Figure 1. The structure of the model was decided by experts. Variables are roughly grouped in ”causes” of failure (seismicity, rain and maintenance), failure modes (over-topping, breaching loss of global stability and piping) and consequences (flooding and costs). Six of the variables are quantified using structured expert judgement, four come from data and the total costs are the sum of economic and human costs. Sixteen (conditional) rank correlations where elicited from experts. In total four experts (different from the experts who decided the structure of the model) participated in the quantification of the model. Every arc in the NPBN is assigned a (conditional) rank correlation between parent and child. To illustrate the protocol of assigning (conditional) rank correlations, consider variable X4 with its parents X1 , X2 and X3 , from Figure 1. A (non-unique) ordering of parents has to be decided. In this particular case the parents were ordered as X3 , X2 , X1 by the analysts. A rank correlation between the first parent and its child is assigned first, in this case, the rank correlation denoted by r3,4 . Then the rank correlation between the child node and its next parent (in the decided order) is considered in the conditional distribution of X4 and X2 given X3 (the previous parent), and is denoted by r4,2|3 . The rank correlation between the child node and its last parent is r4,1|2,3 . The marginal distributions and the dependence structure are quantified separately. Data was available for the quantification of four marginal distributions, namely X1 , X2 , X5 and X7 . For six variables the classical model (Cooke 1991) of structured expert judgement was used. The classical model is a performance based linear pooling (weighted average) model. Loosely speaking, in addition to the variables of interest, experts are queried about calibration variables (see Section 3.1.2). 10

The objective of eliciting calibration variables is to evaluate experts’ performance as uncertainty assessors. Based on their performance they receive a score which is later used as a weight in the weighted combination of their opinions. An example of one question of interest that was asked in order to assign the distribution of variable X3 is: ”What is the current (correct but unknown) number of years between maintenance activities, that would lead any of the seven dams of interest, to an ’as good as new’ condition?. To express the uncertainty associated with the response, please provide the 5th, 50th, and 95th percentiles of your estimate” The format of the questions for the rest of the variables (including calibration variables) is similar. The elicitation of the (conditional) rank correlations required in the model was carried out through structured expert judgement as well. The literature available to guide researchers in the elicitation of a joint distribution is much less than that available for the elicitation of univariate distributions. Methods for eliciting rank correlations from experts have been proposed and used in the past (Clemen et al. 2000). One way is to directly ask experts for an estimate of the rank correlation between pairs of variables. Another way is asking experts for estimates of some other quantity, for example, a conditional probability of exceedence, or probabilities of concordance or discordance, prior to estimating, under certain assumptions (e.g. an underlying copula), the rank correlation of interest. In Morales et al. (2008) the elicitation of rank correlations from domain experts is carried out through conditional probabilities of exceedence, an approach also applied in Hanea et al. (2012). When the (conditional) rank correlations of one child with many parents should be elicited from experts, this approach would require the elicitation of conditional probabilities on a large number of conditioning events, making it un-intuitive for some domain experts. Additionally, while not conclusive, previous results indicate that the most accurate way to obtain a subjective measure of bivariate dependence is simply to ask the expert to estimate the correlation between the two variables in question (Clemen et al. 2000). A total of 20 questions were asked to every expert. For each child node, experts were asked to rank parent variables according to some criterion of their preference, e.g. the largest unconditional rank correlation with the child in absolute value. These questions were meant to help experts in their next assessments. Then for unconditional rank correlations experts would assess P (Xchild > median|Xparent > median). This question, for this application, and for the rank correlation r4,3 is: Consider a situation in which the number of years to undertake maintenance action (X3 ) is above its median (30 years). What is the probability that the loss of global stability (X4 ) is also above its median (1.66)? Observe that, for this question, the median value for the loss of global stability is equal for all experts since it comes from data. The number of years to undertake maintenance actions, however, comes from expert judgement and is different across experts. When combining experts’ assessments it is necessary to take this kind of differences into account. Once the expert has given an assessment for the probability of exceedence, the analyst then calculates the corresponding rank correlation under the copula assumption. After the experts have provided an assessment for the unconditional rank correlation of a given child node to its first parent, then the ratios of the rank correlations of the child with the rest of its parents, to the first rank correlation elicited, are elicited. In this example, the question would be: Given your previous estimates, what is the ratio of r4,2 r4,3 ?. Let us continue with the answers of a particular expert called A. Her answer will depend on the previous estimate which in this case refers to a value of r4,3 = −0.57 computed during the elicitation, and given to the expert. Observe that r4,3 is negative. Hence if the expert believes that the loss of global r stability is positively correlated with rainfall rate, then the ratio r4,2 should be negative. On the other 4,3 r hand, negative correlation between the loss of global stability and rainfall rate will correspond to r4,2 > 0. 4,3 r4,2 Notice that | r4,3 | > 1 corresponds to the expert’s belief that |r4,2 | > |r4,3 | and an analogous situation is r r4,2 observed when | r4,2 | < 1. It is possible to calculate that the ratio r4,3 belongs to the interval (−1.38, 1.38). 4,3 Because X2 and X3 are independent, this ratio is in a symmetric interval around zero. However, if they were correlated, additional constraints would be present for its assessment. The correlation r4,2 will be r4,2 restricted by the experts previous answer. In this example, expert A has stated that r4,3 = 0.25. She then believes that the rank correlation between the loss of global stability and rainfall rate is negative and smaller in absolute value than r4,3 . In particular, r4,2 has a value of −0.14. 11

Loss of global stability has one last parent node for which the conditional rank correlation r4,1|2,3 needs to be specified. Thus the expert would be confronted with the following question: Given your r ?. Similar calculations and restrictions will be applied in order previous estimates, what is the ratio of r4,1 4,3 to obtain a value for r4,1|2,3 . The full quantification of all (conditional) rank correlations follows the same procedure as the one described for the elicitation of r4,3 , r4,2|3 and r4,1|2,3 . For each expert an NPBN is obtained. Once all joint distributions, of all experts are available, these are combined using the weights derived from the calibration questions, as explained previously. Notice that the ordering of variables decided by each expert may be different. This ordering may also be different to the original ordering decided upon by the group of analysts. Since a joint distribution has been obtained for each expert, the analysts may compute the answers that each individual expert would have given, had they been asked questions regarding a fixed ordering of variables in the NPBN. That is, the ordering decided by the group of analysts which, in this case, corresponds to the one in Figure 1. Hence, once the DAG has been agreed upon, the ordering of the nodes is part of a parametrization but not a restriction during the quantification process through expert judgements. After the full quantification of the model, the NPBN is ready to be used for inference. Let us see what the model predicts for one particular dam, namely the Jos´e Antonio Alzate dam (Morales-N´apoles et al. 2014). Because of its geometry and year of construction, there are two variables that can immediately be fixed for the dam under analysis. Firstly, the loss of global stability can be quantified by means of the safety factor of the structure, i.e. 1.95, which is to be used as the observed value in the model. Secondly, the age of the dam is 46 years, which can be associated with the number of years between maintenance activities. Assuming that there have not been any conservation actions since its final construction year, the maintenance observed value can then be assumed to be 46. The left hand side of Figure 2 presents the Jos´e Antonio Alzate dam model. This model was obtained by conditioning on the above mentioned variables. After conditioning, the rest of the probability density distributions are updated. The grey distributions are the unconditional distributions provided for comparison, and the black ones are the conditional distributions of the rest of the variables given the new evidence. To show how the model can be used further, suppose that a hypothetical extraordinary rainfall rate of 15 mm/day takes place in the dam zone and, concurrently, it is known that the seismic frequency in the region corresponds to eight earthquakes ε 5.5 per year. The right hand side of Figure 2 shows the result of this exercise. The anticipated flooding value has increased from 1.71 m/day to 2.22 m/day. Similarly, the predicted human cost moved from 13.9 to 14.7 USD million, and the economic loss changed from 29.4 to 30.1 USD million. This means that the intensification of rain, and the presence of earthquakes at the same time, are expected to produce higher levels of water in the potential flood area and, consequently, superior amounts of both human and economic losses. The application of the model has been discussed at length in Morales-N´apoles et al. (2014) where a risk analysis associated with the seven dams in the region of interest was provided. The use of the model for maintenance decisions and design discharges has also briefly been discussed in Morales-N´apoles et al. (2014). The combination of expert judgements was done in this case with the usual measures for performance evaluation in the classical model. After this research it became clear that measures for evaluating experts performance in dependence assessments would be more suitable. This has led to research in MoralesN´ apoles et al. (2014). 4.1.2

Linking P M2.5 Concentrations to Stationary Source Emissions

Linking a large number of pollutant sources and ambient receptors proves to be a difficult task, given the great multiplicity of sources, the high correlation among nearby monitors, and the inconstancy of relationships between particular sources and receptors amid constantly changing meteorological conditions. We use an NPBN for linking monthly emission sources of sulfur dioxide (SO2 ) to average ambient monthly concentrations of fine particulate matter (P M2.5 ) over the entire eastern U.S. We focus on the relative spatial relationships of sources and receptors, not on their absolute locations. For each receptor we group sources into bins depending on distance and direction from the receptor1 , and these aggregations become the hypothesized causal variables in our NPBN. This study is a joint research of the Resources for the Future Institute in Washington DC, and the 1 The

letters in the names of variables from Figure 3 correspond to distances, and the numbers to compass directions.

12

Delft Institute of Applied Mathematics (Hanea and Harrington 2009). Assembly of the data set was supported by a grant from the Health Effects Institute. The full data set comprises monthly emissions of SO2 gathered from hundreds of electricity generating stations, and monthly mean concentrations of P M2.5 gathered from hundreds of collection sites in the United States over the course of seven years (1999 - 2005). We have 84 monthly averages of emitters and collectors. For each monitoring site we consider 32 variables corresponding to SO2 emission data from its vicinity. The scope is to learn an NPBN from data, on these 33 variables, that captures the dependence present in the data set. An example of such an NPBN is showed in Figure 3. The authors used the algorithm presented in Section 3.1.1. From the 33 variables showed in Figure 3, 17 were chosen. The excluded variables are the variables that correspond to areas where no power plants exist and those further than 500 miles away from the monitor, since emissions close to receptors tend to have stronger effects than emissions further away. Finally another small number of variables is excluded such that the distribution of the locations of remaining power plants is ”uniform” in all directions. The 17 variables left in the model are shown in Figure 4. So are their marginal distributions. No parametric distributions were fitted; the empirical margins were use instead. Before starting to build a model the hypothesis that the data were generated from the joint normal copula is tested. Following the procedure described in Section 3.1.1 one can calculate DER = 0.495E-05 and DNR = 1.326E-05. Based on 1000 simulations, we determine that the 90% central confidence band for DNR is [0.037E-05, 1.695E-05]. The hypothesis that the data were generated from the joint normal copula would not be rejected at the 5% level. One can now start building a model that will reconstruct the dependence structure present in the data set. The 17 variables will be nodes of a DAG. If the DAG has no arcs then DBN = 1. In general, if the DAG is not saturated, then DBN > DNR. Following the general procedure arcs are added between variables whose rank correlations are large. In this case we started with pairs of variables with a correlation greater than 0.7 (in absolute value). By doing so, the value of the DBN decreases. After a number of intermediate steps (where arcs corresponding to decreasing correlations are added), an NPBN with 35 arcs is built. Most of these added arcs correspond to the highest rank correlations. Nevertheless, the main interest is to quantify the relation between the P M2.5 concentration and the rest of the variables, hence arcs that carry information about their direct relationship are added as well. The resultant NPBN has DBN = 1.017E-03 and the 5% precentile of its distribution is 0.582E-04 which is larger than DNR. In consequence, we add 10 more arcs. The resultant structure is shown in Figure 4. The determinant of the rank correlation matrix based on the new NPBN is DBN = 1.261E-04, its 90% central confidence band is [0.059E-04, 1.296E-04] and DNR falls inside this interval. One can conclude that this NPBN is an adequate model of the saturated graph. It is worth stressing that the 17 dimensional joint distribution was described by a dataset with only 84 joint samples. For this reason the empirical distribution of DNR spans two orders of magnitude and the normal copula cannot be rejected. By allowing the dependence structure to be modelled by such copula, the NPBN model excludes the possibility of tail dependence, but it does offer the advantage of obtaining any full conditional distribution in seconds. A Gaussian BN would be a more familiar alternative to the NPBN approach, and equivalent with a standard regression analysis. Nevertheless, the emitters tend to be strongly correlated to each other and weakly correlated to the collectors, hence if we marginalize over a small set of emitters, we have many ”missing covariates” with strong correlations to the included covariates. This will bias the estimates of the regression coeficients. The NPBN method, in contrast, simply models a small set of variables, where other variables have been integrated out. There is no bias; the result of first marginalising then conditionalising is the same as first conditionalising then marginalising. Moreover, in this application the set of regressors have individually weak correlations with the predicted variable, but are collectively important. Because of the size of the dataset the confidence intervals for the regression coefficients may all contain zero. Then their collective importance would be missed. 4.1.3

Causal Models for Air Transport Safety

A project for building a model for air transport safety commissioned by the Dutch Ministry of Transport and Water Management begun on July 2005 (Ale et al. 2008). The project is known as Causal Model for

13

Air Transport Safety (CATS). The motivation for the project is the need for a thorough understanding of the causal factors underlying the risks related to air transport and their relation to different possible consequences. The NPBN representing the CATS model consists of 1504 nodes and 4979 arcs. The quantification of the model used both data and expert judgement. During this project, the need to model rank correlations between discrete variables (or between discrete and continuous variables) came to light. The main theoretical results regarding this problem are presented in Hanea et al. (2007). Moreover, many of the difficulties that arise when expert judgement is needed for both marginal distributions of variables and the dependence between them, were discovered and tackled during this project. 4.1.4

The Benefit-Risk Analysis of Food Consumption

The BENERIS project (Jesionek and Cooke 2007) was initiated in 2006 with the main aim of developing comprehensive methods and tools for the benefit-risk analysis of food consumption. An NPBN model for the benefit-risk analysis of fish consumption was build. The variables involved in the NPBN are fish consumption variables, environmental variables, demographic/personal factors, health effects, exposureresponse factors, baseline responses, and exposure variables. The DAG of the model was constructed based on expertise and known functional relationships between variables. After construction, the model was quantified using field data. The NPBN for the fish case consists of 524 probabilistic nodes, 637 functional nodes (whose corresponding variables are specified as functions of their parents in the DAG) and 1812 arcs. Two limitations of NPBNs were brought to light by this project, and they are both related to the use of functional nodes: functional nodes cannot have probabilistic children, and conditioning on a functional node can only be done when sample based. 4.1.5

The Human Damage in Building Fire

The Dutch National Institute for Public Health and the Environment (RIVM) financed a research project with the goal of showing the opportunities and challenges of developing a decision tool for fire safety in buildings. The Human Damage in Building Fire Model build by Hanea (2009) is a probabilistic decision support tool based on the NPBN methodology. The model is used to estimate the percentage of human deaths caused by fire in a building, and it contains both probabilistic and functional variables. The quantification of the model involved both experts and data. One of the most challenging parts of this project was the quantification using expert judgement. Experts’ input is needed for building the structure, for quantifying marginal distributions and (conditional) rank correlation. The availability of experts that could perform all three tasks in a reliable manner could be a limitation of this approach. Nevertheless the absence of data would jeopardize the quantification of any other model. The advantage of an NPBN model is the availability of structured protocols to elicit its parameters. The following sections sketch more applications of NPBNs in different domains. Appendix 2 offers extra information about them and respective references.

4.2

Reliability of structures

Bayesian Network for the Weigh in Motion System of The Netherlands It is estimated that about 15% of the trucks in The Netherlands are overloaded (Weg-en Waterbouwkunde 2004). Since 2001, the Ministry of Transport, Public Works and Water Management of The Netherlands has implemented the Weigh-in-Motion (WIM) system that measures axle loads, vehicular weight, speed and length of vehicles at eight locations in The Netherlands. The data provided by WIM helps in regulating overloaded traffic and aiding the decision making process for design, maintenance and in general reliability of infrastructure. The dependence structure of the data was modelled with an NPBN. The WIM NPBN consists of 705 nodes and more than 2300 arcs (Morales N´apoles and Steenbergen 2012). The model was fully quantified from WIM data.

14

Through the NPBN it has become clear that the distribution of vehicle loads is necessarily a mixed distribution, and also that the model may be used for the reliable computation of design loads as well as joint design criteria if necessary. In general this model has been used to advise the Dutch Ministry of Infrastructure and Environment mainly with respect to design loads.

4.3

Properties of materials

Technique for Probabilistic Multi-scale Modelling of Materials Investigating mechanical properties of materials is of interest from different perspectives. In particular asphalt’s properties are investigated, for example, in order to improve the quality of asphalt or improve predictions of service life time. An NPBN was used to summarize the properties of such a non homogenous material. The DAG of the NPBN was build based on expert input (about conditional independencies) and data. The model was quantified using field data (Morales N´apoles et al. 2011). It has been shown that the approach used in the model works well for materials with behaviour approximately linear elastic. The adequacy of the model for materials with different type of behaviour needs be investigated further.

4.4 4.4.1

Dynamic NPBNs Permeability Field Estimation

Offshore oil recovery heavily depends on accurate predictions of the flow of fluids like oil, water and gas in an oil reservoir. The ease of fluids’ flow through a porous media is characterized by the permeability of the rock, hence its estimation is of great importance. A reservoir simulator provides synthetic data about the joint behaviour of variables like pressure, saturation, bottom hole pressure, flow rates and permeability, at different snapshots in time. This data was used to learn both the structure and the parameters of a dynamic NPBN. A twin experiment was performed with the aim of recovering a ”true” permeability field using the dynamic NPBN methodology (Hanea et al. 2013). The dependence structure for this application involved 343 variables, and was characterised by many reasonably small pairwise dependences (rather than a few large ones). In this situation, the DAG learned was very close to the saturated DAG (many arcs were needed to represent all the dependence as measured by the determinant). This results in a large number of (conditional) rank correlations to be estimated from data. When learning conditional rank correlations, the loss in significance is a lot more rapid than for learning means and standard deviations, hence larger sample files are imperative. Moreover, the determinant of a correlation matrix of 343 variables is of order 10−200 , and the validation test are inapplicable to such numbers. Despite all these limitations, the NPBN method improved aspects of the state of the art methodology in permeability estimation (which can be easily translated to a fully parametric dynamic BN) 4.4.2

Traffic Prediction in the Netherlands

There have been efforts in the past to investigate traffic through discrete dynamic BNs. The authors of Haak et al. (2010) used a bottleneck in a Dutch highway as case study. Variables included in the model are vehicular speed and vehicular densities at different locations and different snapshots in time. The objective is to predict at a particular time, vehicular speed in a particular location, given the information available at the previous time steps. An NPBN was build in order to add levels of complexity to the previously modelled situation, e.g. more variables, continuous variables, more arcs between the nodes of the DAG. The NPBN was constructed and quantified using data (Worm et al. 2011). The quantification has been done entirely with data from the Monitoring Casco (MONICA) data base system, managed by the Dutch Ministry of Transport, Public Works and Water Management.

4.5

Summary of Completed Applications

For the sake of conciseness we have summarized the completed applications of the NPBN approach in the following table. Both nodes and arcs may be quantified using functions, data, or experts. Sometimes,

15

information from the literature is used in the quantification. This information is considered and listed under the category ”Data”.

4.6 4.6.1

Ongoing Applications Filtration Techniques

An NPBN is used to model the residuals of a deep bed filtration model applied to wastewater effluent filtration2 . Filtration models are used to support design and operation of the filters but the models are rarely validated for the full range of conditions occurring at wastewater treatment plants. Predicting the model residuals for daily and extreme conditions with an NPBN adds the aspect of reliability and uncertainty to the models and makes them more transparent for practical applications. The main challenge of this research is to decide whether the residuals of the model may be explained through an NPBN consisting of environmental variables as well as variables representing different parameters in the physical model. 4.6.2

Flood Defences

The failure probability of a levee section due the failure mechanism piping is modelled with an NPBN3 . The piping mechanism initiates when the upward water pressure in the subsurface exceeds an (uncertain) critical pressure that the soil can withstand. Coupled observations of river water level and evidence of compromised levee integrity (or lack thereof) will be used to update the probability distributions of the soil properties and the critical pressure. The main challenges of this application are computational since numerical instabilities arise when the correlation between variables approaches 1, which is what happens when we densely represent the spatial variability in a levee. This results in a sparser than desired representation of the spatial variability. 4.6.3

Train Disruptions

In the Netherlands, relatively large disruptions occur often, each time leading to a temporary unavailability of the railway system. A probabilistic approach using NPBNs is currently being developed to build a railway disruption length model4 . The NPBN is used as the dependence model between influencing variables and the disruption length. The output of the model will then be used as an input in a decision support system that optimizes the train traffic during and after a disruption. The biggest challenge of this project is modelling the dependence between discrete variables with copulae. Even though the NPBNs are in theory extended to allow hybrid structures from the perspective of the quantification and sampling, when the number of discrete variables exceeds the number of continuous ones, the conditionalisation proves very challenging. In such situations, it seems reasonable to advise the use of discrete BNs (hence the discretization of the continuous part, rather than the use of copulae for the discrete part).

5

Conclusions

We have mentioned a number of approaches to deal with HBNs, amongst which MoTBFs, eBN, and NPBNs. This paper concentrates on the NPBN methodology, hence NPBNs were presented in more detail. The presentation tried to cover three aspects: the quantification of the models, the inference, and the real-life applications. Unlike other approaches, the NPBN methodology does not build on classical BN methods, but rather proposes a new way of dealing with non-parametric continuous distributions in BNs. The nodes of an NPBN represent continuous variables and the influences are quantified with one-parameter (conditional) copulae. The NPBNs are extended to include ordinal discrete variables, but our experience only covers situations where they are outnumbered by continuous variables. When 2 This research is part of the project Integrated Filtration Techniques and it is done in collaboration with Visser & Smit Hanab BV., Grontmij Nederland, TU Delft and Waterschap Veluwe. 3 This work is part of the program Sustainable and Integral Design of Multifunctional Flood Defences, supported by the Dutch Technology Foundation STW. 4 This research is part of the Smart Operational Control Center Rail project, sponsored by ProRAILs ExploRAIL project.

16

modelling a mixed distribution with more discrete than continuous univariate margins, discrete BNs are a more appropriate choice. When data is available, marginal distributions and the parameters of the (conditional) copulae can be easily estimated. Nevertheless empirical data for models that involve substantial uncertainty is often unavailable. In these cases expert judgement is the only alternative. When performed rigorously, expert elicitation and pooling of experts’ opinions is a powerful means for obtaining rational estimates of uncertainty. The structured expert judgement method, briefly discussed in Section 2 and detailed in Section 4.1.1, is used for the quantification of the NPBN model in absence of empirical data. In Weber et al. (2010), the authors touch upon this issue in their conclusions. The structured expert judgement method for uncertainty quantification complements their exposition. Even though structured expert judgement offers an elegant and scientific way of complementing field data, the methods for eliciting dependence are still in their infancy. Much more research is needed in order to gain the same confidence we have when eliciting marginal distributions. A method for learning the structure of an NPBN from data was briefly discussed. This is based on the normal copula assumption which can be validated based on some measure of multivariate dependence. The chosen measure is the determinant of the rank correlation matrix. When the normal copula assumption cannot be validated the procedure becomes unusable. There are several ways of performing inference in NPBNs: sampling (using arbitrary copulae), discretizing and employing the fast updating algorithms for discrete BNs, or assuming the normal copula and performing analytical inference. The possibility of performing fast inference is maybe the most important feature of a HBN model. Nevertheless, fast inference in very large models is usually a challenge. The biggest advantage of the NPBN method is that it can handle hundreds of variables in a very fast manner (when the normal copula is used, inference in a model with hundreds of nodes takes few minutes) while the quantification requires a relatively small number of parameters. It is only fair to emphasize again that the choice of the normal copula might not be appropriate. The price of using different copulae is an increase in computational time and maybe the need of a slight modification of the model.

Appendix 1 Proof of Theorem 3.2 Theorem. Given: 1. a directed acyclic graph with n nodes specifying conditional independence relationships in an NPBN, 2. n variables X1 , ..., Xn , assigned to the nodes, with invertible distribution functions F1 , ..., Fn , 3. the specification (3.2), i = 1, .., n, of (conditional) rank correlations on the arcs of the NPBN, 4. any copula realizing all correlations [−1, 1], 5. the (conditional) independent copula realizing all (conditional) independence relationships encoded by the graph of the NPBN; the joint distribution of the n variables is uniquely determined. This joint distribution satisfies the characteristic factorization (3.1) and the conditional rank correlations in (3.2) are algebraically independent. Proof. Given that all univariate distributions are known, continuous, invertible functions, one can use them to transform the variables to uniforms on (0, 1). Hence, we can assume, without any loss of generality, that all univariate distributions are uniform distributions on (0, 1). The assignment of (conditional) rank correlations is vacuous if {i1 ...ip(i) } = ∅, hence, for the first node (without parents) the first term is assigned vacuously. We assume the joint distribution for {1, ..., i − 1} has been determined. The ith term of the factorization (3.1) involves i − 1 conditioning variables, of which {ip(i)+1 , ...ii−1 } are conditionally independent of i given {i1 , ..., ip(i) }. We assign: ci,ij |i1 ,...,ip(i) {Fi|i1 ...ip(i) (Xi | Xi1 , ..., Xip(i) ), Fij |i1 ...ip(i) (Xij | Xi1 , ..., Xip(i) )} = 1

(5.1)

for ip(i) < ij ≤ i − 1. Then these copulae together with the ones realizing the conditional rank correlations (3.2) are exactly those present on another graphical structure called a D-vine 5 with i variables, 5 The graphical model called vines was introduced in Bedford and Cooke (2002). A vine on n variables is a nested set of trees, where the edges of the j th tree are the nodes of the (j + 1)th tree.

17

denoted Di . The other conditional bivariate distributions on Di are already defined. It follows from well known results for vines that the distribution on {1, ..., i} is uniquely determined. Since the conditional independent copula is used to represent all conditional independence statements of the graph, we have: f1...i (1, ..., i) = fi|1...i−1 (i | 1, ..., i − 1) · f1...i−1 (1, ..., i − 1) = fi|i1 ...ipa(i) (i | i1 , ..., ipa(i) ) · f1...i−1 (1, ..., i − 1) It follows naturally that the factorization (3.1) holds.

Validation tests – essential properties and formulae Testing the normal copula hypothesis depends on the fact that the rank correlation of two standard normal variables Zi and Zj is given by the Pearson transformation of the product moment correlation (Pearson 1907): ρ(Zi , Zj ) 6 arcsin( ) = r(Zi , Zj ) = ρ(Φ(Zi ), Φ(Zj )) π 2 where r denotes the rank correlation, and ρ denotes the product moment correlation. The cumulative distribution function of Zi , applied to Zi , Φ(Zi ), returns the rank of Zi and is uniformly distributed. To determine the empirical rank correlation of a set of bivariate observations of Xi and Xj , we first construct the empirical distribution functions FXi (q) and FXj (q). The empirical rank correlation of two variables Xi and Xj is given by: r(Xi , Xj ) = ρ(FXi (Xi ), FXj (Xj )) (5.2) This involves no modelling hypotheses, and is simply computed from data. The empirical normal version of Xi is defined as: Zi = Φ−1 (FXi (Xi )) (5.3) Applying Pearson’s transformation, we find the empirical normal rank correlation of Xi and Xj to be: ρ(Zi , Zj ) 6 arcsin( ). π 2

(5.4)

This relation is characteristic to the normal distribution. If Xi and Xj are joined by the normal copula, the empirical rank correlation will approach the empirical normal rank correlation as the number of observations becomes large. Otherwise it will not.

Appendix 2: Details of Applications CATS The HBN representing the CATS model is presented in Figure 5. Structured expert judgement was needed in quantifying parts of the model. In total five pilots, five air traffic controllers and one maintenance technician participated in the expert elicitations (Morales-N´apoles 2009). 32 nodes represent discrete random variables and the rest continuous variables. Most of the univariate margins come from different sources of data, however 28 of the univariate margins have a non-parametric distribution elicited through structured expert judgment. All (conditional) conditional rank correlations required by the model were also elicited through structured expert judgment (Morales-N´apoles 2009).

A validation exercise (case validity) was performed in Ale et al. (2010). This consisted in analyzing the crash of the Turkish Airlines flight TK 1951 of February 25, 2009 occurring about 1.5 km before the intended runway at Amsterdam Schiphol Airport. The model shows how vulnerabilities in the aviation system may be detected, by systematically studying the behavior of different parts of the model.

18

BENERIS Partners involved in BENERIS project are the National Institute for Health and Welfare (THL, Finland), TU Delft (Netherlands), Oy Foodfiles Ltd (Finland), Food Safety Authority of Ireland (FSAI, Ireland), Technical University of Denmark (DTU, Denmark), Food Safety Authority of Denmark (FVST, Denmark), Lednac Ltd (Ireland) and Fundacion Privada para la Investigacion Nutricional (FIN, Spain). The work performed by TU Delft, supported by the expertise of THL, has as a product an NPBN model for the benefit-risk analysis of fish consumption. The variables involved in the NPBN are fish consumption variables, environmental variables, demographic/personal factors, health effects, exposure-response factors, baseline responses, and exposure variables. The analysis was performed for the Finnish population. Hence, Finnish databases were used in quantifying the variables of the model. In the fish consumption NPBN the influences between variables are mostly functional. For example, the exposure is calculated as daily fish consumption multiplied by concentration divided by the body weight, and the health effect variables are mathematical functions of exposure variables, exposure-response factors and baseline responses. Nevertheless, several probabilistic influences are also present in the network. They were quantified using data. Quantifying and merging all the components of the study results in the very complex graphical structure presented in Figure 6. TU Delft developed the main core of the model. THL extended it further by adding extra nodes, e.g. the disability adjusted life years (DALYs). DALYs allow expression of the estimated health effects (which have different units, e.g. IQ points, probability, incidence) on a common scale and make meaningful comparisons and conclusions. Unfortunately, the results of these comparisons are not yet available to the public. A case study is also investigated at the moment of writing this paper. The results of the case study will be published in the near future.

Building Fires The model is used to estimate the percentage of human deaths caused by fire in a building, and it is shown in Figure 7. All variables are considered continuous; seven of the probabilistic variables were quantified using information from the literature and ten of them using the classical model for structured expert judgement. All influences/arcs expressed in terms of (conditional) rank correlations were assigned using experts. Four experts participated in the elicitation procedure. Their areas of expertise are: fire prevention, people behavior and evacuation, fire safety in buildings and fire development. Two experts also participated in building the structure of the network. Validation of the model proved difficult due to the lack of available and/or complete data about past fires in buildings that involved human fatalities. The subjectivity of the evacuation exercises and lack of other, similar, integrated models are two more impediments in performing a thorough validation exercise. Nevertheless, in order to investigate and hopefully increase the confidence in the model, as well as to understand its applicability limitations, partial validation studies were carried out. A case validity was performed using the data about the fire that broke out in the K-wing of the detention center of Amsterdam’s Schiphol airport. Based on the available information (Raad 2006), the fire from the Schiphol Cell Complex was analysed using the The Human Damage in Building Fire Model. Most investigations showed that the damage to humans caused by this fire should have been expected.

WIM The dependence structure of the data provided by WIM is also very complex. Thus, the aim would be to be able to summarize the complex multidimensional-distribution with a transparent model. The HBN for the WIM system of The Netherlands has been fully quantified. An earlier and smaller version of the model is described in Morales N´ apoles and Steenbergen (2014). The current model is described in Morales N´ apoles and Steenbergen (2012) and presented in Figure 8. In Morales N´ apoles and Steenbergen (2012) possible uses of the model have been discussed. These include prediction of missing values in the database produced by the WIM system and computation of design vehicular weight in the presence of limited data.

19

The WIM NPBN consists of 705 nodes and more than 2300 arcs. Eight of the nodes correspond to discrete variables. The WIM system measures the variables of interest in four locations/highways, in the right and left directions. Thus, a total of eight sub-models for vehicles from two to eleven axles was built. Within each of the eight locations an axle load distribution for each axle of each vehicle type is considered. The discrete variable corresponds to the number of axles. Speed and length distributions are computed for each of the two to eleven axle vehicle classes. The total weight is a deterministic function of the number of axles and the individual axle loads. The separate locations in the model are connected by a variable called time to rush hour. This variable is unique for all locations. Traffic intensity variables also connect the separate locations. Contrary to the time to rush hour variable, a traffic intensity node is computed for every location. The quantification of the model has been done exclusively with data from the WIM system in the Netherlands. By quantifying a small dynamic BN for vehicular separation and combining it with the model in Figure 8 through simulation, the computation of maximum bending moments for bridges is performed. The approach allows to investigate the traffic configurations leading to the maximum bending moments. That is the probable number of vehicles with their respective number of axles, their individual axle loads, and their relative position, acting on a bridge of a given length in a particular moment of time. One of the possible extensions of the model is the implementation of a virtual WIM system for location(s) where the system is currently not available.

Multi-scale materials modelling One technique for investigating mechanical properties of materials is the use of multi-scale models (Philips 1998). Roughly, some geometrical and mechanical (Poisson’s ratio, Young’s modulus) properties of a sample of material are investigated at a given scale. These are modelled, through Finite Element Methods (FEM) calculations to extract other mechanical properties, for example von Mises stresses. These are used as input at a higher scale to repeat the process. One of the questions raised by this approach is that of the homogeneity of a sample of the material (Philips 1998). Additionally, computations for one sample only might be very expensive in time. In order to partially address these difficulties, statistical approaches using multidimensional copulas have been proposed (Greene et al. 2010). Assembling a multivariate distribution through a BN might have advantages over specifying a multidimensional copula in other ways. With these observations in mind a framework is proposed for material investigation. In 2011 The Netherlands Organization for Applied Scientific Research (TNO) initiated a project called PhysicoChemical modelling of Asphalt. Roughly, the framework proposed in the project consists of the following steps (Morales N´ apoles et al. 2011): 1. Obtain data from the laboratory in the form of images from a material at a given scale. Use these data to statistically characterize the geometrical and physical properties of the material. For example through a distribution of grain size for granular material and the correlation between grain positions in space. 2. Use the characterization done in the previous scale together with random fields in order to generate many realizations of the material of interest (Ostoja-Starzewski 2007). Every realization would be different from all others but they would share the statistics of the geometrical and mechanical properties of the material of interest. Use FEM to investigate the mechanical properties of the material after subjected to some load. 3. Use data obtained in step 2 to learn and quantify a graphical model (NPBN in this case). The model would represent a summary of material properties at different scales and may be used to perform inference. In the left hand side of the Figure 9 we present an image of a sample of asphalt at 24.7 × 24.7µm2 . The model developed in 2011 simplifies the problem by making two assumptions: firstly, asphalt is assumed to be a two phases granular material, some ”grains” (with certain mechanical and geometrical properties) and the continuum (with different mechanical properties than the granular phase). Secondly, the mechanical properties of asphalt are investigated with a linear elastic FEM model.

20

The geometrical properties of this particular sample of asphalt are investigated using a Homogeneous Poisson (HP) random field (Ostoja-Starzewski 2007). The granular phase of the material is idealized by ellipses. The centers of each ellipse conform to the HP random field. Thus we use the HP process first to select a Poisson number of centers on the domain and to position them (uniformly) on the plane. Sizes of the major and minor axes of each ellipse are assigned at random, based on a quantification performed on the original sample. In total 5000 microstructures were generated for this case and each was solved on a 50 × 50 grid with square elements. The NPBN quantified with 5000 microstructures is shown in the right hand side of Figure 9. In total eight variables and 19 direct influences between them are used. The geometrical variables are: the number of ellipses and the ratio of the area of the ellipses to the total domain (variable Volume in the right hand side of Figure 9). Further, the number of ellipse centers falling in the upper-left, upper-right and bottom-right quadrants of the domain are considered. Three mechanical variables are added to the model. These are average von Mises stresses over the whole domain, over the elements in the right hand side of the domain, and over the elements where the two idealized material phases meet. The quantified model may be used in scale transition for example by conditionalising the distribution of mechanical properties on an observed number of grains in the idealized sample, on their size and on their location relative to the four quadrants of interest. The model can easily be extended to account for other geometrical or mechanical properties of interest. All one dimensional distributions and the 19 arcs of the model were quantified with the data generated by the FEM applied to the random microstructures. The NPBN in Figure 9 was validated with the method briefly described in Section 3.1.1. Investigating ways to relax the two main assumptions made in this application are the subject of future research. The framework, represented by the NPBN in Figure 9, is however already valid for materials that show an approximate linear elastic behaviour.

Permeability estimation The problem of estimating parameters like permeability is often referred to as a history matching problem in reservoir simulation. One of the most used methodologies that addresses the history matching problem is data assimilation, more specifically, the ensemble Kalman filter method (EnKF) (Evensen et al. 2007; Aanonsen et al. 2009). Because of the EnKF’s limitations, TNO in collaboration with TU Delft initiated a project with the goal of investigating an alternative approach to the parameter estimation problem (Gheorghe 2010). This approach uses NPBNs and dynamic NPBNs. We consider a synthetic two dimensional squared petroleum reservoir (a 21 X 21 grid) with an injector well in the middle of the field and four producers, one in each corner. One fluid (typically water) is injected through the injector well and oil is pumped out through the producers. Amongst variables of interest we mention: pressure, saturation, bottom hole pressure, flow rate and permeability, corresponding to each grid cell. All variables are continuous and their marginal distributions are taken from synthetic data provided by a reservoir simulator. Using the learning procedure briefly described in Section 3.1.1, one can learn the arcs of an NPBN. Nevertheless, their directionality cannot be extracted from data. In this project, the directionality of the arcs was provided by experts in the field. A twin experiment was performed and various case studies were investigated on parts of the grid. An NPBN on 103 variables and one on 343 variables were considered. In learning the structure of an NPBN we follow the validation steps based on the determinant of the correlation matrix. The determinant of a correlation matrix of 343 variables is of order 10−200 , and the validation tests are unapplicable to such numbers. In this situation, one could use the saturated graph in order to represent the dependence structure present in the synthetic data. However, a saturated graph on 103 variables has 5253 arcs. A picture of such an NPBN would be pointless since it would be practically impossible to visualise anything. The NPBN approach for parameter estimation in reservoir engineering is still in its infancy, yet there are reasons to believe it is a promising approach. More research in this area is definitely needed, and currently undertaken.

Traffic prediction Variables included in the model from Haak et al. (2010) are vehicular speed and vehicular densities at three locations denoted A, B, C. In location A the highway is reduced from three lanes to two lanes.

21

Locations B and C are located (in that order) upstream traffic-wise from location A. The objective was to predict at a particular time t+dt vehicular speed in location A. The discrete model is presented on the left hand side of Figure 10. The conditional probability table for node VA30 (vehicular speed in location A, 30 minutes into the future) requires more than 800000 conditional probabilities. Explanatory variables, speed (V) and density (D) in locations A, B and C, 30 minutes before the desired prediction time, are assumed to be independent. This is mainly due to quantification restrictions. Including other variables in the model would increase the quantification burden and will eventually make the model computationally infeasible. For this reason NPBNs where thought of as one possibility to overcome these difficulties. The continuous and extended version of the model is presented in the right hand side of Figure 10 (Worm et al. 2011). The model in Figure 10 currently includes 12 nodes and 26 arcs. A number of observations are in order. First of all, relaxing the assumption that variables at time t are independent of each other has little effect on the quantification procedure for the NPBN. Variables rain intensity (R) for locations A and B and traffic intensity (I) for locations A, B and C have been included. If the model in right hand side of Figure 10 were discretized similarly to that in the left hand side of the same figure, the conditional probability table for VA30 would have more than 1.3 × 1010 entries. One assumption often made in dynamic discrete BNs is that the model is homogeneous. That is, that the conditional probability tables do not change in time (Murphy 2002a). This would warranty, for example, that the marginal distribution of VA30 in Figure 10 does not change in time. This assumption needs not be present in the NPBN approach. The marginal distribution of VA30 is quantified separately from the dependence structure. Hence, at different times of interest, the marginal distribution of VA30 may stay unchanged while the (conditional) rank correlations with the variables 30 minutes before may be different. In this sense the model is non-homogeneous. The traffic prediction problem in addition to being a dynamic problem requires many explanatory variables. For example knowing what happens in a similar location but on a different highway at the same time might improve our prediction of traffic for a location of interest. Thus models with many more variables than those presented in Figure 10 are expected to be developed for sections of the Dutch road network. The flexibility of NPBNs will become an important asset for extensions of the model in Figure 10 to include multiple locations.

References Aanonsen, S., G. Naedval, D. Oliver, A. Reynolds, and B. Valles (2009). The Ensemble Kalman Filter in Reservoir Engineering. SPE Journal . Aguilera, P., A. Fernndez, F. Reche, and R. Rum (2010). Hybrid Bayesian network classifiers: Application to species distribution models. Environmental Modelling & Software 25 (12), 1630 – 1639. Ale, B., L. Bellamy, R. Cooke, M. Duyvis, D. Kurowicka, P. Lin, O. Morales, A. Roelen, and J. Spouge (2008, May). Causal model for air transport safety. Final Report ISBN 10: 90 369 1724-7, Ministerie van Verkeer en Waterstaat, Rotterdam. Ale, B., L. Bellamy, J. Cooper, D. Ababei, D. Kurowicka, O. Morales, and J. Spouge (2010). Analysis of the crash of tk 1951 using cats. Reliability Engineering and System Safety 95, 469477. Aparicio, A., P. Martnez-Austria, A. Gitrn, and A. Ramrez (2009). Floods in Tabasco, Mexico: a diagnosis and proposal for courses of action. Journal of Flood Risk Management. Bedford, T. and R. Cooke (2002). Vines - A New Graphical Model for Dependent Random Variables. Annals of Statistics 30 (4), 1031–1068. Blauw, S. (2014). Constructing non-parametric Bayesian networks from data. Bachelor Thesis, TU Delft, The Netherlands. Clemen, R., G. Fischer, and R. Winkler (2000). Assessing Dependencies: Some Experimental Results. Management Science 46 (8), 1100–1115. Cobb, B. (2009). Influence diagrams for capacity planning and pricing under uncertainty. Journal of Management Accounting Research 21, 75–97. Cobb, B. (2010). Detecting online credit card fraud with hybrid influence diagrams. Working Paper.

22

Cooke, R. (1991). Experts in Uncertainty : Opinion and Subjective Probability in Science. Environmental Ethics and Science Policy Series. Oxford University Press. Cover, T. and J. Thomas (1991). Elements of information theory. New York: John Wiley & Sons. Delgado-Hern´ andez, D.-J., O. Morales-N´apoles, D. De-Le´on-Escobedo, and J.-C. Arteaga-Arcos (2014). A continuous bayesian network for earth dams’ risk assessment: an application. Structure and Infrastructure Engineering 10 (2), 225–238. Embrechts, P., F. Lindskog, and A. McNeil (2003). Modelling dependence with copulas and applications to risk management. Handbook of heavy tailed distributions in finance 8, 329–384. Evensen, G., J. Hove, H. Meisingset, E. Reiso, K. Seim, and O. Espelid (2007). Using the EnKF for assisted history matching of a North Sea reservoir model. SPE Reservoir Simulation Symposium Huston(Texas), USA. Fernandez, A., I. Perez-Bernabe, R. Rumi, and A. Salmeron (2013). Incorporating Prior Knowledge when Learning Mixtures of Truncated Basis Functions from Data. Frontiers in Artificial Intelligence and Applications 257, 95–104. Gheorghe, M. (2010). Non parametric Bayesian belief nets versus Ensemble Kalman Filter in reservoir simulation. Technical report, MSc Thesis, Delft University of Technology. Greene, M., Y. Liu, W. Chen, and W. Liu (2010). Computational uncertainty analysis in multiresolution materials via stochastic constitutive theory. Comput. Methods Appl. Mech. Engrg. Haak, W. v. d., L. Rothkrantz, P. Wiggers, B. Heijligers, T. Bakri, and D. Vukovic (2010). Modeling Traffic Information using Bayesian Networks. Transaction on transport sciences 3 , 129–136. Hanea, A., D. Kurowicka, and R. Cooke (2006). Hybrid Method for Quantifying and Analyzing Bayesian Belief Nets. Quality and Reliability Engineering International 22 (6), 613–729. Hanea, A., D. Kurowicka, and R. Cooke (2007). The population version of Spearman’s rank correlation coefficient in the case of ordinal discrete random variables. In Proceedings of the Third Brazilian Conference on Statistical Modelling in Insurance and Finance. Hanea, A., D. Kurowicka, R. Cooke, and D. Ababei (2010). Mining and visualising ordinal data with non-parametric continuous BBNs. Computational Statistics and Data Analysis 54 (3), 668–687. Hanea, A. M., M. Gheorghe, R. Hanea, and D. Ababei (2013). Non-parametric Bayesian networks for parameter estimation in reservoir simulation: a graphical take on the ensemble Kalman filter (part I). Computational geosciences 17 (6), 929–949. Hanea, A. M. and W. Harrington (2009). Ordinal PM2.5 Data Mining with Non-Parametric Continuous Bayesian Belief Nets. Information Processes Journal . Hanea, D. (2009). Human Risk of Fire: Building a decision support tool using Bayesian networks. Ph. D. thesis, Technische Universiteit Delft. Hanea, D., H. Jagtman, and B. Ale (2012). Analysis of the Schiphol Cell Complex fire using a Bayesian belief net based model. Reliability Engineering & System Safety 100, 115–124. Huibregtse, E., O. Morales Napoles, L. Hellebrandt, D. Paprotny, and S. de Wit (2015). Climate change in asset management of infrastructure: A risk-based methodology applied to disruption of traffic on road networks due to the flooding of tunnels. Accepted for publication in European Journal of Transport and Infrastructure Research. Jager, W. (2013). Eliciting and Representing Joint Distributions from Experts: Quantification of a Human Performance Model for Risk Analysis. MSc. Thesis, TU Delft, The Netherlands. Jesionek, P. and R. Cooke (2007). Generalized method for modeling dose-response relations application to BENERIS project. Technical report, European Union project. Joe, H. (1997). Multivariate Models and Dependence Concepts. London: Chapman & Hall. Jordan, M., Z. Ghahramani, T. Jaakkola, and L. Saul (1999). An introduction to variational methods for graphical models. Machine Learning Journal 37, 183233. Kozlov, A. and D. Koller (1997). Nonuniform dynamic discretization in hybrid networks. Proceedings of the 13th conference on uncertainty in artificial intelligence. 23

Kurowicka, D. and R. Cooke (2004). Distribution - Free Continuous Bayesian Belief Nets. Proceedings Mathematical Methods in Reliability Conference. Langseth, H., T. Nielsen, R. Rum´ı, and A. Salmer´on (2009). Inference in hybrid Bayesian networks. Reliability Engineering & System Safety 51, 485–498. Langseth, H., T. Nielsen, R. Rum´ı, and A. Salmer´on (2010). Parameter estimation and model selection for mixtures of truncated exponentials. International Journal of Approximate Reasoning 94(10), 1499–1509. Langseth, H., T. Nielsen, R. Rumi, and A. Salmeron (2012). Mixtures of truncated basis functions. International Journal of Approximate Reasoning 53 (2), 212–227. Lauritzen, S. (1992). Propagation of probabilities, means and variances in mixed graphical association models. Journal of the American Statistical Association 87, 10981108. Lerner, U. (2002). Hybrid Bayesian networks for reasoning about complex systems. Ph. D. thesis, Stanford University. Morales, O., D. Kurowicka, and A. Roelen (2008). Eliciting conditional and unconditional rank correlations from conditional probabilities. Reliability Engineering & System Safety 93 (5), 699 – 710. Morales-N´ apoles, O. (2009). Bayesian belief nets in aviation safety and other applications. Ph. D. thesis, Delft University of Technology. Morales-N´ apoles, O., D. Delgado-Hern´andez, D. De-Le´on-Escobedo, and J. Arteaga-Arcos (2014). A continuous bayesian network for earth dams’ risk assessment: Methodology and quantification. Structure and Infrastructure Engineering 10 (5), 589–603. Morales-N´ apoles, O., A. Hanea, and D. Worm (2014). Experimental results about the assessments of conditional rank correlations by experts: Example with air pollution estimates. pp. 1359–1366. Morales N´ apoles, O. and R. Steenbergen (2012). Large-Scale Hybrid Bayesian Network for Traffic Load Modeling from Weigh-in-Motion System Data. accepted in the journal of bridge engineering of the ASCE.. Morales N´ apoles, O. and R. Steenbergen (2014). Analysis of axle and vehicle load properties through bayesian networks based on weigh-in-motion data. Engineering & System Safety 125, 153–164. Morales N´ apoles, O., D. Worm, and B. Dillingh (2011). Framework for probabilistic scale transition in physico-chemical modeling of asphalt. TNO report project number 034.24789 . Murphy, K. (2002a). Dynamic Bayesian Networks: Representation, Inference and Learning. Ph. D. thesis, UC Berkeley, Computer Science Division. Murphy, K. (2002b). An Introduction to Graphical Models. Technical report, Intel Research. Neil, M., M. Tailor, and M. D. (2007). Inference in Bayesian networks using dynamic discretisation. Statistics and Computing 17 (3), 219–33. Nelsen, R. (1999). An Introduction to Copulas. Lecture Notes in Statistics. New York: Springer Verlag. Ostoja-Starzewski, M. (2007). Microstructural Randomness and Scaling in Mechanics of Materials. Modern Mechanics and Mathematics Series.: Chapman & Hall/CRC/Taylor & Francis. Pearl, J. (1988). Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. San Mateo: Morgan Kaufman Publishers. Pearson, K. (1907). Mathematical contributions to the theory of evolution. Biometric Series. VI.Series. Philips, R. (1998). Multiscale modeling in the mechanics of materials. Current Opinion in Solid State and Materials Science 3(6), 526–532. Raad (2006). Brand Cellencomplex Schiphol-Oost: Eindrepport van het onderzoek naar de brand in het detentie-en uitzetcentrum Schiphol-Oost in de nacht van 26 tot 27 oktober 2005. Den Haag: Onderzoeksraad voor veiligheid. Shachter, R. and C. Kenley (1989). Gaussian influence diagrams. Management Science 35 (5), 527–550.

24

Shannon, C. and W. Weaver (1949). The mathematical theory of communication. Urbana, Illinois: University of Illinois Press. Shenoy, P. and J. West (2010). Inference in hybrid Bayesian networks using mixtures of polynomials. International Journal of Approximate Reasoning In Press. Sneller, M. (2013). Implementing quantitative expert judgement in learning the structure of a NPBN. Technical report, Delft University of Technology. Straub, D. and A. Der Kiureghian (2008). An investigation into the combination of Bayesian network with structural reliability methods. In Proc., IFIP WG 7.5 Working Conf. on Reliability and Optimization of Structures, Toluca, Mexico. Straub, D. and A. Der Kiureghian (2009). Bayesian networks as a framework for structural reliability in infrastructure systems. In Proc., ICOSSAR09, CRC Press, Boca Raton. Straub, D. and A. Der Kiureghian (2010). Bayesian Network Enhanced with Structural Reliability Methods: Methodology. Journal of Engineering Mechanics 136 (10), 1248–1258. Weber, P., G. Medina-Oliva, C. Simon, and B. Iung (2010). Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas. Engineering Applications of Artificial Intelligence. Weg-en Waterbouwkunde, D. (2004). Weigh in motion: Het aslastmeetsysteem van Rijkwaterstaat. http://www.rws.nl/rws/dww/home/html/menu5/pdf/DWWwijzer110.pdf. Accesesed on July 28, 2010. Whittaker, J. (1990). Graphical Models in applied multivariate statistics. Chichester: John Wiley and Sons. Worm, D., O. Morales N´ apoles, W. v. d. Haak, and T. Bakri (2011). Continuous Dynamic NonParametric Bayesian Networks: Application to traffic prediction. TNO-report project number 043.01000 . Yule, G. and M. Kendall (1965). An introduction to the theory of statistics. Belmont, California: Charles Griffin & Co. 14th edition. Zhang, K., J. Peters, D. Janzing, and B. Scholkopf (2012). Kernel-based conditional independence test and application in causal discovery. Max Planck Intitute for Intelligent Systems. Zhang, N. and D. Poole (1994). A simple approach to Bayesian network computations. In Proceedings of the Tenth Canadian Conference on Artificial Intelligence.

25

Table 1: Applications of the NPBN method. Application CATS BENERIS Fire Safety P M2.5 Earth Dams Permeability Weigh in Motion Material Modelling Traffic Prediction a These

Continuous Data Experts 618 300 282 0 7 10 17 0 4 6 343 0 705 0 4 0 12 0

Nodes Discrete Data Experts 30 2 242 0 0 0 0 0 0 0 0 0 8 0 4 0 0 0

Functional 586 637 7 0 1 0 8 0 0

Arcs Correlations Data Experts 3034a 26 8 0 0 11 45 0 0 16 89253 0 2344 0 19 0 26 0

Functions 1919 1804 20 0 2 0 528 0 0

correspond to repetitions of the data obtained about (conditional) correlations elicited from experts.

26

Figure 1: The NPBN for Earth Dams Safety.

Figure 2: Jos´e Antonio Alzate dam model.

27

Figure 3: The NPBN for P M2.5 Concentrations.

Figure 4: NPBN for Alabama ambient P M2.5 .

28

Figure 5: The NPBN for CATS.

Figure 6: The NPBN for BENERIS.

29

Figure 7: The NPBN for the Human Damage in Building Fire Model.

Figure 8: The NPBN for the Weigh in Motion System of The Netherlands.

Figure 9: Asphalt sample and fitted NPBN.

30

Figure 10: The discrete and the NPBN model for speed prediction at location A 30 minutes in the future.

31