Robust Learning with Missing Data.
This paper introduces a new method, called the robust Bayesian estimator (RBE), to learn conditional probability distributions from incomplete data sets. The intuition behind the RBE is that, when no information about the pattern of missing data is available, an incomplete database constrains the set of all possible estimates and this paper provides a characterization of these constraints. An experimental comparison with two popular methods to estimate conditional probability distributions from incomplete dataGibbs sampling and the EM algorithmshows a gain in robustness. An application of the RBE to quantify a naive Bayesian classifier from an incomplete data set illustrates its practical relevance.
Introduction
Probabilistic methods play a central role in the development of Artificial
Intelligence (AI) applications because they can deal with for the degrees
of uncertainty so often embedded into real-world problems within a sound
mathematical framework. Unfortunately, the number of constraints needed
to define a probabilistic model, the so called joint probability distribution,
grows exponentially with the number of variables in the domain, thus making
infeasible the straight use of these methods.
However, the assumption that some of the variables are stochastically
independent, given a subset of the remaining variables in the domain, dramatically
reduces the number of constraints needed to specify the probabilistic
model. In order to exploit this property, researchers developed a
new formalism to capture, in graphical terms, the assumptions of independence
in a domain, thus reducing the number of probabilities to be assessed.
This formalism is known as Bayesian Belief Networks (bbns) [Pearl, 1988].
bbns provide a compact representation for encoding probabilistic information
and they may be easily extended into a powerful decision-theoretic
formalism called Influence Diagrams. More technically, a bbn is a direct
acyclic graph in which nodes represent stochastic variables and links represent
conditional dependencies among variables. A variable bears a set of
possible values and may be regarded as a set of mutually exclusive and exhaustive
states, each representing the assignment of a particular value to the
variable. A conditional dependency links a child variable to a set of parent
variables and it is defined by the set of conditional probabilities of each state
of the child variable given each combination of states of the parent variables
in the dependency.
In the original development of bbns, domain experts were supposed to
be their main source of information: their independence assumptions, when
coupled with their subjective assessment of the conditional dependencies
among the variables, produces a sound and compact probabilistic representation
of their domain knowledge. However, bbns have a strong statistical
root and, during the past few years, this root prompted for the development
of methods able to learn bbns directly from databases of cases rather than
from the insight of human domain experts [Cooper and Herskovitz, 1992,
Buntine, 1994, Heckerman and Chickering, 1995]. This choice can be extremely
rewarding when the domain of applications generates large amounts
of statistical information and aspects of the domain knowledge are still unknown
or controversial, or too complex to be encoded as subjective probabilities
2of few domain experts.
Given the nature of bbns, there are two main tasks involved in the learning
process of a bbn from a database: the induction of the graphical model
of conditional independence best fitting the database at hand, and the extraction
of the conditional probabilities defining the dependencies in a given
graphical model from the database. Under certain constraints, this second
task can be accomplished exactly and efficiently when the database is complete
[Cooper and Herskovitz, 1992, Heckerman and Chickering, 1995]. We
can call this assumption the "Database Completeness" assumption. Un-
fortunately, databases are rarely complete: unreported, lost, and corrupted
data are a distinguished feature of real-world databases. In order to move
on applications, methods to learn bbns have to face the challenge of learning
from databases with unreported data. During the past few years, several
methods have been proposed to learn conditional probabilities in bbns
from incomplete databases, either deterministic, such as sequential updating
[Spiegelhalter and Lauritzen, 1990], and [Cowel et al., 1996], or the EM-algorithm
[Dempster et al., 1977], or stochastic such as the Gibbs Sampling
[Neal, 1993]. All these methods share the common assumption that unreported
data are missing at random. Unfortunately this assumption is as
unrealistic as the "Database Completeness" assumption because in the real
world there is often a reason so that data are missing.
This paper introduces a new method to learn conditional probabilities
in bbns from incomplete databases which does not rely on the "Missing at
Random" assumption. The major feature of our method is the ability to
learn bbns which are robust with respect to the pattern of missing data
and to return sets of estimates,instead of point estimates, whose width is a
monotonic increasing function of the information available in the database.
The reminder of this paper is structured as follows. In the next Section we
will start by establishing some notation and reviewing the background and
the motivation of this research. In Section 3 we will outline the theoretical
framework of our robust approach to learn probabilities in bbns and in Section
4 we will describe in details the algorithms we used to implement such
an approach. In Section 5 we will compare the behavior of the system implemented
using our method with an implementation of the Gibbs Sampling
and finally we will draw our conclusions.

Figure

1: The graphical structure of a bbn.
Background
A bbn is defined by a set of variables and a network
structure S defining a graph of conditional dependencies among the elements
of X . We will limit attention to variables which take a finite number
of values and therefore, for each variable X , we can define a set of
mutually exclusive and exhaustive states x representing
the assignments to X i of all its m possible values.

Figure

1 shows a bbn we will use to illustrate our notation. In this case
we have and we assume that each variable may take
two values, true (1) and false (0).
The network structure S defines a set of conditional dependencies among
the variables in X . In the case reported in Figure 1, there are two depen-
dencies: one linking the variable X 3 to its parents X 1 and X 2 , which are
marginally independent, and the other one linking the variable X 4 to its
parent X 3 . Since there are no direct links between
separates X 4 from the nodes we have that X 4 is conditionally independent
of . The conditional dependencies are defined by a
set of twelve conditional probabilities: eight are p(X
and four for p(X
2. These conditional probabilities and the marginal
probabilities of the parent nodes X 1 and X 2 , that is p(X
are then used to write down the joint probability of each
network configuration, for instance fX

Figure

2: The graphical structure of a bbn with the associated parameters.
Using a generic structure S, the joint probability of a particular set of values
of the variables in X , say can be decomposed as
Y
are the parent nodes of X i , and x pa(X i ) denotes the combination
of states of pa(X i ) in x. If X i is a root node, then p(X
In the following we will denote the events X
We are given a database of cases g, each case C
being a set of entries. We will assume that the cases are
mutually independent. Given a network structure S, the task here is to learn
the conditional probabilities defining the dependencies in the bbn from D.
We will consider the conditional probabilities defining the bbn as parameters
so that the joint probability of each case in the database,
say
Y
Table

1: Parameters defining the conditional probabilities of the bbn displayed
in Figure 1.
and ' i parameterizes the probability of x ij given the parent configuration
pa(x i ). For instance, in the network in Figure 1 we need only eight parameters
to describe all the dependencies, two of them are needed to define
the probability distributions of X 1 and X 2 , say
and the other six to define the probability distributions of
is given in Figure 2.
Our task is to assess the set of parameters induced by the database D
over the network structure S. The classical statistical parameter estimation
provides the basis to learn these parameters. When the database is complete,
the common approach is Maximum Likelihood which returns the parameter
values that make the database most likely. Given the values of X in the
database, the joint probability of D is
l
Y
This is a function of the parameters ' only, usually called the likelihood
l('). The Maximum Likelihood estimate of ' is the value which
maximizes l('). For discrete variables the Maximum Likelihood estimates
of the conditional probabilities are the observed relative frequencies of the
relevant cases in the database. Let n(x ij jpa(x i )) be the observed frequency
of cases in the database with x ij , given the parent configuration pa(x i ), and
let n(pa(x i )) be the observed frequency of cases with pa(x i ). The Maximum
Likelihood estimate of the conditional probability of x ij simply
Consider for instance the database given in Table 2. The Maximum
Case

Table

2: An artificial database for the bbn displayed in Figure 1.
Likelihood estimate of ' 3 is 1/3.
The Bayesian approach extend the classical parameter estimation technique
in two ways: (i) the set of parameters ' are regarded as random
variables, and (ii) the likelihood is augmented with a prior, say -('), representing
the observer's belief about the parameters before observing any
data. Given the information in the database, the prior density is updated
in the posterior density using Bayes' theorem, and hence
Z
-(')p(Dj')d':
The Bayesian estimate of ' is then the expectation of ' under the posterior
distribution.
The common assumptions in the Bayesian approach to learn a bbn with
discrete variables are that the parameters (i) are independent from each
other, and (ii) have a Dirichlet distribution, which simplifies to a Beta distribution
for binary variables. Assumption (i) allows to factorise the joint
prior density of ' as
Y
thus allowing "local computations", while (ii) facilitates computations of
the posterior density by taking advantages of conjugate analysis. Full details
are given for instance by [Spiegelhalter and Lauritzen, 1990]. Here we
just outline the standard conjugate analysis with a Dirichlet prior.Consider
for instance the variable X i and let ' be the parameters associated
to the conditional probabilities so that
1. A Dirichlet prior for ' i , denoted as D(ff
a continuous multivariate distribution with density function proportional to
Y
The hyper-parameters ff ij s have the following interpretation: ff
can be regarded as an imaginary sample size needed to formulate this prior
information about ' i , and the mean of ' ij is ff ij =ff i+ , m. Note that
this prior mean is also the probability of x ij given the parent configuration
Hence, a uniform prior with would assign uniform
probabilities to each x ij , given the parent configuration pa(x i ).
With complete data, the posterior distribution of the parameters can
be computed exactly using standard conjugate analysis, yielding a posterior
density and the posterior
means represent their estimates. Thus the Bayes estimate of the conditional
probability of x ij ))g.
Consider again the example in Table 2. If the prior distribution of ' 3 is
so that a priori the conditional probability of X
is 0.5, the posterior distribution would be updated into D(2; 3), yielding a
Bayesian estimate 2/5.
Unfortunately, the situation is quite different when some of the entries
in the database are missing. When a datum is missing, we have to face a
set of possible complete databases, one for each possible value of the variable
for which the datum is missing. Exact analysis would require to compute
the joint posterior distribution of the parameters given each possible
completion of the database, and then to mix these over all possible com-
pletions. This is apparently infeasible. A deterministic method, proposed
by [Spiegelhalter and Lauritzen, 1990] and improved by [Cowel et al., 1996],
provides a way to approximate the exact posterior distribution by processing
data sequentially. However, [Spiegelhalter and Cowel, 1992] show that
this method is not robust enough to cope with systematically missing data:
in this case the estimates rely heavily on the prior distribution.
When deterministic methods fail, current practice has to resort to an
estimate of the posterior means using stochastic methods. Here we will
describe only the most popular stochastic method for Bayesian inference:
the Gibbs Sampling. The basic idea of the Gibbs Sampling algorithm is (i)
for each parameter ' i sample a value from the conditional distribution of ' i
given all the other parameters and the data in the database, (ii) repeat this
for all the parameters, and (iii) iterates these steps several times. It can be
proved that, under broad conditions, this algorithm provides a sample from
the joint posterior distribution of ' given the information in the database.
This sample can then be used to compute empirical estimates of the posterior
means or any other function of the parameters. In practical applications,
the algorithm iterates a number of times and then, when stability seems
to be reached, a final sample from the joint posterior distribution of the
parameters is taken [Buntine, 1996].
When some of the entries in the database are missing, the Gibbs Sampling
treats the missing data as unknown parameters, so that, for each
missing entry, a values is sampled from the conditional distribution of the
corresponding variables, given all the parameters and the available data.
The algorithm is then iterated to reach stability, and then a sample from
the joint posterior distribution is taken which can be used to provide empirical
estimates of the posterior means [Thomas et al., 1992]. This approach
relies on the assumption that the unreported data are missing at random.
When the "Missing at random" assumption is violated, as when data are
systematically missing, such method suffers of a dramatic decrease in ac-
curacy. The completion of the database using the available information in
the database itself leads the learning system to ascribe the missing data
to known values in the database and, in the case of systematically missing
data, to twist the estimates of the probabilities in the database. It is apparent
that this behavior can prevent the applicability of the learning method
to generate bbns because, for the general case, it can produce unreliable
estimates.
Theory
The solution we propose is a robust method to learn parameters in a bbn.
Our method computes the set of possible posterior distributions consistent
with the available information in the database and proceeds by refining this
set as more information becomes available. The set contains all possible
estimates of the parameters in the bbn and can be efficiently computed by
calculating the extreme posterior distributions consistent with the database
D. In this Section, we will first define some basic concepts about convex
sets of probabilities and probability intervals, and then we will describe the
theoretical basis of our learning method.
3.1 Probability Intervals
Traditionally, a probability function p(x) assigns a single real number to the
probability of x. This requirement, called Credal Uniqueness Assumption,
is one of the most controversial points of Bayesian theory. The difficulty
of assessing precise, real-valued quantitative probability measures is a long
standing challenge for probability and decision theory and it motivated the
development of alternative formalisms able to encode probability distributions
as intervals rather than real numbers. The original interest for the notion
of belief functions proposed by Dempster [Dempster, 1967] and further
developed by Shafer [Shafer, 1976] was mainly due to its ability to represent
credal states as intervals rather than point-valued probability measures
[Grosof, 1986].
Several efforts have been addressed to interpret probability intervals
within a coherent Bayesian framework, thus preserving the probabilistic
soundness and the normative character of traditional probability and decision
theory [Levi, 1980, Kyburg, 1983, Stiling and Morrel, 1991]. This approach
regards probability intervals as a convex sets of standard probability
distributions and is therefore referred to as Convex Bayesianism. Convex
Bayesianism subsumes standard Bayesian probability and decision theory
as a special case in which all credal states are evaluated by convex sets of
probability functions containing just one element.
Convex Bayesianism relaxes the Credal Uniqueness Assumption by replacing
the real-valued function p(x) with a convex set of these functions.
A convex set of probability functions P (x) is the set
are two probability functions and 0 - ff - 1 is a real
number. We will use p ffl and p ffl to denote the minimum and the maximum
probability of the set P , respectively.
The intuition behind Convex Bayesianism is that even if an agent cannot
choose a real-valued probability p(x), his information can be enough to constrain
p(x) within a set P (x) of possible values. An appealing feature of this
approach is that each probability function p ff fulfills the requirements
of a standard probability function, thus preserving the probabilistic
soundness and the normative character of traditional Bayesian theory.
The most conservative interpretation of Convex Bayesianism is the so
called Sensitivity Analysis or Robust Bayesian approach [Berger, 1984], within
which the set P (x) is regarded as a set of precise probability functions and
inference is carried out over the possible combinations of these probability
functions. Good [Good, 1962] proposes to model this process as a "black
box" which translates these probability functions into a set of constraints of
an ideal, point-valued probability function. A review of these methods may
be found in [Walley, 1991]. This theoretical framework is especially appealing
when we are faced with the problem of assessing probabilities from an
incomplete database. Indeed, in this case we can safely assume the existence
of such an ideal, point valued probability as the probability value we would
assess if the database was complete.
The suitability of Convex Bayesianism to represent incomplete probabilistic
information motivated the development of computational methods
to reason on the basis of probability intervals [White, 1986, Snow, 1991]. A
natural evolution of this approach leaded to the combination of the computational
advantages provided by conditional independence assumptions with
the expressive power of probability intervals [van der Gaag, 1991]. There-
fore, some efforts have been addressed to extend bbns from real-valued probabilities
to interval probabilities [Breeze and Fertig, 1990, Ramoni, 1995],
thus combining the advantages of conditional independence assumptions
with the explicit representation of ignorance. These efforts provide different
methods to use the bbns we can learn with our method.
3.2 Learning
We will describe the method we propose by using the artificial database
in

Table

2 for the bbn given in Figure 1. Table 3 reports the database
given in Table 2 where some of the entries, denoted by ?, are missing. We
will proceed by analyzing the database sequentially, just to show how we
can derive our results. We stress here that our method does not require a
sequential updating.
We start by assuming total ignorance, and hence the parameters ' i ,
are all given a prior distribution D(1; 1). We also assume that
the are independent. Consider the first entry in the database, this is a
complete case so that we update the distributions of ' 1 , ' 2 , ' 3 and ' 8 in
1). The
distributions of ' 3 are not changed since the entries in C 1 do
not have any information about them. After processing the first case we
thus have that the Bayes estimates of the parameters are:
Case

Table

3: The incomplete database displayed in Table 2.
Consider now case C 2 . The observation on X 3 is missing, and we know that
it can be either 1 or 0. Thus the possible completions of case C 2 are
The two possible completions would yield ' 1
would not be updated, whichever is the com-
pletion, since only the parent configuration 0; 0 is observed. Consider now
the updating of the distributions of ' 3 , ' 7 and ' 8 . We have that ' 3
D(1; 2) and this completion would lead to ' 7
1). Instead C 2 c2
would lead to ' 3
1). Hence the possible Bayes estimates are:
Instead of summarizing this information somehow, we can represent it via
intervals, whose extreme points are the minimum and the maximum Bayes
estimates that we would have from the possible completions of the database.
Thus, for instance, from the two possible posterior distributions of ' 3
we learn
and similarly from the posterior distributions of ' 7 and ' 8 we
learn
3=4. Clearly, we also have
Consider now the third case in the database. Again the observation
on X 3 is missing, and as before we can consider the two possible completions
of this case, and proceed by updating the relevant distributions. Now
however we need to "update intervals", and this can be done by updating
the distributions corresponding to the extreme points of each interval.
This would yield four distributions, and hence 4 Bayes estimates of the relevant
conditional probabilities from which we can extract the extreme ones.
For instance, the completion C 3 c1
would lead to the updating,
among others, of the distribution of ' 7 , so that the two distributions obtained
after processing the first two cases in the database would become
- D(2; 1), and ' 8 is not up-
dated. If we consider the completion C 3 c2
7 is not updated,
1). If we sort
the corresponding estimates in increasing order and find the minimum and
the maximum we obtain p ffl (X
Consider 3=4. The
minimum probability is achieved by assuming that all the completions in
the database assign 1 to X 3 ; the maximum is achieved by assuming that all
the completions in the database assign 0 to X 3 . If we now process the fourth
case, for instance the distribution of ' 7 would be updated by considering the
completions 0; leading to extreme probabilities
4=5. Thus the minimum is
obtained by assigning 0 to X 4 and the maximum by assigning 1 to X 4 .
This can be generalized as follows. Let X i be a binary variable in the
bbn and denote by n ffl (1jpa(x i )) the frequency of cases with
the parent configuration pa(x i ), which have been obtained by completing incomplete
cases. Similarly, let n ffl (0jpa(x i )) denote the frequency of cases with
given the parent configuration pa(x i ), which have been obtained by
completion of incomplete cases. Suppose further that we start from total ig-
3norance, thus the parameter ' i which is associated to p(1jpa(x i )) is assigned
a D(1; 1) prior, before processing the information in the database. Then
and
The minimum and maximum probability of the complementary event are
such that p ffl (1jpa(x i
This result can be easily generalized to discrete variables with k states
leading to
and
Note that in this case the sum of the maximum probability of x ij jpa(x i
and the minimum probabilities of x ih jpa(x i ), h 6= j, is one.
It is worth noting that the bounds depend only on the frequencies of
complete entries in the database, and the "artificial" frequencies of the completed
entries, so that they can be computed in batch mode.
This Section is devoted to the description of the algorithms implementing
the method outlined in Section 3. We will first outline the overall procedure
to extract the conditional probabilities from a database given a network
structure S. Then, we will describe the procedure to store the observations
in the database, and we will analyze the computational complexity of the
algorithms. Finally, we will provide details about the current implementation

4.1

Overview

Section 3 described the process of learning conditional probabilities as a
sequential updating of a prior probability distribution. The nature of our
method allows us to implement the method as a batch procedure which
first parses the database and stores the observations about the variables
and then, as a final step, computes the conditional probabilities needed to
specify a bbn from these observations.
The procedure takes as input a database D, as defined in Section 2,
and a network structure S. The network structure S is identified by a set
of conditional dependencies fd X 1
associated to each variable in
. A dependency dX i
is an ordered tuple (X is a child
variable in X with parent nodes pa(X i ).
The learning procedure takes each case in the database as a statistical
unit and parses it using the dependencies defining the network structure S.
Therefore, for each entry in the case, the procedure recalls the dependency
within which the entry variable appears as a child and identifies the states
of its parent variables in the case. In this way, for each case in the database,
the procedure detects the configuration of states for each dependency in S.
Recall that the probabilities of the states of the child variables given the
states of the parent variables are the parameters ' we want to learn from
the database.
For each configuration the procedure maintains two coun-
ters, say n(x When the detected configuration
does not contain any missing datum, the first counter n(x ij jpa(x i
is increased by one. When the datum for one or more variables in the
combination is missing, the procedure increases by one the second counter
for each configuration of states of the variable of each missing
entry. In other words, the procedure uses the counter n ffl
to ascribe a sort of virtual observation to each possible state of a variable
whose value is missing in the case. Once the database has been parsed, we
just need to collect the counters of each configuration and compute the two
extreme lower bound using formula 3, and the upper bound using formula
4.
Let D be a database and X the set of variables in a bbn. We will denote
as states(X i ) the set of states for the variable X
the set of its immediate predecessors in the bbn. The learning procedure is
defined as follows.
procedure learn(D;X )
while do
while do
while k - jP j do
while l - jX j do
if
The procedure store stores the counters for each configuration of parent
states and it will be described in the next subsection. The procedure collect
simply collects all the counters for each state in each variable in X .
4.2 Storing
It is apparent that the procedure store plays a crucial role for the efficiency
of the procedure. In order to develop an efficient algorithm, we used discrimination
trees to store the parameter counters, following a slightly modified
version of the approach proposed by [Ramoni et al., 1995].
Along this approach, each state of each variable of the network is assigned
to a discrimination tree. Each level of the discrimination tree is defined by
the possible states of a parent variable. Each path in the discrimination
tree represents a possible configuration of parent variables for that state. In
this way, each path is associated to a single parameter in the network. Each
leaf of the discrimination tree holds the pair of counters
Figure

3: The discrimination tree associated to the state X 3 = 1.
For each entry, if there is no missing datum, we just need to
follow a path in the discrimination tree to identify the counters to update.
In order to save memory storage, the discrimination trees are incrementally
built: each branch in the tree is created the first time the procedure needs
to walk through it.
The procedure store uses ordered tuples to identify the variable states
of interest. A state is identified by an ordered tuple (X; x), where X is the
variable and x the reported value. If the value is not reported, x will be
?. The procedure store takes as input three arguments: an ordered tuple
c representing the current state of the child variable, the set A of ordered
tuples for each state of the parent variables of c in the current case, and a
flag dictating whether the update is induced by a missing datum or not.
procedure store(c; A; v)
if c(2) =? then
while do
else
if p(2) =? then
while do
else
counter ffl (p)/counter ffl (p)
else
return
else
The functions counter and counter ffl identify the counters for n and for
respectively, associated to each leaf node in the tree.
In order to illustrate how these procedures work, let's turn back to the
example described in Section 2. The procedure learn starts by parsing the
first line of the database reported in Table 3. It uses the network structure
S depicted in Figure 1 to partition the case into four relevant elements:
corresponding to a parameters in the network structure.
Now, suppose the third entry has to be stored in the discrimination tree
associated to the state X 3 = 1. Figure 3 displays this discrimination tree
associated to the state X 3 = 1. The procedure walks along the solid line in

Figure

3 and updates the counter n(X because the
entry does not include any missing datum.
Suppose now the entry
has to be stored. Figure 3 identifies with a dashed line the paths followed by
the procedure in this case: when the procedure hits a state of the variable
whose datum is not reported, it walks through all the possible branches
following it and updates the n ffl counters at the end of each path.
4.3 Computational Complexity
The learning procedure takes advantage of the modular nature of bbns: it
partions the search space using the dependencies in the network. This partitioning
is the main task performed by the procedure learn. learn starts
by scanning the database D and, for each elements of a row, it scans the
row again to identify the patterns of the dependencies in the bbn. Suppose
the database D contains n columns and m rows, the superior bound
of the execution time of this part of the algorithm is O(gmn 2 ) , where g is
the maximum number of parents for a variable in the bbn. Note that the
number of columns is equal to the number of variables in the bbn, that is
The procedure collect scans the generated discrimination trees and is
applied once during the procedure. The number of discrimination trees generated
during the learning process is
discrimination
tree has a number of leaves
are
the parents of the variable of the state associated to the discrimination tree.
Note that this is the number of conditional distributions we need to learn
to define a conditional dependency.
The main job of the algorithm is left to the procedure store Using discrimination
trees, the time required by the procedure store to ascribe one
entry to the appropriate counter is linear in the number n of parent variables
in the dependency, when the reported entry about the parameter does not
contain any missing datum. When data are missing, the procedure is, in the
worse case, exponential in the number of parent variables with unreported
data in the row. It is worth noting that the main source of complexity in
the algorithm does not depend on the dimension of the database but on the
topology of the bbn.
Our previous example makes clear that the order of variables in the tree
plays a crucial role in the performance of the algorithm: if the positions of
the states of variable X 1 and X 2 were exchanged, the procedure had to walk
through the whole tree rather than just half of it. Efficiency may be gained
by a careful sorting of the variable states in the tree using precompilation
techniques to estimate in advance those variables whose values are missing
more often in the database. This task can be accomplished using common
techniques in the machine learning community to estimate the information
structure of classification trees [Quinlan, 1984].
4.4 Implementation
This method has been implemented in Common Lisp on a Machintosh Performa
6300 under Machintosh Common Lisp. A porting to CLISP running
on a Sun Sparc Station is under development. The system has been implemented
as a module of a general environment for probabilistic inference
called Refinment Architecture (era) [Ramoni, 1995]. era has
been originally developed on a Sun Sparc 10 using the Lucid Common Lisp
development environment. We were very careful to use only standard Common
Lisp resources in order to develop code conforming to the new established
ANSI standard. Therefore, the code should be easily portable on any
Common Lisp development environment.
The implementation of the learning system deeply exploits the modularity
allowed by the Common Lisp Object System (CLOS) protocol included
in the ANSI standard: for instance, the learning algorithm uses the CLOS
classes implemented in the architecture to represent the elements of the net-work
structure. This strategy allows a straightforward integration of the
learning module within the reasoning modules of era. In this way, the
results of the learning process are immediately available to the reasoning
modules of era to draw inferences and make decisions.
5 Experimental Evaluation
Gibbs Sampling is currently the most popular stochastic method for Bayesian
inference in complex problems, such as learning when some of the data are
missing, although its limitations are well-known: the convergence rate is
slow and resource consuming. However, given its popularity, we have compared
the accuracy of our method with one of its implementations. In this
Section, we will report the results of two sets of experimental comparisons,
one using a real-world problem and one using an artificial example. The
aim of these experiments is to compare the accuracy of the parameter estimates
provided by the Gibbs Sampling and our method as the available
information in the database decreases.
Figure

4: The network structure of the bbn used for the first set of experiments

5.1 Materials
In order to experimentally compare our method to the Gibbs Sampling,
we choose the program BUGS [Thomas et al., 1992] which is commonly regarded
as a reliable implementation of such a technique [Buntine, 1996]. In
the following experiments, we used the implementation of BUGS version 0.5
running on a Sun Sparc 5 under SunOS 5.5 and the era implementation of
our method running on a Machintosh PowerBook 5300 under Machintosh
Common Lisp version 3.9. All the experiments reported in this Section share
these materials. Different materials used for each set of experiments will be
illustrated during the description of each experiment.
5.2 Experiment 1: The CHILD Network
In the first set of experiments, we used a well-known medical problem. This
problem has been already used in the bbns literature [Spiegelhalter et al., 1993]
and concerns the early diagnosis of congenital hearth disease in newborn babies

# Name States
1 Birth Asphyxia yes no
Disease PFC TGA Fallot PAIVS TAPVD Lung
3 Age 0-3-days 4-10-days 11-30-days
4 LVH yes no
5 Duct flow Lt-to-Rt None Rt-to-Lt
6 Cardiac mixing None Mild Complete Transp
7 Lung parenchema Normal Congested Abnormal
8 Lung flow Normal Low High
9 Sick yes no
11 Hypoxia in O2 Mild Moderate Severe
Chest X-ray Normal Oligaemic Plethoric Grd-Glass Asy/Patchy
14 Grunting yes no
19 X-ray Report Normal Oligaemic Plethoric Grd-Glass Asy/Patchy
Grunting Report yes no

Table

4: Definition of the variables in the CHILD bbn: the first column
reports the numeric index used in Figure 4, the second the variable name,
and the third its possible states.
5.2.1 Materials

Figure

4 displays the network structure of our medical problem. Table 4
associates each number in the bbn to the name of a variable and reports its
possible values. The clinical problem underlying the bbn depicted in Figure
4 is described in [Frankin et al., 1989]. The task of the bbn is to diagnose
the occurrence of congenital heart disease in pediatric patients using clinical
data reported over the phone by referring pediatricians. Authors report that
the referral process generates a large amount of clinical information but this
information is often incomplete and includes missing and unreported data.
These features, together with the reasonable but realistic size of the
bbn, makes of this problem an ideal testbed for our experiments: the bbn in

Figure

4 is defined by 344 conditional probabilities (although the minimal set
Figure

5: Plots of estimates against amount of information for
the parameters a)
is slightly more than 240 since some of the variables in the bbn are binary)
and the number of cases reported in the original database was more than 100.
Since the aim of our experiment is to test the accuracy of the competing
methods, we cannot use the original database because we need a reliable
measure of the probabilities we want the systems to assess. Therefore, still
using the original network, we generated a complete random sample of 100
cases from a known distribution and we used this sample as the database
to learn from. Using this database and the bbn in Figure 4, we run two
different tests.
Test 1: Missing at Random
The goal of the first test is to characterize the behavior of our method with
respect to the Gibbs Sampling when the "Missing at Random" assumption
holds.
Method. We start with a complete database, where all the parameters are
independent and uniformly distributed, and run both our learning algorithm
and the Gibbs Sampling on it. Then, we proceed by randomly deleting
the 20% of entries in the database, and by running the two methods on
the incomplete database, until the database is empty. For each incomplete
database we run 10,000 iterations of the Gibbs Sampling, which appeared
to be enough to reach stability, and the estimates returned are based on a
final sample of 5,000 values.
Results. Figure 5 shows the estimates of the conditional probabilities
defining the dependency linking the variable n2 (disease) to the variable
inferred by the two learning methods, for 5 different proportions
of completeness. We report this dependency as an example of the overall
behavior of the systems during the test. Stars report the point estimates
given by the Gibbs Sampling, while errorbars indicates 95% confidence interval
about the estimates. Solid lines indicates the lower and upper bounds
of the probability intervals inferred by our method. We choose to represent
the outcomes of the Gibbs Sampling and our system in two different ways
because there is a basic semantic clash between the Gibbs Sampler confidence
intervals and the intervals returned by our system. The estimates
given by the Gibbs Sampling, when the database is incomplete, are based
on the most likely reconstruction of the missing entries, and this relies on
the prior belief about the parameters and the complete data. The 95% confidence
intervals represent the posterior uncertainty about the parameters,
which depend on the inferred database and the prior uncertainty. Thus the
larger the intervals, the less reliable are the estimates. Intervals returned by
our method represent the set of posterior estimates of the parameters that
we would obtain in considering all possible completions of the database. In
particular, both the estimates based on the original database and the estimates
returned by the Gibbs Sampling will be one of them. The width of
our intervals is then a measure of the uncertainty in considering all possible
completions of the database.
When the information in the database is complete, the intervals of our
systems degenerate to a single point, and this point coincides with the exact
estimates. The semantic difference between the intervals returned by the two
systems accounts for some of the slightly tighter intervals returned by the
Gibbs Sampling, together with the fact that, by assuming that the data are
missing at random, the Gibbs Sampling can exploit a piece of information
Figure

Plots of estimates against amount of information for the parameters
a)
not available to our system. Nonetheless, the width of the intervals returned
by the two systems is overall comparable in the case reported in Figure 5,
as well as for all the remaining parameters in the bbn.
The main difference was in the execution time: in the worse case, Gibbs
Sampling took over 37 minutes to run to completion on a Sun Sparc 5, while
our system ran to completion in less than 0.20 seconds on a Machintosh
5.2.3 Test 2: Systematically Missing
What does it happen when the data are not missing at random and therefore
the Gibbs Sampling is using a wrong guess? The aim of the second test is
to compare the behavior of the two systems when the data are not missing
at random, but the value of a variable is systematically removed from the
database.
Method. The procedure used for this test is a slight modified version of
the test for the randomly missing data. In this case, we iteratively deleted
Figure

7: Plots of estimates against amount of information for the parameters
a)
1% of the database by systematically removing the entries reporting the
value normal for the variable n7 (lung parenchema) and we ran the two
learning algorithms. Each run of the Gibbs Sampling is based on 10,000
iterations, and a final sample of 5,000 cases. This procedure was iterated
until no value normal was reported for the variable n7 in the database.
Results. The local independencies induced by the network structure are
such that the modification of values of the variable n7 in the database only
affect the estimates of the conditional probabilities for its immediate predecessors
and successors. Therefore, we will focus on the estimates of the
parameters affected by the changes. Figure 6 shows the estimates provided
by the two systems for the conditional probabilities of the states of variable
n7 given the state of its immediate predecessor in the
bbn. Note that the information percentage reported on the x\Gammaaxis starts
at 98%, meaning that the number of entries for the state n7 = normal account
for the 2% the complete database. At 98%, the database contains no
entry for normal. Results report that the point-valued estimates of
the Gibbs Sampling always fall within the intervals calculated by our sys-
tem. However, plots a and b display the behavior of the Gibbs Sampling
estimates, jumping from the lower bound to the upper bound of the set
calculated by our method as the missing data are replaces with entries n7
normal. The jump testifies a 30% error in the point estimate provided
by the Gibbs Sampling when all the entries normal are missing. The
behavior of both systems does not change when the data are missing on the
parent variables rather than on the child variable. Figure 7 plots the estimates
of some parameters of the dependency linking the parent variables n7
(lung-parench) and n6 (Cardiac Mixing) to the child variable n11 (Hypoxia
in O2). In this case, the initial error of the Gibbs Sampling goes up to the
45%, as shown in the plot b. This error is even more remarkable once we
realize that the 98% of the overall information is still available.
The difference in execution time was comparable to the one of the previous
test: in the worse case, Gibbs Sampling took over 16 minutes to run
to completion on a Sun Sparc 5, while our system ran to completion in less
than 0.20 seconds on a Machintosh PowerBook 5300.
5.3 Experiment 2: An Artificial Network
Results from experiments on the CHILD network show a bias of the point
estimates given by the Gibbs Sampling, although the associated confidence
intervals are always large enough to include the estimates calculated from
the complete database. Since the width of the confidence intervals is a
function of the sample size, we are left with the doubt that a large database
could give tighter intervals around biased estimates. These results prompt
for a more accurate investigation of this behavior as the size of the database
increases. Thus, in this Section, we will use an artificial - and therefore
more controllable - example in order to amplify this bias. The rationale
behind this second experiment is twofold: (i) to display the effect on the
learning process of the bias induced by the "Missing at Random" assumption
when in fact data are systematically missing and the size of the database is
large, and (ii) to show the effect of this bias on the predictive performance
on the bbn.

Figure

8: The simple network used for the second set of experiments.
5.3.1 Materials

Figure

8 shows the simple bbn we used for our second experiment: two
binary variables linked by a dependency. We generated a database of 1000
random cases from the following probability distribution:
The parameters in the bbn were all assumed to be independent and uniformly
distributed.
5.3.2 Learning
Using the bbn displayed in Figure 8, we tried to make clearer the bias effect
detected in the previous set of experiments.
Method. We followed a procedure analogue to that used in the second
test of the previous experiment: we iteratively deleted the 10% of entries
with no value with in the database. We then
run our method and the Gibbs Sampling on each incomplete database. Each
run of the Gibbs Sampling is based on 1,000 iterations to reach stability, and
a final sample of 2,000 values.
Results. Figure 9 shows the parameter estimates given by the two sys-
tems. The bias of the Gibbs Sampling is absolutely clear: when 75% of
the information is available in the database but all the entries of
are missing, the estimate given by the Gibbs Sampling for lies
on the lower extreme of the interval estimated by our method. However,

Figure

9: Plots of estimates agaist amount of information for the parameters
a)
the sample size tight up the confidence interval around the estimate 0.0034
so the "true" estimate 0.4955, which would be computed in the complete
database, is definitely excluded, with an error overpassing the 40%.
The estimate of p(X moves from 0.8919 in the complete
database to 0.662 when no entries with are left in the database. How-
ever, the wide interval associated to the estimate testifies the low confidence
of the Gibbs Sampler in it.
More dramatic is the effect of the Missing at Random assumption on
the estimates of the conditional probability
from 0.1286 in the complete database, to 0.5059 , with a very narrow confidence
interval: 0.4752,0.5367. This narrow interval overestimates the reliability
of the inferred value and excludes the true value inferred from the
complete database.
Execution time was 10 minutes for the Gibbs Sampling and less than
seconds for our system.
5.4 Prediction
The goal of learning bbns is to use them to perform different reasoning tasks,
such as prediction. The aim of this second test is to evaluate the reliability
of the predictions given by the learned bbn.
Materials. We used the bbn learned by the two systems in the previous
test to predict the value of X 2 given an observation about X 1 .
Results. The effect of the strong bias in the estimates returned by the
Gibbs Sampling is remarkable in the predictive performance of the network.
Suppose that X observed and we want to predict the value of X 2 .
Since in this case p(X reduces to p(X Figure 9 c)
plots the marginal probability of p(X Suppose that we use
the estimates learned by the Gibbs Sampling with 75% of the complete
data, when all the entries are missing. The prediction of the Gibbs
Sampling is p(X against the value
0:1286 that we would have inferred from the complete database. Instead,
our method returns the probability interval [0:12; 0:54], thus including the
"true" value.
The results of these experiments match our expectations. The accuracy of
the two systems is overall comparable when data are missing at random.
However, the estimates given by the Gibbs Sampling are prone to bias when
data are systematically missing. The reason of this behavior is easy to
identify. Consider for instance the artificial example in Figure 8. In the
reconstruction of the original database, the Gibbs Sampling exploits the
available information to ascribe the missing entries mainly to X
hence the high confidence on the estimate of the conditional distribution of
These results show also that our method is robust because instead of
"betting" on the most likely complete database that we could infer from the
available information from the incomplete database at hand, our method
returns results which make the problem solver aware of his own ignorance.
This feature is even more important when we consider the remarkable effect
of the strong bias in the estimates returned by the Gibbs Sampling on the
predictive performance of the bbn.
A last word about the execution time. The computational cost of Gibbs
Sampling is a well known issue and the results of our experiments show the
computational advantages of our deterministic method with respect to the
stochastic simulator.
6 Conclusions
Incompleteness is an common feature of real-world databases and the ability
of learning from databases with incomplete data is a basic challenge
researchers have to face in order to move their methods to applications. A
key issue for a method able to learn from incomplete database is the reliability
of the knowledge bases it will generate. The results of our investigation
shows that the common Missing at Random assumption exploited by current
learning methods can dramatically affect the accuracy of their results.
This paper introduced a robust method to learn conditional probabilities
in a bbn which does not rely on this assumption. In order to drop this
assumption, we had to change the overall learning strategy with respect
to traditional Bayesian methods: rather than guessing the value of missing
data on the basis of the available information, our method bounds the set
of all posterior probabilities consistent with the database and proceed by
refining this set as more information becomes available.
The main feature of this method is its robustness with respect to the distribution
of missing data: it does not relies on the assumption that data are
missing at random because it does not try to infer them from the available
information. The basic intuition behind our method is that we are better
off if, rather than trying to complete the database by guessing the value of
missing data, we regard the available information as a set of constraints on
the possible distributions in the database and we reason on the basis of the
set of probability distributions consistent with the database at hand.
An experimental comparison between our method and a powerful stochastic
method shows a remarkable difference in accuracy between the two methods
and the computational advantages of our deterministic method with
respect to the stochastic one.

Acknowledgments

Authors thank Greg Cooper, Pat Langley, Zdenek Zdrahal for their useful
suggestions during the development of this research. Equipment has been
provided by generous donations from Apple Computers and Sun Microsystems



--R

The robust bayesian viewpoint.
Decision making with interval influence diagrams.
Operations for learning with graphical mod- els
A guide to the literature on learning probabilistic networs from data.
A bayesian method for the induction of probabilistic networks from data.
A comparison of sequential learning methods for incomplete data.
Maximum likelihood from incomplete data via the em algorithm.
Upper and lower probabilities induced by multivalued mapping.
Combining clinical judgments and clinical data in expert systems.
Subjective probability as a measure of a non-measurable set
An inequality paradigm for probabilistic knowl- edge
Learning bayesian networks: The combinations of knowledge and statistical data.
Rational belief.
The Enterprise of Knowledge.
Probabilistic inference using markov chain monte carlo methods.
Probabilistic Reasoning in Intelligent Systems: Networks of plausible inference.
Learning efficient classification procedures and their application to chess and games.
An ignorant belief network to forecast glucose concentration from clinical databases.
Ignorant influence diagrams.
A Mathematical Theory of Evidence.
Improved posterior probability estimates from prior and linear constraint system.
Learning in probabilistic expert systems.
Sequential updating of conditional probabilities on directed graphical structures.
Bayesian analysis in expert systems.
Covex bayesian decision theory.
Bugs: A program to perform bayesian inference using gibbs sampling.
Computing probability intervals under independency constraints.
Statistical Reasoning with Imprecise Probabilities.
A posteriori representations based on linear inequality descriptions of a priori conditional probabilities.
--TR
Probabilistic reasoning in intelligent systems: networks of plausible inference
A Bayesian Method for the Induction of Probabilistic Networks from Data
C4.5: programs for machine learning
The EM algorithm for graphical association models with missing data
Learning Bayesian Networks
Irrelevance and parameter learning in Bayesian networks
Bayesian classification (AutoClass)
On the Optimality of the Simple Bayesian Classifier under Zero-One Loss
Bayesian Network Classifiers
Expert Systems and Probabiistic Network Models
Probability Intervals Over Influence Diagrams
Bayesian methods

--CTR
Marco Zaffalon , Marcus Hutter, Robust inference of trees, Annals of Mathematics and Artificial Intelligence, v.45 n.1-2, p.215-239, October   2005
Gert de Cooman , Marco Zaffalon, Updating beliefs with incomplete observations, Artificial Intelligence, v.159 n.1-2, p.75-125, November 2004
Ferat Sahin , M. etin Yavuz , Ziya Arnavut , nder Uluyol, Fault diagnosis for airplane engines using Bayesian networks and distributed particle swarm optimization, Parallel Computing, v.33 n.2, p.124-143, March, 2007
bayesian learning;bayesian networks;bayesian classifiers;missing data;probability intervals