\documentstyle[jair,twoside,11pt,theapa,psfig]{article}
\input epsf %% at top of file

  \jairheading{13}{2000}{155-188}{6/00}{10/00}
  \ShortHeadings{Adaptive Importance Sampling in Bayesian Networks}
  {Cheng \& Druzdzel}
  \firstpageno{155}

\hyphenation{Ba-ye-sian}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{definition}{Definition}

\newcommand{\proof}[1]{\par\noindent{\bf Proof:} {#1} \hfill $\Box$\par}
\newcommand{\Prob}[0]{{\mbox{\rm Pr}}}
\newcommand{\Pa}[0]{{\mbox{\rm Pa}}}
\newcommand{\Paa}[0]{{\mbox{\rm pa}}}
\newcommand{\Anc}[0]{{\mbox{\rm Anc}}}

\newcommand{\ais}[0]{{\sc AIS-BN}}
\newcommand{\sis}[0]{{\sc SIS}}
\newcommand{\lw}[0]{{\sc LW}}
\newcommand{\cpcs}[0]{{\sc CPCS}}
\newcommand{\qmr}[0]{{\sc QMR-DT}}
\newcommand{\andes}[0]{{\sc ANDES}}
\newcommand{\pathfinder}[0]{{\sc PathFinder}}
\newcommand{\pcpt}[0]{ICPT}
\newcommand{\cpt}[0]{CPT}
\newcommand{\mse}[0]{{\it MSE}}

\begin{document}

\title{{\ais}: An Adaptive Importance Sampling Algorithm for Evidential
       Reasoning in Large Bayesian Networks}
\author{\name Jian Cheng \email jcheng@sis.pitt.edu \\
        \name Marek J.\  Druzdzel \email marek@sis.pitt.edu \\
        \addr Decision Systems Laboratory \\
         School of Information Sciences and Intelligent Systems Program \\
         University of Pittsburgh, Pittsburgh, PA 15260 USA}

\maketitle

\begin{abstract}
Stochastic sampling algorithms, while an attractive alternative to
exact algorithms in very large Bayesian network models, have been
observed to perform poorly in evidential reasoning with extremely
unlikely evidence. To address this problem, we propose an adaptive
importance sampling algorithm, \ais, that shows promising
convergence rates even under extreme conditions and seems to
outperform the existing sampling algorithms consistently. Three
sources of this performance improvement are (1) two heuristics for
initialization of the importance function that are based on the
theoretical properties of importance sampling in
finite-dimensional integrals and the structural advantages of
Bayesian networks, (2) a smooth learning method for the importance
function, and (3) a dynamic weighting function for combining
samples from different stages of the algorithm.

We tested the performance of the {\ais} algorithm along with two
state of the art general purpose sampling algorithms, likelihood
weighting \cite{Fung.UAI89,Shachter2.UAI89} and self-importance
sampling \cite{Shachter2.UAI89}. We used in our tests three large
real Bayesian network models available to the scientific
community: the {\cpcs} network \shortcite{Pradhan.UAI94}, the
{\pathfinder} network \cite{heckerman90a}, and the {\andes}
network \cite{conati97}, with evidence as unlikely as $10^{-41}$.
While the {\ais} algorithm always performed better than the other
two algorithms, in the majority of the test cases it achieved
orders of magnitude improvement in precision of the results.
Improvement in speed given a desired precision is even more
dramatic, although we are unable to report numerical results here,
as the other algorithms almost never achieved the precision
reached even by the first few iterations of the {\ais} algorithm.
%
% {\bf Keywords:} Probabilistic reasoning, Bayesian networks,
% stochastic sampling, algorithms, importance sampling
\end{abstract}


\section{Introduction}
\label{sec:intro}

Bayesian networks \cite{pearl88} are increasingly popular tools
for modeling uncertainty in intelligent systems. With practical
models reaching the size of several hundreds of variables
\shortcite<e.g.,>{Pradhan.UAI94,conati97}, it becomes increasingly
important to address the problem of feasibility of probabilistic
inference. Even though several ingenious exact algorithms have
been proposed, in very large models they all stumble on the
theoretically demonstrated NP-hardness of inference
\cite{cooper90}. The significance of this result can be observed
in practice --- exact algorithms applied to large, densely
connected practical networks require either a prohibitive amount
of memory or a prohibitive amount of computation and are unable to
complete. While approximating inference to any desired precision
has been shown to be NP-hard as well \cite{dagum93}, it is for
very complex networks the only alternative that will produce any
result at all. Furthermore, while obtaining the result is crucial
in all applications, precision guarantees may not be
critical for some types of problems and can be traded off against
the speed of computation.
% This, for example, often happens in decision making, where exact
% results may not be needed to prove that one decision option is
% better than another.

A prominent subclass of approximate algorithms is the family of
stochastic sampling algorithms (also called {\em stochastic
simulation} or {\em Monte Carlo} algorithms). The precision
obtained by stochastic sampling generally increases with the
number of samples generated and is fairly unaffected by the
network size. Execution time is fairly independent of the topology
of the network and is linear in the number of samples. Computation
can be interrupted at any time, yielding an anytime property of
the algorithms, important in time-critical applications.

While stochastic sampling performs very well in predictive
inference, diagnostic reasoning, i.e., reasoning from observed
evidence nodes to their ancestors in the network often exhibits
poor convergence. When the number of observations increases,
especially if these observations are unlikely a-priori, stochastic
sampling often fails to converge to reasonable estimates of the
posterior probabilities. Although this problem has been known
since the very first sampling algorithm was proposed by Henrion
\citeyear{Henrion.UAI86}, little has been done to address it
effectively. Furthermore, various sampling algorithms proposed
were tested on simple and small networks, or networks with special
topology, without the presence of extremely unlikely evidence and
the practical significance of this problem has been
underestimated. Given a typical number of samples used in
real-time that are feasible on today's hardware, say $10^6$
samples, the behavior of a stochastic sampling algorithm will be
drastically different for different size networks. While in a
network consisting of 10 nodes and a few observations, it may be
possible to converge to exact probabilities, in very large
networks only a negligibly small fraction of the total sample
space will be probed. One of the practical Bayesian network models
that we used in our tests, a subset of the {\cpcs} network
\cite{Pradhan.UAI94}, consists of 179 nodes. Its total sample
space is larger than $10^{61}$. With $10^6$ samples, we can sample
only $10^{-55}$ fraction of the sample space.

We believe that it is crucial (1) to study the feasibility and
convergence properties of sampling algorithms on very large
practical networks, and (2) to develop sampling algorithms that
will show good convergence under extreme, yet practical
conditions, such as evidential reasoning given extremely unlikely
evidence. After all, small networks can be updated using any of
the existing exact algorithms --- it is precisely the very large
networks where stochastic sampling can be most useful. As to the
likelihood of evidence, we know that stochastic sampling will
generally perform well when it is high \cite{Henrion.UAI86}. So,
it is important to look at those cases in which evidence is very
unlikely. In this paper, we test two existing state of the art
stochastic sampling algorithms for Bayesian networks, likelihood
weighting \cite{Fung.UAI89,Shachter2.UAI89} and self-importance
sampling \cite{Shachter2.UAI89}, on a subset of the {\cpcs}
network with extremely unlikely evidence. We show that they both
exhibit similarly poor convergence rates. We propose a new
sampling algorithm, that we call the {\em adaptive importance
sampling} for Bayesian networks ({\ais}), which is suitable for
evidential reasoning in large multiply-connected Bayesian
networks. The {\ais} algorithm is based on importance sampling,
which is a widely applied method for variance reduction in
simulation that has also been applied in Bayesian networks
\shortcite<e.g.,>{Shachter2.UAI89}. We demonstrate empirically on
three large practical Bayesian network models that the {\ais}
algorithm consistently outperforms the other two algorithms. In
the majority of the test cases, it achieved over two orders of
magnitude improvement in convergence. Improvement in speed given a
desired precision is even more dramatic, although we are unable to
report numerical results here, as the other algorithms never
achieved the precision reached even by the first few iterations of
the {\ais} algorithm. The main sources of improvement are: (1) two
heuristics for the initialization of the importance function that
are based on the theoretical properties of importance sampling in
finite-dimensional integrals and the structural advantages of
Bayesian networks, (2) a smooth learning method for updating the
importance function, and (3) a dynamic weighting function for
combining samples from different stages of the algorithm. We study
the value of the two heuristics used in the {\ais} algorithm: (1)
initialization of the probability distributions of parents of
evidence nodes to uniform distribution and (2) adjusting very
small probabilities in the conditional probability tables, and
show that they both play an important role in the {\ais} algorithm
but only a moderate role in the existing algorithms.

The remainder of this paper is structured as follows.
Section~\ref{sec:importance} first gives a general introduction to
importance sampling in the domain of finite-dimensional integrals,
where it was originally proposed. We show how importance sampling
can be used to compute probabilities in Bayesian networks and how
it can draw additional benefits from the graphical structure of
the network. Then we develop a generalized sampling scheme that
will aid us in reviewing the previously proposed sampling
algorithms and in describing the {\ais} algorithm.
Section~\ref{sec:bnais} describes the {\ais} algorithm. We propose
two heuristics for initialization of the importance function and
discuss their theoretical foundations. We describe a smooth
learning method for the importance function and a dynamic
weighting function for combining samples from different stages of
the algorithm. Section~\ref{sec:experiments} describes the
empirical evaluation of the {\ais} algorithm. Finally,
Section~\ref{sec:discussion} suggests several possible
improvements to the {\ais} algorithm, possible applications of our
learning scheme, and directions for future work.


\section{Importance Sampling Algorithms for Bayesian Networks}
\label{sec:importance}

We feel that it is useful to go back to the theoretical roots of
importance sampling in order to be able to understand the source
of speedup of the {\ais} algorithm relative to the existing state
of the art importance sampling algorithms for Bayesian networks.
We first review the general idea of importance sampling in
finite-dimensional integrals and how it can reduce the sampling
variance. We then discuss the application of importance sampling
to Bayesian networks. Readers interested in more details are
directed to literature on Monte Carlo methods in computation of
finite integrals, such as the excellent exposition by
\shortciteA{rubinstein81} that we are essentially following in the
first section.


\subsection{Mathematical Foundations}
\label{sec:theory}

Let $g({\bf X})$ be a function of $m$ variables ${\bf X} =
(X_1,...,X_m)$ over a domain $\Omega\subset R^m$, such that
computing $g({\bf X})$ for any ${\bf X}$ is feasible. Consider the
problem of approximate computation of the integral
\begin{equation}
  I=\int_\Omega g({\bf X}) \, d{\bf X} \; .
\label{eq:integral}
\end{equation}
Importance sampling approaches this problem by writing the
integral (\ref{eq:integral}) as
\begin{displaymath}
  I=\int_\Omega \frac{g({\bf X})}{f({\bf X})} \, f({\bf X}) \, d{\bf X} \; ,
\end{displaymath}
where $f({\bf X})$, often referred to as the {\em importance
function}, is a probability density function over $\Omega$.
$f({\bf X})$ can be used in importance sampling if there exists an
algorithm for generating samples from $f({\bf X})$ and if the
importance function is zero only when the original function is
zero, i.e., $g({\bf X})$ $\neq 0 \; \Longrightarrow \; f({\bf
X})\neq 0$.

After we have independently sampled $n$ points ${\bf s}_1$, ${\bf
s}_2$, \ldots, ${\bf s}_n$, ${\bf s}_i\in\Omega$, according to the
probability density function $f({\bf X})$, we can estimate the
integral $I$ by
\begin{equation}
 \hat I_n = \frac{1}{n} \sum\limits_{i=1}^n
   \frac{g({\bf s}_i)}{f({\bf s}_i)}
\label{eq:approximate}
\end{equation}
and estimate the variance of $\hat I_n$ by
\begin{equation}
  \widehat{\sigma }^2( \hat I_n)
    = \frac{1}{n\cdot(n-1)}\sum\limits_{i=1}^n
      \left(\frac{g({\bf s}_i)}{f({\bf s}_i)}-\hat I_n\right)^2 \; .
\label{eq:variance}
\end{equation}
It is straightforward to show that this estimator has the
following properties:
\begin{enumerate}
\item
  $E(\hat I_n)=I$
\item
  $\lim_{n\rightarrow\infty} \hat I_n=I$
\item
  $\sqrt{n}\cdot(\hat I_n-I) \;
    \stackrel{n\rightarrow\infty}{\longrightarrow}
      \; \mbox{Normal}(0,\sigma_{f({\bf X})}^2)$, where
\begin{equation}
  \sigma_{f({\bf X})}^2
    = \int_\Omega \left(\frac{g({\bf X})}{f({\bf X})}-I\right)^2
      f({\bf X}) \, d{\bf X}
\label{eq:var1}
\end{equation}
\item
  $E\left(\widehat{\sigma }^2(\hat I_n)\right)
  = \sigma^2(\hat I_n)
  = \sigma_{f({\bf X})}^2/n$
\end{enumerate}
The variance of $\hat I_n$ is proportional to $\sigma_{f({\bf
X})}^2$ and inversely proportional to the number of samples. To
minimize the variance of $\hat I_n$, we can either increase the
number of samples or try to decrease $\sigma_{f({\bf X})}^2$. With
respect to the latter, \shortciteA{rubinstein81} reports the
following useful theorem and corollary.
\begin{theorem}
The minimum of $\sigma_{f({\bf X})}^2$ is equal to
\begin{displaymath}
\sigma_{f({\bf X})}^2
  = \left( \int_\Omega \left| g({\bf X})\right| d{\bf X} \right) ^2-I^2
\end{displaymath}
and occurs when ${\bf X}$ is distributed according to the
following probability density function
\begin{displaymath}
f({\bf X)}=\frac{\left| g({\bf X)}\right| }{\int_\Omega \left|
g({\bf X})\right| d{\bf X}} \; .
\end{displaymath}
\label{th:rubinstein}
\end{theorem}
\begin{corollary}
If $g({\bf X)>}0$, then the optimal probability density function
is
\begin{displaymath}
f({\bf X)}=\frac{g({\bf X)}}I
\end{displaymath}
and $\sigma_{f({\bf X})}^2=0$.
\label{cor:rubinstein}
\end{corollary}
Although in practice sampling from precisely $f({\bf X})=g({\bf
X})/I$ will occur rarely, we expect that functions that are close
enough to it can still reduce the variance effectively. Usually,
the closer the shape of the function $f({\bf X})$ is to the shape
of the function $g({\bf X})$, the smaller is $\sigma_{f({\bf
X})}^2$. In high-dimensional integrals, selection of the
importance function, $f({\bf X})$, is far more critical than
increasing the number of samples, since the former can
dramatically affect $\sigma_{f({\bf X})}^2$. It seems prudent to
put more energy in choosing an importance function whose shape is
as close as possible to that of $g({\bf X})$ than to apply the
brute force method of increasing the number of samples.

It is worth noting here that if $f({\bf X})$ is uniform,
importance sampling becomes a general Monte Carlo sampling.
Another noteworthy property of importance sampling that can be
derived from Equation~\ref{eq:var1} is that we should avoid
$f({\bf X})\ll \left| g({\bf X})-I\cdot f({\bf X})\right|$ in any
part of the domain of sampling, even if $f({\bf X})$ matches well
$g({\bf X})/I$ in important regions. If $f({\bf X})\ll \left|
g({\bf X})-I\cdot f({\bf X})\right|$, the variance can become very
large or even infinite. We can avoid this by adjusting $f({\bf
X})$ to be larger in unimportant regions of the domain of ${\bf
X}$.

While in this section we discussed importance sampling for
continuous variables, the results stated are valid for discrete
variables as well, in which case integration should be substituted
by summation.


\subsection{A Generic Importance Sampling Algorithm for Bayesian Networks}
\label{sec:importancebn}

In the following discussion, all random variables used are
multiple-valued, discrete variables. Capital letters, such as $A$,
$B$, or $C$, denote random variables. Bold capital letters, such
as ${\bf A}$, ${\bf B}$, or ${\bf C}$, denote sets of variables.
Bold capital letter ${\bf E}$ will usually be used to denote the
set of evidence variables. Lower case letters $a$, $b$, $c$ denote
particular instantiations of variables $A$, $B$, and $C$
respectively. Bold lower case letters, such as ${\bf a}$, ${\bf
b}$, ${\bf c}$, denote particular instantiations of sets ${\bf
A}$, ${\bf B}$, and ${\bf C}$ respectively. Bold lower case letter
${\bf e}$, in particular, will be used to denote the observations,
i.e., instantiations of the set of evidence variables ${\bf E}$.
$\Anc(A)$ denotes the set of ancestors of node $A$. $\Pa(A)$
denotes the set of parents (direct ancestors) of node $A$.
$\Paa(A)$ denotes a particular instantiation of $\Pa(A)$.
${\backslash}$ denotes set difference. $\left. {\mbox{\rm
Pa}}(\mathbf{A})\right|  _{\mathbf{E}=\mathbf{e}}$ denotes that we
use the extended vertical bar to indicate substitution of
$\mathbf{e}$ for $\mathbf{E}$ in $\mathbf{A}$.

We know that the joint probability distribution over all variables
of a Bayesian network model, $\Prob({\bf X})$, is the product of
the probability distributions over each of the nodes conditional
on their parents, i.e.,
\begin{equation}
  \Prob({\bf X}) = \prod\limits_{i=1}^n \Prob(X_i|\Pa(X_i)) \; .
\label{eq:joint}
\end{equation}
In order to calculate $\Prob({\bf E}={\bf e})$, we need to sum
over all $\Prob({\bf X}\backslash{\bf E},{\bf E}={\bf e})$.
\begin{equation}
  \Prob({\bf E}={\bf e})
   = \sum\limits_{{\bf X}\backslash{\bf E}}
     \Prob({\bf X}\backslash{\bf E},{\bf E}={\bf e})
\label{eq:sum}
\end{equation}
We can see that Equation~\ref{eq:sum} is almost identical to
Equation~\ref{eq:integral} except that integration is replaced by
summation and the domain $\Omega$ is replaced by ${\bf
X}\backslash{\bf E}$. The theoretical results derived for the
importance sampling that we reviewed in the previous section can
thus be directly applied to computing probabilities in Bayesian
networks.

While there has been previous work on importance sampling-based
algorithms for Bayesian networks, we will postpone the discussion
of this work until the next section. Here we will present a
generic stochastic sampling algorithm that will help us in both
reviewing the prior work and in presenting our algorithm.

\begin{figure}[bp6!]
\centering \fbox{
\parbox{5.3in}{
\renewcommand{\baselinestretch}{1.0}
\begin{minipage}[t]{5.3in}
\begin{enumerate}
\item
  Order the nodes according to their topological order.
\item
  Initialize importance function $\Prob^0({\bf X}\backslash{\bf E})$,
  the desired number of samples $m$, the updating interval $l$, and the score arrays for every node.
\item
  $k\leftarrow 0,$ $T\leftarrow \emptyset$
\item
  {\bf for} $i\leftarrow 1$ {\bf to} $m$ {\bf do}
\item
  \quad {\bf if }($i\;\mbox{mod}\;l==0$) {\bf then}
\item
  \qquad $k\leftarrow k+1$
\item
  \qquad Update importance function $\Prob^k({\bf X}\backslash{\bf E})$
  based on $T$.
  \par
  \quad {\bf end if}
\item
  \quad ${\bf s}_i\leftarrow$
  generate a sample according to $\Prob^k({\bf X}\backslash{\bf E})$
\item
  \quad $T\leftarrow T\cup \{{\bf s}_i\}$
\item
  \quad Calculate $\mbox{Score}({\bf s}_i,\Prob({\bf X}\backslash{\bf E},
    {\bf e}),\Prob^k({\bf X\backslash E}))$ and add it to the corresponding entry of every score array according to the instantiated states.
  \par
  {\bf end for}
\item
  Normalize the score arrays for every node.
\end{enumerate}
\end{minipage}}}
\caption{A generic importance sampling algorithm.} \label{fig:sis}
\end{figure}

The posterior probability $\Prob({\bf a}|{\bf e})$ can be obtained
by first computing $\Prob({\bf a},{\bf e})$ and $\Prob({\bf e})$
and then combining these based on the definition of conditional
probability
\begin{equation}
  \Prob({\bf a}|{\bf e})=\frac{\Prob({\bf a},{\bf e})}{\Prob({\bf e})} \; .
\label{eq:bayes}
\end{equation}
In order to increase the accuracy of results of importance
sampling in computing the posterior probabilities over different
network variables given evidence, we should in general use
different importance functions for $\Prob({\bf a},{\bf e})$ and
for $\Prob({\bf e})$. Doing so increases the computation time only
linearly while the gain in accuracy may be significant given that
obtaining a desired accuracy is exponential in nature. Very often,
it is a common practice to use the same importance function
(usually for $\Prob({\bf e})$) to sample both probabilities. If
the difference between the optimal importance functions for these
two quantities is large, the performance may deteriorate
significantly. Although $\widehat{\Prob}({\bf a,e})$ and
$\widehat{\Prob}({\bf e})$ are unbiased estimators according to
Property~1 (Section~\ref{sec:theory}), $\widehat{\Prob}({\bf
a}|{\bf e})$ obtained by means of Equation~\ref{eq:bayes} is not
an unbiased estimator. However, as the number of samples
increases, the bias decreases and can be ignored altogether when
the sample size is large enough \cite{fishman95}.

Figure~\ref{fig:sis} presents a generic stochastic sampling
algorithm that captures most of the existing sampling algorithms.
Without the loss of generality, we restrict ourselves in our
description to so-called {\em forward sampling}, i.e., generation
of samples in the topological order of the nodes in the network.
The forward sampling order is accomplished by the initialization
performed in Step~1, where parents of each node are placed before
the node itself. In forward sampling, Step~8 of the algorithm, the
actual generation of samples, works as follows. ({\em i}) each
evidence node is instantiated to its observed state and is further
omitted from sample generation; ({\em ii}) each root node is
randomly instantiated to one of its possible states, according to
the importance prior probability of this node, which can be
derived from $\Prob^k({\bf X}\backslash{\bf E})$; ({\em iii}) each
node whose parents are instantiated is randomly instantiated to
one of its possible states, according to the importance
conditional probability distribution of this node given the values
of the parents, which can also be derived from $\Prob^k({\bf
X}\backslash{\bf E})$; ({\em iv}) this procedure is followed until
all nodes are instantiated. A complete instantiation ${\bf s}_i$
of the network based on this method is one sample of the joint
importance probability distribution $\Prob^k({\bf X}\backslash{\bf
E})$ over all variables of the network. The scoring of Step~10
amounts to calculating $\Prob({\bf s}_i,{\bf e})/\Prob^k({\bf
s}_i)$, as required by Equation~\ref{eq:approximate}. The ratio
between the total score sum and the number of samples is an
unbiased estimator of $\Prob({\bf e})$. In Step~10, if we also
count the score sum under the condition ${\bf A}={\bf a}$, i.e.,
that some unobserved variables ${\bf A}$ have the values ${\bf
a}$, the ratio between this score sum and the number of samples is
an unbiased estimator of $\Prob({\bf a},{\bf e})$.

Most existing algorithms focus on the posterior probability
distributions of individual nodes. As we mentioned above, for the
sake of efficiency they count the score sum corresponding to
$\Prob(A=a,{\bf e})$, $A\in {\bf X}\backslash{\bf E}$, and record
it in an score array for node $A$. Each entry of this array
corresponds to a specified state of $A$. This method introduces
additional variance, as opposed to using the importance function
derived from $\Prob^k({\bf X}\backslash{\bf E})$ to sample
$\Prob(A=a,{\bf e})$, $A$ $\in {\bf X}\backslash{\bf E}$,
directly.


\subsection{Existing Importance Sampling Algorithms for Bayesian Networks}
\label{sec:eximportancebn}

The main difference between various stochastic sampling algorithms
is in how they process Steps~2, 7, and~8 in the generic importance
sampling algorithm of Figure~\ref{fig:sis}.

{\em Probabilistic logic sampling} \cite{Henrion.UAI86} is the
simplest and the first proposed sampling algorithm for Bayesian
networks. The importance function is initialized in Step~2 to
$\Prob({\bf X})$ and never updated (Step~7 is null). Without
evidence, $\Prob({\bf X})$ is the optimal importance function for
the evidence set, which is empty anyway. It escapes most authors
that $\Prob({\bf X})$ may be not the optimal importance function
for $\Prob(A=a)$, $A\in{\bf X}$, when $A$ is not a root node. A
mismatch between the optimal and the actually used importance
function may result in a large variance. The sampling process with
evidence is the same as without evidence except that in Step~10 we
do not count the scores for those samples that are inconsistent
with the observed evidence, which amounts to discarding them. When
the evidence is very unlikely, there is a large difference between
$\Prob({\bf X})$ and the optimal importance function. Effectively,
most samples are discarded and the performance of logic sampling
deteriorates badly.

{\em Likelihood weighting} (\lw) \cite{Fung.UAI89,Shachter2.UAI89}
enhances the logic sampling in that it never discards samples. In
likelihood weighting, the importance function in Step~2 is
\begin{displaymath}
  \Prob({\bf X}\backslash{\bf E})
    = \left. \prod_{x_{i} \notin\mathbf{e}}
    \Prob(x_{i}|\Pa(X_{i}))\right| _{\mathbf{E}=\mathbf{e}}.
\end{displaymath}
Likelihood weighting does not update the importance function in
Step~7. Although likelihood weighting is an improvement on logic
sampling, its convergence rate can be still very slow when there
is large difference between the optimal importance function and
$\Prob({\bf X}\backslash{\bf E})$, again especially in situations
when evidence is very unlikely. Because of its simplicity, the
likelihood weighting algorithm has been the most commonly used
simulation method for Bayesian network inference. It often matches
the performance of other, more sophisticated schemes because it is
simple and able to increase its precision by generating more
samples than other algorithms in the same amount of time.

{\em Backward sampling} \cite{Fung.UAI94} changes Step~1 of our
generic algorithm and allows for generating samples from evidence
nodes in the direction that is opposite to the topological order
of nodes in the network. In Step~2, backward sampling uses the
likelihood of some of the observed evidence and some instantiated
nodes to calculate $\Prob^0({\bf X}\backslash{\bf E})$. Although
Fung and del Favero mentioned the possibility of dynamic node
ordering, they did not propose any scheme for updating the
importance function in Step~7. Backward sampling suffers from
problems that are similar to those of likelihood weighting, i.e.,
a possible mismatch between its importance function and the
optimal importance function can lead to poor convergence.

{\em Importance sampling} \cite{Shachter2.UAI89} is the same as
our generic sampling algorithm. Shachter and Peot introduced two
variants of importance sampling: {\em self-importance} (\sis) and
{\em heuristic importance}. The importance function used in the
first step of the self-importance algorithm is
\begin{displaymath}
  \Prob^0({\bf X}\backslash{\bf E})
    = \left. \prod_{x_{i} \notin\mathbf{e}}
    \Prob(x_{i}|\Pa(X_{i}))\right|_{\mathbf{E}=\mathbf{e}} \; .
\end{displaymath}
This function is updated in Step~7. The algorithm tries to revise
the conditional probability tables ({\cpt}s) periodically in order
to make the sampling distribution gradually approach the posterior
distribution. Since the same data are used to update the
importance function and to compute the estimator, this process
introduces bias in the estimator. Heuristic importance first
removes edges from the network until it becomes a polytree, and
then uses a modified version of the polytree algorithm
\cite{pearl86} to compute the likelihood functions for each of the
unobserved nodes. $\Prob^0({\bf X}\backslash{\bf E})$ is a
combination of these likelihood functions with $\Prob({\bf
X}\backslash{\bf E},{\bf e})$. In Step~7 heuristic importance does
not update $\Prob^k({\bf X}\backslash{\bf E})$. As
\shortciteA{Shachter2.UAI89} point out, this heuristic
importance function can still lead to a bad approximation of the
optimal importance function. There exist also other algorithms
such as a combination of self-importance and heuristic importance
\cite{Shachter2.UAI89,shwe91b}. Although some researchers
suggested that this may be a promising direction for the work on
sampling algorithms, we have not seen any results that would
follow up on this.

A separate group of stochastic sampling methods is formed by
so-called {\em Markov Chain Monte Carlo (MCMC)} methods that are
divided into Gibbs sampling, Metropolis sampling, and Hybrid Monte
Carlo sampling \cite{geman84,gilks96,mackay98}. Roughly speaking,
these methods draw random samples from an unknown target
distribution $f({\bf X})$ by biasing the search for this
distribution towards higher probability regions. When applied to
Bayesian networks \cite{pearl87,chavez90} this approach determines
the sampling distribution of a variable from its previous sample
given its Markov blanket \cite{pearl88}. This corresponds to
updating $\Prob^k({\bf X}\backslash{\bf E})$ when sampling every
node. $\Prob^k({\bf X}\backslash{\bf E})$ will converge to the
optimal importance function for $\Prob({\bf e})$ if $\Prob^0({\bf
X}\backslash{\bf E})$ satisfies some ergodic properties
\cite{york92}. Since the convergence to the limiting distribution
is very slow and calculating updates of the sampling distribution
is costly, these algorithms are not used in practice as often as
the simple likelihood weighting scheme.

There are also some other simulation algorithms, such as bounded
variance algorithm \cite{dagum97} and the {\sc AA} algorithm
\shortcite{dagum95}, which are essentially based on the {\lw}
algorithm and the Stopping-Rule Theorem \cite{dagum95}.
\shortciteA{cano96} proposed another importance sampling algorithm
that performed somewhat better than {\lw} in cases with extreme
probability distributions, but, as the authors state, in general
cases it ``produced similar results to the likelihood weighting
algorithm.'' \shortciteA{hernandez98} also applied importance
sampling and reported a moderate improvement on likelihood
weighting.


\subsection{Practical Performance of the Existing Sampling Algorithms}
\label{sec:motivation}

The largest network that has been tested using sampling algorithms
is {\qmr} (Quick Medical Reference --- Decision Theoretic)
\shortcite{shwe91c,shwe91b}, which contains 534 adult diseases and
4,040 findings, with 40,740 arcs depicting disease-to-finding
dependencies. The {\qmr} network belongs to a class of special
bipartite networks and its structure is often referred to as BN2O
\cite{Henrion.UAI91}, because of its two-layer composition:
disease nodes in the top layer and finding nodes in the bottom
layer. Shwe and colleagues used an algorithm combining
self-importance and heuristic importance and tested its
convergence properties on the {\qmr} network. But since the
heuristic method {\em iterative tabular Bayes} (ITB) that makes
use of a version of Bayes' rule is designed for the BN2O networks,
it cannot be generalized to arbitrary networks. Although Shwe and
colleagues concluded that Markov blanket scoring and
self-importance sampling significantly improve the convergence
rate in their model, we cannot extend this conclusion to general
networks. The computation of Markov blanket scoring is more
complex in a general multi-connected network than in a BN2O
network. Also, the experiments conducted lacked a gold-standard
posterior probability distribution that could serve to judge the
convergence rate.

\shortciteA{Pradhan.UAI96} tested an efficient version of the
{\lw} algorithm --- bounded variance algorithm \cite{dagum97} and
the {\sc AA} algorithm \cite{dagum95} on a 146 node, multiply
connected medical diagnostic Bayesian network. One limitation in
their tests is that the probability of evidence in the cases
selected for testing was rather high. Although over 10\% of the
cases had the probability of evidence on the order of $10^{-8}$ or
smaller, a simple calculation based on the reported mean
$\mu=34.5$ number of evidence nodes, shows that the average
probability of an observed state of an evidence node conditional
on its direct predecessors was on the order of
$(10^{-8})^{1/34.5}\approx 0.59$. Given that their algorithm is
essentially based on the {\lw} algorithm, based on our tests we
suspect that the performance will deteriorate on cases where the
evidence is very unlikely. Both algorithms focus on the marginal
probability of one hypothesis node. If there are many queried
nodes, the efficiency may deteriorate.

We have tested the algorithms discussed in
Section~\ref{sec:eximportancebn} on several large networks. Our
experimental results show that in cases with very unlikely
evidence, none of these algorithms converges to reasonable
estimates of the posterior probabilities within a reasonable
amount of time. The convergence becomes worse as the number of
evidence nodes increases. Thus, when using these algorithms in
very large networks, we simply cannot trust the results. We will
present results of tests of the {\lw} and {\sis} algorithms in
more detail in Section~\ref{sec:experiments}.


\section{{\ais}: Adaptive Importance Sampling for Bayesian Networks}
\label{sec:bnais}

The main reason why the existing stochastic sampling algorithms
converge so slowly is that they fail to learn a good importance
function during the sampling process and, effectively, fail to
reduce the sampling variance. When the importance function is
optimal, such as in probabilistic logic sampling without any
evidence, each of the algorithms is capable of converging to
fairly good estimates of the posterior probabilities within
relatively few samples. For example, assuming that the posterior
probabilities are not extreme (i.e., larger than say 0.01), as few
as 1,000 samples may be sufficient to obtain good estimates. In
this section, we present the adaptive importance sampling
algorithm for Bayesian networks (\ais) that, as we will
demonstrate in the next section, performs very well on most tests.
We will first describe the details of the algorithm and prove two
theorems that are useful in learning the optimal importance
sampling function.


\begin{figure}[hbt]
\centering \fbox{
\parbox{5.3in}{
\renewcommand{\baselinestretch}{1.0}
\begin{minipage}[t]{5.3in}
\begin{enumerate}
\item
  Order the nodes according to their topological order.
\item
  Initialize importance function $\Prob^0({\bf X}\backslash{\bf E})$
  using some heuristic methods, initialize weight $w^0$, set the
  desired number of samples $m$ and the updating interval $l$, initialize the score arrays for every node.
\item
  $k\leftarrow 0$, $T\leftarrow \emptyset$, $w_{TScore}\leftarrow
  0$, $w_{\mbox{\it sum}}\leftarrow 0$
\item
  {\bf for} $i\leftarrow 1$ {\bf to} $m$ {\bf do}
\item
  \quad {\bf if} ($i\;\mbox{mod}\;l==0$) {\bf then}
\item
  \qquad $k\leftarrow k+1$
\item
  \qquad Update the importance function $\Prob^k({\bf X}\backslash{\bf E})$
  and $w^k$ based on $T$.
  \par
  \quad {\bf end if}
\item
  \quad ${\bf s}_i\leftarrow$ generate a sample according to
  $\Prob^k({\bf X}\backslash{\bf E})$
\item
  \quad $T\leftarrow T\cup \{{\bf s}_i\}$
\item
  \quad $w_{\mbox{\it iScore}}\leftarrow$
    Score (${\bf s}_i$, $\Prob({\bf X}\backslash{\bf E},{\bf e})$,
    $\Prob^k({\bf X}\backslash{\bf E})$, $w^k$)
\item
  \quad $w_{TScore}\leftarrow w_{TScore}+w_{iScore}$
  \par
  \quad{\footnotesize (Optional: add $w_{iScore}$ to the corresponding entry of every score array)}
\item
  \quad $w_{\mbox{\it sum}}\leftarrow w_{\mbox{\it sum}}+w^k$
  \par
  {\bf end for}
\item
  Output estimate of $\Prob({\bf E})$ as $w_{TScore}/w_{\mbox{\it sum}}$
  \par
  {\footnotesize (Optional: Normalize the score arrays for every node)}
\end{enumerate}
\end{minipage}}}
\caption{The adaptive importance sampling for Bayesian Networks
(\ais) algorithm.} \label{fig:ais}
\end{figure}

\subsection{Basic Algorithm --- \ais}
\label{subsec:ba}

Compared with importance sampling used in normal
finite-dimensional integrals, importance sampling used in Bayesian
networks has several significant advantages. First, the network
joint probability distribution $\Prob({\bf X})$ is decomposable
and can be factored into component parts. Second, the network has
a clear structure, which represents many conditional independence
relationships. These properties are very helpful in estimating the
optimal importance function.

The basic {\ais} algorithm is presented in Figure~\ref{fig:ais}.
The main differences between the {\ais} algorithm and the basic
importance sampling algorithm in Figure~\ref{fig:sis} is that we
introduce a monotonically increasing weight function $w^k$ and two
effective heuristic initialization methods in Step~2. We also
introduce a special learning component in Step~7 to let the
updating process run more smoothly, avoiding oscillation of the
parameters. The score processing in Step~10 is
\begin{displaymath}
  w_{\mbox{\it iScore}}
    = w^k\frac{\Prob({\bf s}_i,{\bf e})}{\Prob^k({\bf s}_i)} \; .
\end{displaymath}
Note that in this respect the algorithm in Figure~\ref{fig:sis}
becomes a special case of {\ais} when $w^k=1$. The reason why we
use $w^k$ is that we want to give different weights to the
sampling results obtained at different stages of the algorithm. As
each stage updates the importance function, they will all have
different distance from the optimal importance function. We
recommend that $w^k\propto 1/\widehat{\sigma}^k$, where
$\widehat{\sigma}^k$ is the standard deviation estimated in
stage~$k$ using Equation~\ref{eq:variance}.\footnote{A similar
weighting scheme based on variance was apparently developed
independently by Ortiz and Kaelbling \citeyear{Ortiz.UAI00}, who
recommend the weight $w^k\propto 1/(\widehat{\sigma}^k)^2$.} In
order to keep $w^k$ monotonically increasing, if $w^k$ is smaller
than $w^{k-1}$, we adjust its value to $w^{k-1}$. This weighting
scheme may introduce bias into the final result. Since the initial
importance sampling functions are often inefficient and introduce
big variance into the results, we also recommend that $w^k=0$ in
the first few stages of the algorithm. We have designed this
weighting scheme to reflect the fact that in practice estimates
with very small estimated variance are usually good estimates.


\subsection{Modifying the Sampling Distribution in \ais}

Based on the theoretical considerations of
Section~\ref{sec:theory}, we know that the crucial element of the
algorithm is converging on a good approximation of the optimal
importance function. In what follows, we first give the optimal
importance function for calculating $\Prob({\bf E}={\bf e})$ and
then discuss how to use the structural advantages of Bayesian
networks to approximate this function. In the sequel, we will use
the symbol $\rho$ to denote the importance sampling function and
$\rho^*$ to denote the optimal importance sampling function.

Since $\Prob({\bf X}\backslash{\bf E},{\bf E}={\bf e})>0$, from
Corollary~\ref{cor:rubinstein} we have
\begin{displaymath}
  \rho({\bf X}\backslash {\bf E})
    = \frac{\Prob({\bf X}\backslash{\bf E},{\bf E}={\bf e})}{\Prob({\bf E}={\bf e})}
    = \Prob({\bf X}|{\bf E}={\bf e}) \; .
\end{displaymath}
The following corollary captures this result.
\begin{corollary}
The optimal importance sampling function $\rho^*({\bf
X}\backslash{\bf E)}$ for calculating $\Prob({\bf E}={\bf e})$ in
{\em Equation~\ref{eq:sum}} is $\Prob({\bf X}|{\bf E}={\bf e})$.
\label{coll2}
\end{corollary}

Although we know the mathematical expression for the optimal
importance sampling function, it is difficult to obtain this
function exactly. In our algorithm, we use the following
importance sampling function
\begin{equation}
  \rho({\bf X}\backslash{\bf E})
    = \prod\limits_{i=1}^n \Prob(X_i|\Pa(X_i),{\bf E}) \; .
\label{eq:ejoint}
\end{equation}
This function partially considers the effect of all the evidence
on every node during the sampling process. When the network
structure is the same as that of the network which has absorbed
the evidence, this function is the optimal importance sampling
function. It is easy to learn and, as our experimental results
show, it is a good approximation to the optimal importance
sampling function. Theoretically, when the posterior structure of
the model changes drastically as the result of observed evidence,
this importance sampling function may perform poorly. We have
tried to find practical networks where this would happen, but to
the day have not encountered a drastic example of this effect.

From Section~\ref{sec:importancebn}, we know that the score sums
corresponding to $\{x_i,\Paa(X_i),{\bf e}\}$ can yield an unbiased
estimator of $\Prob(x_i,\Paa(X_i),{\bf e})$. According to the
definition of conditional probability, we can get an estimator of
$\Prob^{\prime}(x_i|\Paa(X_i),{\bf e})$. This can be achieved by
maintaining an updating table for every node, the structure of
which mimicks the structure of the {\cpt}. Such tables allow us to
decompose the above importance function into components that can
be learned individually. We will call these tables the {\em
importance conditional probability tables} (\pcpt).

\begin{definition}
An {\em importance conditional probability table (\pcpt)} of a
node $X$ is a table of posterior probabilities
$\Prob(X|\Pa(X),{\bf E}={\bf e})$ conditional on the evidence and
indexed by its immediate predecessors, $\Pa(X)$.
\end{definition}
The {\pcpt} tables will be modified during the process of learning
the importance function. Now we will prove a useful theorem that
will lead to considerable savings in the learning process.
\begin{theorem}
\begin{equation}
  X_i\in{\bf X}, X_i\notin\Anc({\bf E})
  \Rightarrow \Prob(X_i|\Pa(X_i),{\bf E})=\Prob(X_i|\Pa(X_i)) \; .
\label{eq:sim}
\end{equation}
\label{th3}
\end{theorem}
\proof{Suppose we have set the values of all the parents of node
$X_i$ to $\Paa(X_i)$. Node $X_i$ is dependent on evidence ${\bf
E}$ given $\Paa(X_i)$ only when $X_i$ is $d$-connecting with ${\bf
E}$ given $\Paa(X_i)$ \cite{pearl88}. According to the definition
of $d$-connectivity, this happens only when there exists a member
of $X_i$'s descendants that belongs to the set of evidence nodes
${\bf E}$. In other words $X_i\notin\Anc({\bf E})$.}

Theorem~\ref{th3} is very important for the {\ais} algorithm. It
states essentially that the {\pcpt} tables of those nodes that are
not ancestors of the evidence nodes are equal to the {\cpt} tables
throughout the learning process. We only need to learn the {\pcpt}
tables for the ancestors of the evidence nodes. Very often this
can lead to significant savings in computation. If, for example,
all evidence nodes are root nodes, we have our {\pcpt} tables for
every node already and the {\ais} algorithm becomes identical to
the likelihood weighting algorithm. Without evidence, the {\ais}
algorithm becomes identical to the probabilistic logic sampling
algorithm.

\begin{figure}[bp13!]
\centering \fbox{
\parbox{5.3in}{
\renewcommand{\baselinestretch}{1.0}
\begin{minipage}[t]{5.3in}
{\em Input:} Initialized importance function $\Prob^0({\bf
X}\backslash{\bf E})$, learning rate $\eta (k)$.
\\
{\em Output:} An estimated importance function $\Prob^S({\bf
X}\backslash{\bf E})$.
\par
{\bf for} stage $k\leftarrow 0$ {\bf to} $S$ {\bf do}
\begin{enumerate}
\item
  Sample $l$ points ${\bf s}_1^k$, ${\bf s}_2^k$, \ldots,
  ${\bf s}_l^k$ independently according to the current importance
  function $\Prob^k({\bf X}\backslash{\bf E})$.
\item
  For every node $X_i$ such that $X_i\in{\bf X}\backslash{\bf E}$
  and $X_i\notin\Anc({\bf E})$ count score sums corresponding
  to $\{x_i,\Paa(X_i),{\bf e}\}$ and estimate
  $\Prob^{\prime}(x_i|\Paa(X_i),{\bf e)}$ based on
  ${\bf s}_1^k$, ${\bf s}_2^k$, \ldots, ${\bf s}_l^k$.
\item
  Update $\Prob^k({\bf X}\backslash{\bf E})$ according to
  the following formula:
  \begin{eqnarray*}
  \lefteqn{\Prob^{k+1}(x_i|\Paa(X_i),{\bf e}) =} \\
    && \Prob^k(x_i|\Paa(X_i),{\bf e})+\eta(k)
    \cdot\left(\Prob^{\prime}(x_i|\Paa(X_i),
    {\bf e})-\Prob^k(x_i|\Paa(X_i),{\bf e})\right)
  \end{eqnarray*}
\end{enumerate}
{\bf end for}
\end{minipage}}}
\caption{The {\ais} algorithm for learning the optimal importance
function.} \label{fig:learn}
\end{figure}

It is worth pointing out that for some $X_i$,
$\Prob(X_i|P_a(X_i),{\bf E})$ (i.e., the {\pcpt} table for
$X_i$), can be easily calculated using exact methods. For example,
when $X_i$ is the only parent of an evidence node $E_j$ and $E_j$
is the only child of $X_i$, the posterior probability distribution
of $X_i$ is straightforward to compute exactly. Since the focus of
the current paper is on sampling, the test results reported in
this paper do not include this improvement of the {\ais}
algorithm.

Figure~\ref{fig:learn} lists an algorithm that implements Step~7
of the basic {\ais} algorithm listed in Figure~\ref{fig:ais}. When
we estimate $\Prob^{\prime}(x_i|\Paa(X_i),{\bf e})$, we only use
the samples obtained at the current stage. One reason for this is
that the information obtained in previous stages has been absorbed
by $\Prob^k({\bf X}\backslash{\bf E})$. The other reason is that
in principle, each successive iteration is more accurate than the
previous one and the importance function is closer to the optimal
importance function. Thus, the samples generated by
$\Prob^{k+1}({\bf X}\backslash{\bf E})$ are better than those
generated by $\Prob^k({\bf X}\backslash{\bf E})$.
$\Prob^{\prime}(X_i|\Paa(X_i),{\bf e})-\Prob^k(X_i|\Paa(X_i),{\bf
e})$ corresponds to the vector of first partial derivatives in the
direction of the maximum decrease in the error. $\eta(k)$ is a
positive function that determines the {\em learning rate}. When
$\eta(k)=0$ (lower bound), we do not update our importance
function. When $\eta(k)=1$ (upper bound), at each stage we discard
the old function. The convergence speed is directly related to
$\eta(k)$. If it is small, the convergence will be very slow due
to the large number of updating steps needed to reach a local
minimum. On the other hand, if it is large, convergence rate will
be initially very fast, but the algorithm will eventually start to
oscillate and thus may not reach a minimum. There are many papers
in the field of neural network learning that discuss how to choose
the learning rate and let estimated importance function converge
quickly to the destination function. Any method that can improve
learning rate should be applicable to this algorithm. Currently,
we use the following function proposed by Ritter et al.\
\citeyear{ritter91}
\begin{equation}
  \eta(k)=a\left(\frac{b}{a}\right)^{k/k_{\max }} ,
\label{eq:learning}
\end{equation}
where $a$ is the initial learning rate and $b$ is the learning
rate in the last step. This function has been reported to perform
well in neural network learning \cite{ritter91}.


\subsection{Heuristic Initialization in {\ais}}
\label{sec:heuristic}

The dimensionality of the problem of Bayesian network inference is
equal to the number of variables in a network, which in the
networks considered in this paper can be very high. As a result,
the learning space of the optimal importance function is very
large. Choice of the initial importance function $\Prob^0({\bf
X}\backslash{\bf E})$ is an important factor affecting the
learning --- an initial value of the importance function that is
close to the optimal importance function can greatly affect the
speed of convergence. In this section, we present two heuristics
that help to achieve this goal.

Due to their explicit encoding of the structure of a decomposable
joint probability distribution, Bayesian networks offer
computational advantages compared to finite-dimensional integrals.
A possible first approximation of the optimal importance function
is the prior probability distribution over the network variables,
$\Prob({\bf X)}$. We propose an improvement on this
initialization. We know that the effect of evidence nodes on a
node will be attenuated as the path length of that node to
evidence nodes is increased \cite{henrion89a} and the most
affected nodes are the direct ancestors of the evidence nodes.
Initializing the {\pcpt} tables of the parents of the evidence
nodes to uniform distributions in our experience improves the
convergence rate. Furthermore, the {\cpt} tables of the parents of
an evidence node $E$ may be not favorable to the observed state
$e$ if the probability of $E=e$ without any condition is less than
a small value, such as $\Pr (E=e)<1/(2\cdot n_{E})$, where $n_{E}$
is the number of outcomes of node $E$. Based on this observation,
we change the {\cpt} tables of the parents of an evidence node $E$
to uniform distributions in our experiment only when
$\Pr(E=e)<1/(2\cdot n_{E})$, otherwise we leave them unchanged.
This kind of initialization involves the knowledge of $\Pr (E=e)$,
the marginal probability without evidence. Probabilistic logic
sampling \cite{Henrion.UAI86} enhanced by Latin hypercube sampling
\cite{cheng00a} or quasi-Monte Carlo methods \cite{Cheng.UAI00}
will produce a very good estimate of $\Pr (E=e)$. This is an
one-time effort that can be made at the model building stage and
is worth pursuing to any desired precision.

Another serious problem related to sampling are extremely small
probabilities. Suppose there exists a root node with a state $s$
that has the prior probability $\Prob(s)=0.0001$. Let the
posterior probability of this state given evidence be
$\Prob(s|{\bf E})=0.8$. A simple calculation shows that if we
update the importance function every $1,000$ samples, we can
expect to hit $s$ only once every $10$ updates. Thus $s$'s
convergence rate will be very slow. We can overcome this problem
by setting a threshold $\theta$ and replacing every probability
$p<\theta$ in the network by $\theta$.\footnote{This
initialization heuristic was apparently developed independently by
Ortiz and Kaelbling \citeyear{Ortiz.UAI00}.} At the same time, we
subtract ($\theta -p$) from the largest probability in the same
conditional probability distribution. For example, the value of
$\theta =10/l$, where $l$ is the updating interval, will allow us
to sample $10$ times more often in the first stage of the
algorithm. If this state turns out to be more likely (having a
large weight), we can increase its probability even more in order
to converge to the correct answer faster. Considering that we
should avoid $f({\bf X})\ll \left| g({\bf X})-I\cdot f({\bf
X})\right|$ in an unimportant region as discussed in
Section~\ref{sec:theory}, we need to make this threshold larger.
We have found that the convergence rate is quite sensitive to this
threshold. Based on our empirical tests, we suggest to use $\theta
=0.04$ in networks whose maximum number of outcomes per node does
not exceed five. A smaller threshold might lead to fast
convergence in some cases but slow convergence in others. If one
threshold does not work, changing it in a specific network will
usually improve convergence rate.


\subsection{Selection of Parameters}

There are several tunable parameters in the {\ais} algorithm. We
base the choice of these parameters on the Central Limit Theorem
(CLT). According to CLT, if $Z_1$, $Z_2$, \ldots, $Z_n$ are
independent and identically distributed random variables with
$E(Z_{i})=\mu_{Z}$ and Var$(Z_{i})=\sigma_{Z}^{2}$, $i=1$, ...,
$n$, then $\overline{Z} =(Z_{1}+...+Z_{n})/n$ is approximately
normally distributed when $n$ is sufficiently large. Thus,
\begin{equation}
{\lim_{n\rightarrow \infty }}P(\frac{\left|
  \overline{Z}-\mu_{z}\right| }{ \mu_{z}}
  \geq \frac{\sigma_{Z}/\sqrt{n}}{\mu_z}\cdot t)
  =\frac{2}{\sqrt{ 2\pi }}\int_{t}^{\infty}e^{-x^{2}/2}dx \; .
\label{eq:CLT}
\end{equation}
Although this approximation holds when $n$ approaches infinity,
CLT is known to be very robust and lead to excellent
approximations even for small $n$. The formula of
Equation~\ref{eq:CLT} is an {\em (}$\varepsilon_{r}${\em , }$
\delta ${\em ) Relative Approximation,} which is an estimate
$\overline{\mu}$ of $\mu$ that satisfies
\begin{displaymath}
P(\frac{\left|\overline{\mu}-\mu\right|}{\mu}
  \geq \varepsilon_{r})\leq \delta \; .
\end{displaymath}
If $\delta$ has been fixed,
\begin{displaymath}
\varepsilon_{r}=\frac{\sigma_{Z}/\sqrt{n}}{\mu_{z}}\cdot \Phi
_{Z}^{-1}(\frac{\delta }{2}) \; ,
\end{displaymath}
where
$\Phi_{Z}(z)=\frac{1}{\sqrt{2\pi}}\int_{z}^{\infty}e^{-x^{2}/2}dx$.
Since in our sampling problem, $\mu_{z}$ (corresponding to
$\Prob({\bf E})$ in Figure~\ref{fig:ais}) has been fixed, setting
$\varepsilon_{r}$ to a smaller value amounts to letting
$\sigma_{Z}/\sqrt{n}$ be smaller. So, we can adjust the parameters
based on $\sigma_{Z}/\sqrt{n}$, which can be estimated using
Equation~\ref{eq:variance}. It is also the theoretical intuition
behind our recommendation $w^{k}\propto 1/\widehat{\sigma }^{k}$ in
Section~\ref{subsec:ba}. While we expect that this should work well
in most networks, no guarantees can be given here --- there exist
always some extreme cases in sampling algorithms in which no good
estimate of variance can be obtained.


\subsection{A Generalization of {\ais}:
  The Problem of Estimating $\Prob({\bf a}|{\bf e})$}
\label{sec:estimation}

A typical focus of systems based on Bayesian networks is the
posterior probability of various outcomes of individual variables
given evidence, $\Prob(a|{\bf e})$. This can be generalized to the
computation of the posterior probability of a particular
instantiation of a set of variables given evidence, i.e.,
$\Prob({\bf A}={\bf a}|{\bf e})$. There are two methods that are
capable of performing this computation. The first method is very
efficient at the expense of precision. The second method is less
efficient, but offers in general better convergence rates. Both
methods are based on Equation~\ref{eq:bayes}.

The first method reuses the samples generated to estimate
$\Prob({\bf e})$ in estimating $\Prob({\bf a},{\bf e})$.
Estimation of $\Prob({\bf a},{\bf e})$ amounts to counting the
scored sum under the condition ${\bf A}={\bf a}$. The main
advantage of this method is its efficiency --- we can use the same
set of samples to estimate the posterior probability of any state
of a subset of the network given evidence. Its main disadvantage
is that the variance of the estimated $\Prob({\bf a}, {\bf e})$
can be large, especially when the numerical value of $\Prob({\bf
a}|{\bf e})$ is extreme. This method is the most widely used
approach in the existing stochastic sampling algorithms.

The second method, used much more rarely
\shortcite<e.g.,>{cano96,Pradhan.UAI96,dagum97}, calls for
estimating $\Prob({\bf e})$ and $\Prob({\bf a},{\bf e})$
separately. After estimating $\Prob({\bf e})$, an additional call
to the algorithm is made for each instantiation ${\bf a}$ of the
set of variables of interest ${\bf A}$. $\Prob({\bf a},{\bf e})$
is estimated by sampling the network with the set of observations
${\bf e}$ extended by ${\bf A}={\bf a}$. The main advantage of
this method is that it is much better at reducing variance than
the first method. Its main disadvantage is the computational cost
associated with sampling for possibly many combinations of states
of nodes of interest.

\shortciteA{cano96} suggested a modified version of the second
method. Suppose that we are interested in the posterior
distribution $\Prob({\bf a}_i|{\bf e})$ for all possible values
${\bf a_i}$ of ${\bf A}$, $i=1$, 2, \ldots, $k$. We can estimate
$\Prob({\bf a}_i,{\bf e})$ for each $i=1$, \ldots, $k$ separately,
and use the value $\sum_{i=1}^k \Prob({\bf a}_i,{\bf e})$ as an
estimate for $\Prob({\bf e})$. The assumption behind this approach
is that the estimate of $\Prob({\bf e})$ will be very accurate
because of the large sample from which it is drawn. However, even
if we can guarantee small variance in every $\Prob({\bf a}_i,{\bf
e})$, we cannot guarantee that their sum will also have a small
variance. So, in the {\ais} algorithm we only use the pure form of
each of the methods. The algorithm listed in Figure~\ref{fig:ais}
is based on the first method when the optional computations in
Steps~12 and~13 are performed. An algorithm corresponding to the
second method skips the optional steps and calls the basic {\ais}
algorithm twice to estimate $\Prob({\bf e})$ and $\Prob({\bf
a},{\bf e})$ separately.

The first method is very attractive because of its simplicity and
possible computational efficiency. However, as we have shown in
Section~\ref {sec:importancebn}, the performance of a sampling
algorithm that uses just one set of samples (as in the first
method above) to estimate $\Prob({\bf a}|{\bf e})$ will
deteriorate if the difference between the optimal importance
functions for $\Prob({\bf a}$,${\bf e})$ and $\Prob({\bf e})$ is
large. If the main focus of the computation is high accuracy of
the posterior probability distribution of a small number of nodes,
we strongly recommend to use the algorithm based on the second
method. Also, this algorithm can be easily used to estimate
confidence intervals of the solution.

\section{Experimental Results}
\label{sec:experiments}

In this section, we first describe the experimental method used in
our tests. Our tests focus on the {\cpcs} network, which is one of
the largest and most realistic networks available and for which we
know precisely which nodes are observable. We were, therefore,
able to generate very realistic test cases. Since the {\ais}
algorithm uses two initialization heuristics, we designed an
experiment that studies the contribution of each of these two
heuristics to the performance of the algorithm. To probe the
extent of {\ais} algorithm's excellent performance, we test it on
several real and large networks.


\subsection{Experimental Method}

We performed empirical tests comparing the {\ais} algorithm to the
likelihood weighting ({\lw}) and the self-importance sampling
({\sis}) algorithms. The two algorithms are basically the state of
the art general purpose belief updating algorithms. The AA
\cite{dagum95} and the bounded variance \cite{dagum97} algorithms,
which were suggested by a reviewer, are essentially enhanced
special purpose versions of the basic {\lw} algorithm. Our
implementation of the three algorithms relied on essentially the
same code with separate functions only when the algorithms
differed. It is fair to assume, therefore, that the observed
differences are purely due to the theoretical differences among
the algorithms and not due to the efficiency of implementation. In
order to make the comparison of the {\ais} algorithm to {\lw} and
{\sis} fair, we used the first method of computation
(Section~\ref{sec:estimation}), i.e., one that relies on single
sampling rather than calling the basic {\ais} algorithm twice.

We measured the accuracy of approximation achieved by the
simulation in terms of the Mean Square Error ({\mse}), i.e.,
square root of the sum of square differences between
$\Prob^{\prime }(x_{ij})$ and $\Prob(x_{ij})$, the sampled and the
exact marginal probabilities of state $j$ ($j=1,2,\ldots ,n_{i}$)
of node $i$, such that $X_{i}\notin {\bf E}$. More precisely,
\begin{displaymath}
{\mbox{\mse}}=\sqrt{\frac{1}{\sum_{X_{i}\in {\bf N}\backslash {\bf
E}}n_{i}} \sum_{X_{i}\in {\bf N}\backslash{\bf
E}}\sum_{j=1}^{n_{i}}(\Prob ^{\prime
}(x_{ij})-\Prob(x_{ij}))^{2}}\;\;,
\end{displaymath}
where ${\bf N}$ is the set of all nodes, ${\bf E}$ is the set of
evidence nodes, and $n_{i}$ is the number of outcomes of node $i$.
In all diagrams, the reported {\mse} is averaged over 10 runs. We
used the clustering algorithm \cite{lauritzen88} to compute the
gold standard results for our comparisons of the mean square
error. We performed all experiments on a Pentium~II, 333 MHz
Windows computer.

While {\mse} is not perfect, it is the simplest way of capturing
error that lends itself to further theoretical analysis. For
example, it is possible to derive analytically the idealized
convergence rate in terms of {\mse}, which, in turn, can be used
to judge the quality of the algorithm. {\mse} has been used in
virtually all previous tests of sampling algorithms, which allows
interested readers to tie the current results to the past studies.
A reviewer offered an interesting suggestion of using
cross-entropy or some other technique that weights small changes
near zero much more strongly than the equivalent size change in
the middle of the $[0,1]$ interval. Such measure would penalize
the algorithm for imprecisions of possibly several orders of
magnitude in very small probabilities. While this idea is
interesting, we are not aware of any theoretical reasons as to why
this measure would make a difference in comparisons between
{\ais}, {\lw} and {\sis} algorithms. The {\mse}, as we mentioned
above, will allow us to compare the empirically determined
convergence rate to the theoretically derived ideal convergence
rate. Theoretically, the {\mse} is inversely proportional to the
square root of the sample size.

Since there are several tunable parameters used in the {\ais}
algorithm, we list the values of the parameters used in our test:
$l=2,500$; $w^{k}=0$ for $k\leq 9$ and $w^{k}=1$ otherwise. We
stopped the updating process in Step~7 of Figure~\ref{fig:ais}
after $k\geq 10$. In other words, we used only the samples
collected in the last step of the algorithm. The learning
parameters used in our algorithm are $k_{\max }=10$, $a=0.4$, and
$b=0.14$ (see Equation~\ref{eq:learning}). We used an empirically
determined value of the threshold $\theta =0.04$
(Section~\ref{sec:heuristic}). We only change the {\cpt} tables of
the parents of a special evidence node $A$ to uniform
distributions when $\Pr (A=a)<1/(2\cdot n_{A})$. Some of the
parameters are a matter of design decision (e.g., the number of
samples in our tests), others were chosen empirically. Although we
have found that these parameters may have different optimal
values for different Bayesian networks, we used the above
values in all our tests of the {\ais} algorithm described in this
paper. Since the same set of parameters led to spectacular
improvement in accuracy in all tested networks, it is fair to say
that the superiority of the {\ais} algorithm to the other
algorithms is not too sensitive to the values of the parameters.

For the {\sis} algorithm, $w^k=1$ by the design of the algorithm.
We used $l=2,500$. The updating function in Step~7 of
Figure~\ref{fig:sis} is that of \cite{shwe91c,cousins93}:
\begin{displaymath}
\Prob_{\mbox{\it new}}^k(x_i|{\Paa}(X_i),{\bf e})
  = \frac{\Prob(x_i|{\Paa}(X_i))
   + k\cdot \widehat{\Prob}_{\mbox{\it current}}
     (x_i|{\Paa}(X_i),{\bf e})} {1+k} \; ,
\end{displaymath}
where $\Prob(x_i|{\Paa}(X_i))$ is the original sampling
distribution, $\widehat{\Prob}_{\mbox{\it
current}}(x_i|{\Paa}(X_i),{\bf e})$ is an equivalent of our
{\pcpt} tables estimator based on all currently available
information, and $k$ is the updating step.


\subsection{Results for the {\cpcs} Network}

The main network used in our tests is a subset of the {\cpcs}
(Computer-based Patient Case Study) model \cite{Pradhan.UAI94}, a
large multiply-connected multi-layer network consisting of 422
multi-valued nodes and covering a subset of the domain of internal
medicine. Among the 422 nodes, 14 nodes describe diseases, 33
nodes describe history and risk factors, and the remaining 375
nodes describe various findings related to the diseases. The
{\cpcs} network is among the largest real networks available to
the research community at the present time. The {\cpcs} network
contains many extreme probabilities, typically on the order of
$10^{-4}$. Our analysis is based on a subset of 179 nodes of the
{\cpcs} network, created by Max Henrion and Malcolm Pradhan. We
used this smaller version in order to be able to compute the exact
solution for the purpose of measuring approximation error in the
sampling algorithms.

The {\ais} algorithm has some learning overhead. The following
comparison of execution time vs.\ number of samples may give the
reader an idea of this overhead. Updating the {\cpcs} network with
20 evidence nodes on our system takes the {\ais} algorithm a total
of 8.4 seconds to learn. It generates subsequently 3,640 samples
per second, while the {\sis} algorithm generates 2,631 samples per
second, and the {\lw} algorithm generates 4,167 samples per
second. In order to remain conservative towards the {\ais}
algorithm, in all our experiments we fixed the execution time of
the algorithms (our limit was 60 seconds) rather than the number
of samples. In the {\cpcs} network with 20 evidence nodes, in 60
seconds, {\ais} generates about 188,000 samples, {\sis} generates
about 158,000 samples and {\lw} generates about 250,000 samples.

\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=evidencedistribute.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{The probability distribution of evidence $\Prob({\bf
E}={\bf e})$ in our experiments.} \label{fig:distribution}
\end{figure}

We generated a total of 75 test cases consisting of five sequences
of 15 test cases each. We ran each test case 10 times, each time
with a different setting of the random number seed. Each sequence
had a progressively higher number of evidence nodes: 15, 20, 25,
30, and 35 evidence nodes respectively. The evidence nodes were
chosen randomly (equiprobable sampling without replacement) from
those nodes that described various plausible medical findings.
Almost all of these nodes were leaf nodes in the network. We
believe that this constituted very realistic test cases for the
algorithms. The distribution of the prior probability of evidence,
$\Prob({\bf E}={\bf e})$, across all test runs of our experiments
is shown in Figure~\ref{fig:distribution}. The least likely
evidence was $5.54\times 10^{-42}$, the most likely evidence was
$1.37\times 10^{-9}$, and the median was $7\times 10^{-24}$.
\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=cpcs20e.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{A typical plot of convergence of the tested sampling
algorithms in our experiments --- Mean Square Error as a function
of the execution time for a subset of the {\cpcs} network with 20
evidence nodes chosen randomly among plausible medical
observations ($\Prob({\bf E}={\bf e})=3.33\times 10^{-26}$ in this
particular case) for the {\ais}, the {\sis}, and the {\lw}
algorithms. The curve for the {\ais} algorithm is very close to
the horizontal axis.} \label{fig:cpcs20e}
\end{figure}
\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=cpcs20ais.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{The lower part of the plot of Figure~\ref{fig:cpcs20e}
showing the convergence of the {\ais} algorithm to correct
posterior probabilities. } \label{fig:cpcs20ais}
\end{figure}

Figures~\ref{fig:cpcs20e} and~\ref{fig:cpcs20ais} show a typical
plot of convergence of the tested sampling algorithms in our
experiments. The case illustrated involves updating the {\cpcs}
network with 20 evidence nodes. We plot the {\mse} after the
initial 15 seconds during which the algorithms start converging.
In particular, the learning step of the {\ais} algorithm is
usually completed within the first 9 seconds. We ran the three
algorithms in this case for 150 seconds rather than the 60 seconds
in the actual experiment in order to be able to observe a wider
range of convergence. The plot of the {\mse} for the {\ais}
algorithm almost touches the $X$ axis in Figure~\ref{fig:cpcs20e}.
Figure~\ref{fig:cpcs20ais} shows the same plot in a finer scale in
order to show more detail in the {\ais} convergence curve. It is
clear that the {\ais} algorithm dramatically improves the
convergence rate. We can also see that the results of {\ais}
converge to exact results very fast as the sampling time
increases. In the case captured in Figures~\ref{fig:cpcs20e}
and~\ref {fig:cpcs20ais}, a tenfold increase in the sampling time
(after subtracting the overhead for the {\ais} algorithm, it
corresponds to a 21.5-fold increase in the number of samples)
results in a 4.55-fold decrease of the {\mse} (to {\mse}$\leq
0.00048$). The observed convergence of both {\sis} and {\lw}
algorithms was poor. A tenfold increase in sampling time had
practically no effect on accuracy. Please note that this is a very
typical case observed in our experiments.

\begin{figure}[ht]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=learning.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{Convergence of the conditional probabilities during the
example run of the {\ais} algorithm captured in
Figure~\ref{fig:cpcs20ais}. The displayed fragment of the
conditional probability table belongs to node {\em gasAcute} which
is a parent of one of the evidence nodes.} \label{fig:learning}
\end{figure}
\begin{table}[tbh]
\centering
\begin{tabular}{|l|c|c|c|}
\hline & Original {\cpt} & Exact {\pcpt} & Learned {\pcpt} \\
\hline ``Absent'' & 0.99631 & 0.0037 & 0.015 \\ ``Mild'' & 0.00183
& 0.1560 & 0.164 \\ ``Moderate'' & 0.00093 & 0.1190 & 0.131 \\
``Severe'' & 0.00093 & 0.7213 & 0.690 \\ \hline
\end{tabular}
\caption{A fragment of the conditional probability table of a node
of the {\cpcs} network (node {\em gasAcute}, parents {\em
hepAcute}={\em Mild} and {\em wbcTotTho}={\em False}) in
Figure~\ref{fig:cpcs20ais}.} \label{tab:learn}
\end{table}

Figure~\ref{fig:learning} illustrates the {\pcpt} learning process
of the {\ais} algorithm for the sample case shown in Figure~\ref
{fig:cpcs20ais}. The displayed conditional probabilities belong to
the node {\em gasAcute} which is a parent of two evidence nodes,
{\em difInfGasMuc} and {\em abdPaiExaMea}. The node {\em gasAcute}
has four states: ``absent,'' ``mild,'' ``moderate,'' and
``severe'', and two parents. We randomly chose a combination of
its parents' states as our displayed configuration. The original
{\cpt} for this configuration without evidence, the exact {\pcpt}
with evidence and the learned {\pcpt} with evidence are summarized
numerically in Table~\ref{tab:learn}. Figure~\ref{fig:learning}
illustrates that the learned importance conditional probabilities
begin to converge to the exact results stably after three updating
steps. The learned probabilities in Step~10 are close to the exact
results. In this example, the difference between
$\Prob(x_i|{\Paa}(X_i),{\bf e})$ and $\Prob(x_i|{\Paa}(X_i))$ is
very large. Sampling from $\Prob(x_i|{\Paa}(X_i))$ instead of
$\Prob(x_i|{ \Paa}(X_i),{\bf e})$ would introduce large variance
into our results.

\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=indivprob.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{Performance of the {\ais}, {\sis}, and {\lw} algorithms:
Mean Square Error for each of the 75 individual test cases plotted
against the probability of evidence. The sampling time is 60 seconds.}
\label{fig:individual}
\end{figure}
\begin{table}[tbh]
\centering
\begin{tabular}{|r|c|c|c|}
\hline & {\ais} & {\sis} & {\lw} \\ \hline $\mu$ & 0.00082 & 0.110
& 0.148 \\ $\sigma$ & 0.00022 & 0.076 & 0.093 \\ $\min$ & 0.00049
& 0.0016 & 0.0031 \\ $\mbox{median}$ & 0.00078 & 0.105 & 0.154 \\
$\max$ & 0.00184 & 0.316 & 0.343 \\ \hline
\end{tabular}
\caption{Summary of the simulation results for all of the 75
simulation cases on the {\cpcs} network.
Figure~\ref{fig:individual} shows each of the 75 cases
graphically.} \label{tab:sum}
\end{table}

Figure~\ref{fig:individual} shows the {\mse} for all 75 test cases
in our experiments with the summary statistics in
Table~\ref{tab:sum}. A paired one-tailed $t$-test resulted in
statistically highly significant differences between the {\ais}
and {\sis} algorithms ($p<3.1\times 10^{-20}$), and also between
the {\sis} and {\lw} algorithms ($p<1.7\times 10^{-8}$). As far as
the magnitude of difference is concerned, {\ais} was significantly
better than {\sis}. {\sis} was better than {\lw}, but the
difference was small. The mean {\mse}s of {\sis} and {\lw}
algorithms were both greater than $0.1$, which suggests that
neither of these algorithms is suitable for large Bayesian
networks.

\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=ratio.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{The ratio of {\mse} between {\sis} and {\ais} versus
percentage.} \label{fig:ratio}
\end{figure}

The graph in Figure~\ref{fig:ratio} shows the {\mse} ratio between
the {\ais} and {\sis} algorithms. We can see that the percentage
of the cases whose ratio was greater than 100 (two orders of
magnitude improvement!) is 60\%. In other words, we obtained two
orders of magnitude improvement in {\mse} in more than half of the
cases. In 80\% cases, the ratio was greater than 50. The smallest
ratio in our experiments was 2.67, which happened when posterior
probabilities were dominated by the prior probabilities. In that
case, even though the {\lw} and {\sis} algorithms converged very
fast, their {\mse} was still far larger than that of {\ais}.

Our next experiment aimed at showing how close the {\ais}
algorithm can approach the best possible sampling results. If we
know the optimal importance sampling function, the convergence of
the {\ais} algorithm should be the same as that of forward
sampling without evidence. In other words, the results of the
probabilistic logic sampling algorithm without evidence approach
the limit of how well stochastic sampling can perform. We ran the
logic sampling algorithm on the {\cpcs} network without evidence
mimicking the test runs of the {\ais} algorithm, i.e., 5 blocks of
15 runs, each repeated 10 times with a different random number
seed. The number of samples generated was equal to the average
number of samples generated by the {\ais} algorithm for each
series of 15 test runs. We obtained the average {\mse} $\mu
=0.00057$, with $\sigma =0.000025$, $\min =0.00052$, and $\max
=0.00065$. The best results should be around this range. From
Table~\ref{tab:sum}, we can see that the minimum {\mse} for the
{\ais} algorithm was $0.00049$, within the range of the optimal
result. The mean {\mse} in {\ais} is $0.00082$, not too far from
the optimal results. The standard deviation, $\sigma $, is
significantly larger in the {\ais} algorithm, but this is
understandable given that the process of learning the optimal
importance function is heuristic in nature. It is not difficult to
understand that there exist a difference between the {\ais}
results and the optimal results. First, the {\ais} algorithm in
our tests updated the sampling distribution only 10 times, which
may be too few times to let it converge to the optimal importance
distribution. Second, even if the algorithm has converged to the
optimal importance distribution, the sampling algorithm will still
let the parameter oscillate around this distribution and there
will be always small differences between the two distributions.

Figure~\ref{fig:convergence} shows the convergence rate for all
tested cases for a four-fold increase in sampling time (between 15
and 60 seconds).
\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=converge.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{The distribution of the convergence ratio of the {\ais},
{\sis}, and {\lw} algorithms when the number of samples increases
four times.} \label{fig:convergence}
\end{figure}
We adjusted the convergence ratio of the {\ais} algorithm by
dividing it by a constant. According to
Equation~\ref{eq:variance}, the theoretically expected convergence
ratio for a four-fold increase in the number of samples should be
around two. There are about 96\% cases among the {\ais} runs whose
ratio lays in the interval (1.75, 2.25], in a sharp contrast to
11\% and 13\% cases in the {\sis} and {\lw} algorithms. The ratios
of the remaining 4\% cases in {\ais} lay in the interval [2.25,
2.5]. In the {\sis} and {\lw} algorithms, the percentage of cases
whose ratio were smaller than 1.5 was 71\% and 77\% respectively.
Less than 1.5 means that the number of samples was too small to
estimate variance and the results cannot be trusted. The ratio
greater than 2.25 means possibly that 60 seconds was long enough
to estimate the variance, but 15 seconds was too short.


\subsection{The Role of {\ais} Heuristics in Performance Improvement}

From the above experimental results we can see that the {\ais}
algorithm can improve the sampling performance significantly. Our
next series of tests focused on studying the role of the two
{\ais} initialization heuristics. The first is initializing the
{\pcpt} tables of the parents of evidence to uniform
distributions, denoted by {\em U}. The second is adjusting small
probabilities, denoted by {\em S}. We denote {\ais} without any
heuristic initialization method to be the {\sc AIS} algorithm.
{\sc AIS}+{\em U}+{\em S} equals {\ais}. We compared the following
versions of the algorithms: {\sis}, {\sc AIS}, {\sis}+{\em U},
{\sc AIS}+{\em U}, {\sis}+{\em S}, {\sc AIS}+{\em S}, {\sis}+{\em
U}+{\em S}, {\sc AIS}+{\em U}+{\em S}. All algorithms with {\sis}
used the same number of samples as {\sis}. All algorithms with
{\sc AIS} used the same number of samples as {\ais}. We tested
these algorithms on the same 75 test cases used in the previous
experiment. Figure~\ref{fig:compare} shows the {\mse} for each of
the sampling algorithms with the summary statistics in
Table~\ref{tab:sumCompare}. Even though the {\sc AIS} algorithm is
better than the {\sis} algorithm, the difference is not as large
as in case of the {\sc AIS}+{\em U}, {\sc AIS}+{\em S}, and {\ais}
algorithms. It seems that heuristic initialization methods help
much. The results for the {\sis}+{\em S}$,$ {\sis}+{\em U},
{\sis}+{\em U}+{\em S} algorithms suggest that although heuristic
initialization methods can improve performance, they alone cannot
improve too much. It is fair to say that significant performance
improvement in the {\ais} algorithm is coming from the combination
of {\sc AIS} with heuristic methods, not any method alone. It is
not difficult to understand that, as only with good heuristic
initialization methods is it possible to let the learning process
quickly exit oscillation areas. Although both {\em S} and {\em U}
methods alone can improve the performance, the improvement is
moderate compared to the combination of the two.
\begin{figure}[tbh]
\begin{center}
\begin{minipage}{7cm}
\psfig{file=compare.eps,width=9.5cm}
\end{minipage}
\end{center}
\caption{A comparison of different algorithms in the {\cpcs}
network. Each bar is based on 75 test cases. The dotted bar shows
the {\mse} for the {\sis} algorithm while the gray bar shows the
{\mse} for the {\sc AIS} algorithm. } \label{fig:compare}
\end{figure}
\begin{table}[tbh]
\centering
\begin{tabular}{|r|c|c|c|c|c|c|c|c|}
\hline & SIS & AIS & SIS+U & AIS+U & SIS+S & AIS+S & SIS+U+S &
{\ais}
\\
 \hline $\mu$ & 0.110 & 0.060 & 0.050 & 0.0084 & 0.075
   & 0.0015 & 0.050 & 0.00082 \\
 $\sigma$ & 0.076 & 0.049 & 0.052 & 0.025 & 0.074
   & 0.0016 & 0.059 & 0.00022 \\
 $\min$ & 0.0016 & 0.00074 & 0.0011 & 0.00058 & 0.00072
   & 0.00056 & 0.00086 & 0.00049 \\
 $\mbox{median}$ & 0.105 & 0.045 & 0.031 & 0.0014 & 0.052
   & 0.00087 & 0.028 & 0.00078 \\
 $\max$ & 0.316 & 0.207 & 0.212 & 0.208 & 0.279 & 0.0085
   & 0.265 & 0.0018 \\ \hline
\end{tabular}
\caption{Summary of the simulation results for different
algorithms in the {\cpcs} network.} \label{tab:sumCompare}
\end{table}

\subsection{Results for Other Networks}

In order to make sure that the {\ais} algorithm performs well in
general, we tested it on two other large networks.

The first network that we used in our tests is the {\pathfinder}
network \cite{heckerman90a}, which is the core element of an
expert system that assists surgical pathologists with the
diagnosis of lymph-node diseases. There are two versions of this
network. We used the larger version, consisting of 135 nodes. In
contrast to the {\cpcs} network, {\pathfinder} contains many
conditional probabilities that are equal to 1, which reflects
deterministic relationships in certain settings. To make the
sampling challenging, we randomly selected 20 evidence nodes from
among the leaf nodes. Each of these was an observable node (David
Heckerman, personal communication). We verified in each case that
the probability of so selected evidence was not equal to zero.

We fixed the execution time of the algorithms to be 60 seconds.
The learning overhead for the {\ais} algorithm in the
{\pathfinder} network was about 3.5 seconds. In 60 seconds, {\ais}
generated about 366,000 samples, {\sis} generated about 250,000
samples and {\lw} generated about 2,700,000 samples. The reason
why {\lw} could generate more than 10 times as many samples as
{\sis} within the same amount of time is that the {\lw} algorithm
terminates sample generation at a very early stage in many
samples, when the weight of a sample becomes zero. This is a
result of determinism in the probability tables, mentioned above.
We will see that {\lw} benefits greatly from generating more
samples. The other parameters used in {\ais} were the same as
those used in the {\cpcs} network.

We tested 20 cases, each with randomly selected 20 evidence nodes. The
reported {\mse} for each case was averaged over 10 runs. Some of the
runs of the {\sis} and {\lw} algorithms did not
manage to generate any effective samples (the weight score sum was
equal to zero). {\sis} had only 75\% effective runs and {\lw} had
only 89\% effective runs, which means that in some runs {\sis} and
{\lw} were unable to yield any information about the posterior
distributions. In all those cases, we discarded the run and only
averaged over the effective runs. All runs in the {\ais} algorithm
were effective. We report our experimental
results with the summary statistics in Table~\ref{tab:pathfinder}.
From these data, we can see that the {\ais} algorithm is still
significantly better than the {\sis} and {\lw} algorithms. Since
the {\lw} algorithm can generate more than ten times the number of
samples than the {\sis} algorithm, its performance is better than
that of the {\sis} algorithm.
\begin{table}[tbh]
\centering
\begin{tabular}{|r|c|c|c|}
 \hline & {\ais} & {\sis} & {\lw} \\
 \hline $\mu$ & 0.00050 & 0.166 & 0.089 \\
 $\sigma$ & 0.00037 & 0.107 & 0.0707 \\
 $\min$ & 0.00025 & 0.00116 & 0.00080 \\
 $\mbox{median}$ & 0.00037 & 0.184 & 0.0866 \\
 $\max$ & 0.0017 & 0.467 & 0.294 \\
 effective runs & 200 & 150 & 178 \\
 \hline
\end{tabular}
\caption{Summary of the simulation results for all of the 20
simulation cases on the {\pathfinder} network. }
\label{tab:pathfinder}
\end{table}

The second network that we tested was one of the {\andes} networks
\cite{conati97}. {\andes} is an intelligent tutoring system for
classical Newtonian physics that is being developed by a team of
researchers at the Learning Research and Development Center at the
University of Pittsburgh and researchers at the United States
Naval Academy. The student model in {\andes} uses a Bayesian
network to do long--term knowledge assessment, plan recognition,
and prediction of students' actions during problem solving. We
selected the largest {\andes} network that was available to us,
consisting of 223 nodes.

In contrast to the previous two networks, the depth of the
{\andes} network was significantly larger and so was its
connectivity. There were only 22 leaf nodes. It is quite
predictable that this kind of networks will pose difficulties to
learning. We selected 20 evidence nodes randomly from the
potential evidence nodes and tested 20 cases. All parameters were
the same as those used in the {\cpcs} network. We fixed the
execution time of the algorithms to be 60 seconds. The learning
overhead for the {\ais} algorithm in the {\andes} network was 13.4
seconds. In 60 seconds, {\ais} generated about 114,000 samples,
{\sis} generated about 98,000 samples and {\lw} generated about
180,000 samples. In this network, {\lw} still can generate almost
two times the number of samples generated by the {\sis} algorithm.

We report our experimental results with the summary statistics in
Table~\ref{tab:andes}. The results show that also in the {\andes}
network the {\ais} algorithm was significantly better than the
{\sis} and {\lw} algorithms. Since {\lw} generated almost two
times the number of samples that were generated by the {\sis}
algorithm, its performance was better than that of the {\sis}
algorithm.

\begin{table}[bbh]
\centering
\begin{tabular}{|r|c|c|c|}
\hline & {\ais} & {\sis} & {\lw} \\ \hline $\mu$ & 0.0059 & 0.0628
& 0.0404 \\ $\sigma$ & 0.0049 & 0.102 & 0.0539 \\ $\min$ & 0.0023
& 0.0028 & 0.0028 \\ $\mbox{median}$ & 0.0045 & 0.0190 & 0.0198 \\
$\max$ & 0.0237 & 0.321 & 0.221 \\ \hline
\end{tabular}
\caption{Summary of the simulation results for all of the 20
simulation cases on the {\andes} network. } \label{tab:andes}
\end{table}

While the {\ais} algorithm is on the average an order of magnitude
more precise than the other two algorithms, the performance
improvement is smaller than in the other two networks. The reason
why the performance improvement of the {\ais} algorithm over the
{\sis} and {\lw} algorithms in the {\andes} network is smaller
compared to that in the {\cpcs} and {\pathfinder} networks is
that: (1) The {\andes} network used in our tests was apparently
not challenging enough for sampling algorithms in general. In the
{\andes} network, {\sis} and {\lw} also can perform well in some
cases. The minimum {\mse} of {\sis} and {\lw} in our tested cases
is almost the same as that of {\ais}. (2) The number of samples
generated by {\ais} in this network is significantly smaller than
that in the previous two networks and {\ais} needs more time to
learn. Although increasing the number of samples will improve the
performance of all three algorithms, it improves the performance
of {\ais} more since the convergence ratio of the {\ais} algorithm
is usually larger than that of {\sis} and {\lw} (see
Figure~\ref{fig:convergence}). (3) The parameters that we used in
this network were tuned for the {\cpcs} network. (4) The large
depth and fewer leaf nodes of the {\andes} network pose some
difficulties to learning.

\section{Discussion}
\label{sec:discussion}

There is a fundamental trade-off in the {\ais} algorithm between
the time spent on learning the importance function and the time
spent on sampling. Our current approach, which we believe to be
reasonable, is to stop learning at the point when the importance
function is good enough. In our experiments we stopped learning
after 10 iterations.

There are several ways of improving the initialization of the
conditional probability tables at the outset of the {\ais}
algorithm. In the current version of the algorithm, we initialize
the {\pcpt} table of every parent $N$ of an evidence node $E$
($N\in {\Pa}(E)$, $E\in {\bf E}$) to the uniform distribution when
$\Pr (E=e)<1/(2\cdot n_{E})$. This can be improved further. We can
extend the initialization to those nodes that are severely
affected by the evidence. They can be identified by examining the
network structure and local {\cpt}s.

We can view the learning process of the {\ais} algorithm as a
network rebuilding process. The algorithm constructs a new network
whose structure is the same as the original network (except that
we delete the evidence nodes and corresponding arcs). The
constructed network models the joint probability distribution
$\rho({\bf X}\backslash{\bf E})$ in Equation~\ref{eq:ejoint},
which approaches the optimal importance function. We use the
learned $\rho^\prime$ to approximate this distribution. If
$\rho^\prime$ approximates $\Prob({\bf X}|{\bf E})$ accurately
enough, we can use this new network to solve other approximate
tasks, such as the problem of computing the Maximum A-Posterior
assignment (MAP) \cite{pearl88}, finding $k$ most likely scenarios
\cite{seroussi94}, etc. A large advantage of this approach is that
we can solve each of these problems as if the network had no
evidence nodes.

We know that Markov blanket scoring can improve convergence rates
in some sampling algorithms \cite{shwe91b}. It may also be applied
to the {\ais} algorithm to improve its convergence rate. According
to Property~4 (Section~\ref{sec:theory}), any technique that can
reduce the variance $\sigma_{\Prob({\bf e})}^2$ will reduce the
variance of $\widehat{\Prob}({\bf e})$ and correspondingly improve
the sampling performance. Since the variance of stratified
sampling \cite{rubinstein81} is never much worse than that
of random sampling, and can be much better, it can improve the
convergence rate. We expect some other variance reduction methods
in statistics, such as: ({\it i}) the expected value of a random
variable; ({\it ii}) antithetic variants correlations (stratified
sampling, Latin hypercube sampling, etc.); and ({\it iii})
systematic sampling, will also improve the sampling performance.

Current learning algorithm used a simple approach. Some heuristic
learning methods, such as adjusting learning rates according to
changes of the error \cite{jacobs88}, should also be applicable to
our algorithm. There are several tunable parameters in the {\ais}
algorithm. Finding the optimal values of these parameters for any
given network is another interesting research topic.

It is worth observing that the plots presented in
Figure~\ref{fig:individual} are fairly flat. In other words, in
our tests the convergence of the sampling algorithms did not
depend too strongly on the probability of evidence. This seems to
contradict the common belief that forward sampling schemes suffer
from unlikely evidence. {\ais} for one shows a fairly flat plot.
The convergence of the {\sis} and {\lw} algorithms seems to
decrease slightly with unlikely evidence. It is possible that all
three algorithms will perform much worse when the probability of
evidence drops below some threshold value, which our tests failed
to approach. Until this relationship has been studied carefully,
we conjecture that the probability of evidence is not a good
measure of difficulty of approximate inference.

Given that the problem of approximating probabilistic inference is
NP-hard, there exist networks that will be challenging for any
algorithm and we have no doubt that even the {\ais} algorithm will
perform poorly on them. To the day, we have not found such
networks. There is one characteristic of networks that may be
challenging to the {\ais} algorithm. In general, when the number
of parameters that need to be learned by the {\ais} algorithm
increases, its performance will deteriorate. Nodes with many
parents, for example, are challenging to the {\ais} learning
algorithm, as it has to update the {\pcpt} tables under all
combinations of the parent nodes. It is possible that conditional
probability distributions with causal independence properties,
such as Noisy-OR distributions
\cite{pearl88,henrion89a,Diez.UAI93,Srinivas.UAI93,heckerman94},
common in very large practical networks, can be treated
differently and lead to considerable savings in the learning time.

One direction of testing approximate algorithms, suggested to us
by a reviewer, is to use very large networks for which exact
solution cannot be computed at all. In this case, one can try to
infer from the difference in variance at various stages of the
algorithm whether it is converging or not. This is a very
interesting idea that is worth exploring, especially when combined
with theoretical work on stopping criteria in the line of the work
of Dagum and Luby \citeyear{dagum97}.


\section{Conclusion}
\label{sec:conclusion}

Computational complexity remains a major problem in application of
probability theory and decision theory in knowledge-based systems.
It is important to develop schemes that improve the performance of
updating algorithms --- even though the theoretically demonstrated
worst case will remain NP--hard, many practical cases may become
tractable.

In this paper, we studied importance sampling in Bayesian
networks. After reviewing the most important theoretical results
related to importance sampling in finite-dimensional integrals, we
proposed a new algorithm for importance sampling in Bayesian
networks that we call {\em adaptive importance sampling} ({\ais}).
While the process of learning the optimal importance function for
the {\ais} algorithm is computationally intractable, based on the
theory of importance sampling in finite-dimensional integrals we
proposed several heuristics that seem to work very well in
practice. We proposed heuristic methods for initializing the
importance function that we have shown to accelerate the learning
process, a smooth learning method for updating importance function
using the structural advantages of Bayesian networks, and a
dynamic weighting function for combining samples from different
stages of the algorithm. All these methods help the {\ais}
algorithm to get fairly accurate estimates of the posterior
probabilities in a limited time. Of the two applied heuristics,
adjustment of small probabilities, seems to lead to the largest
improvement in performance, although the largest decrease in
{\mse} is achieved by a combination of the two heuristics with the
{\ais} algorithm.

The {\ais} algorithm can lead to a dramatic improvement in the
convergence rates in large Bayesian networks with evidence
compared to the existing state of the art algorithms. We compared
the performance of the {\ais} algorithm to the performance of
likelihood weighting and self-importance sampling on a large
practical model, the {\cpcs} network, with evidence as unlikely as
$5.54\times 10^{-42}$ and typically $7\times 1.0^{-24}$. In our
experiments, we observed that the {\ais} algorithm was always
better than likelihood weighting and self-importance sampling and
in over $60\%$ of the cases it reached over two orders of
magnitude improvement in accuracy. Tests performed on the other
two networks, {\pathfinder} and {\andes}, yielded similar results.

Although there may exist other approximate algorithms that will
prove superior to {\ais} in networks with special structure or
distribution, the {\ais} algorithm is simple and robust for
general evidential reasoning problems in large multiply-connected
Bayesian networks.

\acks{We thank anonymous referees for several insightful comments
that led to a substantial improvement of the paper. This research
was supported by the National Science Foundation under Faculty
Early Career Development (CAREER) Program, grant IRI--9624629, and
by the Air Force Office of Scientific Research grants
F49620--97--1--0225 and F49620--00--1--0112. An earlier version of
this paper has received the 2000 School of Information Sciences
Robert R. Korfhage Award, University of Pittsburgh. Malcolm
Pradhan and Max Henrion of the Institute for Decision Systems
Research shared with us the {\cpcs} network with a kind permission
from the developers of the {\sc Internist} system at the
University of Pittsburgh. We thank David Heckerman for the
{\pathfinder} network and Abigail Gerner for the {\andes} network
used in our tests. All experimental data have been obtained using
{\sc SMILE}, a Bayesian inference engine developed at the Decision
Systems Laboratory and available at {\tt
http://www2.sis.pitt.edu/$\sim$genie.}}


\bibliographystyle{theapa}

\begin{thebibliography}{}

\bibitem[\protect\BCAY{Cano, Hernandez, \BBA\ Moral}{Cano
  et~al.}{1996}]{cano96}
Cano, J.~E., Hernandez, L.~D., \BBA\ Moral, S. \BBOP1996\BBCP.
\newblock \BBOQ Importance sampling algorithms for the propagation of
  probabilities in belief networks\BBCQ\
\newblock {\Bem International Journal of Approximate Reasoning}, {\Bem 15},
  77--92.

\bibitem[\protect\BCAY{Chavez \BBA\ Cooper}{Chavez \BBA\
  Cooper}{1990}]{chavez90}
Chavez, M.~R.\BBACOMMA\  \BBA\ Cooper, G.~F. \BBOP1990\BBCP.
\newblock \BBOQ A randomized approximation algorithm for probabilistic
  inference on {B}ayesian belief networks\BBCQ\
\newblock {\Bem Networks}, {\Bem 20\/}(5), 661--685.

\bibitem[\protect\BCAY{Cheng \BBA\ Druzdzel}{Cheng \BBA\
  Druzdzel}{2000a}]{Cheng.UAI00}
Cheng, J.\BBACOMMA\  \BBA\ Druzdzel, M.~J. \BBOP2000a\BBCP.
\newblock \BBOQ Computational investigations of low�-discrepancy sequences in
  simulation algorithms for {B}ayesian networks\BBCQ\
\newblock In {\Bem Proceedings of the Sixteenth Annual Conference on
  Uncertainty in Artificial Intelligence (UAI--2000)}, \BPGS\ 72--81\ San
  Francisco, CA. Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{Cheng \BBA\ Druzdzel}{Cheng \BBA\
  Druzdzel}{2000b}]{cheng00a}
Cheng, J.\BBACOMMA\  \BBA\ Druzdzel, M.~J. \BBOP2000b\BBCP.
\newblock \BBOQ Latin hypercube sampling in {B}ayesian networks\BBCQ\
\newblock In {\Bem Proceedings of the 13th International Florida Artificial
  Intelligence Research Symposium Conference (FLAIRS-2000)}, \BPGS\ 287--292\
  Orlando, Florida.

\bibitem[\protect\BCAY{Conati, Gertner, VanLehn, \BBA\ Druzdzel}{Conati
  et~al.}{1997}]{conati97}
Conati, C., Gertner, A.~S., VanLehn, K., \BBA\ Druzdzel, M.~J.
\BBOP1997\BBCP.
\newblock \BBOQ On-line student modeling for coached problem solving using
  {B}ayesian networks\BBCQ\
\newblock In {\Bem Proceedings of the Sixth International Conference on User
  Modeling (UM--96)}, \BPGS\ 231--242\ Vienna, New York. Springer Verlag.

\bibitem[\protect\BCAY{Cooper}{Cooper}{1990}]{cooper90}
Cooper, G.~F. \BBOP1990\BBCP.
\newblock \BBOQ The computational complexity of probabilistic inference using
  {B}ayesian belief networks\BBCQ\
\newblock {\Bem Artificial Intelligence}, {\Bem 42\/}(2--3), 393--405.

\bibitem[\protect\BCAY{Cousins, Chen, \BBA\ Frisse}{Cousins
  et~al.}{1993}]{cousins93}
Cousins, S.~B., Chen, W., \BBA\ Frisse, M.~E. \BBOP1993\BBCP.
\newblock \BBOQ A tutorial introduction to stochastic simulation algorithm for
  belief networks\BBCQ\
\newblock In {\Bem Artificial Intelligence in Medicine}, \BCH~5, \BPGS\
  315--340. Elsevier Science Publishers B.V.

\bibitem[\protect\BCAY{Dagum, Karp, Luby, \BBA\ Ross}{Dagum
  et~al.}{1995}]{dagum95}
Dagum, P., Karp, R., Luby, M., \BBA\ Ross, S. \BBOP1995\BBCP.
\newblock \BBOQ An optimal algorithm for {M}onte {C}arlo estimation (extended
  abstract)\BBCQ\
\newblock In {\Bem Proceedings of the 36th IEEE Symposium on Foundations of
  Computer Science}, \BPGS\ 142--149\ Portland, Oregon.

\bibitem[\protect\BCAY{Dagum \BBA\ Luby}{Dagum \BBA\ Luby}{1993}]{dagum93}
Dagum, P.\BBACOMMA\  \BBA\ Luby, M. \BBOP1993\BBCP.
\newblock \BBOQ Approximating probabilistic inference in {B}ayesian belief
  networks is {NP}-hard\BBCQ\
\newblock {\Bem Artificial Intelligence}, {\Bem 60\/}(1), 141--153.

\bibitem[\protect\BCAY{Dagum \BBA\ Luby}{Dagum \BBA\ Luby}{1997}]{dagum97}
Dagum, P.\BBACOMMA\  \BBA\ Luby, M. \BBOP1997\BBCP.
\newblock \BBOQ An optimal approximation algorithm for {B}ayesian
  inference\BBCQ\
\newblock {\Bem Artificial Intelligence}, {\Bem 93}, 1--27.

\bibitem[\protect\BCAY{Diez}{Diez}{1993}]{Diez.UAI93}
Diez, F.~J. \BBOP1993\BBCP.
\newblock \BBOQ Parameter adjustment in {B}ayes networks. {T}he generalized
  noisy {OR-gate}\BBCQ\
\newblock In {\Bem Proceedings of the Ninth Annual Conference on Uncertainty in
  Artificial Intelligence (UAI--93)}, \BPGS\ 99--105\ San Francisco, CA. Morgan
  Kaufmann Publishers.

\bibitem[\protect\BCAY{Fishman}{Fishman}{1995}]{fishman95}
Fishman, G.~S. \BBOP1995\BBCP.
\newblock {\Bem {M}onte {C}arlo: concepts, algorithms, and applications}.
\newblock Springer-Verlag.

\bibitem[\protect\BCAY{Fung \BBA\ Chang}{Fung \BBA\ Chang}{1989}]{Fung.UAI89}
Fung, R.\BBACOMMA\  \BBA\ Chang, K.-C. \BBOP1989\BBCP.
\newblock \BBOQ Weighing and integrating evidence for stochastic simulation in
  {B}ayesian networks\BBCQ\
\newblock In {\Bem Uncertainty in Artificial Intelligence~5}, \BPGS\ 209--219\
  New York, N. Y. Elsevier Science Publishing Company, Inc.

\bibitem[\protect\BCAY{Fung \BBA\ del Favero}{Fung \BBA\ del
  Favero}{1994}]{Fung.UAI94}
Fung, R.\BBACOMMA\  \BBA\ del Favero, B. \BBOP1994\BBCP.
\newblock \BBOQ Backward simulation in {B}ayesian networks\BBCQ\
\newblock In {\Bem Proceedings of the Tenth Annual Conference on Uncertainty in
  Artificial Intelligence (UAI--94)}, \BPGS\ 227--234\ San Francisco, CA.
  Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{Geman \BBA\ Geman}{Geman \BBA\ Geman}{1984}]{geman84}
Geman, S.\BBACOMMA\  \BBA\ Geman, D. \BBOP1984\BBCP.
\newblock \BBOQ Stochastic relaxations, {G}ibbs distributions and the
  {B}ayesian restoration of images\BBCQ\
\newblock {\Bem IEEE Transactions on Pattern Analysis and Machine
  Intelligence}, {\Bem 6\/}(6), 721--742.

\bibitem[\protect\BCAY{Gilks, Richardson, \BBA\ Spiegelhalter}{Gilks
  et~al.}{1996}]{gilks96}
Gilks, W., Richardson, S., \BBA\ Spiegelhalter, D. \BBOP1996\BBCP.
\newblock {\Bem Markov chain {M}onte {C}arlo in practice}.
\newblock Chapman and Hall.

\bibitem[\protect\BCAY{Heckerman \BBA\ Breese}{Heckerman \BBA\
  Breese}{1994}]{heckerman94}
Heckerman, D.\BBACOMMA\  \BBA\ Breese, J.~S. \BBOP1994\BBCP.
\newblock \BBOQ A new look at causal independence\BBCQ\
\newblock In {\Bem Proceedings of the Tenth Annual Conference on Uncertainty in
  Artificial Intelligence (UAI--94)}, \BPGS\ 286--292\ San Mateo, CA. Morgan
  Kaufmann Publishers, Inc.

\bibitem[\protect\BCAY{Heckerman, Horvitz, \BBA\ Nathwani}{Heckerman
  et~al.}{1990}]{heckerman90a}
Heckerman, D.~E., Horvitz, E.~J., \BBA\ Nathwani, B.~N.
\BBOP1990\BBCP.
\newblock \BBOQ Toward normative expert systems: The {P}athfinder project\BBCQ\
\newblock \BTR\ KSL--90--08, Medical Computer Science Group, Section on Medical
  Informatics, Stanford University, Stanford, CA.

\bibitem[\protect\BCAY{Henrion}{Henrion}{1988}]{Henrion.UAI86}
Henrion, M. \BBOP1988\BBCP.
\newblock \BBOQ Propagating uncertainty in {B}ayesian networks by probabilistic
  logic sampling\BBCQ\
\newblock In {\Bem Uncertainty in Artificial Intellgience~2}, \BPGS\ 149--163\
  New York, N. Y. Elsevier Science Publishing Company, Inc.

\bibitem[\protect\BCAY{Henrion}{Henrion}{1989}]{henrion89a}
Henrion, M. \BBOP1989\BBCP.
\newblock \BBOQ Some practical issues in constructing belief networks\BBCQ\
\newblock In Kanal, L., Levitt, T., \BBA\ Lemmer, J.\BEDS, {\Bem Uncertainty in
  Artificial Intelligence~3}, \BPGS\ 161--173. Elsevier Science Publishers
  B.V., North Holland.

\bibitem[\protect\BCAY{Henrion}{Henrion}{1991}]{Henrion.UAI91}
Henrion, M. \BBOP1991\BBCP.
\newblock \BBOQ Search-based methods to bound diagnostic probabilities in very
  large belief nets\BBCQ\
\newblock In {\Bem Proceedings of the Seventh Annual Conference on Uncertainty
  in Artificial Intelligence (UAI--91)}, \BPGS\ 142--150\ San Mateo,
  California. Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{Hernandez, Moral, \BBA\ Antonio}{Hernandez
  et~al.}{1998}]{hernandez98}
Hernandez, L.~D., Moral, S., \BBA\ Antonio, S. \BBOP1998\BBCP.
\newblock \BBOQ A {M}onte {C}arlo algorithm for probabilistic propagation in
  belief networks based on importance sampling and stratified simulation
  techniques\BBCQ\
\newblock {\Bem International Journal of Approximate Reasoning}, {\Bem 18},
  53--91.

\bibitem[\protect\BCAY{Jacobs}{Jacobs}{1988}]{jacobs88}
Jacobs, R.~A. \BBOP1988\BBCP.
\newblock \BBOQ Increased rates of convergence through learning rate
  adaptation\BBCQ\
\newblock {\Bem Neural Networks}, {\Bem 1}, 295--307.

\bibitem[\protect\BCAY{Lauritzen \BBA\ Spiegelhalter}{Lauritzen \BBA\
  Spiegelhalter}{1988}]{lauritzen88}
Lauritzen, S.~L.\BBACOMMA\  \BBA\ Spiegelhalter, D.~J.
\BBOP1988\BBCP.
\newblock \BBOQ Local computations with probabilities on graphical structures
  and their application to expert systems\BBCQ\
\newblock {\Bem Journal of the Royal Statistical Society, Series B
  (Methodological)}, {\Bem 50\/}(2), 157--224.

\bibitem[\protect\BCAY{MacKay}{MacKay}{1998}]{mackay98}
MacKay, D. \BBOP1998\BBCP.
\newblock \BBOQ Intro to {M}onte {C}arlo methods\BBCQ\
\newblock In Jordan, M.~I.\BED, {\Bem Learning in Graphical Models}. The MIT
  Press, Cambridge, Massachusetts.

\bibitem[\protect\BCAY{Ortiz \BBA\ Kaelbling}{Ortiz \BBA\
  Kaelbling}{2000}]{Ortiz.UAI00}
Ortiz, L.~E.\BBACOMMA\  \BBA\ Kaelbling, L.~P. \BBOP2000\BBCP.
\newblock \BBOQ Adaptive importance sampling for estimation in structured
  domains\BBCQ\
\newblock In {\Bem Proceedings of the Sixteenth Annual Conference on
  Uncertainty in Artificial Intelligence (UAI--2000)}, \BPGS\ 446--454\ San
  Francisco, CA. Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{Pearl}{Pearl}{1986}]{pearl86}
Pearl, J. \BBOP1986\BBCP.
\newblock \BBOQ Fusion, propagation, and structuring in belief networks\BBCQ\
\newblock {\Bem Artificial Intelligence}, {\Bem 29\/}(3), 241--288.

\bibitem[\protect\BCAY{Pearl}{Pearl}{1987}]{pearl87}
Pearl, J. \BBOP1987\BBCP.
\newblock \BBOQ Evidential reasoning using stochastic simulation of causal
  models\BBCQ\
\newblock {\Bem Artifical Intelligence}, {\Bem 32}, 245--257.

\bibitem[\protect\BCAY{Pearl}{Pearl}{1988}]{pearl88}
Pearl, J. \BBOP1988\BBCP.
\newblock {\Bem Probabilistic Reasoning in Intelligent Systems: Networks of
  Plausible Inference}.
\newblock Morgan Kaufmann Publishers, Inc., San Mateo, CA.

\bibitem[\protect\BCAY{Pradhan \BBA\ Dagum}{Pradhan \BBA\
  Dagum}{1996}]{Pradhan.UAI96}
Pradhan, M.\BBACOMMA\  \BBA\ Dagum, P. \BBOP1996\BBCP.
\newblock \BBOQ Optimal {M}onte {C}arlo inference\BBCQ\
\newblock In {\Bem Proceedings of the Twelfth Annual Conference on Uncertainty
  in Artificial Intelligence (UAI--96)}, \BPGS\ 446--453\ San Francisco, CA.
  Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{Pradhan, Provan, Middleton, \BBA\ Henrion}{Pradhan
  et~al.}{1994}]{Pradhan.UAI94}
Pradhan, M., Provan, G., Middleton, B., \BBA\ Henrion, M.
\BBOP1994\BBCP.
\newblock \BBOQ Knowledge engineering for large belief networks\BBCQ\
\newblock In {\Bem Proceedings of the Tenth Annual Conference on Uncertainty in
  Artificial Intelligence (UAI--94)}, \BPGS\ 484--490\ San Francisco, CA.
  Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{Ritter, Martinetz, \BBA\ Schulten}{Ritter
  et~al.}{1991}]{ritter91}
Ritter, H., Martinetz, T., \BBA\ Schulten, K. \BBOP1991\BBCP.
\newblock {\Bem Neuronale Netze}.
\newblock Addison-Wesley, M\"unchen.

\bibitem[\protect\BCAY{Rubinstein}{Rubinstein}{1981}]{rubinstein81}
Rubinstein, R.~Y. \BBOP1981\BBCP.
\newblock {\Bem Simulation and the {M}onte {C}arlo Method}.
\newblock John Wiley \& Sons.

\bibitem[\protect\BCAY{Seroussi \BBA\ Golmard}{Seroussi \BBA\
  Golmard}{1994}]{seroussi94}
Seroussi, B.\BBACOMMA\  \BBA\ Golmard, J.~L. \BBOP1994\BBCP.
\newblock \BBOQ An algorithm directly finding the {K} most probable
  configurations in {B}ayesian networks\BBCQ\
\newblock {\Bem International Journal of Approximate Reasoning}, {\Bem 11},
  205--233.

\bibitem[\protect\BCAY{Shachter \BBA\ Peot}{Shachter \BBA\
  Peot}{1989}]{Shachter2.UAI89}
Shachter, R.~D.\BBACOMMA\  \BBA\ Peot, M.~A. \BBOP1989\BBCP.
\newblock \BBOQ Simulation approaches to general probabilistic inference on
  belief networks\BBCQ\
\newblock In {\Bem Uncertainty in Artificial Intelligence~5}, \BPGS\ 221--231\
  New York, N. Y. Elsevier Science Publishing Company, Inc.

\bibitem[\protect\BCAY{Shwe \BBA\ Cooper}{Shwe \BBA\ Cooper}{1991}]{shwe91b}
Shwe, M.~A.\BBACOMMA\  \BBA\ Cooper, G.~F. \BBOP1991\BBCP.
\newblock \BBOQ An empirical analysis of likelihood-weighting simulation on a
  large, multiply-connected medical belief network\BBCQ\
\newblock {\Bem Computers and Biomedical Research}, {\Bem 24\/}(5), 453--475.

\bibitem[\protect\BCAY{Shwe, Middleton, Heckerman, Henrion, Horvitz, \BBA\
  Lehmann}{Shwe et~al.}{1991}]{shwe91c}
Shwe, M., Middleton, B., Heckerman, D., Henrion, M., Horvitz, E.,
\BBA\
  Lehmann, H. \BBOP1991\BBCP.
\newblock \BBOQ Probabilistic diagnosis using a reformulation of the
  {INTERNIST--1/QMR} knowledge base: {I}. {T}he probabilistic model and
  inference algorithms\BBCQ\
\newblock {\Bem Methods of Information in Medicine}, {\Bem 30\/}(4), 241--255.

\bibitem[\protect\BCAY{Srinivas}{Srinivas}{1993}]{Srinivas.UAI93}
Srinivas, S. \BBOP1993\BBCP.
\newblock \BBOQ A generalization of the noisy-{OR} model\BBCQ\
\newblock In {\Bem Proceedings of the Ninth Annual Conference on Uncertainty in
  Artificial Intelligence (UAI--93)}, \BPGS\ 208--215\ San Francisco, CA.
  Morgan Kaufmann Publishers.

\bibitem[\protect\BCAY{York}{York}{1992}]{york92}
York, J. \BBOP1992\BBCP.
\newblock \BBOQ Use of the {G}ibbs sampler in expert systems\BBCQ\
\newblock {\Bem Artificial Intelligence}, {\Bem 56}, 115--130.
\end{thebibliography}

\end{document}
