Decision Analysis by Augmented Probability Simulation.
We provide a generic Monte Carlo method to find the alternative of maximum expected utility in a decision analysis. We define an artificial distribution on the product space of alternatives and states, and show that the optimal alternative is the mode of the implied marginal distribution on the alternatives. After drawing a sample from the artificial distribution, we may use exploratory data analysis tools to approximately identify the optimal alternative. We illustrate our method for some important types of influence diagrams.
Introduction
1.1 Decision Analysis by Simulation
Decision Analysis provides a framework for solving decision making problems under uncer-
tainty, based on finding an alternative with maximum expected utility. While conceptually
simple, the actual solution of the maximization problem may be extremely involved, e.g.,
when the probability model is complex, the set of alternatives is continuous, or when a sequence
of decisions is included. Therefore, only particular probability models are studied,
such as the multivariate Gaussian in Shachter and Kenley (1989). Inclusion of continuous
variables in simple problems is carried out through discretization (Miller and Rice 1983,
Smith 1991), through summaries of the first few moments and derivatives (Smith 1993),
or through approximations by means of Gaussian mixtures (Poland 1994). In complicated
problems, there may be no hope for an exact solution method and we may have to turn to
approximate methods, specifically simulation.
As observed in Pearl (1988, p311) and Cooper (1989), in principle any simulation method
to solve Bayesian networks (BN) may be used to solve decision problems represented by influence
diagrams (ID) by means of sequentially instantiating decision nodes and computing
expected values. Cooper notes that, for a given instantiation of the decision nodes, the
computation of the expected value at the value node can be reformulated as a computation
of a posterior distribution in an artificially created additional random node. The problem
of solving BNs is summarized, for example, in Shachter and Peot (1990). Exact algorithms,
e.g. using clique join trees (Lauritzen and Spiegelhalter 1988), cutset conditioning (Pearl
1986) or arc reversal (Shachter 1986, 1988) proved to be intractable in many real-world
networks, leading to approximate inference algorithms based on simulation methods. These
include short run algorithms, such as Logic Sampling (Henrion 1988); Likelihood Weighting
(Shachter and Peot 1990) and its improved modifications, Bounded Variance and AA algorithms
(Pradhan and Dagum 1996); and long run algorithms, using Markov chain Monte
Carlo methods like Gibbs sampling (Pearl 1987, Hrycej 1990, York 1992) or hybrid strategies
(Brewer et al. 1996).
However, as Matzkevich and Abramson (1995) note, we only have a couple of outlines of
simulation methods specifically for IDs in Jenzarli (1995) and Charnes and Shenoy (1996).
Whereas the first one combines stochastic dynamic programming and Gibbs sampling, the
latter simulates iid observations from only a small set of chance variables for each decision
node instead of using the entire distribution. Both become intractable when continuous
decision spaces are included.
In recent statistical literature the same problem, i.e., that of finding the optimal action
in a decision problem, has been considered in M-uller and Parmigiani (1996) and Carlin,
Kadane and Gelfand (1998), among others. Again, all these approaches use Monte Carlo
simulation to evaluate the expected utility of given instantiations of nodes.
1.2 Augmented Probability Simulation
In this paper we propose a scheme which differs in important ways from the above mentioned
approaches. Since they use simulation to evaluate expected utilities (losses) for given
instantiations of the decision nodes, they do not accomodate continuous variables, especially
decision variables, unless a discretization is carried out or the probability distributions are in
a conjugate framework. In contrast, we go a step further and define an artificial distribution
on all nodes, including the decision nodes. We show that simulation from this artificial augmented
probability model is equivalent to solving the original decision problem. The specific
strength of the proposed method is its generality. The algorithm can, in principle, accomodate
arbitrary probability models and utility functions, as long as it is possible to pointwise
evaluate the probability density and the utility function for any chosen value of all involved
nodes. Evaluation of the probability density up to a constant factor suffices. The idea of
augmenting the probability model to transform the optimization problem into a simulation
problem is not entirely new. For example, Shachter and Peot (1992) have proposed a similar
approach which involves augmenting the probability model to include the decision nodes and
thus transforms the original optimization problem into a simulation problem. But to the best
of our knowledge the approach described here is the first to solve this simulation problem by
systematically exploiting Markov chain Monte Carlo simulation methods recently developed
in the statistical literature.
The method starts by considering an artificial distribution on the space of alternatives
and states. The distribution is defined in such a way that its marginal on the space of
alternatives is proportional to the expected utility of the alternative and, consequently, the
optimal alternative coincides with the mode of the marginal. Then, the proposed simulation
based strategy follows these steps: (i) draw a sample from the artificial distribution; (ii)
marginalise it to the space of alternatives; and, (iii) find the mode of the sample as a way
of approximating the optimal alternative. A key issue is how to sample from the artificial
distribution. For that we introduce Markov chain Monte Carlo (MCMC) algorithms. See,
for example, Smith and Roberts (1993), Tierney (1994) or Tanner (1994) for a review of
MCMC methods.
Section 2 describes the basic strategy with a simple example. Section 3 is of a more
technical nature and provides generic methods to sample approximately from the artificial
distribution and identify the mode of the sample. Section 4 discusses application examples.
Section 5 compares our method with alternative schemes and identifies situations which call
for different approaches.
2 Basic Approach
Here we outline the basic approach. Assume we have to choose under uncertainty an alternative
d from a set A. The set of states ' is \Theta. We propose as optimal the alternative
d   with maximum expected utility:
u(d; ')p d (')d'] ; where u(d; ') is the
utility function modeling preferences over consequences and p d (') is the probability distribution
modeling beliefs, possibly influenced by actions. When the problem is structurally
complicated, say a heavily asymmetric and dense, large influence diagram with continuous
non-Gaussian random variables, non quadratic utility functions and/or continuous sets of
alternatives at decision nodes, finding the exact solution might be analytically and computationally
intractable, and we might need an approximate solution method. We shall provide
such an approximation based on simulation.
Assume that p d (') ? 0, for all pairs (d; '), and u(d; ') is positive and integrable. Define an
artificial distribution over the product space A\Theta\Theta with density h proportional to the product
of utility and probability, specifically h(d; ') / u(d; ') \Delta p d ('): Note that the artificial distribution
h is chosen so that the marginal on the alternatives is h(d)
Hence, the optimal alternative d   coincides with the mode of the marginal of the artificial
distribution h in the space of alternatives. As a consequence, we can solve the expected
utility maximization problem approximately with the following simulation based
(i) draw a random sample from the distribution h(d; '); (ii) convert it to a random sample
from the marginal h(d); and (iii) find the mode of this sample.
This augmented probability model simulation is conceptually different from other simulation
algorithms reviewed earlier. Simulation is not used to pointwise evaluate expected
utilities for each decision alternative. Instead, simulation generates the artificial probability
model h(\Delta) on the augmented state vector (d; ').
The key steps are (i) and (iii). For (ii), since we use simulation to generate from h(d; '),
we can get a marginal sample from h(d) by simply discarding the simulated ' values. For
(iii) we rely mainly on tools from exploratory data analysis, as we describe in Section 3.3.
For (i), we shall introduce generic Markov chain simulation methods. Their underlying idea
is simple. We wish to generate a sample from a distribution over a certain space, but cannot
do this directly. Suppose, however, that we can construct a Markov chain with the same
state space, which is straightforward to simulate from and whose equilibrium distribution is
the desired distribution. If we simulate sufficiently many iterations, after dropping an initial
transient phase, we may use the simulated values of the chain as an approximate sample
from the desired distribution. We shall provide several algorithms for constructing chains
with the desired equilibrium distribution, in our case the artificial distribution h, in Section
3.2. In the rest of this section we shall provide an algorithm and a simple example, so that
readers may grasp the basic idea, without entering into technical details. Readers familiar
with MCMC simulation may skip directly to Section 3.
The strategy we propose now is very simple, but may be only undertaken in limited cases.
Suppose the conditional distributions h(dj') and h('jd) are available for efficient random
variate generation. Then, we suggest the following scheme, which is known as the Gibbs
sampler in the statistical literature (Gelfand and Smith, 1990): (i) Start at an arbitrary
value
steps (ii) and (iii) until convergence is judged.
As a consequence of results in Tierney (1994) and Roberts and Smith (1994) we have:
Proposition 1 If the utility function is positive and integrable, p d (') ? 0, for all pairs (d; '),
and A and \Theta are intervals in IR n , the above scheme defines a Markov chain with stationary
distribution h.
It is impossible to give generally applicable results about when to terminate iterations in
Markov chain Monte Carlo simulations. It is well known that this is a difficult theoretical
problem, see, e.g., Robert (1995) and Polson (1996), who discuss approaches to find the
number of iterations that will ensure convergence in total variation norm within a given
distance to the true stationary distribution. However, practical convergence may be judged
with a number of criteria, see, e.g., Cowles and Carlin (1996) or Brooks and Roberts (1999).
Most of these methods have been implemented in CODA (Best et al, 1995), which we have
used in our examples. Once practical convergence has been judged, say after
we may record the next N iterations of the simulation output (d use
as an approximate sample from h(d). From that we may try to assess the mode.
We illustrate the above approach with an artificial example, adapted from Shenoy (1994).
Example 1. A physician has to determine a policy for treating patients suspected of
suffering from a disease D. D causes a pathological state P that, in turn, causes symptom
S to be exhibited. The physician observes whether a patient is
exhibiting the symptom. Based on this information, she either treats T the patient
(for P and D) or not 0). The physician's utility function depends on T; P and D,
as shown in Table 1. The value 0.001 was changed from the original value (0) to adapt
to the general result in Proposition 1. The probability of disease D (D = 1) is 0.1. For
patients known to suffer from D, 80% suffer from P 1). On the other hand, for

Table

1: The probability model p(D; the physician's utility function u(T;
h(d; D; S). The probabilities used in steps (ii) and (iii) of the Markov chain Monte Carlo
scheme described in the text are proportional to the entries in the appropriate column and
row, respectively, of the h section at the right of the table.
D
patients known not to suffer from D (D = 0), 15% suffer from P. For patients known to
suffer from P patients known not to suffer from
We assume that D and S are probabilistically independent given
P . To implement the proposed algorithm, we need to find the conditional distributions
h('jd) and h(dj'). In this case, is the decision
taken if the symptom is exhibited, and d 0 , if it is not exhibited. d means to treat
(not to treat) the patient. Let p(D; the probabilities
given in the above description. With h(d; D;
Our proposed method goes as follows: (i) Start at an arbitrary decision (d 0
steps (ii) and (iii) until convergence is judged.
Once convergence is judged, we record the next N iterations of the algorithm
use (d as an approximate sample from the marginal in d of the artificial
distribution. We leave out some values between those recorded to avoid serial correlation.
Since alternatives are finite in number, we just need to inspect the histogram to approximate
the mode. From a simulated sample of size 1000, we find that the optimal decision is
that is treat if symptom is present and not treat if symptom is absent.
Note in the example how the proposed augmented probability model simulation differs
from other simulation methods proposed for the solution of IDs. We use one simulation over
the joint ('; d) space to simulate from h(\Delta) instead of many small simulations to evaluate
expected utilities for each possible decision one at a time. Of course, the previous example is
extremely simple in that we are able to sample from h(dj') and h('jd), and, by inspection of
the histogram, we may approximate the modes. The next sections deal with more complex
cases.
3 Sampling from the Artificial Distribution
We shall provide here a generic method to sample from the artificial distribution h(\Delta). Typ-
ically, this distribution will not be straightforward to simulate from, requiring generation
from possibly high dimensional models, including complex probability and utility functions,
continuous decision and chance nodes, and possibly conditioning on observed data. MCMC
simulation schemes are the most commonly used methods known to accomodate such gen-
hence we choose them.
Given the enormous interest in IDs as a tool for structuring and solving decision problems,
see, e.g., Matzkevich and Abramson (1995), we concentrate on such structures. An ID is a
directed graph representation of a decision problem as a probability network with additional
nodes representing decisions and values. For notational purposes, we shall partition the set
of nodes into five subsets, differentiating three types of chance nodes: (i) Decision nodes
d, representing decisions to be made. (ii) Chance nodes, including random variables x
observed prior to making the decision, i.e., data available at the time of decision making;
not yet observed random variables y, i.e., data which will only be observed after making the
decisions; and unobservable random variables ', i.e., unknown parameters; (iii) One value
node u representing the utility function u(d; x; '; y). Figure 1 provides a simple generic ID for
our scheme. An ID is solved by determining a decision d   with maximum expected utility.
y
x
d
\Phi-
oe
@
@ @R
'i
'i
\Deltaff
'/
'i
@
@
@ @

Figure

1: A generic influence diagram for our scheme.
This requires marginalizing over all chance nodes (y; '), conditioning on x, and maximizing
over d. See Shachter (1986) for a complete description and an algorithm to solve IDs.
The method we propose here is applicable to IDs with non-sequential structure, i.e.,
decision nodes must not have any chance nodes as predecessors which have distributions
depending, in turn, on other decision nodes. Except for some technical conditions there will
be no further requirements.
3.1 The Probability Model Defined by Influence Diagrams
An ID defines the conditional distributions p(xj'), p(') and p d (yj'), a joint distribution
on (x; '; y) via p d (x; ';
p(')p(xj')p d (yj'); for ('; y) given the observed nodes x. Typically, x and y are independent
given ', allowing the given factorization, and p(') does not depend on d. If a particular
problem does not fit this setup, modifications of the proposed algorithm are straightforward.
In the context of this probability model, solving the ID amounts to maximizing the expected
utility over d, where p d ('; yjx) is the relevant distribution to compute this expectation.
In summary, solving the ID amounts to finding
d
Z
We shall solve this problem approximately by simulation. Augment the probability measure
to a probability model for ('; d) by defining a joint p.d.f.
The mode of the implied marginal distribution h(d) /
R
corresponds to the optimal decision d   . The underlying rationale of our method is to simulate
a Markov chain in ('; defined to have h('; d) as its asymptotic distribution. For
big enough t, the simulated values successive states of the simulated process
provide, approximately, a Monte Carlo sample from h('; d). Note that the simulation is
defined on an augmented probability model h(d; '; y) rather than on p d ('; y) for each possible
instantiation of the actions d, as traditional methods do. By considering the marginal distribution
of d t in this Monte Carlo sample, we can infer the optimal decision using methods
such as those discussed in Section 3.3.
The key issue is the definition of a Markov chain with the desired limiting distribution
h(:). For that, we capitalise on recent work in numerical Bayesian inference concerning the
application of Markov chain Monte Carlo methods to explore high dimensional distributions
which do not allow analytic solutions for expectations, marginal distributions, etc.
3.2 Markov Chain Monte Carlo Simulation
We shall provide a general algorithm which will be valid for all IDs satisfying the structural
conditions specified above and some minor technical conditions discussed below. The algorithm
we describe is of the Metropolis type (Tierney 1994): we generate a new candidate for
the states from a probing distribution, and then move to that new state or stay at the old
one according to certain probabilities. We do this transition in three steps, for d, ' and y.
We only require to be able to evaluate the utility function u(d; x; '; y) and the probability
distributions p d (yj'), p('), p(xj'), for any relevant d; x; '; y. This will typically be possible,
since the definition of the ID includes explicit specification of these distributions, i.e., the
modeler is likely to specify well-known distributions.
The scheme requires specification of probing distributions g 1 , g 2 and g 3 . The choice
of probing distributions g j (:j:) is conceptually arbitrary, with the only constraint that the
resulting Markov chain should be irreducible and aperiodic. As we shall argue, whenever
possible, we assume symmetric probing distributions, i.e., satisfying
example, g(ajb) could be a multivariate normal distribution N(b; \Sigma) for some \Sigma. Details
about the choice of probing distribution are discussed in the appendix. We then have:
Algorithm 1.
1. Start at values (d parameters and outcomes, and set
2. Let u
Generate a "proposal" ~
djd
d;
Compute
d;'
h(d
d (y
With probability a 1 , set d
d; otherwise, keep d
3. Let u
Generate a "proposal" ~
Compute
';y
With probability a 2 , set '
4. Let u
Generate a proposal ~
y).
Compute
u3
With probability a 3 set y
5. 1. Repeat steps 2 through 4 until chain is judged to have practically
converged.
This algorithm defines a Markov chain, with h('; d) as stationary distribution. The generality
of this algorithm comes at a price, namely possible slow convergence. Depending
on the application, long simulation runs might be required to attain practical convergence.
However, this fully general algorithm is rarely required.
Many problems allow simpler algorithms based on using p('jx) and p d (yj') to generate
proposals. Algorithm 2, given below, only requires a probing distribution g( ~
djd) for d, evaluation
of the utility function and algorithms to generate from p('jx) and p d (yj'). While
simulating from p d (yj') is typically straightforward, simulating from p('jx) is not. In gen-
eral, this distribution will not be explicitly specified in the ID, but needs to be computed
through repeated applications of Bayes formula, or several arc reversals in the language of
IDs. However, note that simulating from p('jx) amounts to solving the statistical inference
problem of generating from the posterior distribution on ' given the data x. Hence, we can
appeal to versions of posterior simulation schemes appropriate for a variety of important inference
problems recently discussed in the Bayesian literature (see, e.g., Smith and Roberts
1993; Tanner 1994; and Tierney 1994). Before starting the algorithm described below, we
generate a sufficiently large Monte Carlo sample from p('jx) by whatever simulation method
is most appropriate.
Algorithm 2.
1. Start at values (d
2. Evaluate u
3. Generate ( ~
d; ~
djd (i\Gamma1) )p ~
djd (i\Gamma1) )p( ~
'jx)p ~
d (~yj ~
4. Evaluate ~
d; x; ~
5. Compute
d; ~
h(d
d; ~
7. convergence is practically judged.
In step 3, generation of ' - p('jx) is done using the simulated Monte Carlo sample generated
Algorithm 3. The algorithm simplifies if x is missing in the ID, i.e., if no data is given at
the time of the decision. The associated Algorithm 3 would be stated as Algorithm 2, with
the proposal distribution in step 3 replaced by ( ~
d; ~
djd)p ~
d (~yj ~
Sampling from p d ('; will be feasible in general, since these distributions are
defined explicitly in the ID.
3.3 Finding the Optimal Solution
The MCMC simulation provides us with an approximate simulated sample f(d
y), from which we deduce an approximate sample (d
the marginal h(d). The mode of h(d) is an approximation of the optimal alternative.
In the case of discrete alternatives, the problem is simple since we only have to count the
number of times each element has appeared, and choose the one with the highest frequency.
It may be worthwhile retaining not one but several of the most frequent decisions, and study
them in further detail, as a way of conducting sensitivity analysis.
In the case of continuous alternatives, as a first approach we may use graphical exploratory
data analysis tools, especially with low dimensional decision vectors. When the
decision vector d is one- or two dimensional, we may produce the histogram (or a smooth
version) and inspect it to identify modes. For higher dimensional decision vectors d, we
propose to consider the problem as one of cluster analysis. Modes of h(d) correspond to d's
with higher density, which suggests looking for regions with higher concentration of sampled
d's. This leads us to compute a hierarchical cluster tree for the simulated points d t . Since
we are assuming h to be a density with respect to Lebesgue measure in IR n , and we are
interested in identifying regions where the optimal alternative might lie, we suggest using
complete linkage with Euclidean distance. Once we have a classification tree, we cut at a
certain height and obtain the corresponding clusters. The location of the largest cluster
indicates the area of the best decision. Again, as before, it may be useful to keep several
larger clusters and explore the corresponding regions. The result of course would depend
on the cutting height, but by exploring several heights we may be able to identify several
decisions of interest. We illustrate the approach in Section 4.2.
4.1 Example 2: A Medical Decision Making Problem
We illustrate the algorithm with a case study concerning the determination of optimal aphere-
sis designs for cancer patients undergoing chemotherapy. Palmer and M-uller (1998) describe
the clinical background and solve the problem by large scale Monte Carlo integration.
Between a pre-treatment and start of chemotherapy, stem cells (CD34) are collected to
allow later reconstitution of white blood cell components. Depending on the pre-treatment,
the first stem cell collection process (apheresis) is scheduled on the fifth or seventh day after
pre-treatment. A decision is to be made on which days between pre-treatment and treatment
we should schedule stem cell collections so as to (i) collect some target number of cells; and
(ii) minimize the number of aphereses. We have data on I = 22 past patients, and for the
first day of the new patient.
denote the observed CD34 count for patient i on
day t ij . Also, y
shall designate the i-th's patient data and
the combined data vector. Palmer and M-uller (1998) specify the following probability model
for this process. The likelihood is based on the observation that the typical profile of stem
cell counts over days shows first a rise after pre-treatment, reaches a maximum, and then
slowly declines back towards the base level, as shown in Figure 2. To model such shapes we
use a nonlinear regression model. Let g(t; a Gamma probability
density function with parameters chosen to imply mean and variance
matching e and s 2 and rescaled by
We use g(\Delta; e; s) to parametrize a nonlinear regression for the profiles through time of each
The prior model on the patient specific parameters is hierarchical. Patient i undergoes
one of two possible pre-treatments x i 2 f1; 2g, which serves as a covariate to specify the
first level prior: ' i - N(j x i
The hyperprior at the second level is common for both
cases: 2. The model is completed with a prior on V and oe

Figure

shows observed counts y ij and fitted profiles
typical patients. For a new patient I +1 denote with
) the (unknown)
stem cell counts on days t
. For a first day t 0 , we already have a count y h0 . Using
the notation introduced at the beginning of Section 3, is the observed
Figure

2: Three typical patients. The dashed lines connect the data points. The solid curve
plots the fitted profile, using the described probability model.
data vector, is the future data vector, and are
the unobservable parameters in the model. Given the typical profile, the optimal decision
will schedule aphereses for all days between some initial day d 0 and a final day d 1 , i.e., the
decision parameter is
Let A be the event of failing to collect a target number y   of stem cells,
y   g, where L h is the volume of blood processed at each stem cell collection for the new pa-
tient. Let n a the number of scheduled stem cell collections. The
utility function is u(d; x; '; c is the sampling cost and p a
penalty for underachievement of the target. We need to maximize over d the expected utility
R
Note that the probability model p(')p(x; yj') does not
depend on the decision nodes, but there is data x influencing the belief model. Since p('jx)
may be actually sampled with a Markov chain Monte Carlo method described in Palmer
and M-uller (1998), we use Algorithm 2 to solve the problem. To ensure a positive utility
function we add a constant offset to u(\Delta).
We found the optimal design d   for a future patient with the above belief and preference
model when 10:0. For a patient undergoing treatment x first observation
5, the optimal apheresis schedule for the remaining six days was
found to be given by d the decision space is two dimensional, we can
do this by a simple inspection of the histogram. Figure 3 plots the estimated distribution
Figure

3: The grey shades show a histogram of the simulated d t for the medical problem.
Inspection of h(d) reveals the optimal decision at
4.2 Example 3: A Water Reservoir Management Problem
In R'ios Insua et al (1997), we describe a complex multiperiod decision analysis problem
concerning the management of two reservoirs: Lake Kariba (K) and Cahora Bassa (C).
Here we solve a simplified version using the proposed MCMC approach to simulate from the
augmented probability model.
We want to find, for a given month, optimal values to be announced for release from
K and C through turbines and spillgates, d k
respectively. The actual amounts of
water released depend on the water available, which is uncertain, since there is uncertainty
about inflows i k and i c to the reservoirs. There is a forecasting model for both i k and i c ,
the latter being dependent on the water released from K and the incremental inflows (inc),
which, in turn, depend on a parameter fi. The preference model combines utilities for both
K and C. Those for K depend on the energy deficit (def), the final storage (sto k ) and the
amount of water spilled (spi). Those for C depend on the energy produced (ene) and the
final storage (sto c ). Initial storages s k and s c have influence as well over actual releases.

Figure

4 shows the influence diagram representing the problem. Nodes with double border
are either known values or deterministic functions of their predecessors. They are required
to compute the value node u, but will not show up in the probability model. In terms of
rel
sto k
I c ene
sto c
'i
@
@
@
@
@
@
\Gamma\Psi
\Deltaff
\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi*
-:
A
A
A
A
A
A
A
A
A
A
AU
\Deltaff
J-
-BBBBBBBN
'i
\Gamma\Psi
\Theta
\Theta
\Theta
\Theta
\Theta
\Theta
\Theta
J-
\Omega \Omega \Omega \Omega \Omega \Omega \Omega OE
\Phi \Phi \Phi \Phi*
Z
Z
Z
Z Z~
ae
ae
ae
ae
ae=
@
@
@
@
@
@
@ @R
-A
A
A
A
A
A
A AU b
'i
'i
'i
'i
'i
'i
'i
ae-
---Figure

4: Influence diagram for the reservoir problem.
our notation, the problem includes four decision nodes
nodes i k and fi.

Figure

5 shows some profiles of the histogram of the simulated d t - h(d), generated by
Algorithm 3. The decision parameter is four dimensional. Hence we used a four dimensional
cells) to record a four dimensional histogram of the simulated
states. Simple inspection of the empirical distribution allows to read off the optimal release
at d
200g. The solution is based on 100,000
simulated values from the Markov chain Monte Carlo scheme. Figure 5 illustrates also
another feature of our method which is a simple sensitivity analysis procedure at no extra
cost. Darkness of Figure 5b suggests that expected utility is rather flat when releases through
turbines are fixed at their optimal values, hence suggesting insensitivity with respect to
changes in spill. On the other hand, Figure 5a, with just one dark area where the estimated
optimum is, suggests that expected utility is fairly peaked in release through turbines, and
hence very sensitive to changes in energy releases.
Alternatively, as discussed in Section 3.3, we consider a hierarchical cluster tree of the
simulation output. The dots in Figure 5 show the solution based on cutting a hierarchical
cluster tree of 1000 simulated values d - h(d) at height 2000 and finding the cluster with
the most members. The optimum is found at d
476g. This comes reasonably close to the optimum estimated earlier.
5.1 Comparison with Alternative Schemes
The scheme described in Algorithms 1, 2 and 3 transforms the original expected utility maximization
problem (1) into a simulation problem. Our scheme is very generic, in the sense
of accomodating arbitrary probability models, be they discrete or continuous, and utility
(a) f(d k

Figure

5: (a) Expected utility as a function of release for energy (d k
spill fixed at
the optimal levels; and (b) as a function of spill (d k
release through turbines fixed at
the optima. The diamond indicates the optimal decision. The dots indicate the simulations
in the largest cluster of the hierarchical clustering tree cut at height 2000.
functions, as long as the probability density (or probability mass function) and the utility
function are pointwise evaluable. The main difference with other simulation schemes earlier
considered in the literature is that instead of using simulation to evaluate expected utilities
(losses) for each possible instantiation of decisions, we use simulation from an artificial
auxiliary model which augments the original probability model to include an artificial distribution
over decision nodes. Whether one or the other approach is more efficient depends on
the specifics of the considered decision problem. No general comparisons are possible. Even
in specific examples, performance will depend heavily on arbitrary choices like the amount
of discretization, which is necessary for many methods; run length of the involved Monte
Carlo simulations; chosen MCMC scheme, etc. However, some general observations about
the relative efficiency of the methods are possible.
In problems with few alternatives, analytic solutions using methods like arc reversal
(Shachter 1986), and simulation methods which use simulation to pointwise evaluate expected
utilities, like Likelihood Weighting (Shachter and Peot 1990), are typically more
efficient than simulation over the auxiliary probability model. Bielza and Shenoy (1998) discuss
a decision problem (the "reactor problem") with 6 possible actions, and chance nodes
with less than 10 possible joint outcomes. An exact solution using Shachter's (1986) algorithm
requires one arc reversal and the largest state space used during the solution phase
contains 4 variables. By comparison, we implemented the same example using augmented
probability simulation, following Algorithm 3. We used 100,000 iterations in the MCMC
simulation. The computational effort of one iteration is approximately comparable to one
arc reversal. Thus the exact solution is clearly far more efficient in terms of computing time.
Alternatively, consider simulation to compute the expected utility of each of the six possible
actions, using, for example, Likelihood Weighting. Considering the involved numerical standard
errors, Monte Carlo simulation sizes of around 1000 simulations for each alternative
decision would be adequate. Thus, also Likelihood Weighting dominates simulation from the
augmented probability model.
In problems where the optimal decision is to be computed conditional on some already
available data x the comparison changes, especially if the posterior distribution of the unknown
parameters is significantly different from the initial prior distribution, i.e., under low
prior probability for the evidence x. Consider, for example, the application reported in Section
4.1, which is not amenable to exact methods. Using Monte Carlo simulation to compute
expected utilities for alternative decisions, we can no longer generate independent samples.
Following Jenzarli's (1995) proposal we could use Gibbs sampling to compute expected util-
ities. Depending on the specific choices of the implemented MCMC scheme and termination
criteria, one would typically use on the order of 10,000 iterations (Palmer and M-uller 1998).
Discretizing the sample space, one could, in principle, also use Logic Sampling (Henrion,
1988). However, Logic Sampling would not be advisable for this problem since the fraction
of simulated experiments which generate variables corresponding to the actual observations
would be close to zero (i.e., p(x
in the notation of Shachter and Peot, 1990). For
similar reasons, Likelihood Weighting (Shachter and Peot 1990) would fail. Since only leaf
nodes are observed, the sample scores would be proportional to the likelihood function, i.e.,
the scheme would amount to importance sampling using the prior probability model as importance
sampling function. This can, however, be addressed using bounded variance type
algorithms as discussed, for example, in Pradhan and Dagum (1996).
Finally, many decision problems involve continuous decision variables, like the example
considered in Section 4.2. Continuous decision parameters create no problem for simulation
from the augmented probability model, but would not allow a straightforward application
of any scheme based on evaluating expected utilities for one decision at a time. Even if
discretization was used, say on grid, the resulting number of alternative
actions renders such schemes difficult to use.
5.2 Conclusion
Complex decision problems may render impossible the application of exact methods to obtain
optimal decisions. As a consequence, we should look for approximation methods, including
simulation.
We have proposed a simulation based strategy for approximating optimal decisions in a
decision analysis. Our experiments and examples suggest that this approach may be very
powerful. Implementation of the algorithms is fairly straightforward based on the schemes
provided. Specific cases may require simple modifications such as the ones suggested in
Section 3.2. The exploration of the sample in search for modes may be done with standard
statistical software. As we mentioned in the discussion of Example 3, one feature of our
method is the provision of simple sensitivity analysis features, at no extra cost.
A number of challenging problems remain, particularly perhaps, the extension of our
scheme to sequential decisions. The straightforward approach of expanding the model to
non-sequential normal form may only be applied when the number of decision nodes is small.
Another challenging problem would be to develop a computational environment based on
our approach. It would be also interesting to develop further methods to look for modes in
multivariate settings.
Similar ideas may be pursued to solve traditional statistical optimal design problems.
From a formal point of view, an optimal design problem can be described as a stochastic
optimization problem (1). This is explored in Clyde, M-uller and Parmigiani (1995) for the
special case of Algorithm 3 with continuous sample spaces and non-sequential setup.


Appendix

Implementation
The choice of the probing distributions g j (:j:) in Algorithm 1 is conceptually arbitrary, with
the only constraint that the resulting Markov chain be irreducible and aperiodic.
In the statement and proofs of the proposed algorithms, we assumed g to be symmetric
in its arguments, i.e., is a continuous parameter, we propose to
use a normal kernel g( ~
appropriately chosen covariance matrix \Sigma,
for example, a diagonal matrix with diagonal entries corresponding to reasonable step sizes
in each of the decision parameters. Good values for the step size can be found by trial and
error with a few values. In a particular setup, Gelman, Roberts and Gilks (1996) show that
the optimal choice of step size should result in average acceptance probabilities around 25%,
and similarly, for other parameters.
If d is discrete, a simple choice for g( ~
djd) could generate
0.5. Of course, many other problem specific choices are possible. In Example 2, e.g., we
define choosing with probability 1/6 one of six possible moves: (i) increase
d 0 and d 1 by 1 day; (ii) decrease d 0 and d 1 by 1; (iii) increase d 0 by 1; (iv) decrease d 0 by
etc.
Should symmetry of g be violated, an additional factor g(dj ~
djd) would be added in
the expressions for acceptance probabilities. This would correspond to Metropolis-Hastings
steps, rather than Metropolis steps. Convergence proofs for the proposed scheme are simple,
based on results in Tierney (1994) and Roberts and Smith (1994).

Acknowledgments

Research suported by grants from the National Science Foundation,
CAM, CICYT and the Iberdrola Foundation. Parts of it took place while Peter M-uller was
visiting UPM and David R'ios Insua was visiting CNR-IAMI. We are grateful to discussions
with Mike Wiper.



--R

"A comparison of graphical techniques for asymmetric decision problems"
"A comparison of hybrid strategies for Gibbs sampling in mixed graphical models,"
"Assesing convergence of Markov chain Monte Carlo algo- rithms,"
"Approaches for optimal sequential decision analysis in clinical trials,"
"A forward Monte Carlo method for solving influence diagrams using local computation,"
"Exploring expected utility surfaces by Markov changes,"
"A method for using belief networks as influence diagrams,"
"Markov chain Monte Carlo convergence diagnostics: a comparative review,"
"Efficient Metropolis jumping rules,"
Bayesian Statistics 5
"Sampling based approaches to calculating marginal densities,"
"Propagating uncertainty in bayesian networks by probabilistic logic sampling"
"Gibbs sampling in Bayesian networks,"
Solving Influence Diagrams using Gibbs sampling
"Local computations with probabilities on graphical structures and their application to expert systems,"
"Decision analytic networks in Artificial Intelligence,"
"Discrete approximations of probability distributions,"
"Optimal design via curve fitting of Monte Carlo experi- ments,"
"Bayesian Optimal Design in Population Models of Hematologic Data,"
"Fusion, propagation and structuring in belief networks,"
Probabilistic reasoning in intelligent systems
Decision Analysis with continuous and discrete variables: a mixture distribution approach
"Convergence of Markov chain Monte Carlo algorithms,"
"Optimal Monte Carlo estimation of belief network inference,"
"Bayesian methods in reservoir opera- tions: the Zambezi river case,"
"Convergence control methods for Markov chain Monte Carlo algorithms,"
"Simple conditions for the convergence of the Gibbs sampler and Metropolis-Hastings algorithms,"
"Probabilistic inference and influence diagrams,"
"Gaussian Influence Diagrams,"
"Simulation approaches to general probabilistic inference on belief net- works,"
"Decision making using probabilistic inference methods,"
"A comparison of graphical techniques for decision analysis,"
"Bayesian computational methods,"
"Moment methods for Decision Analysis,"
"Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods,"
Tools for Statistical Inference
"Markov chains for exploring posterior distributions (with discussion),"
"Use of the Gibbs sampler in expert systems,"
--TR
simulation;decision analysis;influence diagrams;optimal design;Markov Chain Monte Carlo