 %Latex2e: Documentclass and packages
\documentclass{article}
\usepackage{chicago}

%%% Text dimensions and basic formatting
\setlength{\textheight}{20cm}
\setlength{\textwidth}{15.7cm}
\setlength{\topmargin}{0.5cm}
\setlength{\oddsidemargin}{0.5cm}
\setlength{\evensidemargin}{0.5cm}
\setlength{\parskip}{10pt}
\setlength{\unitlength}{1mm}

%%% Double spacing, as indicated by the LaTeX Companion book.
\renewcommand{\baselinestretch}{1.67}

%%% Newcommands
\newcommand{\ex}[1]{E\!\left[ #1 \right]}
\newcommand{\lexo}{\underline{E}}
\newcommand{\uexo}{\overline{E}}
\newcommand{\lex}[1]{\underline{E}\!\left[ #1 \right]}
\newcommand{\uex}[1]{\overline{E}\!\left[ #1 \right]}

\newcommand{\pr}[1]{P\!\left( #1 \right)}
\newcommand{\pro}{P}
\newcommand{\lpro}{\underline{P}}
\newcommand{\upro}{\overline{P}}
\newcommand{\lpr}[1]{\underline{P}\!\left( #1 \right)}
\newcommand{\upr}[1]{\overline{P}\!\left( #1 \right)}

\newcommand{\distr}[1]{P\!\left( #1 \right)}
\newcommand{\distro}{P}
\newcommand{\ldistro}{\underline{P}}
\newcommand{\udistro}{\overline{P}}
\newcommand{\ldistr}[1]{\underline{P}\!\left( #1 \right)}
\newcommand{\udistr}[1]{\overline{P}\!\left( #1 \right)}

\newcommand{\den}[1]{p\!\left( #1 \right)}
\newcommand{\deno}{p}
\newcommand{\ldeno}{\underline{p}}
\newcommand{\udeno}{\overline{p}}
\newcommand{\lden}[1]{\underline{p}\!\left( #1 \right)}
\newcommand{\uden}[1]{\overline{p}\!\left( #1 \right)}

%%% Newtheorems
\newtheorem{Example}{Example}
\newtheorem{Theorem}{Theorem}

\title{\vspace{-1.8cm}Calculation of Posterior Bounds Given Convex Sets of Prior
	 Probability Measures and Likelihood Functions}

%\author{Fabio Gagliardi Cozman\thanks{Fabio G. Cozman is Assistant Professor,
%	Escola Polit\'ecnica, Universidade de S\~ao Paulo, Av. Prof.
%	Mello Moraes 2231, Cidade Universit\'aria 05508-900, S\~ao Paulo, Brazil
%	(fgcozman@usp.br, http://www.cs.cmu.edu/\~{}fgcozman/home.html).
%	This research was partially conducted while the
%	author was with the Robotics Institute, Carnegie Mellon University.
%	The project was partially supported by NASA under Grant NAGW-1175;
%	the author was supported under a scholarship from CNPq, Brazil.} \\
%	University of S\~ao Paulo, S\~ao Paulo, Brazil}

\date{}

\begin{document}

\maketitle

\begin{abstract}
This paper presents alternatives and improvements to Lavine's algorithm,
currently the most popular method for calculation of posterior expectation
bounds induced by sets of probability measures. Firstly, methods from
probabilistic logic and Walley's and  White-Snow's algorithms are reviewed
and compared to Lavine's algorithm. Secondly, the calculation of posterior
bounds is reduced to a fractional programming problem: from the unifying
perspective of fractional programming, Lavine's algorithm is derived from
Dinkelbach's algorithm, and the White-Snow algorithm is similar to
the Charnes-Cooper transformation. From this analysis, a novel algorithm
for expectation bounds is derived. This algorithm provides a complete solution
for the calculation of expectation bounds from priors and likelihood functions
specified as convex sets of measures. This novel algorithm is then
extended to handle the situation where several independent identically
distributed measurements are available. Examples are analyzed
through a software package that performs robust inferences and that
is publicly available.
\end{abstract}

{\flushleft
{\em Keywords:} Bayesian robustness analysis, lower and upper envelopes,
		Lavine's algorithm, White-Snow algorithm, fractional
		programming. (AMS classification: Primary 62G35; Secondary 90C90,
		90C32, 28B20.)}
% 62G35  Robustness
% 90C90  Applications of mathematical programming
% 90C32  Fractional programming
% 90C05  Linear programming
% 28B20  Set-valued set functions and measures


\section{INTRODUCTION}

This paper presents alternatives and improvements to Lavine's algorithm,
currently the most popular method for calculation of posterior expectation
bounds induced by sets of probability measures. Models based on convex sets
of measures have been advocated as more realistic, meaningful and
flexible than standard probability models (e.g., as discussed
by \citeN{Giron80}, \citeN{Levi80} and \citeN{Walley91}).

The most important practical application of convex sets
of probability measures is {\em robustness analysis}.
The idea, detailed by \citeN{Berger90},
is to employ sets of probability measures to represent
perturbations and variations in a probabilistic model.
The goal of robustness analysis is to produce bounds on expected values.
The interval between the upper and lower bounds induced by a
convex set of probability measures reflects the quality of model and data;
small intervals indicate {\em robustness} to perturbations.

\citeN{Lavine91} has proposed a bracketing algorithm, aimed at
robustness analysis, that has captured great attention in recent years.
The main goal of this paper is to present algorithms that conduct
robustness analysis more generally than
Lavine's algorithm. A second important goal of the paper is to
indicate the connections between the research on convex sets
of measures that has been developed in Artificial Intelligence
and Statistics. Note that the word {\em algorithm} is used here
informally to mean any sequence of operations that produces probability
and expectation bounds.

Firstly, algorithms previously developed in Artificial Intelligence
are compared to Lavine's algorithm. Methods from probabilistic logic
are reviewed, and the White-Snow algorithm is analyzed and presented
in a uniform notation. Walley's algorithm is also compared to Lavine's
algorithm. Secondly, the calculation of posterior bounds
is reduced to a fractional programming problem: from the unifying
perspective of fractional programming, Lavine's algorithm is identical
to Dinkelbach's algorithm, and the White-Snow algorithm is similar to
the Charnes-Cooper transformation. From this analysis, a novel algorithm
for expectation bounds is derived. This algorithm provides a complete
solution for the calculation of expectation bounds from priors and
likelihood functions specified as convex sets of measures. This novel
algorithm is then extended to handle the situation where several independent
identically distributed experiments are conducted, given that the experiments
are modeled by convex sets of likelihood functions.

The novel algorithms developed in this paper have been implemented
in a Matlab$^{TM}$ package that is publicly available. To illustrate
the algorithms in this package, two examples are analyzed
in Section~\ref{section:Examples}.

\section{THE THEORY OF CONVEX SETS OF PROBABILITY MEASURES}
\label{QuasiBayesianTheory}

This section presents the basic technical results used in
this paper. Several theories of inference advocate sets of
probability measures as representations for statements of uncertainty;
for example, the quasi-Bayesian theory of \citeN{Giron80},
inner/outer
measures \cite{Suppes74,Good83,Ruspini87,Halpern92}, lower
probability theory \cite{Smith61,Fine88,Breese91,Chrisman96},
the intervalism theory of \citeN{Kyburg87}, the
convex Bayesian theory of \citeN{Levi80}, the theory of coherent
lower previsions of \citeN{Walley91}, and the
very general theory of probability/utility sets of \citeN{Seidenfeld95}.
Some theories, like Dempster-Shafer theory,
use representations that can be recast in terms of
convex sets of measures \cite{Dempster67,Kyburg87,Wasserman90}.
This article adopts the ideas advocated by Walley's theory of coherent
lower previsions, but emphasizes an interpretation of the concepts
that is based on convex sets of probability
measures, much in the spirit of the quasi-Bayesian theory of \citeN{Giron80}.

A closed convex set of probability measures is called a {\em credal set}
by \citeN{Levi80}. A credal set containing joint probability measures for
a collection of variables is called a {\em joint credal set}. Denote by
$\pr{\cdot}$ either the probability $\pr{A}$ of an event $A$ or
the distribution $\distr{X}$ of a random variable $X$.
Use the symbol $\den{X}$ to denote either the
probability mass function of $X$ (when $X$ is discrete-valued) or the
density function of $X$ (when $X$ is continuous-valued). A credal set
defined by a set of distributions $\distr{X}$ is denoted by $K(X)$.
To simplify the discussion, assume that every
distribution $\distr{X}$ is defined on the power set of values of $X$.

\subsection{Lower and Upper Expectations, Envelopes and Densities}

{\em Lower and upper expectations} for a bounded function $f(X)$,
given a credal set $K$, are defined as:
\[
\lex{f(X)} = \min_{\pro \in K} E_p[f(X)],
\hspace{1cm}
\uex{f(X)} = \max_{\pro \in K} E_p[f(X)],
\]
where $E_p[f(X)]$ is the standard expectation of $f(X)$.
Lower expectations can be obtained from upper expectations through the
expression $\lex{f(X)} = - \uex{-f(X)}$.
There is a one-to-one correspondence between credal sets and collections
of coherent lower (or upper) expectations (the definition of coherence
is given by \citeN{Walley91}). Any credal set generates a unique
collection of coherent lower (or upper) expectations and vice-versa.

Given a credal set $K$, a probability interval is induced for every event $A$:
\[
\lpr{A} = \min_{\pro \in K}  \pr{A},
\hspace{1cm}
\upr{A} = \max_{\pro \in K} \pr{A}.
\]
The set-functions $\lpr{A}$ and $\upr{A}$ are called {\em lower and upper
envelopes} respectively. Similarly, the functions
$\lden{X} = \min_{\pro \in K} \den{X}$ and
$\uden{X} = \max_{\pro \in K} \den{X}$ are called {\em lower and upper densities}
respectively. Lower envelopes can be obtained from upper envelopes
through the expression $\lpr{A} = 1 - \upr{A^c}$.
For any event $A$, the lower envelope
of $A$ is obtained by taking the lower expectation of the indicator
function $\delta_A(X)$, which is one if $X \in A$ and zero otherwise.

\subsection{Conditioning}

Convex sets of conditional probability measures are used to represent
the beliefs held by a decision maker {\em given} an event. A credal
set defined by conditional distributions $\distr{X|A}$ (for variable $X$
and event $A$) is denoted by $K(X|A)$. A collection of credal sets $K(X|y)$,
indexed by a variable $Y$, is denoted by $K(X|Y)$. To simplify terminology,
$K(X|Y)$ is also called a credal set.

A credal set $K(X|Y)$ is {\em separately specified} when the constraints
that define $K(X|y_1)$ do not interfere with the constraints that define
$K(X|y_2)$ for $y_1 \neq y_2$. More generally, a collection of credal
sets is separately specified when measures can be independently selected
from the credal sets. A conditional credal set $K(X|Y)$ is separately 
specified when each credal set $K(X|y)$ is specified by coherent lower 
expectations given all values of the conditioning 
variable \cite[Chapter~6]{Walley91}.

Inference is performed by applying Bayes rule to
each measure in a credal set; the posterior credal set
is the union of all posterior probability measures (more details
can be found in the Internet at http://www.cs.cmu.edu/\~{}qBayes/Tutorial).
To obtain a posterior credal set, one has to apply Bayes
rule only to the extreme points of a joint credal set and then take the convex hull
of the resulting posterior probability measures \cite{Giron80,Levi80}.
To obtain maximum and minimum values of posterior probabilities,
one must look only at the extreme points of the posterior credal
set \cite[Section~6.4.2]{Walley91}.

Given a joint credal set $K(X, Y)$ and a bounded function $f(X)$,
the functionals $\uex{f(X)|Y}$ and $\uex{g(X, Y)}$ (for arbitrary
bounded functions $f(X)$ and $g(X, Y)$) are closely related
by the {\em generalized Bayes rule} (first proposed by
\citeN[Section~6.4.1]{Walley91}). For an event $A$ such that
$\lpr{A} > 0$, $\uex{f(X)|A}$ is the unique value of $\mu$ that
solves the equation:
\begin{equation}
\label{equation:GeneralizedBayesRule}
\uexo \left[
	\left( f(X) - \mu \right) \delta_A(Y) \right] = 0.
\end{equation}
Note that the generalized Bayes rule uses operations on the
joint credal set $K(X, Y)$ to generate conditional values; the same
technique is used in Lavine's algorithm (Expression (\ref{equation:GeneralLavine})).

Suppose that a credal set $K(X)$ and a separately specified credal
set $K(Y|X)$ are given. To obtain the posterior expectation $\uex{f(X)|y}$,
it is necessary to apply the generalized Bayes rule to the joint credal
set $K(X,Y)$ that has the correct marginal credal set $K(X)$ and the correct
conditional credal set $K(Y|X)$. The following theorem investigates this
situation.

\begin{Theorem}
\label{theorem:LowerUpperLikelihoods}
Consider a bounded function $f(X)$ and suppose that $K(X)$ and
$K(Y|X)$ are separately specified. For a given $y$,
define the {\em lower likelihood}
$L_y(x) = \lden{y|x}$ and the {\em upper likelihood} $U_y(x) = \uden{y|x}$.
If $\lden{y} > 0$, $\uex{f(X)|y}$ is equal to the unique value of $\mu$
that satisfies the equation
\[
\uexo \left[ \left( f(X) - \mu \right) p_{\mu}(y| X) \right] = 0,
\; \mbox{ where } \;
p_{\mu}(y|x) = \left\{
\begin{array}{ccc}
U_y(x) & \mbox{ if } & f(x) \geq \mu, \\
L_y(x) & \mbox{ if } & f(x) < \mu.
\end{array}
\right.
\]
$\lex{f(X)|y}$ is obtained by solving a similar equation (obtained by
replacing $\lexo$ for $\uexo$, greater than by smaller than and vice-versa).
\end{Theorem}

{\bf Proof.} Direct from \citeN[Section~8.5.3]{Walley91}.

\subsection{Inferences with Convex Sets of Probability Measures in Probabilistic
Logic and Robust Statistics}

A credal set can be understood as a collection
of constraints on probability measures.
This viewpoint has been emphasized in Artificial Intelligence,
particularly after the discussion of probabilistic logic by \citeN{Nilsson86}.
Several variants of probabilistic logic exist
\cite{Bacchus90,Halpern90,Thone92,Frisch94,Lukasiewicz95,DuboisSMC96};
the purpose of the work is always to start from a collection of linear
constraints on the probability of propositions, and obtain probability
bounds through linear programming. Conditional and posterior constraints
can be handled to a limited extent by \citeN{Halpern90}
and \citeN{Lukasiewicz96}.

A different line of research has focused on the properties of
2-monotone and infinite monotone Choquet capacities. These
credal sets admit closed-form expressions for posterior quantities when a
single measurement is performed \cite{Halpern92,WassermanKadane90,Walley81}.
\citeN{Chrisman95} has modified these closed-form expressions to handle
sequences of measurements.

Robust Statistics focuses mostly on inferences involving a
conditional distribution $\distr{X|\theta}$ and a prior credal set $K(\Theta)$
\cite{Berger85Book,Kadane84,Berger90,Wasserman92}. Inference is equated
to calculation of the posterior
credal set $K(\Theta|x)$. A number of results are available for
important cases, such as density ratio classes \cite{Berger90} and
unimodal and symmetric constraints \cite{Berger90}.

There are three common ways to specify credal sets:
either by specifying a finite collection of extreme points, or by
specifying a finite collection of linear inequalities, or by
generating a neighborhood of measures
with convenient properties (for example, credal sets generated by 2-monotone
capacities).\footnote{This classification of problems, and the fact
that Lavine's algorithm can use $f(X)\delta_A(X)$ rather than
$f(X)$, to compute its starting point, were suggested to me by Peter Walley.}
The first situation is relatively simple:
If a credal set $K(X)$ is specified by a collection of extreme points
and $A$ is an event such that $\lpr{A} > 0$, 
the conditional expectation $E_p[f(X)|A]$ can be computed
for each extreme point of $K(X)$, and $\uex{f(X)|A}$
is the maximum of these conditional expectations. 
The second situation (collection of linear inequalities) can be
handled through Lavine's and White-Snow's algorithms.
The third situation (neighborhoods of measures) can be
handled through Lavine's and Walley's algorithms.
All algorithms assume that the lower envelope
of the conditioning event is larger than zero, an
assumption that is taken for granted in this paper.

Only a few authors consider the possibility that prior {\em and}
conditional credal sets be specified \cite{Lavine91,Pericchi94,Walley91}).
An algorithm for prior and conditional credal sets,
that completely solves this problem, is derived through linear fractional
programming methods in Section \ref{section:LinearFractionalProgramming}.

\section{LAVINE'S AND WALLEY'S ALGORITHMS}
\label{section:LavineAlgorithm}

The posterior upper expectation for a bounded function $f(X)$ conditional
on event~$A$ with positive lower envelope is:
\begin{equation}
\label{equation:FirstLavine}
\uex{f(X)|A} = \max \left[
	\frac{ E_p[f(X)\delta_A(X)] }{ E_p[\delta_A(X)] } \right],
\end{equation}
where the maximization is with respect to
all the measures in the joint credal set.
To simplify the presentation, only upper expectations
are considered; lower expectations can be obtained
through similar computations.

Lavine's algorithm can be informally stated as follows.
Pick a real number $\mu$ in the interval
$[\inf f(X)\delta_A(X),$ $\sup f(X)\delta_A(X)]$ and
check whether $\uex{f(X)|A}$ is larger, smaller or equal to $\mu$, and
respectively increase $\mu$, decrease $\mu$ or stop. When $f(X)$
is bounded, repetition of this procedure is certain to produce an interval
containing $\uex{f(X)|A}$; the algorithm stops when this interval is small
enough, or the number of repetitions exceeds some threshold.

Note that \cite{Lavine91}:
\begin{equation}
\label{equation:GeneralLavine}
\uex{f(X)|A} > \mu
\Leftrightarrow
\max \left[ E_p[f(X)\delta_A(X)] - \mu E_p[\delta_A(X)] \right] > 0
\Leftrightarrow
\uex{ (f(X) - \mu) \delta_A(X) } > 0,
\end{equation}
demonstrating that Lavine's algorithm is simply a bracketing
solution for the generalized Bayes rule
(Expression (\ref{equation:GeneralizedBayesRule})).

Lavine's algorithm is a viable choice if the calculation of
$\uex{ (f(X) - \mu) \delta_A(X) }$ is in fact easier than the
nonlinear maximization in Expression (\ref{equation:FirstLavine}).
For example, Lavine's algorithm is effective when applied to
density ratio classes and 2-monotone Choquet capacities.
Apart from these special cases, Lavine's algorithm is practical
when a single conditional distribution $\distr{X|\theta}$ and a prior
credal $K(\theta)$ are specified, because in this case each iteration
of the algorithm (Expression (\ref{equation:GeneralLavine}))
demands the maximization of a linear functional. Such a maximization
is straightforward when $\theta$ has finitely many values or is
discretized.

\citeN[Note~6.4.1]{Walley91} has proposed an iterative algorithm
to obtain posterior quantities by calculating a series of upper expectations.
In Walley's algorithm, $\uex{f(X)|A}$ is obtained by iterating
\[
\mu_{i+1} = \mu_i + 2 \uex{(f(X) - \mu_i) \delta_A(X)}/(\upr{A} + \lpr{A}).
\]
Walley also proved that the error at step $n$ is bounded by $c \Delta^n$,
where $c$ is a constant 
and $\Delta = (\upr{A} - \lpr{A})/(\upr{A} + \lpr{A})$; consequently,
Walley's algorithm has the property that $\epsilon_{n+1} = \Delta \epsilon_n$,
where $\epsilon_n$ is a bound on the error at step $n$. 
Walley's algorithm is attractive when $\uex{ (f(X) - \mu_i) \delta_A(X) }$ 
can be easily computed (for example, models generated by 2-monotone
Choquet capacities). To compare Walley's algorithm to Lavine's,
note first that Lavine's algorithm has linear convergence. If bisection
is used in Lavine's algorithm, then $\epsilon_{n+1} = (1/2) \epsilon_n$,
and Walley's algorithm converges faster than Lavine's when $\Delta < (1/2)$.

In the remainder of this paper, it is assumed that $\theta$
has a finite number of values and $K(\theta)$ is finitely
generated, so that techniques of linear programming
can be applied and compared to Lavine's algorithm.

\section{WHITE-SNOW ALGORITHM}
\label{section:WhiteSnowAlgorithm}

Lavine's algorithm requires the solution of a parametric linear
program\footnote{The use of parametric linear programming in Lavine's
algorithm was pointed out by one of the refrees.} in the value of $\mu$ when
(i) a single conditional distribution $\distr{X|\theta}$ is specified; and
(ii) a prior credal $K(\theta)$ is specified through linear constraints
$\mathbf{A} [ \den{\theta_1} \dots \den{\theta_n} ]^T \leq \mathbf{B}$,
where $\mathbf{A}$ and $\mathbf{B}$ are matrices of appropriate dimensions.
An algorithm that produces upper and lower expectations by direct
application of linear programming has been developed by \citeN{White86}
and improved by \citeN{Snow91SMC} for these situations. The White-Snow
algorithm was developed in the context of
expert systems and no connection to Lavine's algorithm has been mentioned
in the literature so far. The White-Snow algorithm is summarized
below in uniform notation.

For a given $x$, define:
\begin{itemize}
\item the vectors $\alpha_i = \den{\theta_i}$, $\beta_i = \den{x|\theta_i}$ and
      $f_i = f(\theta_i)$;
\item the matrix $\mathbf{C} = \mathbf{A} - \mathbf{B} \mathbf{1}$, where
$\mathbf{1}$ is a row vector of ones;
\item the matrix
$\mathbf{D} = \mathbf{C} \times
	\mbox{ diag} \left[ \beta_1^{-1}, \dots , \beta_n^{-1} \right]$.
\end{itemize}

The calculation of a posterior upper expectation is:
\[
\uex{f(\theta)|x} = \max_{\alpha} \left[ \frac{ \sum_i f_i \alpha_i \beta_i }
                             { \sum_j \alpha_j \beta_j } \right],
\]
subject to:
\[
\mathbf{C} \alpha \leq 0,
\hspace{1cm}
\sum_i \alpha_i = 1,
\hspace{1cm}
\alpha_i \geq 0.
\]

\citeN{White86} proposed the following change of variables:
\[
\gamma_i = \frac{\alpha_i \beta_i}{ \sum_j \alpha_j \beta_j },
\]
which reduces the calculation of the posterior upper expectation to
a linear program:
\[
\uex{f(\theta)|x} = \max_{\gamma} \left[ \sum_i f_i \gamma_i \right],
\]
subject to:
\[
\mathbf{D} \gamma \leq 0,
\hspace{1cm}
\sum_i \gamma_i = 1,
\hspace{1cm}
\gamma_i \geq 0.
\]

To obtain the matrix $\mathbf{D}$, all $\beta_i$ must be larger
than zero. When $\beta_i = 0$, set the variable $\gamma_i$
to zero and discard it from the equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Snow's solution is actually to create an auxiliary variable that
%%% selects which probability values are zero; this has the disadvantage
%%% of allowing unbounded solutions for $\gamma_i$ in the linear program.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\citeN{Snow91SMC} has proved that the solution of this linear program yields
the correct posterior upper expectation.

\section{POSTERIOR BOUNDS THROUGH LINEAR FRACTIONAL PROGRAMMING}
\label{section:LinearFractionalProgramming}

The algorithms presented in Sections \ref{section:LavineAlgorithm}
and \ref{section:WhiteSnowAlgorithm} were derived through special
properties of upper expectations. Lavine's and White-Snow's algorithms 
can be derived
in a more direct and general way as algorithms for linear fractional
programming. Only recent references associate linear fractional programs
to upper expectations \cite{Betro96,Jaumard91,Luo96,Zaffalon97};
only \citeN{Pacifico94} compare linear fractional programming to
Lavine's algorithm.

Linear fractional programming studies the maximization of ratios
of linear functions \cite{Ibaraki81,Schaible81}. Consider the
linear fractional program:
\[
\max_{\alpha} \left[ \frac{ \sum_i f_i \alpha_i \beta_i }
                 { \sum_j \alpha_j \beta_j } \right],
\]
where the $\beta_i$, $f_i$ are given and the $\alpha_i$ are
subject to linear constraints.
There are two main algorithms for linear fractional
programming:
\begin{itemize}
\item Create a ``parameterized'' problem for a parameter $\mu$:
\[
\max_{\alpha} \left[ \sum_i (f_i - \mu) \alpha_i \beta_i  \right],
\]
subject to the same constraints on the original problem; and search
for the value of $\mu$ such that the maximum is zero. This method
is variously called Dinkelbach or Jagannatham algorithm, and dates
back to the sixties. Note that this method
is identical to Lavine's algorithm.

\item Transform the problem by a change of variables:
\[
\gamma'_i = \frac{\alpha_i}{ \sum_j \alpha_j \beta_j },
\]
which reduces the calculation of the posterior upper expectation to
a linear program:
\[
\max_{\gamma'} \left[ \sum_i f_i \beta_i \gamma'_i \right],
\]
subject to:
\[
\mathbf{C} \gamma' \leq 0,
\hspace{1cm}
\sum_i \beta_i \gamma'_i = 1,
\hspace{1cm}
\gamma'_i \geq 0.
\]
This is called the Charnes-Cooper method, and also dates back to the
sixties. Notice that this method is similar to White-Snow's
algorithm; the only difference is that $\gamma_i = \gamma'_i \beta_i$.
On one hand, the
Charnes-Cooper method has the advantage of automatically handling the case
$\beta_i = 0$; on the other hand, the White-Snow algorithm 
has the advantage that the posterior
credal set is directly represented by the variables $\gamma_i$.
\end{itemize}

\section{POSTERIOR UPPER EXPECTATIONS FOR PRIOR AND CONDITIONAL CREDAL SETS}
\label{section:PosteriorUpperExpectationPriorConditional}

The main results of the paper are described in this section. A novel
method for the calculation of posterior bounds given separately
specified prior and conditional credal sets ($K(\theta)$ and $K(X|\theta)$
respectively) is proposed and then extended to the case of multiple
independent measurements.
A simplified version of this problem, where an interval is associated
with each value of $\distr{X|\theta}$, is studied by \citeN{Snow96}; the same
techniques are used in the following derivation. Approximate bounds for
the same problem are provided by \citeN{Salo96SMC} using different techniques.
\citeN[Page~18]{Walley96} presents closed-form solutions for a narrower version
of the problem where both prior and conditional credal sets are defined by
density bounded families.

First note that Theorem \ref{theorem:LowerUpperLikelihoods} demonstrates
that separately specified prior and conditional credal sets can be
handled through the function
$h(\mu) = \uex{ (f(\theta) - \mu) p_{\mu}(x|\theta)}$.
The function $h(\mu)$
is strictly decreasing with $\mu$, so the solution of $h(\mu) = 0$ can be
obtained by bracketing $\mu$ in the 
interval $[\inf f(\theta), \sup f(\theta)]$.
This generates a sequence of linear programs, and closely mimics Lavine's
bracketing algorithm.

It is actually possible to compute posterior upper expectations given sets
of priors and likelihood functions through a single linear program.
By Theorem \ref{theorem:LowerUpperLikelihoods}, the upper expectation
is attained only by likelihood functions that assume either $U_x(\theta_i)$
or $L_x(\theta_i)$ for each $\theta_i$.
Take the vectors $\alpha$ and $f$ defined in Section
\ref{section:WhiteSnowAlgorithm}. 
Define two new vectors, $\alpha'$ and $\alpha''$,
each with the same length as $\alpha$.
Consider the following linear fractional program:
\[
\max_{\alpha', \alpha''}
	\left[
		\frac{ \sum_i \left( f_i L_x(\theta_i) \alpha'_i +
			   f_i U_x(\theta_i) \alpha''_i \right) }
		     { \sum_j \left( L_x(\theta_j) \alpha'_j +
			   U_x(\theta_j) \alpha''_j \right) } \right],
\]
subject to:
\[
\mathbf{C} ( \alpha' + \alpha'') \leq 0,
\hspace{1cm}
\sum_i ( \alpha'_i + \alpha''_i) = 1,
\hspace{1cm}
\alpha'_i \geq 0,
\hspace{1cm}
\alpha''_i \geq 0.
\]
Now the Charnes-Cooper transformation can be applied and
the lower expectation can be obtained through a linear program.
For each $\theta_i$, a solution of this linear
program has either $\alpha'_i = 0$ or $\alpha''_i = 0$.
This automatically selects the correct likelihood value.
This method is more general than the procedure
derived by \citeN{Snow96}, and the derivation is shorter and
simpler; note that the algorithm automatically
handles situations where the likelihoods are zero, while Snow
has to consider a variety of special cases.

Suppose now that a sequence of measurements $X_1, \dots, X_n$, is given,
and the measurements are all taken to be independent and modeled by
identical sets $K(X_i|\theta)$ of likelihood functions. The objective
of robust inference is to calculate the bounds
$\lex{f(\theta) | x_1, \dots, x_n}$ and $\uex{f(\theta) | x_1, \dots, x_n}$.

There are several different interpretations to the statement
that ``measurements are independent'' (as discussed by
\citeN{Chrisman96Independence} and \citeN{deCamposUAI95}).
Results derived in this section are valid for several
definitions in the literature, but, for definiteness,
Walley's definition of independence \cite[Chapter~9]{Walley91}
is adopted throughout the discussion.

Variables $X$ and $Y$ are independent given a value of $Z$ when
the lower expectation $\lex{f(X)|y, z}$ is equal to the lower
expectation $\lex{f(X)|z}$ for any bounded function $f(X)$, {\em and}
the lower expectation $\lex{g(Y)|x, z}$ is equal to the lower
expectation $\lex{g(Y)|z}$ for any bounded function $g(Y)$.
For example, the statement ``measurements $X_1, \dots, X_n$ are independent
given $\theta$'' indicates that the credal set
$K(X_i|X_1, \dots, X_{i-1}, X_{i+1}, \dots,  X_n, \theta)$
contains the same functions as the credal set  $K(X_i|\theta)$.

The central problem is how to calculate the lower and upper
likelihoods $L_{x_1, \dots, x_n}(\theta)$ and $U_{x_1, \dots, x_n}(\theta)$;
once these functions are calculated, the same algorithm
developed in Section \ref{section:PosteriorUpperExpectationPriorConditional}
can be used to obtain posterior bounds. The following simple result
yields the likelihoods:

\begin{Theorem}
The upper and lower likelihoods are given by:
\[
U_{x_1, \dots, x_n}(\theta) = \prod_i \left[ \uden{x_i|\theta} \right],
\hspace{1cm}
L_{x_1, \dots, x_n}(\theta) = \prod_i \left[ \lden{x_i|\theta} \right].
\]
\end{Theorem}

{\bf Proof.} Note the inequalities:
\begin{equation}
\label{equation:LowerProductInequality}
L_{x_1, \dots, x_n}(\theta) = \min \den{x_1, \dots, x_n|\theta} =
 \min \left[ \prod_i \den{x_i|x_{i+1}, \dots, x_n, \theta} \right] \geq
		\prod_i \left[ \lden{x_i|\theta} \right],
\end{equation}
\begin{equation}
\label{equation:UpperProductInequality}
U_{x_1, \dots, x_n}(\theta) = \max \den{x_1, \dots, x_n|\theta} =
 \max \left[ \prod_i \den{x_i|x_{i+1}, \dots, x_n, \theta} \right] \leq
		\prod_i \left[ \uden{x_i|\theta} \right].
\end{equation}
To obtain the lower and upper likelihoods, construct distributions
$p'$ and $p''$. Construct $p'(x_1, \dots, x_n|\theta)$ by multiplying
together the likelihood functions that yield the minimum likelihood
value for the measurements. Construct $p''(x_1, \dots, x_n|\theta)$
by multiplying together the likelihood functions that yield the maximum
likelihood value for the measurements.
From that,
$p'(x_1, \dots, x_n|\theta) = \prod_i \left[ \lden{x_i|\theta} \right]$
and $p''(x_1, \dots, x_n|\theta) = \prod_i \left[ \uden{x_i|\theta} \right]$.
These equalities demonstrate that the inequalities
(\ref{equation:LowerProductInequality}) and
(\ref{equation:UpperProductInequality}) are in fact tight.

\section{EXAMPLES}
\label{section:Examples}

The methods developed in this paper can be implemented on top of
standard linear programming packages. The following examples were
solved through procedures implemented with the optimization
facilities available in the {\em Matlab}$^{TM}$ system
(Appendix \ref{appendix:Software}).

\begin{Example}[\citeN{White86}]
\label{WhiteSnowExample}
Consider a variable $\theta$ with four values
$\{ \theta_1, \theta_2, \theta_3, \theta_4 \}$, and the following
constraints on the marginal prior measure of $\theta$:
\[
2.5 \den{\theta_1} \geq \den{\theta_4} \geq 2 \den{\theta_1},
\hspace{1cm}
10 \den{\theta_3} \geq \den{\theta_2} \geq 9 \den{\theta_3},
\hspace{1cm}
\den{\theta_2} = 5 \den{\theta_4}.
\]
Suppose the following lower and upper likelihoods are given
for a measurement $x$:

\begin{center}
\begin{tabular}{cccc}
$L_x(\theta_1) = 0.9$, &
$L_x(\theta_2) = 0.1125$, &
$L_x(\theta_3) = 0.05625$, &
$L_x(\theta_4) = 0.1125$, \\
$U_x(\theta_1) = 0.95$, &
$U_x(\theta_2) = 0.1357$, &
$U_x(\theta_3) = 0.1357$, &
$U_x(\theta_4) = 0.1357$.
\end{tabular}
\end{center}
\end{Example}

Consider the calculation of the lower density
$\lden{\theta_1|x} =
	\min_{\alpha', \alpha''}
		\left( 0.9 \alpha'_1 + 0.95 \alpha''_1 \right)$,
where $\alpha'$ and $\alpha''$ are vectors with four elements,
subject to:
\begin{itemize}
\item $\mathbf{C} [ \alpha' + \alpha'' ] \leq 0$
where the matrix $\mathbf{C}$ is obtained as described in Section
\ref{section:PosteriorUpperExpectationPriorConditional};
%\[
%\mathbf{C} = \left[ \begin{array}{cccccccc}
%      -2.5 & 0 & 0 & 1 & -2.5 &  0 & 0 &  1 \\
%      2 & 0 & 0 & -1 & 2 & 0 & 0 & -1 \\
%      0 & -1 & 0 & 5 & 0 & -1 & 0 & 5 \\
%      0 & 1 & 0 & -5 & 0 & 1 & 0 & -5 \\
%      0 & -1 & 9 & 0 & 0 & -1 & 9 & 0 \\
%      0 & 1 & -10 & 0 & 0 & 1 & -10 & 0 \end{array} \right].
%\]
\item $[0.9, 0.1125, 0.0562, 0.1125] \alpha' +
	[0.95, 0.1357, 0.1357, 0.1357] \alpha'' = 1$;
\item $\alpha'_i \geq 0$ and $\alpha''_i \geq 0$.
\end{itemize}
The lower density $\lden{\theta_1|X} = 0.2881$ is obtained through
linear programming. The minimizing $\alpha'$ is $[0.3201, 0, 0, 0]$
and the minimizing $\alpha''$ is $[0, 4.0013, 0.4446, 0.8003]$.

The bounds obtained through linear fractional programming
are only valid if the conditional credal sets are specified
separately for each value of $\theta$. White's original
example specified the conditional probability measures through
linear inequalities:
\[
\den{x|\theta_2} = \den{x|\theta_4},
\hspace{1cm}
\den{x|\theta_3} \leq \den{x|\theta_2} \leq 2 \den{x|\theta_3},
\]
\[
7 \den{x|\theta_2} \leq \den{x|\theta_1} \leq 8 \den{x|\theta_2},
\hspace{1cm}
\den{x|\theta_3} \geq 0.01,
\hspace{1cm}
0.9 \leq \den{x|\theta_1} \leq 0.95.
\]
In this case the bounds produced by linear fractional
programming are not tight, because
Theorem \ref{theorem:LowerUpperLikelihoods} does not apply.

\begin{Example}
\label{GaussianExample}
Consider a discrete variable $\theta$ ranging from 120 to 180
in unitary increments. A beta distribution is used to specify
a prior probability mass function:
\[
r(\theta_i) = \gamma \mathbf{Beta}
		\left( \frac{ \theta_i - 120 }{60}, 15, 11 \right),
\]
where $\gamma$ is a constant to ensure that $\sum_i r(\theta_i) = 1$.
Suppose that $r(\theta)$ is not deemed a reliable model and an
$\epsilon$-contaminated model is taken with $\epsilon = 0.15$.
The prior credal set $K(\theta)$ is defined by all densities of the form
$\den{\theta_i} = (1-\epsilon) r(\theta_i) + \epsilon q(\theta_i)$,
where $q(\theta)$ is an arbitrary probability mass function.

Several measurements $X_j$ of the variable $\theta$ are taken.
All measurements are independent and identically modeled by
a set of likelihoods $K(X_j|\theta)$ defined by all densities
of the form:
\[
\den{X_j|\theta_i} =
    \frac{1}{\sqrt{2 \pi \sigma}} \exp
		\left( - \frac{(X_j - \theta_i)^2}{2 \sigma} \right), \;
		\mbox{ where } \sigma \in [ 10, 11 ].
\]
\end{Example}

Consider the calculation of $\lex{\theta|x_1, \dots, x_n}$ and
$\uex{\theta|x_1, \dots, x_n}$. Two steps must be taken: (i) formulation
of constraints for prior measures, and (ii) calculation of lower and
upper likelihoods.

The constraints that must be satisfied by any prior
measure in $K(\theta)$ can be formulated as:
\[
\den{\theta_i} \geq (1-\epsilon) r(\theta_i),
\hspace{1cm}
\den{\theta_i} \geq 0,
\hspace{1cm}
\sum_i \den{\theta_i} = 1.
\]

Any likelihood function
$\den{x_j|\theta_i} \propto \exp(-(x_j-\theta_i)^2/(2\sigma))/\sqrt{\sigma}$
has a single maximum at $\sigma = (x_j - \theta_i)^2$.
Consequently, to find the maximum and minimum of each likelihood
$\den{x_j|\theta_i}$, it is necessary to check the values of the likelihood
for $\sigma = (x_j - \theta_i)^2$, $\sigma = 10$ (minimum value of
variance) and $\sigma = 11$ (maximum value of variance).
The full lower likelihood $L_{x_1, \dots, x_n}(\theta)$ is
the product of all lower bounds on the individual
likelihoods (Section \ref{section:PosteriorUpperExpectationPriorConditional});
likewise, the full upper likelihod $U_{x_1, \dots, x_n}(\theta)$ is
the product of all upper bounds on the individual likelihoods.

Consider, for example, the sequence of 20 measurements in
Table \ref{MeasurementsBounds}. The table also shows the lower
and upper posterior expectations, $\lex{\theta|x_1, \dots, x_n}$ and
$\uex{\theta|x_1, \dots, x_n}$, obtained with these measurements
through linear fractional programming.

\section{CONCLUSION}
\label{section:Discussion}

Two new technical results are noteworthy in this paper. First,
a unified perspective for Lavine's algorithm and the White-Snow
algorithm is derived, based on linear fractional programming.
Related algorithms derived in Artificial Intelligence and robust
Statistics are brought together in this manner, despite the
fact that their motivations have stemmed from disparate applications.

Second, the main contribution is the algorithm for calculation
of posterior quantities given separately specified
prior and conditional credal sets. This algorithm and its extension
to independent identically distributed measurements provides a complete
solution for robust Bayesian inferences based on separately specified
credal sets.

\appendix

\section*{APPENDIX: SOFTWARE FOR INFERENCES}
\label{appendix:Software}

The {\em Matlab}$^{TM}$ procedures used in the paper are
publicly available in the Internet
at the address http://www.cs.cmu.edu/\~{}qbayes/RobustInferences/Matlab/;
information about {\em Matlab}$^{TM}$ can be found in the Internet
at http://www.matlab.com.
The following files are available:
(i) {\tt q1.m} is a procedure that generates bounds from inequalities
in the prior but uses a single likelihood function; (ii) {\tt q2.m}
is a procedure  that generates bounds from inequalities in the
prior and bounds in the likelihood functions; (iii) {\tt cond.m}
is a procedure that generates likelihood bounds for a
set of Gaussian distributions with a range of variances;
(iv) {\tt ws.m} contains the matrices and operations that
solve Example \ref{WhiteSnowExample}; and (v) {\tt ex.m} contains
the matrices and operations that solve Example \ref{GaussianExample}.

\section*{Acknowledgements}

I thank Eric Krotkov for substantial support during the research
that led to this work, Teddy Seidenfeld for introducing me to
the theory of convex sets of probability measures, Peter Walley
for many important suggestions, and
Jay Kadane and Paul Snow for their insightful comments regarding the work.
One of the refrees pointed out the work of \citeN{Betro96},
and another refree indicated the relationship between Lavine's algorithm
and parametric linear programming.

\bibliographystyle{chicago}

\bibliography{../Isipta/consulted}

\clearpage

\begin{table}[h]
{\small
\begin{center}
\begin{tabular}{cc}
\multicolumn{2}{c}{Measurements} \\ \hline
156.2684 & 161.2516 \\
150.7508 & 157.2864 \\
153.5161 & 126.2255 \\
143.0349 & 147.2622 \\
166.9614 & 146.7706 \\
150.5906 & 153.1799 \\
157.5622 & 144.8883 \\
154.0049 & 149.9796 \\
136.5862 & 166.0651 \\
153.7504 & 158.4765
\end{tabular}
\hspace{1.5cm}
\begin{tabular}{lc}
       & Lower and upper posterior expectations \\ \hline
$\theta$ &                     [ 141.5769, 150.5769 ] \\
$\theta|x_1$ &                  [ 147.0183, 151.1553 ] \\
$\theta|x_1, x_2$ &             [ 147.5942, 150.8717 ] \\
$\theta|x_1, \dots, x_3$ &      [ 148.3935, 151.7279 ] \\
$\theta|x_1, \dots, x_4$ &      [ 147.4245, 150.2526 ] \\
$\theta|x_1, \dots, x_5$ &      [ 149.4921, 153.6234 ] \\
$\theta|x_1, \dots, x_{10}$ &   [ 149.1507, 153.1745 ] \\
$\theta|x_1, \dots, x_{15}$ &   [ 147.9065, 152.5462 ] \\
$\theta|x_1, \dots, x_{20}$ &   [ 148.7107, 153.8154 ]
\end{tabular}
\end{center}}
\caption{Twenty measurements and lower and upper posterior expectations.}
\label{MeasurementsBounds}
\end{table}

\end{document}


