\documentstyle[probs,jair,twoside,11pt,myepsf,theapa]{article}

%
%
%  NOTE:  the TR version has an index and table of contents
%     %TR - precedes parts for the TR version
%     %JAIR - preceds parts for the JAIR version
%

\newcommand{\incomp}[2]{ #1 \,:\, #2 }
\newcommand{\Tau}{T}
%TR \newcommand{\myindex}[1]{{\em #1\index{#1}}}
%TR \makeindex
\newcommand{\myindex}[1]{{\em #1}}  %JAIR
\newcommand{\hlinewithspace}
    {\rule[-5pt]{0cm}{1pt}\hline\rule[12pt]{0cm}{1pt}}
\newcommand{\newlinewithspace}
    {\rule[-5pt]{0cm}{1pt}\\\rule[12pt]{0cm}{1pt}}
\newcommand{\newlinewithbigspace}
    {\rule[-8pt]{0cm}{1pt}\\\rule[17pt]{0cm}{1pt}}
\jairheading{2}{1994}{159-225}{4/94}{12/94}  %JAIR
\ShortHeadings{Learning with Graphical Models}{Buntine}
\firstpageno{159}  %JAIR

\begin{document}
%\setstretch{1.2}

\ShortHeadings{Learning with Graphical Models}{Buntine}

\title{\bf Operations for Learning with Graphical Models}
\author{\name Wray L. Buntine \email wray@kronos.arc.nasa.gov \\
        \addr  RIACS \& NASA Ames Research Center, Mail Stop 269--2 \\
         Moffett Field, CA 94035--1000, USA } 
       

% \date{\vskip 15mm{\normalsize (This version: \today)}}
\maketitle

\begin{abstract}
This paper is a multidisciplinary review of empirical, statistical
learning from a graphical model perspective.
Well-known examples of graphical models include Bayesian networks,
directed graphs representing a Markov chain,
and undirected networks representing a Markov field.
These graphical models are extended to model
data analysis and empirical learning using the notation of plates.
Graphical operations for simplifying and manipulating
a problem are provided including decomposition,
differentiation,  and the manipulation of probability
models from the exponential family.
Two standard algorithm schemas for learning are reviewed 
in a graphical framework: Gibbs sampling and the expectation
maximization algorithm.
Using these operations and schemas, some popular algorithms
can be synthesized from their graphical specification.
This includes versions of linear regression,
techniques for feed-forward networks, and learning
Gaussian and discrete Bayesian networks from data.
The paper concludes by sketching
some implications for data analysis
and summarizing how some popular algorithms
fall within the framework presented.

The main original contributions here are the decomposition
techniques and the demonstration that
graphical models provide a framework
for understanding and developing complex learning algorithms.
\end{abstract}

%TR \newpage
%TR \tableofcontents

%TR \newpage

\section{Introduction}

A probabilistic graphical model
is graph where the nodes represent variables
and the arcs (directed or undirected)
represent dependencies between variables.
They are used to define a mathematical
form for the joint or conditional probability 
distribution between variables.
Graphical models come in various forms:
 Bayesian networks used to represent
causal and probabilistic processes, 
data-flow diagrams used to represent
deterministic computation,
influence diagrams used to represent decision processes, 
and undirected Markov networks (random fields) used to represent
correlation for images and hidden causes.

Graphical models 
are used in domains such as diagnosis,
probabilistic expert systems,
and, more recently, in 
planning and control \cite{Dean.wellman.91,Chan.shachter.UAI92},
dynamic systems and time-series \cite{Kjaerulff.UAI92,Dagum.galper.ksl93},
and general data analysis 
\shortcite{Gilks.jrssb} and statistics
\cite{Whittaker.nets}.
This paper shows the task of learning can also be modeled
with graphical models.
This meta-level use of graphical models
was first suggested by Spiegelhalter and Lauritzen 
\citeyear{Spiegelhalter.laur.net} 
in the context of learning probabilities for Bayesian networks.

Graphical models provide a representation for the decomposition
of complex problems.  They also have an associated
set of mathematics and algorithms
for their manipulation.
When graphical models are discussed,
both the graphical formalism and the associated algorithms 
and mathematics are implicitly included.
In fact, the graphical formalism is unnecessary for the technical 
development of the approach, but its use conveys
the important structural information of a problem in
a natural visual manner.
Graphical operations manipulate the underlying structure of
a problem
unhindered by the fine detail of the connecting functional
and distributional equations.
This structuring process is
important in the same way that a high-level programming language
leads to higher productivity over assembly language. 

A graphical model can be developed to 
represent the basic prediction done by linear regression,
a Bayesian network for an expert system, a hidden Markov model,
or a connectionist feed-forward network
\cite{Buntine.rep}.
A graphical model can also be used
to represent and reason about the task of 
{\em learning} the parameters,
weights, and structure of each of these representations.
An extension of the standard graphical model
 that allows this kind of
learning to be represented  is used here.
The extension is the notion of a {\em plate} introduced by 
Spiegelhalter\footnote{The notion of a ``replicated node''
was my version of this developed independently.
I have adopted the notation of Spiegelhalter and colleagues for uniformity.}
\citeyear{Spiegelhalter.comm}.
Plates allow samples to be explicitly represented on the
graphical model, and thus reasoned about and manipulated.
This makes data analysis problems explicit in much the same
way that utility and decision nodes are used for decision
analysis problems \cite{Shachter.opres}.

This paper develops a framework in which the basic computational 
techniques for learning can be directly
applied to graphical models.
This forms the basis of a computational theory of Bayesian learning
using the language of graphical models.
By a {\em computational theory} we mean that the approach 
shows how 
a wide variety of learning algorithms can be created from
graphical specifications and a few simple algorithmic criteria.
The basic computational techniques of probabilistic 
(Bayesian) inference used in this computational theory
of learning are widely reviewed
\cite{Tanner.93,Press.89,Kass.Raftery.571,Neal.review,Bretthorst.MAXENT93}.
These include various exact methods,
Markov chain Monte Carlo methods such as Gibbs sampling,
the expectation maximization (EM) algorithm,
and the Laplace approximation.  
More specialized computational techniques also exist for
handling missing values \cite{Little.Rubin.87},
making a batch algorithm incremental, and
adapting an algorithm to handle large samples.
With creative combination, these techniques are able to address
a wide range of data analysis problems.

The paper provides the blueprint
for a software toolkit that can be used to
construct many data analysis and learning algorithms
based on a graphical specification.
The conceptual  architecture for such a toolkit
is given in Figure~\ref{fig-compile}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{compilealg.ps}{5in}
\end{center}
\vskip -0.2in
\caption{A software generator}
\label{fig-compile}
\end{figure}
Probability and decision theory are used to decompose
a problem into a computational prescription,
and then search and optimization techniques are used to
fill the prescription.
A version of this toolkit already exists
using Gibbs sampling as the general computational
scheme \shortcite{Gilks.etal.stat}.
The list of algorithms that can be constructed in one form
or another by the scheme in Figure~\ref{fig-compile} is impressive.
But the real gain from the scheme
does not arise from the potential re-implementation of existing
software, but from understanding gained by putting
these in a common language,
the ability to create novel hybrid algorithms,
and the ability to tailor special purpose algorithms for
specific problems.

This paper is tutorial in the sense that
it collects material from different communities
and presents it in the language of graphical models.
This paper introduces graphical models, to represent
first-order inference and learning.
Second, this paper develops and reviews a number of operations
on graphical models.
Finally, this paper gives some examples of developing
learning algorithms using combinations
of the basic operations.
The main original contribution is
the demonstration that
graphical models provide a framework
for understanding and developing complex learning algorithms.

In more detail, the paper covers the following topics.
\begin{description}
\item[Introduction to graphical models: ]
Graphical models are used in two ways.
\begin{description}
\item[Graphical models for representing inference tasks: ]
Section~\ref{sct-graph} reviews some
basics of graphical models.
Graphical models provide a means of representing
patterns of inference.
\item[Adapting graphical representations to represent learning: ]
Section~\ref{sct-bl} discusses the representation
of the problem of learning within graphical models
using the notion of a plate.
\end{description}
\item[Operations on graphical models: ]
Operations take a graphical representation
of a learning problem and simplify it
or perform an important calculation required to solve the
problem.
\begin{description}
\item[Operations using closed-form solutions: ]
Section~\ref{sct-basicop} covers those classes of learning
problems where closed-form solutions to learning are known.
This section adapts standard statistical methods to graphical models.
\item[Other basic operations: ]
Other operations on graphs are required to be able to
handle more complex problems.  These are covered in
Section~\ref{sct-oper}, including:
\begin{description}
\item[Decomposition:]  Breaking a learning
problem into independent components and evaluating
each component.
\item[Differentiation:]   Computing the derivative of a probability
or log probability with respect to variables on the graph.
Differentiation can be decomposed into
operations local to  groups of nodes in the graph
as is popular in neural networks.
\end{description}
\item[Some approximate operations: ]
Some approximate algorithms follow naturally
from the above methods.
Section~\ref{sct-approx} reviews Gibbs sampling and 
its deterministic cousin, the EM algorithm.
\end{description}
\item[Some example algorithms: ]
The closed-form solutions to learning can sometimes be used to
form a fast inner loop of more complex algorithms.
Section~\ref{ssct-partial} illustrates how graphical 
models help here.
\end{description}
The conclusion lists some common algorithms and their
derivation within the above framework.
Proofs of lemmas and theorems are collected in 
Appendix~\ref{app-proofs}.
%TR , and
%TR an index to major definitions appears on the last page.

\section{Introduction to graphical models}
\label{sct-graph}

This section introduces graphical models.
The brief tour is necessary before introducing
the operations for learning.

Graphical models offer  a unified
qualitative and quantitative
framework for representing and reasoning with probabilities
and independencies.  
They combine a representation for uncertain problems 
with techniques for performing inference.
Flexible toolkits and systems exist for applying these techniques
\cite{Srinivas.breeze.UAI90,Andersen.etal.IJCAI89,Cowell.BayesS4}.
Graphical models are based on the notion of independence,
which is worth repeating here.
\begin{defn}
$A$ is independent of $B$ given $C$
if 
$p(A,B|C) = p(A|C) p(B|C)$ whenever $p(C)\neq 0$,
for all $A,B,C$.
\end{defn}
The theory of independence as a basic tool for
knowledge structuring is developed by Dawid
\citeyear{Dawid79} and Pearl \citeyear{Pearl.88}.
A graphical model can be equated with the
set of probability distributions that satisfy its implied
constraints.
Two graphical models are \myindex{equivalent probability models}
if their corresponding sets of satisfying probability distributions
are equivalent.


\subsection{Directed graphical models}
\label{sct-dags}

The basic kind of graphical model is the Bayesian network,
also called belief net,
which is most popular in artificial intelligence.
See Charniak \citeyear{Charniak.aimag},
Shachter and Heckerman \citeyear{Shachter.aimag},
and Pearl \citeyear{Pearl.88} 
for an introduction.
This is also a graphical representation for
a Markov chain.
A Bayesian network is a graphical model that uses
directed arcs exclusively to form a directed acyclic graph (DAG),
(i.e., a directed graph without directed cycles).
Figure~\ref{fig-smp}, adapted from \cite{Shachter.aimag}
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{disease.ps}{3in}
\end{center}
\vskip -0.2in
\caption{A simplified medical problem}
\label{fig-smp}
\end{figure}
shows a simple Bayesian network
for a simplified medical problem.
The graphical model represents a
conditional decomposition
of the joint probability
(see \cite{Laur.dawid.networks}
for more details and interpretations).
This decomposition works as follows
(full variable names have been abbreviated).
\begin{eqnarray}
\label{eq-bayesde}
\lefteqn{p(Age, Occ, Clim, Dis, Symp|M)  ~=} \\
&&      
\nonumber p(Age|M) \, p(Occ|M) \, p(Clim|M) \, p(Dis | Age, Occ, Clim,M) \,
        p(Symp | Dis,M) ~,
\end{eqnarray}
where $M$ is the conditioning context, for instance the expert's prior
knowledge and the choice of the graphical model in Figure~\ref{fig-smp}.
Each variable is written conditioned on its \myindex{parents},
where $parents(x)$ is the set of variables with a directed arc into $x$.
The general form for this equation for a set of variables $X$ is:
\begin{equation}
\label{eq-dags}
p(X|M) ~=~ \prod_{x \in X} p(x|parents(x),M)~.
\end{equation}
This equation is the \myindex{interpretation} of a Bayesian
network used in this paper.

\subsection{Undirected graphical models}

Another popular form of graphical model is an undirected
graph, sometimes called a Markov network \cite{Pearl.88}.
This is a graphical model for a Markov random field.
Markov random fields became used in statistics with
the advent of the Hammersley-Clifford theorem
\cite{Besag.aism}.
A variant of the theorem is given later in Theorem~\ref{th-ugs}.
Markov random fields are used in imaging and spatial reasoning
\cite{Ripley.1981,Geman.geman.pami,Besag.aism}
and various stochastic models in neural networks
\cite{Hertz.krogh.palmer}.
Undirected graphs are also important because they simplify
the theory of Bayesian networks \cite{Laur.dawid.networks}.

Figure~\ref{fig-imag} shows a simple $4\times4$  image
and an undirected model for the image.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{image.ps}{4in}
\end{center}
\vskip -0.2in
\caption{A simple $4\times4$ mage and its graphical model}
\label{fig-imag}
\end{figure}
This model is based on the  first degree Markov assumption;
that is, the current pixel is only directly influenced by
pixels positioned next to it, as indicated by the undirected arcs
between variables
$p_{i,j}$ and $p_{i,j+1}$, $p_{i,j}$ and $p_{i+1,j}$, etc.
Each node $x$ (corresponding to a pixel)
has its set of \myindex{neighbors}---those nodes
 it is directly connected to by an
undirected arc.
For instance, the neighbors of $p_{1,1}$ are $p_{1,2}$, $p_{2,2}$
and $p_{2,1}$.
For the variable/node $x$, denote these by $neighbors(x)$.

In general, there is no formula for undirected graphs
in terms of conditional probabilities
corresponding to Equation~(\ref{eq-bayesde}) for the Bayesian network of
Figure~\ref{fig-smp}.
However, a functional decomposition does exist 
in another form based on the maximal cliques
in Figure~\ref{fig-imag}.
\myindex{Maximal cliques} are subgraphs 
that are fully connected but are not strictly contained
in other fully connected subgraphs.
These are the 9 sets of $2\times 2$ cliques
such as $\{ p_{1,2}, p_{1,3}, p_{2,2}, p_{2,3} \}$.
The interpretation of the graph is that the joint probability
is a product over functions of the maximal cliques.
\begin{eqnarray}
\label{eq-imagef}
\lefteqn{p(p_{1,1}, \ldots, p_{4,4}) ~=} \\
&& \nonumber
   f_1(p_{1,1}, p_{1,2}, p_{2,1}, p_{2,2}) \,
        f_2(p_{1,2}, p_{1,3}, p_{2,2}, p_{2,3}) \,
                f_3(p_{1,3}, p_{1,4}, p_{2,3}, p_{2,4}) \\
&& \nonumber
  f_4(p_{2,1}, p_{2,2}, p_{3,1}, p_{3,2}) \,
        f_5(p_{2,2}, p_{2,3}, p_{3,2}, p_{3,3}) \,
                f_6(p_{2,3}, p_{2,4}, p_{3,3}, p_{3,4}) \\
&& \nonumber
  f_7(p_{3,1}, p_{3,2}, p_{4,1}, p_{4,2})  \,
        f_8(p_{3,2}, p_{3,3}, p_{4,2}, p_{4,3})  \,
                f_9(p_{3,3}, p_{3,4}, p_{4,3}, p_{4,4}) ~,
\end{eqnarray}
for some functions $f_1, \ldots, f_9$ defined up to
a constant.
From this formula it follows 
that $p_{1,3}$ is 
conditionally independent of 
its non-neighbors given its neighbors
$p_{1,2},  p_{2,2}, p_{2,3}, p_{1,4}, p_{2,4}$.

The general form for Equation~(\ref{eq-imagef})
for a set of variables $X$ is given 
in the next theorem.
Compare this with Equation~\ref{eq-dags}.
\begin{theorem}
\label{th-ugs}
An undirected graph $G$ is on variables in the set $X$.
The set of maximal cliques on $G$ is $Cliques(G)\subset 2^X$.
The distribution $p(X)$ (probability or probability density)
is strictly positive in the domain $\times_{x\in X} domain(x)$.
Then under the distribution $p(X)$, 
$x$ is independent of
$X-\{x\}-neighbors(x)$ given $neighbors(x)$
for all $x \in X$
(Frydenberg \citeyear{Frydenberg.scand} refers to this
 condition as local G-Markovian)
if, and only if,
$p(X)$ has the functional representation
\begin{equation}
\label{eq-ugs}
p(X) ~=~  \prod_{C \in  Cliques(G)} f_C (C) ~,
\end{equation}
for some functions $f_C>0$.
\end{theorem}
The general form of this theorem for finite discrete
domains is called the Hammersley-Clifford Theorem 
\cite{Geman.book,Besag.aism}.  Again, this equation is
used as the \myindex{interpretation} of a Markov network.

\subsection{Conditional probability models}
\label{ssct-cond}

Consider the conditional probability
$p(Dis | Age, Occ, Clim)$ found in
the simple medical problem
 from Figure~\ref{fig-smp} and Equation~(\ref{eq-bayesde}).
This conditional probability models how the disease should
vary for given values of age, occupation, and climate.
Class probability trees \cite{Cart,Quinlan.92},
graphs and rules \cite{Rivest.mlj,Oliver.93,kohavi-oodg-aaai}, 
and feed-forward networks are representations
devised to express conditional models in different ways.
In statistics, the conditional distributions are also
represented as regression
models and generalized linear models 
\cite{Mccullagh.nelder.glm}.

The models of Figure~\ref{fig-smp} and Equation~(\ref{eq-bayesde}) and
Figure~\ref{fig-imag} and Equation~(\ref{eq-imagef})
show how the joint distribution is composed from simpler components.
That is, they give a global model of variables
in a problem.
The conditional probability models, in contrast,
give a model for a subset of variables conditioned
on knowing the values of another subset.
In diagnosis the concern may be
a particular direction for reasoning,
such as predicting the disease given
patient details and symptoms,
so the full joint model provides unnecessary detail.
The full joint model may
require extra parameters and thus more data to learn.
In supervised learning applications,
the general view is that conditional models are superior unless
prior knowledge dictates a full joint model is more appropriate.
This distinction is sometimes referred to as the
diagnostic versus the discriminant approach to 
classification \cite{Dawid.70s}.

There are a number of ways for explicitly
representing conditional probability models.
Any joint distribution
implicitly gives the conditional distribution for any subset of variables,
by definition of conditional probability.
For instance, if there is a model for $p(Age, Occ, Clim, Dis, Symp)$,
then by the definition of conditional
probability  a conditional model follows:
\begin{eqnarray*}
p(Dis|Age, Occ, Clim, Symp)  &=& \frac{p(Age, Occ, Clim, Dis, Symp)}{
                \sum_{Dis} p(Age, Occ, Clim, Dis, Symp)} ~ \\
& \propto& p(Age, Occ, Clim, Dis, Symp) ~.
\end{eqnarray*}
Conditional distributions can also be represented
by a single node that is labeled
to identify which functional form the node takes.
For instance, in the graphs to follow,
labeled Gaussian nodes, linear nodes, and other
standard forms are all used.
Conditional models such as
rule sets and feed-forward networks can be constructed
by the use of special deterministic nodes.
For instance, Figure~\ref{fig-rules} shows four model
constructs.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{rules.ps}{6in}
\end{center}
\vskip -0.2in
\caption{Graphical conditional models}
\label{fig-rules}
\end{figure}
In each case, the input variables to the conditional models
represented have been shaded.
This shading means the values of these variables are 
{\em known} or given.
In Figure~\ref{fig-rules}(b), the shading indicates
that the value for $x$
is known, but the value for $c$ is unknown.
Presumably $c$ will be predicted using $x$.
Figure~\ref{fig-rules}(a) represents a rule set.
The nodes with double ovals are \myindex{deterministic functions}
of their inputs in contrast to the usual nodes, which
are probabilistic functions.
This means that the value for $rule_1$
is a deterministic function of $x_1, \ldots, x_n$.
Notice that this implies that the
        values for $rule_1$ and the others are known as well.
The conditional probability for a deterministic node,
as required for Equation~(\ref{eq-dags}),
is treated as a delta function.
Figure~\ref{fig-rules}(c) corresponds to the statement:
\[
p(unit|x_1,\ldots,x_n) ~=~ \delta_{unit~=~f(x_1,\ldots,x_n)}
        ~=~\left\{ \begin{array}{lr}
                1 & ~~~~\mbox{if }unit~=~f(x_1,\ldots,x_n)~, \\
                0 & ~~~~\mbox{otherwise~} \end{array} \right.
\]
for some function $f$ not specified.
A Bayesian network constructed entirely of double ovals
is equivalent to a data flow graph where the inputs 
 are shaded.
The analysis of deterministic nodes in Bayesian networks
and, more generally, in influence diagrams
is considered by Shachter \citeyear{Shachter.networks}.
For some purposes, deterministic nodes are best treated
as intermediate variables and removed from the problem.
The method for doing this,
variable elimination,
is given later in Lemma~\ref{lem-nodet}.

The logical or conjunctive form of each rule 
in Figure~\ref{fig-rules}(a)
is not expressed in the graph, and presumably
would be given in the formulas 
accompanying the graph, however the basic functional structure
of the rule set exists.
In Figure~\ref{fig-rules}(b), a node has been labeled with
its functional type.
The functional type for this node with a Boolean variable $c$
 is the function,  
\begin{equation}
\label{eq-log}
p( c =1| x ) ~=~ \frac{ 1 }{ 1 + e^x } ~=~ Sigmoid(x) ~=~Logistic^{-1}(x)
\end{equation}
which maps a real value $x$ onto a probability in $(0,1)$
  for the binary variable $c$.
This function is the inverse of the logistic or logit function
used in generalized linear models \cite{Mccullagh.nelder.glm},
and is also the sigmoid function used in 
feed-forward neural networks.
Figure~\ref{fig-rules}(c) uses a deterministic node to
reproduce a single unit from a connectionist feed-forward network,
where the unit's activation is computed via a sigmoid function.
Figure~\ref{fig-rules}(d) is a simple univariate Gaussian,
which makes $y$ normally distributed with mean $\mu$
and standard deviation $\sigma$.
Here the node is labeled in italics to indicate
its conditional type.

At a more general level, networks can be conditional.
%?  conditional trees ??
Figure~\ref{fig-smpcond}
shows two conditional versions of the simple medical problem.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{diseasecond.ps}{6in}
\end{center}
\vskip -0.2in
\caption{Two equivalent conditional models of the medical problem}
\label{fig-smpcond}
\end{figure}
If the shading of nodes is ignored,
the joint probability, 
$p(Age, Occ, Clim, Dis, Symp)$ for the two graphs (a) and (b) is:
\begin{eqnarray*}
&& p(Age) \, p(Occ|Age) \, p(Clim| Age, Occ ) \, p(Dis | Age, Occ, Clim) 
        \, p(Symp | Age, Dis )  \\
&& p(Age) \, p(Occ) \, p(Clim ) \, p(Dis | Age, Occ, Clim) 
        \, p(Symp | Age, Dis ) ~.
\end{eqnarray*}
However, because four of the five nodes are shaded, this means 
their values are known.
The conditional distributions computed from the above are identical:
\begin{eqnarray*}
\lefteqn{p(Dis | Age, Occ, Clim, Symp)} \\
&=& \frac{ p(Dis | Age, Occ, Clim) \, p(Symp | Age, Dis ) }
         {\sum_{Dis} p(Dis | Age, Occ, Clim) \, p(Symp |Age,  Dis ) } ~. 
\end{eqnarray*}
Why do these distinct graphs become identical when
viewed from the conditional perspective?
Because conditional components of the model corresponding to
age, occupation and climate cancel out when the conditional distribution
is formed.
However, the symptoms node has the unknown variable
disease as a parent, so the arc from age to symptoms is kept.

More generally, the following simple lemma applies and is
derived directly from Equation~\ref{eq-dags}.
\begin{lemma}
\label{lem-remgiven}
Given a Bayesian network $G$ with some nodes shaded representing
a conditional probability distribution,
if a node $X$ and all its parents 
have their values given,
then the Bayesian network $G'$ created by deleting all
the arcs into $X$ represents an equivalent 
probability model to the Bayesian network $G$.
\end{lemma}
This does not mean, for instance in the
graphs just discussed,
that there is no causal or influential links between
the variables age, occupation, and climate,
rather that their effects become irrelevant in the conditional model
considered because their values are already known.
A corresponding result holds for undirected graphs,
and follows directly from Theorem~\ref{th-ugs}.
\begin{lemma}
\label{lem-remgivenug}
Given an undirected graph $G$ with some nodes shaded representing
a conditional probability distribution,
delete an arc between nodes $A$ and $B$ if all their common neighbors
are given.
The resultant graph $G'$ represents an equivalent 
probability model to the graph $G$.
\end{lemma}


\subsection{Mixed graphical models}
\label{ssct-cg}

Undirected and directed graphs can also be mixed in a sequence.
These mixed graphs are called \myindex{chain graphs}
\cite{Wermuth.laur.jrssb,Frydenberg.scand}.
These chain graphs are sometimes used here,
However, a precise understanding of them is
not required for this paper.
A simple chain graph is given in Figure~\ref{fig-smpmx}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{diseasemix.ps}{3.5in}
\end{center}
\vskip -0.2in
\caption{An expanded medical problem}
\label{fig-smpmx}
\end{figure}
In this case, the single disease
node and single symptom node of Figure~\ref{fig-smp}
are expanded 
to represent the case where there are two possibly co-occurring
diseases and three possibly co-occurring symptoms.
The medical specialist may have said something like:
``Lung disease and heart disease can influence
        each other, or may have some hidden common cause;
        however, it is often difficult to tell which is
        the cause of which, if at all.''
In the causal model, join the two disease
nodes by an undirected arc to represent direct influence.
Likewise for the symptom nodes.
The resultant joint distribution takes the form:
\begin{eqnarray}
\label{eq-bayescg}
\lefteqn{p(Age, Occ, Clim, Heart\mbox{-}Dis, Lung\mbox{-}Dis,
        Symp\mbox{-}A, Symp\mbox{-}B, Symp\mbox{-}C)  ~=} \\
&&      
\nonumber p(Age) \, p(Occ) \, p(Clim) \, 
        p(Heart\mbox{-}Dis, Lung\mbox{-}Dis | Age, Occ, Clim) \\
&&      \nonumber p(Symp\mbox{-}A, Symp\mbox{-}B, Symp\mbox{-}C | 
        Heart\mbox{-}Dis, Lung\mbox{-}Dis) 
\end{eqnarray}
where the last two conditional probabilities
can take on an arbitrary form.
Notice that the probabilities now have more than one
variable on the left side.

In general,
a  chain graph consists of a chain of 
undirected graphs connected  by directed arcs.
Any cycle through the graph cannot have directed arcs going in opposite
directions.
Chain graphs can be interpreted as  Bayesian networks defined
over the components of the chain instead of the original 
variables.
This goes as follows:
\begin{defn}
Given a subgraph $G$ over some variables $X$,
the \myindex{chain components} are subsets of $X$ that are
maximal undirected connected subgraphs in a chain graph $G$
\cite{Frydenberg.scand}.
Furthermore, let
$chain\mbox{-}components(A)$ denote all nodes in the same
chain component as at least one variable in $A$.
\end{defn}
The chain components for the graph above,
ordered consistently with the directed arcs,
are $\{Age\}$, $\{Occ\}$, $\{Clim\}$, 
$\{Heart\mbox{-}Dis, Lung\mbox{-}Dis\}$,
     and $\{Symp\mbox{-}A, Symp\mbox{-}B, Symp\mbox{-}C \}$.
Informally, a chain graph over variables $X$ with chain components
given by the set $\Tau$
is interpreted first as the decomposition
corresponding to the decomposition of Bayesian
networks in Equation~(\ref{eq-dags}):
\begin{equation}
\label{eq-chdag}
p(X|M) ~=~ \prod_{\tau \in \Tau} p(\tau|parents(\tau),M)~
\end{equation}
where\index{parents}
\[
parents(A) ~=~ \bigcup_{a \in A} parents(a) - A~.
\]
Each component probability $p(\tau|parents(\tau),M)$
has a form similar to Equation~(\ref{eq-ugs}).

Sometimes, to process graphs of this form without having
to consider the mathematics of chain graphs,
the following device is used.
\begin{comment}
\label{comm-cg}
When a set of nodes $U$ in a chain graph
form a clique (a fully connected subgraph), 
and all have identical children and parents otherwise,
then the set of nodes can be replaced a single node
representing the cross product of the variables.
\end{comment}
This operation for Figure~\ref{fig-smpmx}
is done to get Figure~\ref{fig-smpmxjn}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{diseasemixjn.ps}{3in}
\end{center}
\vskip -0.2in
\caption{An expanded medical problem}
\label{fig-smpmxjn}
\end{figure}
Furthermore, chain graphs are sometimes used
here where Bayesian networks can
also be used.  In this case:
\begin{comment}
\label{comm-cgtodag}
The chain components of a Bayesian network are the singleton sets
of individual variables in the graph.
Furthermore, $chain\mbox{-}components(A)=A$.
\end{comment}

Chain graphs can be decomposed into a chain
of directed and undirected graphs.
An example is given in Figure~\ref{fig-chaing}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{chaing.ps}{5in}
\end{center}
\vskip -0.2in
\caption{Decomposing a chain graph}
\label{fig-chaing}
\end{figure}
Figure~\ref{fig-chaing}(a) shows the original chain graph.
Figure~\ref{fig-chaing}(b) shows its directed and undirected
components together with the Bayesian network on the right
showing how they are pieced together.
Having done this decomposition,
the components are analyzed
using all the machinery of directed and undirected graphs.
The interpretation of these graphs in terms of
independence statements and the implied functional
form of the joint probability is a combination
of the previous two forms given in
Equation~(\ref{eq-dags}) and Theorem~\ref{th-ugs}, based on
 \cite[Theorem~4.1]{Frydenberg.scand},
and on the interpretation of conditional graphical models 
in Section~\ref{ssct-cond}.

\section{Introduction to learning with graphical models}
\label{sct-bl}

A simplified inference problem is represented in
Figure~\ref{fig-sbn}.
Here, the nodes $var_1$, $var_2$ and $var_3$ are shaded.
This represents that the value of these nodes 
is given, so the inference task is to predict
the value of the remaining variable $class$.
This graph matches the so-called ``idiot's'' Bayes classifier 
\cite{Duda.hart:73,langley:aaai92}
used in supervised learning for its speed and simplicity.
The probabilities on this network are easily learned from data
about the three input variables  $var_1$, $var_2$ and $var_3$,
and $class$.
This graph also matches an unsupervised
learning problem where the class $class$ is not in the data
but is hidden.   An unsupervised learning algorithm
learns hidden classes \cite{Cheeseman.AAAI88,McLachlan.basford}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{smplbayes.ps}{1in}
\end{center}
\vskip -0.2in
\caption{A simple classification problem}
\label{fig-sbn}
\end{figure}

The implied joint for these variables read from the graph is:
\begin{equation}
\label{eq-scp}
p(class,var_1,var_2,var_3) ~=~
p(class)\, p(var_1 | class)\, p(var_2 | class)\, p(var_3 | class) ~.
\end{equation}
The Bayesian classifier gets its name because it is derived
by applying Bayes theorem to this joint to get the
conditional formula:
\begin{equation}
\label{eq-sbc}
p(class | var_1,var_2,var_3 ) ~=~
\frac{p(class) p(var_1 | class)\, p(var_2 | class)\, p(var_3 | class)}
{\sum_{class}p(class)  
        p(var_1 | class)\, p(var_2 | class)\, p(var_3 | class) }~.
\end{equation}
The same formula is used to predict the hidden class for
objects in the simple unsupervised learning framework.
Again, this formula, and corresponding formula for
 more general classifiers,
 can be found automatically by using
exact methods for inference on Bayesian networks.

Consider the simple model given in Figure~\ref{fig-sbn}.
If the matching unsupervised learning problem for
this model was represented, a sample of $N$ cases
of the variables would be observed, 
with the first case being $var_{1,1}$, $var_{2,1}$, $var_{3,1}$,
and the $N$-th case being $var_{1,N}$, $var_{2,N}$, $var_{3,N}$.
The corresponding hidden classes, 
$class_1$ to $class_N$, would  not be observed,
but interest would be in performing inference about the parameters
needed to specify the hidden classes.
The learning problem is represented in Figure~\ref{fig-sbnl}.
This includes two added features:
an explicit representation of the model parameters
$\phi$ and $\theta$, and
a representation of the sample as $N$ repeated subgraphs.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{learnbayes.ps}{5in}
\end{center}
\vskip -0.2in
\caption{Learning the simple classification}
\label{fig-sbnl}
\end{figure}
The parameter $\phi$ (a vector of class probabilities)
gives the proportions for the
hidden classes, and the three parameters $\theta_1$, $\theta_2$
and $\theta_3$ give how the variables are distributed within
each hidden class.
For instance,
if there are 10 classes, then
$\phi$ is a vector of 10 class probabilities
such that the prior probability of a case being in class $c$ is
$\phi_c$.
If $var_1$ is a binary variable, then $\theta_1$ would be
10 probabilities, one for each class, such that
if the case is known to be in class $c$,
then the probability $var_1$ is true is given by $\theta_{1,c}$
and the probability $var_1$ is false is given by $1-\theta_{1,c}$.
This yields the following equations:
\begin{eqnarray*}
p(class=c|\phi,M) &=& \phi_c ~,\\
p(var_j=true|class=c,\theta_j,M) &=&  \theta_{j,c} ~.
\end{eqnarray*}
The unknown model parameters
$\phi$, $\theta_1$, $\theta_2$ and $\theta_3$
are included in the
graphical model to explicitly represent
all unknown variables
and parameters in the learning problem.

\subsection{Introduction to Bayesian learning}
\label{ssct-introlearn}

Now is a useful time to introduce the basic terminology of
Bayesian learning theory.
This is not an introduction to the field.
Introductions are given in
\cite{Bretthorst.MAXENT93,Press.89,Loredo.92,Bernardo.Smith.94,Cheeseman.cmdtf}.
This section reviews notions such as the
sample likelihood and Bayes factor, important
for subsequent results.

For the above unsupervised learning 
problem there is the {\em model},
$M$, which is the use of the hidden class
and the particular graphical structure of Figure~\ref{fig-sbnl}.
There are  data assumed to be independently sampled,
and there are the \myindex{parameters} of the model 
($\phi$, $\theta_1$, etc.).
In order to use the theory,
it must be assumed that the model is correct.  That is, the ``true''
distribution for the data can be assumed to come from this model
with some parameters.
In practice, hopefully the model assumptions are sufficiently
close to the truth.  Different sets of model assumptions may
be tried.
Typically, the ``true'' model parameters are unknown, although
there may be some rough idea about their values.
Sometimes, several models are considered
(for instance different kinds of Bayesian networks),
but it is assumed that just one of them is correct.
Model selection or model averaging
methods are used to deal with them.

For the Bayesian classifier above, 
a subjective probability is placed over the 
model parameters, in the form 
$p(\phi, \theta_1, \theta_2, \theta_3|M)$.
This is called the \myindex{prior probability}.
Bayesian
statistics and decision theory is distinguished from all other
statistical approaches in that it places 
initial probability distributions, the prior probability,
over unknown model parameters.  
If the model is a feed-forward neural network,
then a prior probability needs to be placed over the
network weights and the standard deviation of the error.
If the model is linear regression with Gaussian error, then it
is over the linear parameters $\theta$ 
and the standard deviation of the error.
Prior probabilities are an active area of research
and are discussed in most introductions to Bayesian methods.

The next important component is the sample \myindex{likelihood},
which, on the basis of the model assumptions
$M$ and given a set of parameters 
$\phi$, $\theta_1$, $\theta_2$ and $\theta_3$,
says how likely the sample of data was.
This is
$p( sample |\phi, \theta_1, \theta_2, \theta_3,M)$.
The model needs to completely determine
the sample likelihood.
The sample likelihood is the basis of
the maximum likelihood principle
and many hypothesis testing methods \cite{Casella.berger.90}.
This combines with the prior to form the
\myindex{posterior probability}:
\[
p(\phi, \theta_1, \theta_2, \theta_3 | sample, M) ~=~
        \frac{ p( sample |\phi, \theta_1, \theta_2, \theta_3,M) 
        p(\phi, \theta_1, \theta_2, \theta_3|M) }
        { p(sample | M) } ~.
\]
This equation is Bayes theorem and the term
$p(sample | M)$ is derived from the prior and sample likelihood
using an integration or sum that is often difficult
to do:
\begin{equation}
\label{eq-evid}
p(sample | M)  ~=~\int_{\phi, \theta_1, \theta_2, \theta_3}
        p( sample |\phi, \theta_1, \theta_2, \theta_3,M) 
        p(\phi, \theta_1, \theta_2, \theta_3|M) \,
        \diff (\phi, \theta_1, \theta_2, \theta_3) ~.
\end{equation}
This term is called the \myindex{evidence} for model $M$,
or \myindex{model likelihood},
and is the basis for most Bayesian model selection,
model averaging methods,
and Bayesian hypothesis testing
methods using Bayes factors 
\cite{Smith.etal.jrssb,Kass.Raftery.571}.
The \myindex{Bayes factor} is a relative quantity used to compare one model
$M_1$ with another $M_2$:
\[
Bayes\mbox{-}factor(M_2,M_1) ~=~ 
        \frac{p(sample |M_2)}{p(sample | M_1)} ~.
\]
Kass and Raftery \citeyear{Kass.Raftery.571}
review the large variety of methods available
for computing or estimating the evidence for a model
including numerical integration,
importance sampling, and the Laplace approximation.
In implementation, the log of the Bayes factor is used
to keep the arithmetic within reasonable bounds.
The log of the evidence can still produce large numbers,
and since rounding errors in floating point arithmetic scales with
the order of magnitude, the log Bayes factor is the preferred
quantity to consider in implementation.
The evidence is often simpler in mathematical analysis.
The Bayes factor is the Bayesian equivalent to the 
likelihood ratio test used in orthodox statistics
and developed by Wilks.  
See Casella and Berger \citeyear{Casella.berger.90} 
for an introduction
and Vuong \citeyear{Vuong.econo} for a recent review.

The evidence and Bayes factors
are fundamental to Bayesian methods.
It is often the case that a complex
``non-parametric'' model (a statistical term that
loosely translates as ``many and varied parameter'' model) be used
for a problem, rather than a simple
model with some fixed number of parameters.
Examples of such models are decision trees, most neural networks, and
Bayesian networks.
For instance, suppose two models are proposed with
$M_1$ and $M_2$ being two Bayesian networks suggested by the domain expert.
These are given in Figure~\ref{fig-mods1and2}.
They are over two multinomial variables
$var_1$ and $var_2$ and two Gaussian variables $x_1$ and $x_2$.
Model $M_2$ has an additional arc going from the
discrete variable $var_2$
to the real valued variable $x_1$.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{mods1and2.ps}{5in}
\end{center}
\vskip -0.2in
\caption{Two graphical models.  Which should learning select?}
\label{fig-mods1and2}
\end{figure}

The parameters $\theta_{1},\theta_2,\mu_1,\sigma_1,\mu_2,\sigma_2$
 for model $M_1$ parameterize probability distributions
 for the first Bayesian network, 
and the parameters $\theta_{1},\theta_2,\mu_1',\sigma_1',\mu_2,\sigma_2$
for the second.  The task is to learn not only a set of
parameters, but also to select a Bayesian network from the two.
The Bayes factor gives the comparative worth of the two models.
This simple example extends in principle to selecting a single
decision tree, rule set, or Bayesian network from the
huge number available from attributes in the domain.
In this case compare the posterior probabilities
of the two models $p(M_1|sample)$ and $p(M_2|sample)$.
Assuming the truth falls in one or other model,
the first is computed using Bayes theorem as:
\begin{eqnarray*}
p(M_1|sample) &=&  \frac{ p(sample | M_1) p(M_1) }
                        {p(sample | M_1) p(M_1) +p(sample | M_2) p(M_2) } ~,\\
        &=& \frac{ 1 }{ 1 + Bayes\mbox{-}factor(M_2,M_1) 
                \frac{p(M_2)}{p(M_1) } } ~.
\end{eqnarray*}
More generally, when multiple models exist, it still 
holds that:
\[
\frac{p(M_2|sample)}{p(M_1|sample)} ~=~Bayes\mbox{-}factor(M_2,M_1) 
                \frac{p(M_2)}{p(M_1) }.
\]
Notice that the computation requires of each model its prior
and its evidence.
The second form reduces the computation to the relative quantity 
being the Bayes factor, and a ratio of the priors.
Bayesian hypothesis testing corresponds to checking if the
Bayes factor of the null hypothesis compared to the
alternative hypothesis is very small or very large.
Bayesian model building corresponds to searching for a model
with a high value of $p(M|sample) \propto p(sample | M) p(M)$,
which usually involves comparing Bayes factors of this model
with alternative models during the search.

When making an estimate about a new case $x$,
the estimate becomes:
\begin{eqnarray*}
\lefteqn{p(x|sample,\{M_1,M_2\}) }\\
&=&  p(M_1|sample) \, p(x|sample,M_1) +
p(M_2|sample) \,p(x|sample,M_2) \\
&=&
\frac{p(sample|M_1)\,p(M_1) \,p(x|sample,M_1)+
p(sample|M_2)\,p(M_2)\, p(x|sample,M_2)}
  {p(sample|M_1)\,p(M_1) + p(sample|M_2)\,p(M_2)} ~.
\end{eqnarray*}
The predictions of the individual models is averaged
according to the model posteriors $p(M_1|sample)$
and $p(M_2|sample) = 1 - p(M_1|sample)$.
The general components used in this calculation are
the model priors, the evidence for each model or the Bayes factors,
and the prediction for the new case made for each model.

This process of \myindex{model averaging}
happens in general.
A typical non-parametric problem would be to learn
class probability trees from data.
The number of class probability tree models is super-exponential 
in the number of features.
Even when learning Bayesian networks from data the
number of Bayesian networks is at best quadratic  in the number of features.
Doing an exhaustive search of these spaces
and doing the full averaging implied by the equation above
is computationally infeasible in general.
It may be the case that 15 models have posterior
probabilities $p(M|sample)$ between 0.1 and 0.01,
and several thousand more models have posteriors from
0.001 to 0.0000001.
Rather than select a single model,
a representative set of several models might be chosen
and averaged using the identity:
\[
p(x|sample) ~=~ \sum_{i}  p(M_i|sample) \, p(x|sample,M_i) ~.
\]
The general averaging process is depicted in Figure~\ref{fig-multnets}
where a Gibbs sampler is used to generate a representative
subset of models with high posterior.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{multnetsM.ps}{5in}
\end{center}
\vskip -0.2in
\caption{Averaging over multiple Bayesian networks}
\label{fig-multnets}
\end{figure}
This kind of computation is done for class probability
trees where representative sets of trees are found using
a heuristic branch and bound algorithm \cite{Buntine.AISIII},
and for learning Bayesian networks \cite{Madigan.raftery.jasa}.
A sampling scheme for Bayesian networks  is presented in
Section~\ref{ssct-lbn}.

\subsection{Plates:  representing learning with graphical models}
\label{ssct-plates}

In their current form, graphical models do not
allow the convenient representation of a learning problem.
There are four important points to be observed regarding the
use of graphical models to improve their
suitability for learning.
Consider again the unsupervised learning system
described in the introduction of Section~\ref{sct-bl}.
\begin{itemize}
\item
The unknown model parameters
$\phi$, $\theta_1$, $\theta_2$, and $\theta_3$
are included in the
graphical model to explicitly represent
all variables in the problem, even model parameters.
By including these in the probabilistic
model, an explicitly Bayesian model is constructed.
Every variable in a graphical model,
even unknown model parameters,
has a defined prior probability.
\item
The learning sample is a repeated set
of measured variables so the
basic model of Figure~\ref{fig-sbn}
appears duplicated as many times as there are
cases in the sample,
as shown in Figure~\ref{fig-sbnl}.
Clearly, this awkward repetition will occur whenever homogeneous
data is being modeled (typical in learning).
Techniques for handling this repetition form a major
part of this paper.
\item
Neither graph in
 Figures~\ref{fig-sbn} and~\ref{fig-sbnl} represents the
goal of learning.
For learning to be goal directed,
additional information needs to be included in the
graph:    how is learned knowledge evaluated or
how can subsequent performance be measured?
This is the role of decision theory  and it
is modeled in graphical form using influence diagrams
\cite{Shachter.opres}.
This is not discussed here,
but is covered in \cite{Buntine.rep}.
\item
Finally, it must be possible to take a graphical representation
of a learning problem and the goal of learning and
construct an algorithm to solve the problem.
Subsequent sections discuss techniques for this.
\end{itemize}

Consider a simplified version of the same unsupervised problem.
In fact, the simplest possible learning problem 
containing uncertainty goes as follows:
there is a biased coin with an unknown bias for heads $\theta$.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{coin1.ps}{4in}
\end{center}
\vskip -0.2in
\caption{Tossing a coin:   model without and with a plate}
\label{fig-coin}
\end{figure}
That is, the long-run frequency of getting heads for this coin
        on a fair toss is $\theta$.
The coin is tossed $N$ times and each time the binary
variable $heads_i$ is recorded.
The graphical model for this is in Figure~\ref{fig-coin}(a).
The $heads_i$ nodes are shaded because 
their values are  given, but the $\theta$ node is not.
The $\theta$ node has a $Beta(1.5,1.5)$ prior.
This assumes $\theta$ is distributed according to the
Beta distribution with parameters $\alpha_1=1.5$ and $\alpha_2=1.5$,
\begin{equation}
\label{eq-bino}
p(\theta|\alpha_1,\alpha_2) ~=~
        \frac{\theta^{\alpha_1-1} \, (1-\theta)^{\alpha_2-1}}
        {Beta(\alpha_1,\alpha_2)} 
\end{equation}
where $Beta(,)$ is the standard beta function 
given in many mathematical tables.
This prior is plotted in Figure~\ref{fig-betaprior}.  
A $Beta(1.0,1.0)$ prior, for instance, is uniform in $\theta$,
whereas $Beta(1.5,1.5)$ slightly favors values closer to 0.5---a fairer coin.
\begin{figure}[htbp]
\hspace*{-0.5in}\myepsf{mybeta.ps}   % ???
\vskip -0.2in
\caption{The $Beta(1.5,1.5)$ prior ($\alpha=1.5$)
  and other priors on $\theta$}
\label{fig-betaprior}
\end{figure}
Figure~\ref{fig-coin}(b) is an equivalent graphical model
using the notation of plates.
The repeated group, in this case the $heads_i$ nodes,
is replaced by a single node with a box around it.
The box is referred to as a  \myindex{plate}.
and implies that
\begin{itemize}
\item
 the enclosed subgraph is duplicated $N$ times
(into a ``stack'' of plates),
\item
the enclosed variables are indexed, and
\item
any exterior-interior links are duplicated.
\end{itemize}

In Section~\ref{sct-dags} it was shown that any Bayesian network
has a corresponding form for the joint probability
of variables in the Bayesian network.
The same applies to plates.
The plate indicates that a product ($\prod$) will appear
in the corresponding form.
The probability equation for
Figure~\ref{fig-sbnl}, read directly from the graph, is:
\begin{eqnarray*}
\lefteqn{p(\phi, \theta_1, \theta_2, \theta_3, class_1, 
        var_{1,1}, var_{2,1}, var_{3,1}, \ldots, 
        class_N, 
        var_{1,N}, var_{2,N}, var_{3,N}) ~=~} \\
&& p(\phi) \, p( \theta_1) \, p( \theta_2) \, p( \theta_3) \,
        p(class_1| \phi) \, p(var_{1,1}|class_1,\theta_1)
        \, p(var_{2,1}|class_1,\theta_2)\, p(var_{3,1}|class_1,\theta_3)
         \\
&&    ~ \ldots ~ p(class_N| \phi) \, p(var_{1,N}|class_N,\theta_1)
        \, p(var_{2,N}|class_N,\theta_2)\, p(var_{3,N}|class_N,\theta_3) ~.
\end{eqnarray*}
The corresponding equation using product notation is:
\begin{eqnarray*}
\lefteqn  {p(\phi, \theta_1, \theta_2, \theta_3, 
        class_i, 
        var_{1,i}, var_{2,i}, var_{3,i} \, :\, i=1,\ldots, N) 
        ~=~ p(\phi) \, p( \theta_1) \, p( \theta_2) \, p( \theta_3) } \\
&&         \prod_{i=1}^N
        p(class_i| \phi) \, p(var_{1,i}|class_i,\theta_1)
        \, p(var_{2,i}|class_i,\theta_2)\, p(var_{3,i}|class_i,\theta_3) ~.
\end{eqnarray*}
These two equations are equivalent.
However, the differences
in their written form corresponds to the differences
in their graphical form.
Each plate is converted into a product:
the joint probability ignoring the plates is written,
a product ($\prod$) is added for each plate to index the
variables inside it.
If there are two disjoint plates then there are two disjoint products.
 Overlapping plates yield overlapping products.
The corresponding transformation for the unsupervised
learning problem of Figure~\ref{fig-sbnl}
is given in Figure~\ref{fig-sbnlp}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{auto3p.ps}{2in}
\end{center}
\vskip -0.2in
\caption{Simple unsupervised learning, with a plate}
\label{fig-sbnlp}
\end{figure}
Notice, the hidden $class$ variable is not shaded
so it is not given.
The corresponding transformation for the supervised
learning problem of Figure~\ref{fig-sbnl},
where the classes are given,
and thus corresponds to the idiot's Bayes classifier,
is identical to Figure~\ref{fig-sbnlp} except that the
$class$ variable is shaded because the classes
are now part of the training sample.

Many learning problems 
can be similarly modeled with plates.
Write down the graphical model for the
full learning problem with only a single case provided.
Put a box around the data part of the model,
pull out the model parameters
(for instance, the weights of the network or the classification parameters),
and ensure they are unshaded because they are unknown.
Now add the data set size ($N$) to the bottom left corner.

The notion of a plate is formalized below.
This formalization is included for
use in subsequent proofs.
\begin{defn}
A chain graph $G$ with plates on variable set $X$ consists of a
chain graph $G'$ on variables $X$ with additional boxes 
called plates placed around groups
of variables.  Only directed arcs can cross plate boundaries,
and plates can be overlapping.
Each plate $P$ has an integer $N_P$ in the bottom left corner indicating
its cardinality.
Each plate indexes the variables inside it with values $i=1,\ldots,N_P$.
Each variable $V\in X$ occurs in some subset of the
plates.  Let $indval(V)$\index{$indval$} denote the set of values for
indices corresponding to these plates.
That is,  $indval(V)$ is
the cross product of index sets $\{1,\ldots,N_P\}$ for plates $P$ containing $V$.
\end{defn}

A graph with plates can be expanded to remove the plates.
Figure~\ref{fig-sbnl} is the expanded form of Figure~\ref{fig-sbnlp}.
Given a chain graph with plates $G$ on variables
$X$, construct the \myindex{expanded graph} as follows:
\begin{itemize}
\item
For each variable $V\in X$,
add a node for $V_i$ for each $i\in indvar(V)$.
\item
For each undirected arc between variables $U$ and $V$,
add an undirected arc between $U_i$ and $V_i$ 
for $i\in indvar(V)=indvar(U)$.
\item
For each directed arc between variables $U$ and $V$,
add a directed arc between $U_i$ and $V_j$ for $i\in indvar(V)$
and $j\in indvar(V)$ where $i$ and $j$ have identical values
for index components from the same plate.
\end{itemize}
The parents \index{parents}
for indexed variables in a graph with plates are the parents
in the expanded graph. 
\[
parents(U) ~=~ \bigcup_{i\in indval(U)} parents(U_i)~.
\]

A graph with plates is interpreted using the following product form.
If the product form for the chain graph $G'$ without plates 
with chain components $\Tau$ is
\[
p(X|M(G')) ~=~ \prod_{\tau \in \Tau} p(\tau|parents(\tau),M)~,
\]
then the product form for the chain graph $G$ with plates
has a product for each plate:
\begin{equation}
\label{eq-chgpf}
p(X|M(G)) ~=~ \prod_{\tau \in \Tau} 
      \prod_{i\in indval(\tau)}   p(\tau_i|parents(\tau_i),M)~.
\end{equation}
This is given by the expanded version of the graph.
Testing for independence on chain graphs with
plates involves expanding 
the plates.  In some cases, this can be simplified.


\section{Exact operations on graphical models}
\label{sct-basicop}

This section introduces basic inference methods on graphs
without plates and exact inference methods on
graphs with plates.  While there are no common
machine learning algorithms explained in this section,
the operations explained 
are the mathematical basis of most fast learning
algorithms.  Therefore, the importance of these
basic operations should not be underestimated.
Their use within more well-known learning algorithms
is explained in later sections.

Once a graphical model is developed to represent a problem,
the graph can be manipulated using various
exact or approximate transformations to
simplify the problem.
This section reviews basic exact transformations available:
arc reversal, node removal,
and exact removal of plates by recursive arc reversal.
The summary of operations  emphasizes the
computational aspects.
A graphical model has an associated
set of definitions or tables for the basic functions and conditional
probabilities implied by the graph,
the operations given below effect both the graphical
structure and these underlying mathematical specifications.
In both cases,  the process of making these
transformations should be constructive
so that a graphical specification for a learning problem
can be converted into an algorithm.  

There are several generic approaches for performing
inference on directed and undirected networks without plates.
These approaches are mentioned, but will not be covered in detail.
The first approach is exact and corresponds to removing
independent or irrelevant information from the graph,
then attempting to optimize an exact probabilistic computation
by finding a reordering of the variables.
The second approach to performing inference
is approximate and corresponds to
approximate algorithms such as Gibbs sampling,
and other Markov chain Monte Carlo methods
\cite{Hrycej.aij,Hertz.krogh.palmer,Neal.review}.
In some cases, the complexity of the first approach is inherently
exponential in the number of variables,
so the second can be more efficient.
The two approaches can be combined in some cases
after appropriate reformulation of the problem
\cite{Dagum.horvitz.UAI92}.


\subsection{Exact inference without plates}
\label{ssct-without}

The exact inference approach has been highly refined for
the case where all variables are discrete.
It is not surprising that
 available algorithms have strong similarities
\cite{Shachter.AS.UAI94} since the major choice points
involve the ordering of the summation and
whether this ordering is selected dynamically or statically.
Other special classes of inference algorithms include the cases
where the model is a 
multivariate Gaussian \cite{Shachter.kenley.89,Whittaker.nets},
or corresponds to some specific diagnostic structure,
such as two-level believe networks
with a level of symptoms connected to a level of diseases
\cite{Henrion.idbn}.
This subsection reviews some simple, exact 
transformations on graphical models without plates.
Two representative methods are covered
but are by no means optimal: 
arc reversal and arc removal.
They are important, however, because they are the
building blocks on which methods for graphs with plates are based.
Many more sophisticated variations and combinations
of these algorithms exist in the literature,
including the handling of deterministic nodes \cite{Shachter.networks}
and chain graphs and undirected graphs
\cite{Frydenberg.scand}.

\subsubsection{Arc reversal}

Two basic steps for inference are to marginalize nuisance parameters
or to condition on new evidence.
This may require evaluating probability variables in a different order.
The \myindex{arc reversal} operator 
interchanges the order of two nodes connected
by a directed arc \cite{Shachter.opres}.
This operator corresponds to Bayes theorem and
is used, for instance, to automate the derivation of 
Equation~(\ref{eq-sbc}) from Equation~(\ref{eq-scp}).
The operator applies to directed acyclic graphs
and to chain graphs
where  $a$ and $b$  are  adjacent chain components.

\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{arcrev.ps}{5in}
\end{center}
\vskip -0.2in
\caption{Arc reversal: reversing nodes $a$ and $b$}
\label{fig-arcrev}
\end{figure}
Consider a fragment of a graphical model as given
in the left of Figure~\ref{fig-arcrev}.
The equation for this fragment is:
\[
p( a,b | A) ~=~ p(b|parents(b)) \, p(a|b,parents(a)) ~.
\]
Suppose nodes $a$ and $b$ need to be reordered.
Assume that between $a$ and $b$ there is no directed path 
of length greater than one.
If there was, then reversing the arc between $a$ and $b$ would create
a cycle, which is forbidden in a Bayesian network and chain graph.
The formula for the variable reordering can be found by applying Bayes theorem
to the above equation.
\begin{eqnarray*}
p(a|A) &=& \sum_b p(b|parents(b)) \, p(a|parents(a)) ~, \\
p(b|a,A) &=& \frac{p(b|parents(b)) \, p(a|parents(a))}
                         {\sum_b p(b|parents(b)) \, p(a|parents(a))} ~.
\end{eqnarray*}
The corresponding graph is given in the right
of Figure~\ref{fig-arcrev}.
Notice that the effect on the graph is that nodes for $a$ and $b$ now
share their parents.
This is an important point.
If all of $a$'s parents were also $b$'s, and vice versa,
excepting $b$ ($parents(a) = parents(b) \cup \{b\}$),
then the graph would be unchanged
except for the direction of the arc between $a$ and $b$.
Regardless, the probability tables or formula associated
with the graph also need to be updated.
If the variables are discrete and full conditional probability
tables are maintained, then this operation
requires instantiating the set
$\{a,b\} \cup parents(a) \cup parents(b)$ in all ways,
which is exponential in the number of variables.

\subsubsection{Arc and node removal}
\label{ssct-arcrem}

Some variables in a graph are part of the model, but are not
important for the goal of the data analysis.
These are called \myindex{nuisance parameters}.
An unshaded node $y$ without children (no outward going arcs) that is neither an action
node nor a utility node can always be
removed from a Bayesian network.
This corresponds to leaving out the term $p(y|parents(y))$
in the product of Equation~\ref{eq-dags}. 
Given
\[
p(a,b,y) ~=~ p(a)\,p(b|a)\, p(y|a,b) 
\]
then $y$ can be \myindex{marginalized}  out trivially
to yield:
\[
p(a,b) ~=~ p(a)\,p(b|a)~.
\]
More generally this applies to chain graphs---a chain component 
whose nodes are all unshaded and have no children
can be  removed.
If $y$ is a node without children,
then remove the node with $y$ from the graph
and the arcs to it;  ignore the factor $p(y|parents(y))$ in the
full joint form.
Consider that the $i$-th case in Figure~\ref{fig-sbnl} 
(nodes $class_i, var_{1,i}, var_{2,i}, var_{3,i}$)
can be removed from the model
without affecting the rest of the graph.

\subsection{Removal of plates by exact methods}
\label{ssct-suff}

Consider the simple coins problem of Figure~\ref{fig-coin} again.
The graph represents the joint probability for
$p(\theta, heads_1, \ldots, heads_N)$.
The main question of interest here is the conditional
probability of $\theta$ given the data $heads_1, \ldots, heads_N$.
This could be obtained through repeated arc reversals
between $\theta$ and $heads_1$,
then between $\theta$ and $heads_2$, and so on, 
until all the data appears
before $\theta$ in the directed graph.   % ???? show
Doing this repeated series of applications of Bayes theorem
yields  a fully connected graph
with $(N+1)N/2$ arcs.
The corresponding formula for the posterior simplified with
Lemma~\ref{lem-remgiven} is also simple:
\begin{equation}
\label{eq-binop}
p(\theta |  heads_1, \ldots, heads_N, \alpha_1=1.5,\alpha_2=1.5) ~=~
        \frac{\theta^{\alpha_1-1+p} \, (1-\theta)^{\alpha_2-1+n}}
        {Beta(\alpha_1+p,\alpha_2+n)} 
\end{equation}
where $p$ is the number of heads in the sequence and $n=N-p$
is the number of tails.
This is a worthwhile introductory exercise in 
Bayesian decision theory \cite{howard70}
that should be familiar to most students of statistics.
Compare this with Equation~(\ref{eq-bino}).
There are several important points to notice about this result:
\begin{itemize}
\item
Effectively, this does a parameter update,
$\alpha_1' = \alpha_1+p$ and $\alpha_2' = \alpha_2+n$,
requiring no search or numerical optimization.
The whole sequence of tosses, irrespective
of its length and the ordering of the heads and tails,
can be summed up with two numbers.
These summary statistics are called \myindex{sufficient statistics}
because,  assuming the model
used is correct,
they are sufficient to explain all that is
important about $\theta$ in the data,.
\item
The corresponding graph can be simplified as shown in 
Figure~\ref{fig-coin2}.
The plate is efficiently removed and replaced
by the sufficient statistics (two numbers)
irrespective of the size of the sample.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{coin2.ps}{2in}
\end{center}
\vskip -0.2in
\caption{Removing the plate in the coin problem}
\label{fig-coin2}
\end{figure}
\item
The posterior distribution has a simple form.
Furthermore, all the moments of $\theta$,
$\log \theta$, and
$\log (1-\theta)$
for the distribution can be computed
as simple functions 
of the normalizing constant, $ Beta(\alpha_1',\alpha_2')$.
For instance:
\begin{eqnarray*}
\expec{\theta|heads_1,\ldots,heads_N,\alpha_1,\alpha_2}{\log \theta} &=&
        \frac{ \partial \log Beta(\alpha_1', \alpha_2') }{\partial \alpha_1}
                ~,\\
\expec{\theta|heads_1,\ldots,heads_N,\alpha_1,\alpha_2}{\theta} &=&
                \frac{Beta(\alpha_1'+1, \alpha_2') }
                {Beta(\alpha_1', \alpha_2')} ~,\\
\expec{\theta|heads_1,\ldots,heads_N,\alpha_1,\alpha_2}
        {\left(\theta-\overline{\theta}\right)^2} &=&
        \frac{Beta(\alpha_1'+2, \alpha_2') }
                {Beta(\alpha_1', \alpha_2')} -
        \frac{Beta^2(\alpha_1'+1, \alpha_2') }
                {Beta^2(\alpha_1', \alpha_2')}~.
\end{eqnarray*}
This result might seem somewhat obscure,
but it is a general property holding for a large
class of distributions that allows some
averages to be calculated by symbolic manipulation
of the normalizing constant.
\end{itemize}

\subsection{The exponential family}
\label{ssct-expon}

This result generalizes to
a much larger class of distributions referred to as the exponential
family \cite{Casella.berger.90,Degroot.70}.  This includes 
standard undergraduate distributions such as Gaussians,
Chi squared, and Gamma,
and many more complex distributions constructed from
simple components including
class probability trees over discrete input domains \cite{Buntine.AISIII}, 
simple discrete and Gaussian versions of a Bayesian network
\cite{Whittaker.nets},
and linear regression with a Gaussian error.
Thus, these are a broad and not insignificant class of distributions
that are given in the definition below.
Their general form has a linear combination of parameters
and data in the exponential.
\begin{defn}
\label{defn-expfam}
A space $X$ is \myindex{independent} of the
parameter $\theta$ if the space remains the same when just
$\theta$ is changed.
If the domains of $x,y$
 are independent of $\theta$,
then the conditional distribution for $x$ given $y$,
$p(x|y,\theta,M)$,
is in the \myindex{exponential family}
when
\begin{equation}
\label{eq-exp}
p(x|y,\theta,M) ~=~ \frac{h(x,y)}{Z(\theta)} \exp \left( \sum_{i=1}^k
                w_i(\theta) t_i(x,y) \right) 
\end{equation}
for some functions $w_i$, $t_i$, $h$ and $Z$ and some integer $k$,
for $h(x,y)>0$.
The \myindex{normalization constant}
$Z(\theta)$ is known as the \myindex{partition function}.
\end{defn}
Notice the functional form of Equation~(\ref{eq-exp}) is
similar to the functional form for an undirected
graph of Equation~(\ref{eq-ugs}), as holds
in many cases for a Markov random field.
For the previous coin tossing example, both the coin tossing distribution
(a binomial on $heads_i$) and the posterior distribution on the model
parameters ($\theta$) are in the exponential family.
To see this, notice the following rewrites of the original
probabilities.
These make the components $w_i$, $t_i$ and $Z$ explicit.
\begin{eqnarray*}
&& p(heads|\theta) ~=~ \exp \left(  1_{heads=true}\,\log \theta +
                1_{heads=false}\,\log (1-\theta) \right) ~, \\
\lefteqn{p(\theta | heads_1, \ldots, heads_N, \alpha_1, \alpha_2) ~=~}\\
     &&   \frac{1}{Beta(\alpha_1+p,\alpha_2+n)} \exp \left(
                (\alpha_1+p-1) \, \log \theta +
                (\alpha_2+n-1) \, \log (1-\theta) \right) ~.
\end{eqnarray*}
Table~\ref{tab-ef} in  Appendix~\ref{app-ef} gives a selection of
distributions, and their functional form.
Further details can be found in most textbooks on probability
distributions \cite{Degroot.70,Bernardo.Smith.94}.

The following is a simple graphical reinterpretation
 of the Pitman-Koopman Theorem
from statistics \cite{Jeffreys.61,Degroot.70}.
In Figure~\ref{fig-suff}(a),
$T(x_*,y_*)$ is a statistic of fixed dimension
independent of the sample size $N$
(corresponding to $n,p$ in the coin tossing example).
The theorem says  that the sample 
in Figure~\ref{fig-suff}(a) can be
summarized in statistics,  as shown in 
Figure~\ref{fig-suff}(b),
if and only if
the probability distribution for $x|y,\theta$ is in the exponential family.
In this case,  $T(x_*,y_*)$
is a sufficient statistic.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{suffstatc.ps}{4in}
\end{center}
\vskip -0.2in
\caption{The generalized graph for plate removal}
\label{fig-suff}
\end{figure}
\begin{theorem}
\label{thm-plate}
{\em (}\myindex{Recursive arc-reversal}{\em )}.
Consider the model $M$ represented by the
graphical model for a sample of size $N$
given in Figure~\ref{fig-suff}(a).
Have $x$ in the domain $X$ and
$y$ in the domain $Y$, both domains
are independent  of $\theta$,
and both domains have components
that are real valued or finite discrete.
Let the conditional distribution for
$x$ given $y,\theta$ be $f(x|y,\theta)$,
which is positive for all $x \in X$.
If first derivatives exist with respect to
all real valued components
of $x$ and $y$,
the plate removal operation applies for all samples
$x_* = x_1,\ldots,x_N$, $y_* = y_1,\ldots,y_N$, and $\theta$,
as given in Figure~\ref{fig-suff}(b),
for some sufficient statistics $T(x_*,y_*)$ of dimension independent of $N$
if and only if the conditional distribution for
$x$ given $y,\theta$ is in the exponential family,
given by Equation~(\ref{eq-exp}).
In this case, $T(x_*,y_*)$ is an
invertible function of the $k$ averages:
\[
        \frac{1}{N} \sum_{j=1}^N t_i(x_j) ~~:~~i=1,\ldots,k  ~.
\]
\end{theorem}
In some cases, this extends to domains $X$ and $Y$ dependent 
on $\theta$ \cite{Jeffreys.61}.

Sufficient statistics for a distribution
from the exponential family are easily read from the functional form
by taking a logarithm.
For instance, for the multivariate Gaussian,
the sufficient statistics are
$x_i$ for $i=1,\ldots,d$ and $x_ix_j$ for $0\leq i\leq j \leq d$,
and the normalizing constant $Z(\mu,\Sigma)$ is given by:
\[
Z(\mu,\Sigma) ~=~ \frac{(2\pi)^{d/2}}{\det^{1/2}\Sigma}
                \exp\left(\frac{1}{2}
        \mu^\dagger \Sigma \mu \right) ~.
\]

As for coin tossing, it generally holds that
if a sampling distribution (a binomial on $heads_i$)
is in the exponential family,
then the posterior distribution for the model
parameters ($\theta$) can also be cast as exponential family.
This is only useful when the normalizing
constant 
(this is $Beta(\alpha_1+p,\alpha_2+n)$ in the coin tossing example) 
and its derivatives are  readily computed.
\begin{lemma}
\label{lem-conj}
{\em (}{\em The conjugacy property}{\em )}.
In the context in Theorem~\ref{thm-plate}, assume the distribution for
$x$ given $y,\theta$ can be represented by the exponential family.
Factor the normalizing constant $Z(\theta)$ into two components,
$Z(\theta) = Z_1(\theta)Z_2$,
where the second is the constant part independent of $\theta$.
Assume the prior on $\theta$ takes the form:
\begin{equation}
\label{eq-conj}
p(\theta|\tau,M) ~=~ \frac{f(\theta)}{Z_\theta(\tau)}
         \exp \left( \tau_{k+1}(\log 1/Z_1(\theta)) + 
                \sum_{i=1}^k \tau_i w_i(\theta) \right) 
\end{equation}
for some $k+1$ dimensional parameter $\tau$,
where $Z_\theta(\tau)$ is the appropriate normalizing constant
and $f(\theta)$ is any function.
Then the posterior distribution for $\theta$,
$p(\theta|\tau,x_1,\ldots,x_N,M)$,
is also represented by Equation~(\ref{eq-conj})
with the parameters
\begin{eqnarray*}
\tau_{k+1} ' &=& \tau_{k+1} + N  \\
\tau_i' &=& \tau_{i} + \sum_{j=1}^N t_i(x_j,y_j)  ~~~~~~~~~~\mbox{$i=1,\ldots,k$}~.
\end{eqnarray*}
\end{lemma}
When the function $f(\theta)$ is trivial, for instance uniformly
equal to 1,
then the distribution in Equation~(\ref{eq-conj})
is referred to as the \myindex{conjugate} distribution,
which means it has a mathematical form mirroring
that of the sample likelihood.
The prior parameters $\tau$, by looking at the
update equations in the lemma, can be thought of as corresponding
to the sufficient statistics from some ``prior sample'' and
given by $\tau_{k+1}$.

This property is useful for analytic and computational purposes.
Once the posterior distribution is found,
and assuming it is one of the standard distributions,
the property can easily be established.
Table~\ref{tab-conj} in Appendix~\ref{app-ef}
gives some standard conjugate prior distributions
for those in Table~\ref{tab-ef},
and Table~\ref{tab-conjpost} gives their matching posteriors.
More extensive summaries of this are given by DeGroot
\citeyear{Degroot.70} and Bernardo and Smith 
\citeyear{Bernardo.Smith.94}.
The parameters for these priors can be set using standard
reference priors \cite{Box.tiao,Bernardo.Smith.94}
or elicited from a domain expert.

There are several other important consequences of the Pitman-Koopman
Theorem or recursive arc reversal that should not go unnoticed.
\begin{comment}
\label{comm-disc}
If $x,y$ are discrete and finite valued,
then the distribution $p(x|y,\theta)$ can be represented
as a member of the exponential family.
This holds because a positive finite discrete distribution 
can always be represented
as an extended case statement
in the form
\[
p(x|y,\theta) ~=~  \exp\left( \sum_{i=1,k} 1_{t_i(x,y)} 
f_i(\theta) \right)  
\]
where the boolean functions $t_i(x,y)$ are a set of mutually
exclusive and exhaustive conditions.
The indicator function $1_A$ has the
value 1 if the boolean $A$ is true and 0 otherwise.
The main importance of the exponential
family is in continuous or integer domains.
Of course, since a large class of
functions $\log p(x|y,\theta)$ can always be
approximated arbitrarily well by a polynomial in $x,y$ and $\theta$
with sufficiently many terms, the exponential family covers a broad
class of distributions.
\end{comment}

The application of the exponential family
to learning is perhaps the earliest published result on
 computational learning theory.
The following two interpretations of the recursive arc reversal theorem are
relevant mainly for distributions involving continuous variables.
\begin{comment}
An incremental learning algorithm with finite memory
must compress the information it has seen so far in the training
sample into a smaller set of statistics.
This can only be done
without sacrificing information in the sample, in a context
where all probabilities are positive,
if the hypothesis or search space of 
learning is a distribution from the exponential family.
\end{comment}
\begin{comment}
The computational requirements for
learning an exponential family distribution
are guaranteed to be linear in the sample size:
first compute the sufficient statistics and then learning
proceeds independently of the sample size.
This could be exponential in the dimension of the feature space, however.
\end{comment}

Furthermore, in the case where the functions $w_i$ are {\em full rank} in
$\theta$
(dimension of $\theta$ is $k$, same as $w$, and the
Jacobian of $w$ with respect to $\theta$ is invertible, 
$\det\left(\frac{\diff w(\theta)}{\diff \theta}\right)\neq 0$),
various moments of the distribution can easily be found.
For this situation, the function $w^{-1}$, when it exists, is called the
\myindex{link function} \cite{Mccullagh.nelder.glm}.
\begin{lemma}
\label{lem-efsolve}
Consider the notation of Definition~\ref{defn-expfam}.
If the link function $w^{-1}$ for an exponential family distribution exists,
then moments of functions of $t_i(x,y)$ and $\exp \left( t_i(x,y) \right)$
can be expressed in terms of derivatives and direct applications
 of the functions $Z$, $t_i$, $w_i$, and $w^{-1}$.
If the normalizing constant $Z$ and the link function are in closed form,
then so will the moments.
\end{lemma}
Techniques for doing these symbolic calculations are given
in Appendix~\ref{app-ef}.
Exponential family distributions then fall into two
groups.  There are those where the normalizing constant
and link function are known,
such as the Gaussian.  One can efficiently compute
their moments and determine the functional form
of their conjugate distributions up to the normalizing constant.
For others this is not the case.
For others, such as a Markov random field used in image processing,
moments can generally only be computed by an approximation
process like  Gibbs sampling given in Section~\ref{ssct-gibbs}.

\subsection{Linear regression: an example}
\label{ssct-linreg}

As an example, consider the problem of linear regression
with Gaussian error described in Figure~\ref{fig-ln}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{linreg.ps}{3in}
\end{center}
\vskip -0.2in
\caption{Linear regression with Gaussian error}
\label{fig-ln}
\end{figure}
This is an instance of a generalized linear model and
has a linear construction
\[
m ~=~\sum_{j=1}^M \theta_j basis_j(x_1,\ldots,x_n)
\]
 at its core.
The $M$ basis functions are known deterministic functions of the
input variables $x_1, \ldots, x_n$.
These would typically be nonlinear orthogonal functions
such as Legendre polynomials.
These combine linearly with the parameters $\theta$ to produce
the mean $m$ for the Gaussian.

The corresponding learning problem represented
as plates is expressed in Figure~\ref{fig-linregl}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{linregl.ps}{3.5in}
\end{center}
\vskip -0.2in
\caption{The linear regression problem}
\label{fig-linregl}
\end{figure}
The joint probability for this model is as follows:
\[
p(\sigma)\,p(\theta|\sigma)
        \frac{1}{\left(\sqrt{2\pi}\sigma\right)^N } 
                   e^{ - \frac{1}{2\sigma^2}
                     \sum_{i=1}^N
              \left(y_i-\sum_{j=1}^M \theta_j basis_j(x_{.,i})\right)^2} 
\]
where $x_{.,i}$ denotes the vector of values for the
$i$-th datum, 
$x_{1,i},\ldots,x_{n,i}$.
Linear regression with Gaussian error falls in the exponential
family because a Gaussian is in the exponential family
        and the mean of the simple Gaussian is a linear function
        of the regression parameters (see Lemma~\ref{lem-exp}).

For this case, the correspondence to the exponential family
is drawn as follows.
The individual data likelihoods, 
$p(y| x_1,\ldots, x_n,\theta,\sigma)$, need only be considered.
Expand the probability to show it is a linear
sum of data terms and parameter terms.
\begin{eqnarray*}
\lefteqn{p(y| x_1,\ldots, x_n,\theta,\sigma) 
         } \\
&=&
        \frac{1}{\sqrt{2\pi}\sigma}
        \exp \left( -\frac{1}{2\sigma^2} \left( y - 
                \sum_{j=1}^M basis_j(x_{.})\, \theta_j \right)^2 \right)
                 ~, \\
        &=& 
        \frac{1}{\sqrt{2\pi}\sigma} 
        \exp \left( -\frac{1}{2\sigma^2} y^2 
                - \sum_{j,k=1}^M basis_j(x_{.}) \, basis_k(x_{.}) \,
                        \frac{\theta_j\theta_k}{2\sigma^2}
                + \sum_{j=1}^M basis_j(x_{.})\,y\, 
                        \frac{\theta_j}{2\sigma^2} \right) ~.
\end{eqnarray*}
The data likelihood in this last line
can be seen to be in the same form as the general exponential
family where the sufficient statistics are the various data terms
in the exponential.
Also, the link function does not exist because there
are $M$ parameters and $M(M+1)/2$ sufficient statistics.

The model of Figure~\ref{fig-linregl}
can therefore be simplified to the graph in Figure~\ref{fig-linregls},
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{linregls.ps}{2in}
\end{center}
\vskip -0.2in
\caption{The linear regression problem with the plate removed}
\label{fig-linregls}
\end{figure}
where $q$ and $S$ are the usual
sample means and covariances
obtained from the so-called normal equations
of linear regression.
$S$ is a matrix of dimension $M$ (the number of basis functions)
and $q$ is a vector of dimension $M$.
\begin{eqnarray*}
S_{j,k} &=&  \frac{1}{N} \sum_{i=1}^N basis_j(x_{.,i}) \,basis_k(x_{.,i}) 
                \\
q_j &=& \frac{1}{N} \sum_{i=1}^N basis_j(x_{.,i}) \, y_i \\
ysq &=& \frac{1}{N} \sum_{i=1}^N y_i^2 ~.
\end{eqnarray*}
These three sufficient statistics
can be read directly from the data likelihood above.

Consider the formula:
\[
\int_y p(y| x_1,\ldots, x_n,\theta,\sigma)  \diff y ~.
\]
Differentiating with respect to $\theta_i$
shows that the expected value of $y$ given $x_1,\ldots, x_n$ and
$\theta,\sigma$ is the mean (as expected):
\[
\expec{y|x_1,\ldots, x_n,\theta,\sigma}{y} ~=~ m ~=~
        \sum_{j=1}^M  basis_j(x_{.})\, \theta_j ~.
\]
Differentiating with respect to $\sigma$ shows that
 the expected error from the mean is $\sigma^2$
(again expected):
\[
\expec{y|x_1,\ldots, x_n,\theta,\sigma}{(y-m)^2} ~=~
        \sigma^2 ~.
\]
Higher-order derivatives give formula for higher-order moments
such as skewness and kurtosis \cite{Casella.berger.90},
which are functions of the second, third and fourth central moments.
While these are well known for the Gaussian,
the interesting point is that
these formula are constructed 
by differentiating the component functions
in Equation~(\ref{eq-exp})
without recourse to integration.
Finally, the conjugate distribution for the
parameters $\theta$ in this linear regression
problem is the multivariate Gaussian distribution
when $\sigma$ is known,
and for $\sigma^2$ is the inverted Gamma.


\section{Recognizing and using the exponential family}

How can recursive arc reversal be applied automatically
to a graphical model?
First, when a graphical model
or some subset of a graphical model falls in the exponential family
needs to be identified.
If each conditional distribution in a Bayesian network or 
chain component in a chain graph is
exponential family, then the full joint is exponential family.
The following lemma gives this with some additional
conditions for deterministic nodes.
This applies to Bayesian networks using Comment~\ref{comm-cgtodag}.
\begin{lemma}
\label{lem-exp}
A chain graph has a single plate.
Let the non-deterministic variables inside the plate
be $X$, and the deterministic variables be $Y$.
Let the variables outside the plate be $\theta$.
If:
\begin{enumerate}
\item
All arcs crossing a plate boundary
are directed into the plate.
\item
For all  chain components $\tau$, the  conditional distribution
$p(\tau | parents(\tau ))$ is from the
exponential family with data variables from $(X,Y)$ and
model parameters from $\theta$;
furthermore
$\log p(\tau  |parents(\tau ),\theta)$ is a polynomial function of 
variables in $Y$.
\item
Each variable  $y\in Y$ 
can be expressed as a deterministic function of the form
\[
y ~=~ \sum_{i=1}^l u_i(X)\,v_i(\theta) 
\]
for some functions $u_i,v_i$.
\end{enumerate}
Then the conditional distribution 
$p(X,Y|\theta)$ is from the exponential family.
\end{lemma}

Second, how can these results be used when
the model does not fall in the exponential family?
There are two categories of techniques available in this
context.
In both cases, the algorithms concerned can be
constructed from the graphical specifications.
The two new classes together with the 
recursive arc reversal case
are given in Figure~\ref{fig-generic}.
In each case, (I) denotes the graphical configuration
and (II) denotes the operations and simplifications
 performed by the algorithm.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{generic.ps}{6in}
\end{center}
\vskip -0.2in
\caption{Three categories of algorithms using the exponential family}
\label{fig-generic}
\end{figure}
When the various normalization
constants are known in closed form and appropriate moments
and Bayes factors can be computed quickly, all three algorithm schemas
have reasonable computational properties.  

The first category is where a useful subset of the model
does fall into the exponential family.
This is represented by the \myindex{partial exponential family}
in Figure~\ref{fig-generic}.
The part of the problem that is exponential family is simplified
using the recursive arc reversal of Theorem~\ref{thm-plate},
and the remaining part of the problem is typically
handled approximately.
Decision trees and Bayesian networks over multinomial or Gaussian variables
also fall into this category.
This happens because when the structure of the 
tree or Bayesian network is given the
remaining problem is composed of a product of multinomials or Gaussians.
This is the basis of various Bayesian algorithms developed
for these problems 
\cite{Buntine.AISIII,Madigan.raftery.jasa,Buntine.UAI91,Spiegel.dawid.ss,HGC.TR9409}.
Strictly speaking, decision trees and Bayesian networks over multinomial 
or Gaussian variables are in the exponential family
(see Comment~\ref{comm-disc}).
However, it is more computationally convenient to treat them
this way.  This category is discussed more in Section~\ref{ssct-partial}.

The second category is where, if  some hidden
variables are introduced  into the data, the problem becomes exponential
family if  the hidden values are known.
This is represented by the \myindex{mixture model} in 
Figure~\ref{fig-generic}.
Mixture models \cite{Titter.smith.85,Poland.thesis}
are used to model unsupervised learning,
incomplete data in the classification problems, robust regression,
and general density estimation.
Mixture models extend the exponential family to a rich class
of distributions, so this second category is an important
one in practice.
General methods for handling these problems correspond
to Gibbs sampling (and other Markov chain Monte Carlo methods)
discussed in Section~\ref{ssct-gibbsplates}
and its deterministic counterpart the expectation maximization 
algorithm, discussed in Section~\ref{ssct-EM}.
As shown in Figure~\ref{fig-generic},
these  algorithms cycle back and forth between
a process that re-estimates $c$ given $\theta$
using first-order inference
and a process that uses the fast exponential family algorithms
to re-estimate $\theta$ given $c$.


\section{Other operations on graphical models}
\label{sct-oper}

The recursive arc reversal theorem of Section~\ref{ssct-suff}
characterizes when plates can be readily removed and the
sample summarized in some statistics.
Outside of these cases, more
general classes of approximate algorithms exist.
Several of these are introduced in Section~\ref{sct-approx}
and more detail is given, for instance, by Tanner
\citeyear{Tanner.93}.
These more general algorithms require a number of
basic operations be performed on graphs:
\begin{description}
\item[Decomposition: ]
A learning problem can sometimes be decomposed into
simpler sub-problems with each yielding to separate analysis.
One form of decomposition of learning problems is considered in
Section~\ref{ssct-decomp}.
Another related form that applies to undirected graphs
is developed by Dawid and Lauritzen 
\citeyear{Dawid.lauritzen.anstat}.
Other forms of decomposition can be done at the modeling
level, where the initial model is constructed in
a manner requiring fewer parameters,
as is Heckerman's similarity networks
\citeyear{Heckerman.91}.
\item[Exact Bayes factors: ]
Model selection and averaging methods are used to
deal with multiple models 
\cite{Kass.Raftery.571,Buntine.AISIII,Stewart.stat,Madigan.raftery.jasa}.
These require the computation of Bayes factors
for models constructed during search.
Exact methods for computing Bayes factors are considered
in Section~\ref{ssct-select}.
\item[Derivatives: ]
Various approximation and search algorithms require derivatives
be calculated, as discussed next.
\end{description}

\subsection{Derivatives}
\label{ssct-deriv}

An important operation on graphs is the
calculation of derivatives of parameters.
This is useful after conditioning on the known data
to do approximate inference.
Numerical optimization using derivatives can be
done to search for MAP values of parameters,
or to apply  the Laplace approximation 
to estimate moments.
This section shows how to compute derivatives using
operations local to each node.
The computation is therefore easily parallelized,
as is popular, for instance, in neural networks.

Suppose a graph is used to compile a function
that searches for the MAP values of
parameters in the graph conditioned on the known data.
In general, this requires use of numerical optimization
methods \cite{Gill.murray.etal.81}.
To use a gradient descent, conjugate gradient
or Levenberg-Marquardt
approach requires calculation of first derivatives.
To use a Newton-Raphson
approach requires calculation of second derivatives, as well.
While this could be done numerically by difference
approximations,  %  ????  explain these
more accurate calculations exist.
Methods for symbolically differentiating networks
of functions, and piecing together the results
to produce global derivatives 
are well understood \cite{ADoA91}.
For instance, 
software is available for taking a function
defined in Fortran, C++ code, or some other language, 
to produce a second function
that computes the exact derivative.
These problems are also well understood for feed-forward networks
\cite{WerbosMcAvoySu:92,Buntine.weigend.2nd},
and graphical models with plates only add some additional complexity.
The basic results are reproduced in this section
and some simple examples given to highlight
special characteristics arising from their use with chain graphs.

Consider the problem of learning a feed-forward network.
A simple feed-forward network is given in Figure~\ref{fig-ffnnl}(a).
The corresponding learning problem is given in Figure~\ref{fig-ffnnl}(b),
representing the feed-forward network as a Bayesian network.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{feedforwlearnplus.ps}{5.5in}
\end{center}
\vskip -0.2in
\caption{Learning a feed-forward network}
\label{fig-ffnnl}
\end{figure}
Here the sigmoid units of the network are modeled with deterministic
nodes, and the network output represents the mean of a
bivariate Gaussian with inverse variance matrix $\Sigma$.
Because of the nonlinear sigmoid function making 
the deterministic mapping from inputs $x_1,x_2,x_3$
to the means $m_1,m_2$, this learning problem has no
reasonable component falling in the exponential family.
A rough fallback method is to calculate a
MAP value for the weight parameters.
This would be the method used for the
Laplace approximation \cite{Buntine.weigend.bbp,mackay:bp}
covered in \cite{Tanner.93,Tierney.kadane.jasa}.
The setting of priors for feed-forward networks is 
difficult \cite{Mackay.cl93,Nowlan.NIPS4,Wolpert.NIPS6},
and it will not be considered here other than assuming
a prior is used, $p(w)$.
The graph implies the following posterior  probability:
\begin{eqnarray}
\label{eq-ffp}
\lefteqn{
p(\Sigma,w_1,\ldots,w_5~|~o_{1,i},o_{2,i},x_{1,i},x_{2,i},x_{3,i}
        ~:~i=1,\ldots,N)  } && \\
\nonumber
&\propto& p(\Sigma)\,p(w_1,\ldots,w_5)\,\prod_{i=1}^N \frac{\det^{1/2}\Sigma}{2\pi}
                \exp\left( \frac{1}{2}  (o_i-m_i)^\dagger 
                        \Sigma (o_i-m_i) \right) \\
\nonumber
m_i &=&  Sigmoid(w_i^\dagger h) \\
\nonumber
h_i &=& Sigmoid(w_{i+2}^\dagger x) ~.
\end{eqnarray}
The undirected clique on the parameters $w$ indicates
the prior has a term $p(w)$.
Suppose the posterior is differentiated with respect to the 
parameters $w_4$.
The result is well known to the neural network community
since this kind of calculation
yields the standard back-propagation equations.

Rather than work through this calculation, instead look at the 
general case.
To develop the general formula for
differentiating a graphical model,
a few more concepts are needed.
Deterministic nodes form islands of determinism
within the uncertainty represented by the graph.
Partial derivatives within each island can be calculated
via recursive use of the chain rule, for instance,
by forward or backward propagation of derivatives
through the equations.
For instance,
\myindex{forward propagation} for the above network gives:
\[
\frac{\partial m_1}{\partial w_4} ~=~ 
     \sum_{i=1}^3 \frac{\partial m_1}{\partial h_i}
           \frac{\partial h_i}{\partial w_4}  ~.
\]
This is called forward propagation because the derivatives
with respect to $w_4$ are propagated forward in the network.
In contrast, \myindex{backward propagation} would propagate derivatives of
$m_1$ with respect to different variables backwards.
For each island of determinism, the important variables are the
output variables, and  their derivatives are required.
So for the feed-forward network above,
partial derivatives of $m_1$,  $m_2$, and
$m_3$ with respect to $w_4$ are required.
\begin{defn}
\label{def-nd}
The \myindex{non-deterministic children} of a node $x$,
denoted $ndchildren(x)$,
are the set of non-deterministic 
variables $y$ such that there exists a directed
path from $x$ to $y$ given by $x, y_1, \ldots, y_n, y$,
with all intermediate variables ($y_1, \ldots, y_n$) being
deterministic.
The \myindex{non-deterministic parents} of a node $x$,
denoted $ndparents(x)$,
are the set of non-deterministic
variables $y$ such that there exists a directed
path from $y$ to $x$ given by $y, y_1, \ldots, y_n, x$,
with all intermediate variables ($y_1, \ldots, y_n$) being
deterministic.
The \myindex{deterministic children} of a node $x$,
denoted  $detchildren(x)$,
are the set of deterministic 
variables $y$ that are children of $x$.
The \myindex{deterministic parents} of a node $x$,
denoted  $detparents(x)$,
are the set of deterministic 
variables $y$ that are parents of $x$.
\end{defn}
For instance, in the model in Figure~\ref{fig-ffnnl},
the non-deterministic children of $w_3$ are $o_1$ and $o_2$.
Deterministic nodes can be removed from a graph by rewriting
the equations represented into the remaining variables of the graph.
Because some graphical operations do not apply to deterministic nodes, 
this removal is often done
implicitly within a theorem.
This goes as follows:
\begin{lemma}
\label{lem-nodet}
A chain graph $G$ with nodes $X$ has deterministic nodes $Y\subset X$.
The chain graph $G'$ is created by adding to $G$
 a directed arc from every node to its non-deterministic children,
and by deleting the deterministic nodes $Y$.
The graphs $G$ and $G'$ are equivalent probability models on the
nodes $X-Y$.
\end{lemma} %  OK now apply this to the figure

The general formula for differentiating
Bayesian networks with plates and deterministic nodes
is given below in Lemma~\ref{lem-differ}.
This is nothing more than the chain rule for differentiation,
but it is important to notice the network structure of the computation.
When partial derivatives are computed over networks,
there are local and global partial derivatives that can be
different.
Consider the feed-forward network of Figure~\ref{fig-ffnnl} again.
On this figure,
place an extra arc from $w_4$ to $m_2$.
Now consider the partial derivative 
of $m_2$ with respect to $w_4$.
The value of $m_2$ is influenced by $w_4$
directly, as the new arc shows, and indirectly via $h_2$.
When computing a partial derivative involving indirect
influences, we need to differentiate between the
direct and indirect effects.
Various notations are used for this
\cite{WerbosMcAvoySu:92,Buntine.weigend.2nd}.
Here the notation of a local versus global derivative is used.
The \myindex{local partial derivative} is subscripted with an
$l$, $\partial/\partial_l$
and represents the partial derivative computed at the
node using only the direct influences---the parents.
For the example of the partial derivative 
of $m_2$ with respect to $w_4$, the various local
partial derivatives combine to produce the global
partial derivative:
\[
\frac{\partial m_2}{\partial w_4} ~=~
        \frac{\partial m_2}{\partial_l w_4} + 
    \frac{\partial m_2}{\partial_l h_2}\frac{\partial h_2}{\partial_l w_4} ~.
\]
This is equivalent to:
\[
\frac{\partial m_2}{\partial w_4} ~=~
        \frac{\partial m_2}{\partial_l w_4} + 
    \frac{\partial m_2}{\partial_l h_2}\frac{\partial h_2}{\partial_l w_4}
 +  \frac{\partial m_2}{\partial_l h_1}\frac{\partial h_1}{\partial_l w_4}~,
\]
since $\frac{\partial h_1}{\partial_l w_4}=0$.

In general,
the (global) partial derivative for an index variable $\theta_i$
is the sum of the
\begin{itemize}
\item
local partial derivative 
at the node containing $\theta_i$,
\item
the
partial derivatives for each child of $\theta_i$
that is also a non-deterministic 
child, and
\item
combinations of (global) partial derivatives for deterministic
children found by backward or forward
propagation of derivatives.
\end{itemize}
\begin{lemma}
\label{lem-differ}
{\em (Differentiation)}.
A model $M$ is represented by a Bayesian network $G$ with plates and
deterministic nodes on variables $X$.
Denote the known variables in $X$ by $K$ and the unknown variables 
by $U = X-K$.
Let the conditional probability represented by the graph $G$
be $p(U|K,M)$.
Let $\theta$ be some unknown variable in the graph,
and let $1_{nd(\theta)}$ be 1 if $\theta$ is non-deterministic
        and 0 otherwise.
If $\theta$ occurs inside a plate then let $i$ be some
arbitrary valid index ($i \in indval(\theta)$), otherwise let $i$ be null.
Then:
\begin{eqnarray}
\label{eq-differ} 
\frac{\partial \log p(U|K,M)}{\partial \theta_i}  &=&
         1_{nd(\theta_i)}
           \frac{\partial \log p(\theta_i|parents(\theta_i))}
          {\partial_l \theta_i}   \\
\nonumber
&&  ~~~~+ \sum_{x \in ndchildren(\theta_i)\cap children(\theta_i)  } 
             \frac{\partial \log p(x|parents(x))}{\partial_l \theta_i} \\
\nonumber
&&  ~~~~+ \sum_{x \in ndchildren(\theta_i) } 
              \sum_{y\in detparents(x), y\neq \theta_i} 
                        \frac{\partial \log p(x|parents(x))}{\partial_l y}
                        \frac{\partial y}{\partial \theta_i}  ~.
\end{eqnarray}
Furthermore, if $Y\subset U$ is some subset of the unknown variables,
then the partial derivative of the probability
of $Y$ given the known variables, $p(Y|K,M)$ is an expected
value of the above probabilities:
\begin{equation}
\label{eq-margdiff}
\frac{\partial \log p(Y|K,M)}{\partial \theta_i} ~=~
        \expec{U-Y|Y,K,M}{\frac{\partial \log p(U|K,M)}{\partial \theta_i}}~.
\end{equation}
\end{lemma}
Equation~(\ref{eq-differ})
contains only one global partial derivative which is
 inside the double sum on the right  side. 
This is the partial derivative
${\partial y}/{\partial \theta_i }$
and can be computed
from its local island of determinism using the 
chain rule of differentiation,
for instance, using  forward propagation from $\theta_i$'s
deterministic children, or backward propagation
from $\theta_i$'s non-deterministic parents.
%   ??? really need to give example for Figure~\ref{fig-ffnnl}

To apply the Differentiation Lemma
on problems like feed-forward networks or unsupervised learning,
the lemma needs to be extended to chain graphs.  This means
differentiating Markov networks as well as Bayesian networks,
and handling the expected value in Equation~(\ref{eq-margdiff}).
These extensions are explained below after first giving
two examples.

As a first example, consider the feed-forward network problem of 
Figure~\ref{fig-ffnnl}.
By treating the two output units in the
feed-forward network as a single variable,
a Cartesian product
$(o_1,o_2)$, the above Differentiation Lemma can now
be applied directly to the feed-forward network model of 
Figure~\ref{fig-ffnnl}.   
This uses the simplification given in
Section~\ref{ssct-cg} with Comment~\ref{comm-cg}. 
Let $P$ be the joint probability for the
feed-forward network model, given in Equation~(\ref{eq-ffp}).
The non-deterministic children of $w_4$ are the single chain
component consisting of the 
two variables $o_1$ and $o_2$. 
Its parents are the set $\{m_1,m_2\}$.
Consider the Differentiation Lemma.
There are no children of $w_4$ that are also
non-deterministic, so the middle sum
in Equation~(\ref{eq-differ}) of the Differentiation Lemma is empty.
Then the lemma yields, after expanding out the inner most sum:
\begin{eqnarray*}
\frac{\partial \log P}{\partial w} &=& 
\frac{\partial \log p(w)}{\partial_l w_4}  + \\ &&
\sum_{i=1}^N   
\left(
   \frac{\partial \log p(o_{1,i},o_{2,i} | \Sigma, m_1,m_2)  }
       {\partial_l m_1}
        \frac{\partial m_1}{\partial_l w_4} +
\frac{\partial \log p(o_{1,i},o_{2,i} | \Sigma, m_1,m_2)  }{\partial_l m_2}
        \frac{\partial m_2}{\partial_l w_4}  \right)~,
\end{eqnarray*}
where $p(o_{1,i},o_{2,i} | \Sigma, m_1,m_2) $ is the 
two-dimensional Gaussian,
and $\frac{\partial m_i}{\partial w_4}$ is from the
global derivative but evaluates to a local derivative.

As a second example, reconsider the simple unsupervised learning
problem given in the introduction to Section~\ref{sct-bl}.
The likelihood for a single datum given the 
model parameters is a marginal of the form:
\[
p(var_1=1,var_2=0,var_3=1|\phi,\theta) ~=~ \sum_{c=1}^{10}
 \phi_c \, \theta_{1,c}\, (1-\theta_{2,c}) \, \theta_{3,c}  ~.
\]
Taking the logarithm of the full case probability
$p(class,var_1,var_2,var_3|\phi,\theta)$
reveals the vectors of components $w$ and $t$ of the exponential distribution:
\begin{eqnarray*}
\lefteqn{\log p(class,var_1,var_2,var_3|\phi,\theta) }\\
&=&   \sum_{c=1}^{10} 1_{class=c}\,\log \phi_c +
          \sum_{j=1}^3 \sum_{c=1}^{10} 
                                \left(1_{class=c,var_j=true}\,\log\theta_{j,c} + 
                                   1_{class=c,var_j=false}\,\log (1-\theta_{j,c})\right)  ~. 
\end{eqnarray*}
Notice that the normalizing constant  $Z(\phi,\theta)$ is 1 in this case.
Consider finding the partial derivative
$\partial \log p(var_1,var_2,var_3|\phi,\theta) /\partial \theta_{2,5}$.
This is done for each case when differentiating the posterior
or the likelihood of the unsupervised learning model.
Applying Equation~(\ref{eq-margdiff}) to this yields:
\begin{eqnarray*}
\lefteqn{\frac{\partial \log p(var_1,var_2,var_3|\phi,\theta)}
        {\partial \theta_{2,5}}}  \\
&=&  \sum_{d=1}^{10} 
        \frac{\partial \log \theta_{2,d}}{\partial \theta_{2,5}}
      \expec{class=c|var_1,var_2,var_3,\phi,\theta}{    
            1_{class=d,var_2=true}  }  \\
&&~~~~~~ +   \frac{\partial \log (1-\theta_{2,d})}{\partial \theta_{2,5}}
          \expec{class=c|var_1,var_2,var_3,\phi,\theta}{        
       1_{class=d,var_2=true}  }   \\
&=&  \frac{1}{\theta_{2,5}} 1_{var_2=true}\,
          p(class=5| var_1,var_2,var_3,\phi,\theta) \\
&&  ~~~~~~        + \frac{1}{1-\theta_{2,5}} 1_{var_2=false} \,
         p(class=5| var_1,var_2,var_3,\phi,\theta) ~.
\end{eqnarray*}
Notice that the derivative is computed  by  doing 
first-order inference to find 
$p(class=5| var_1,var_2,var_3,\phi,\theta)$,
as noted by Russell, Binder, and Koller
\citeyear{Russell.TR824}.
This property holds in general for exponential
family models with missing or unknown
variables.   Derivatives are calculated
by some first-order inference followed by a combination
with derivatives of the $w$ functions.
Consider the notation for the exponential family introduced
previously in Definition~\ref{defn-expfam},
where the functional form is:
\[
p(x|y,\theta,M) ~=~ \frac{h(x,y)}{Z(\theta)} \exp \left( \sum_{i=1}^k
                w_i(\theta) t_i(x,y) \right) ~.
\]
Consider the partial derivative of a marginal of this
probability,
$p(x-u|y,\theta,M)$, for $u \subset x$.
Using
Equation~(\ref{eq-margdiff}) in the Differentiation Lemma,
the partial derivative becomes:
\begin{equation}
\label{eq-margdiffef}
\frac{ \partial \log p(x-u|y,\theta,M)}{\partial \theta}
   ~=~   
       \sum_{i=1}^k \frac{ \partial w_i(\theta)}{\partial \theta} 
              \expec{u|x-u,y,\theta}{t_i(x,y) }
       - \frac{ \partial Z(\theta)}{\partial \theta }
      ~.
\end{equation}
If the partition function is not known in closed form
(the case with the Boltzmann machine) then
the final derivative $\partial Z(\theta)/\partial \theta $
is approximated
(the key formula for doing this is Equation~(\ref{eq-expeqs1})
in Appendix~\ref{app-ef}).

To extend the Differentiation Lemma
to chain graphs, use the trick illustrated with the feed-forward
network.
First, interpret the chain graph as a Bayesian network on
chain components, as done in Equation~(\ref{eq-chdag}),
then apply the Differentiation Lemma.
Finally, evaluate necessary local partial
derivatives with respect to $\theta_i$ of each individual chain component.
Since undirected graphs are not necessarily normalized,
this may present a problem.
In general, there is an undirected graph $G'$ on
variables $X\cup Y$.
Following Theorem~\ref{th-ugs}, the general form is:  %????
\[
p(X|Y) ~=~  \frac{\prod_{C \in  Cliques(G')} f_C (C)}
            {\sum_X \prod_{C \in  Cliques(G')} f_C (C)}  ~.
\]
The local partial derivative with respect to $x$ becomes:
\begin{equation}
\label{eq-differug}
\frac{ \partial \log p(X|Y)}{\partial_l x}
   ~=~  \left( \sum_{C \in  Cliques(G'), x\in C} 
        \frac{ \partial \log f_C (C)}{\partial_l x} \right) -    
    \expec{X|Y} 
       {\sum_{C \in  Cliques(G'), x\in C} \frac{ \partial \log f_C(C)} 
     {\partial_l x} } ~.
\end{equation}
The difficulty here is computing the expected value
in the formula, which comes from the normalizing constant.
Indeed, this computation forms the core of 
the early Boltzmann machine
algorithm \cite{Hertz.krogh.palmer}.  
In general, this must be done using
something like Gibbs sampling and the techniques of
Section~\ref{ssct-gibbs} can be applied directly.

\subsection{Decomposing learning problems}
\label{ssct-decomp}

Learning problems can be decomposed into sub-problems
in some cases.
While the material in this section applies generally to these sorts
of decompositions,
this section considers one simple example and then
proves some general
results on problems decomposition.
Problem decompositions can also be recomputed on the fly
to create a search through a space of models
that takes advantage of decompositions that exist.
A general result is also presented on incremental decomposition.
These results are simple applications of known methods for
testing independence \cite{Frydenberg.scand,Laur.dawid.networks},
with some added complication
because of the use of plates.

Consider the simple learning problem given in
Section~\ref{sct-bl},
Figure~\ref{fig-mods1and2} over two multinomial variables
$var_1$ and $var_2$, and two Gaussian variables $x_1$ and $x_2$.
For this problem we have specified two alternative
models, model $M_1$ and model $M_2$.
Model $M_2$ has an additional arc going from the
discrete variable $var_2$
to the real valued variable $x_1$.
We will use this subsequently to discuss local search of these models
evaluated by their evidence.

A manipulation of the conditional 
distribution for this model,
making use of Lemma~\ref{lem-remgiven},
yields, for model $M_1$, the conditional distribution
given in Figure~\ref{fig-mod1simp}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{mods1simp.ps}{4.5in}
\end{center}
\vskip -0.2in
\caption{A simplification of model $M_1$}
\label{fig-mod1simp}
\end{figure}
When parameters, $\theta_1$, $\theta_2$ are {\it a priori}
independent, and their data likelihoods
do not introduce cross terms between them,
the parameters become {\it a posteriori} independent as well.
This occurs for
$\theta_1$, $\theta_2$, and the set $\{\mu_1,\sigma_1\}$.
This model simplification also implies the evidence for
model $M_1$ decomposes similarly.
Denote the sample of the variable $x_1$ as
$x_{1,*} = x_{1,1},\ldots,x_{1,N}$, 
and likewise for $var_1$ and $var_2$.
In this case, the result is:
\begin{equation}
\label{eq-m1fac}
evidence(M_1) ~=~ 
p(var_{1,*}|M_1) \,
p(var_{2,*}|var_{1,*},M_1) \,
p(x_{1,*}|var_{1,*},M_1) \,
p(x_{2,*}|x_{1,*},var_{1,*},M_1) ~.
\end{equation}
The evidence for model $M_2$ is similar except that the posterior
distribution of $\mu_1$ and $\sigma_1$ is replaced by
the posterior distribution for $\mu_1'$ and $\sigma_1'$.

This result is general, and applies to
Bayesian networks, undirected graphs, and more generally to chain
graphs.  Similar results are covered
by  Dawid and Lauritzen 
\citeyear{Dawid.lauritzen.anstat} for a family
of models they call hyper-Markov.
The general result described above
is an application of the rules of independence applied
to plates.
This uses the notion of non-deterministic children and parents
introduced in Definition~\ref{def-nd}.
It also requires a notion of local dependence,
which is called the Markov blanket, following Pearl
\citeyear{Pearl.88}, since it is a generalization
of the equivalent set for Bayesian networks.
\begin{defn}
We have a chain graph $G$ without plates.
The \myindex{Markov blanket} of a node $u$ is
all neighbors, non-deterministic parents,
non-deterministic children, and non-deterministic parents
of the children and their chain components:
\begin{eqnarray}
\label{eq-markovblanket}
Markov\mbox{-}blanket(u) &=& 
neighbors(u)\cup ndparents(u)\cup ndchildren(u) \\
\nonumber
        && \cup ~
        ndparents(chain\mbox{-}components(ndchildren(u))) ~.
\end{eqnarray}
\end{defn}
From Frydenberg \citeyear{Frydenberg.scand} it follows that
$u$ is independent of the other non-deterministic
variables in the graph $G$
given the Markov blanket.

To perform the simplification depicted in Figure~\ref{fig-mod1simp},
it is sufficient then to find the finest partitioning
of the model parameters such that they are independent.
The decomposition in Figure~\ref{fig-mod1simp}
represents the finest such partition of model $M_1$.
The evidence for the model will then factor according to
the partition, as given for model $M_1$ in Equation~(\ref{eq-m1fac}).
For this task there is the following theorem,
depicted graphically in Figure~\ref{fig-decomp2}.
\begin{theorem}
\label{thm-decomp1}
{\em (Decomposition)}.
A model $M$ is represented by a chain graph $G$ with 
plates.
Let the variables in the graph be $X$.
There are $P$ possibly empty subsets of the variables $X$,
$X_i$ for $i=1,\ldots,P$ such that $unknown(X_i)$ is a partition
of $unknown(X)$.
This induces a decomposition of the graph $G$
into $P$ subgraphs $G_i$ where:
\begin{itemize}
\item
the graph $G_i$ contains the nodes $X_i$ and any
arcs and plates occurring on these nodes,  and
\item
the potential functions for cliques in $G_i$ 
are equivalent to those in $G$.
\end{itemize}
The induced decomposition represents the unique finest
 equivalent independence model to the original graph if and only if
$X_i$ for $i=1,\ldots,P$ is the finest collection of sets such that,
when ignoring plates,
for every unknown node $u$ in  $X_i$,
its Markov blanket is also in $X_i$.
This finest decomposition takes $O(|X|^2)$ to compute.
Furthermore,
the evidence for $M$ now becomes a product over each subgraph:
\begin{equation}
\label{eq-modfac}
evidence(M) ~=~ p(known(X_*) | M) ~= ~ f_0\, 
        \prod_i f_i( known(X_{i,*}) ) 
\end{equation}
for some functions $f_i$ (given in the proof).
\end{theorem}

Figure~\ref{fig-decomp2} shows how this decomposition works
when there are unknown nodes.
Figure~\ref{fig-decomp2}(a)
shows the basic problem
and Figure~\ref{fig-decomp2}(b) shows the finest decomposition.
Notice the bottom component cannot be further decomposed
because the variable $x_1$ is unknown.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{decomp2.ps}{5in}
\end{center}
\vskip -0.2in
\caption{The incremental decomposition of a model}
\label{fig-decomp2}
\end{figure}

In some cases, the functions $f_i$ given in 
the Decomposition Theorem in Equation~(\ref{eq-modfac})
have a clean interpretation:
they are equal to the evidence for the subgraphs.
This result can be obtained from the following corollary.
\begin{coroll}
\label{coroll-decomp2}
{\em (Local Evidence)}.
In the context of Theorem~\ref{thm-decomp1},
suppose there exists a set of chain components
$\tau_j$ from the graph ignoring plates
such that $X_j ~=~ \tau_j \cup ndparents(\tau_j)$,
where $unknown(ndparents(\tau_j))=\emptyset$.
Then
\[
f_j( known(X_{j,*}) )  ~=~ 
                p(known(\tau_j)_*| ndparents(\tau_j)_*,M)   ~.
\] 
\end{coroll}
If we denote the $j$-th subgraph by model $M^S_j$,
then this term is the conditional evidence for model $M^S_j$
given $ndparents(\tau_j)_*$.
Denote by $M^S_0$ the maximal subgraph on known variables only
(induced by $cliques_0$ as given in the proof of the 
Decomposition Theorem).
If the condition of Corollary~\ref{coroll-decomp2} 
holds for $M^S_j$ for $j=0,1,\ldots,P$,
then it follows that the evidence for the model $M$ is 
equal to the product of the evidence for each subgraph:
\begin{equation}
\label{eq-modfacfull}
evidence(M)  ~= ~
        \prod_{i=0}^P  evidence(M^S_i) ~.
\end{equation}
This holds in general if the original graph $G$ is a Bayesian network,
as used in learning Bayesian networks 
\cite{Buntine.UAI91,Cooper.mlj}.
\begin{coroll}
\label{cor-modfacfull}
Equation~(\ref{eq-modfacfull}) holds if
the parent graph $G$ is a Bayesian network with plates.
\end{coroll}

In general, we might consider searching through a family of
graphical models.
To do this local search \cite{johnson.etal:focs85}
or numerical optimization can be used
to find high posterior models, or Markov chain Monte Carlo
methods to select a sample of representative models,
as discussed in Section~\ref{ssct-gibbsplates}.
To do this, how to represent a family of models must be shown.
Figure~\ref{fig-modsgen}, for instance,
is similar to the models of Figure~\ref{fig-mods1and2}
except that some arcs are hatched.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{modsgen.ps}{2.5in}
\end{center}
\vskip -0.2in
\caption{A family of models (optional arcs hatched)}
\label{fig-modsgen}
\end{figure}
This is used to indicate that these arcs are optional.
To instantiate a hatched arc they can either be
removed or replaced with a full arc.
This graphical model then represents many different
models, for all $2^4$ possible instantiations
of the arcs.
Prior probabilities for these models could be generated
using a scheme such as in \cite[p54]{Buntine.UAI91}
or \cite{HGC.TR9409}, where
a prior probability is assigned by a domain expert for
different parts of the model, arcs and parameters,
and the prior for a full model found by multiplication.
The family of models 
given by Figure~\ref{fig-modsgen}
includes those of Figure~\ref{fig-mods1and2}
as instances.   %  ??? too many x-refs
During search or sampling, an important property is the
Bayes factor for the two models,
$Bayes\mbox{-}factor(M_2,M_1)$,
as described in Section~\ref{ssct-introlearn}.
Because of the Decomposition Theorem and its corollary, 
the Bayes factor for $M_2$ versus $M_1$ can be found by looking
at local Bayes factors.
The difference between models $M_1$ and $M_2$
is the parent for the variable $x_1$,
\[
Bayes\mbox{-}factor(M_2,M_1) ~=~ 
  \frac{ p(x_{1,*}|var_{1,*}, var_{2,*},M_2) }
  { p(x_{1,*}|var_{1,*},M_1) }  ~.
\]
That is, the Bayes factor can be computed from only considering
the models involving $\mu_1,\sigma_1$ and $\mu_1',\sigma_1'$.

This incremental modification of evidence, Bayes factors,
 and finest decompositions
is also general, and follows directly from the independence test.
A similar property
for undirected graphs is given in \cite{Dawid.lauritzen.anstat}.
This  is developed below for the case of directed arcs
and non-deterministic variables.
Handling deterministic variables will require
repeated application of these results, because several non-deterministic 
variables may be effected when adding a single arc between
deterministic variables.
\begin{lemma}
\label{lem-localbf}
{\em (Incremental decomposition)}.
For a graph $G$ in the context of Theorem~\ref{thm-decomp1},
we have two non-deterministic variables $U$ and $V$
such that $U$ is given.
Consider adding or removing
 a directed arc from $U$ to $V$.
To update the finest
decomposition of $G$,
there is a unique subgraph containing
the unknown variables in $ndparents(chain\mbox{-}component(V))$.
To this subgraph add or delete an arc from $U$ to $V$,
and add or delete $U$ to the subgraph if required.
\end{lemma}

Shaded non-deterministic parents can be added at will to nodes
in a graph, and the finest decomposition remains unchanged
except for a few additional arcs.  
 The use of hatched arcs in these contexts
causes no additional trouble to the decomposition
process.  That is,
the finest decomposition for a graph with plates and hatched
directed arcs is formed  as if the arcs where unhatched directed arcs.
The evidence is adjusted during the search by adding the different
parents as required.

\subsection{Bayes factors for the exponential family}
\label{ssct-select}

To make use of the decomposition results
in a learning system it is necessary to be able to generate Bayes factors or 
evidence
for the component models.
For models in the exponential family, whose normalization constant
is known in closed form,
this turns out to be easy. 
If these exact computations are not available,
various approximation methods can be used 
to compute the evidence or Bayes factors \cite{Kass.Raftery.571};
some are discussed in Section~\ref{ssct-cfa}.

If the conjugate distribution for an exponential
family model and its derivatives can be readily computed,
then the Bayes factor for the model can be found
in closed form.
Along with the above decomposition methods,
this result is an important
basis of many fast Bayesian algorithms considering multiple
models.
It is used explicitly or implicitly in all Bayesian
methods for learning decision trees, directed graphical models
(with discrete or Gaussian variables), and linear regression
\cite{Buntine.AISIII,Spiegel.dawid.ss}.

For instance, if the normalizing constant 
$Z_\theta(\tau)$ in Lemma~\ref{lem-conj}
was known in closed form,
then the Bayes factor can be readily computed.
\begin{lemma}
\label{lem-exactbf}
Consider the context of Lemma~\ref{lem-conj}.
Then the {\em model likelihood} or {\em evidence},
given by $evidence(M) = p(x_1,\ldots, x_N|y_1,\ldots,y_N,M)$,
can be computed as:
\begin{eqnarray*}
evidence(M) &=& \frac{p(\theta|\tau) \prod_{j=1}^N p(x_j|y_j,\theta)}
                        {p(\theta|\tau')}  \\
        &=& \frac{Z_\theta(\tau')}{Z_\theta(\tau)Z_2^N} ~.
\end{eqnarray*}
\end{lemma}
For $y|x \sim$ Gaussian this involves multiplying out the
two sets of normalizing constants for the Gaussian and Gamma
distributions.
The evidence for some common exponential family
distributions is given in
Appendix~\ref{app-ef} in Table~\ref{tab-bf} 

For instance, consider the learning problem given
in Figure~\ref{fig-mod1simp}.
Assume that the variables $var_1$ and $var_2$ are both
binary (0 or 1) and that 
the parameters $\theta_1$ and $\theta_2$ are interpreted as follows:
\begin{eqnarray*}
p(var_1=0|\theta_1) &=& \theta_1 ~, \\
p(var_2=0|var_1=0,\theta_2) &=& \theta_{2,0|0} ~, \\
p(var_2=0|var_1=1,\theta_2) &=& \theta_{2,0|1} ~.
\end{eqnarray*}
If we use Dirichlet priors for these parameters,
as shown in Table~\ref{tab-conj},
then the priors are:
\begin{eqnarray*}
(\theta_1,1-\theta_1) & \sim &
        \mbox{Dirichlet}(\alpha_{1,0},\alpha_{1,1})~,\\
(\theta_{2,0|j},1-\theta_{2,0|j})& \sim &
        \mbox{Dirichlet}(\alpha_{2,0|j},\alpha_{2,1|j})~~~~
        ~~~~~~~~\mbox{for }j=0,1~,
\end{eqnarray*}
where $\theta_{2,0|0}$ is {\em a priori} independent of $\theta_{2,0|1}$.
The choice of priors for these distributions
is discussed in \cite{Box.tiao,Bernardo.Smith.94}.
Denote the corresponding sufficient statistics as
$n_{1,j}$ (equal to the number of data where $var_1=j$)
and $n_{2,j|i}$ (equal to the number of data where
        $var_2=j$ and $var_1=i$).
Then the first two terms of the evidence for model $M_1$,
read directly from Table~\ref{tab-bf}, can be written as:
\begin{eqnarray*}
p(var_{1,*}|M_1) &=& 
        \frac{Beta(n_{1,0}+\alpha_{1,0}, n_{1,1} +\alpha_{1,1})}
           {Beta(\alpha_{1,0},  \alpha_{1,1})}  ~,\\
p(var_{2,*}|var_{1,*},M_1) &=&
        \frac{Beta(n_{2,0|0}+\alpha_{2,0|0}, n_{2,1|0} +\alpha_{2,1|0})}
           {Beta(\alpha_{2,0|0}, \alpha_{2,1|0})}
        \frac{Beta(n_{2,0|1}+\alpha_{2,0|1}, n_{2,1|1} +\alpha_{2,1|1})}
           {Beta(\alpha_{2,0|1}, \alpha_{2,1|1})} ~.
\end{eqnarray*}

Assume the variables $x_1$ and $x_2$ are
Gaussian with means given by
\begin{eqnarray*}
 && \mu_{1|0} ~~~~~~~~~~~~\mbox{when $var_1=0$}, \\
 && \mu_{1|1} ~~~~~~~~~~~~\mbox{when $var_1=1$},\\
 && \mu_{2|0,1}  + \mu_{2|0,2} x_1 ~~~~~~~~~~~~\mbox{when $var_1=0$}~, \\
 && \mu_{2|1,1}  + \mu_{2|1,2} x_1 ~~~~~~~~~~~~\mbox{when $var_1=1$} 
\end{eqnarray*}
and variances $\sigma_{1|j}$ and $\sigma_{2|j}$ respectively.
In this case, we split the data set into two parts,
those when $var_1=0$, and those when $var_1=1$.
Each get their own parameters, sufficient statistics,
and contribution to the evidence.
Conjugate priors from Table~\ref{tab-conj} in Appendix~\ref{app-ef} 
(using $y|x \sim$ Gaussian)
are indexed accordingly as:
\begin{eqnarray*}
\mu_{i|j}|\sigma_{i|j}  & \sim &
        \mbox{Gaussian}(\mu_{0,i|j},\frac{\Sigma_{0,i|j}}{\sigma_{i|j}^2}) 
        ~~~~~~~~\mbox{for }i=1,2\mbox{ and }j=0,1~, \\
\sigma_{i|j}^{-2}  & \sim & 
        \mbox{Gamma}(\delta_{0,i|j}/2,\beta_{0,i|j}) 
        ~~~~~~~~\mbox{for }i=1,2\mbox{ and }j=0,1~. \\
\end{eqnarray*}
Notice that $\Sigma_{0,i|j}$ is one-dimensional when $i=0$ and
two-dimensional when $i=2$.
Suitable sufficient statistics for this situation are read from Table~\ref{tab-conjpost}
by looking at the data summaries used there.
This can be simplified for $x_1$  because $d=1$ and $y_1$ for the Gaussian is
uniformly 1.
Thus the sufficient statistics for $x_1$ become
the means and variances
for the different values of $var_1$.
Denote
$\overline{x_{1|0}}$ and $\overline{x_{1|1}}$ as the sample means of
$x_1$ when $var_1=0,1$, respectively,
and $s^2_{1|0}$ and $s^2_{1|1}$ their corresponding sample variances.
This cannot be done for the second case, so we
use the notation from Table~\ref{tab-conjpost},
where $\Sigma,\overline{\mu},\beta$ from
Table~\ref{tab-conjpost} become,
respectively, $\Sigma_{2|j},\overline{\mu_{2|j}},\beta_{2|j}$.
Change the vector $y$ to $(1,x_1)$ when making the calculations
indicated here.
The sufficient statistics are, for each case of $var_1=j$:
\begin{eqnarray*}
S_{2|j} &=&  \sum_{i=1}^N 1_{var_{1,i}=j} \, y_iy_i^\dagger ~, \\
m_{2|j} &=& \sum_{i=1}^N 1_{var_{1,i}=j} \, x_iy_i ~, \\
s^2_{2|j} &=&  \sum_{i=1}^N 1_{var_{1,i}=j} \,
                (x_i-\overline{\mu_{2|j}}^\dagger y_i)^2 ~.
\end{eqnarray*}
The evidence for the last two terms can now be read from Table~\ref{tab-bf}.
This becomes:
\begin{eqnarray*}
p(x_{1,*}|var_{1,*},M_1) &=&
        \prod_{j=0,1}
    \frac{\sqrt{\Sigma_{0,1|j}} }{\pi^{n_{1,j}/2}\sqrt{\Sigma_{0,1|j}+n_{1,j}}}
                \frac{\Gamma((\delta_{0,1|j}+n_{1,0})/2)}
                  {\Gamma(\delta_{0,1|j}/2)\beta_{0,1|j}^{\delta_{0,1|j}/2}  }
                         \\
        && ~~~~~~~~
        \left(\beta_{0,1|j}+s^2_{1|j}+
                \frac{\Sigma_0 N}{\Sigma_0+N}(\overline{x}-\mu_{0,1|j})^2
                 \right)^{(\delta_{0,1|j}+n_{1,j})/2} ~,\\
p(x_{2,*}|x_{1,*},var_{1,*},M_1) &=&      
          \prod_{j=0,1}
    \frac{\det^{1/2}{\Sigma_{0,2|j}} }
        {\pi^{n_{1,j}/2}\det^{1/2}\Sigma_{2|j}}
                \frac{\Gamma((\delta_{0,2|j}+n_{1,0})/2)
                        \beta_{2|j}^{(\delta_{0,2|j}+n_{1,j})/2} }
                  {\Gamma(\delta_{0,2|j}/2)\beta_{0,2|j}^{\delta_{0,2|j}/2}  } ~.
\end{eqnarray*}
The final simplification of the model is given in Figure~\ref{fig-mod1simpf}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{mods1simpsuff.ps}{5in}
\end{center}
\vskip -0.2in
\caption{The full simplification of model $M_1$}
\label{fig-mod1simpf}
\end{figure}

\section{Approximate methods on graphical models}
\label{sct-approx}

Exact algorithms for learning of any reasonable size
invariably involve the recursive arc reversal theorem of 
Section~\ref{ssct-suff}.
Most learning methods, however, use approximate algorithms at some level.
The most common uses of the exponential family
within approximation algorithms were summed up in
Figure~\ref{fig-generic}.  
Various other methods for inference on plates
can be applied either at the
model level or the parameter level:
Gibbs sampling, first described in Section~\ref{ssct-gibbs},
other more general Markov chain Monte Carlo algorithms,
EM style algorithms  \cite{Dempster.etal.EM},
and various closed form approximations
such as the mean field approximation, 
and the Laplace approximation
\cite{Berger.85,Azevedo.UAI94}.
This section summarizes the main families of these
approximate methods.

\subsection{Gibbs sampling}
\label{ssct-gibbs}

Gibbs sampling is the basic tool of simulation
and can be applied to most probability distributions
\cite{Geman.geman.pami,Gilks.jrssb,Ripley.ssim}
as long as the full joint has no zeros
(all variable instantiations are possible).
It is a special case of the general Markov chain Monte Carlo
methods for approximate inference \cite{Ripley.ssim,Neal.review}.
Gibbs sampling can be applied to virtually any graphical model
whether there are plates, undirected or directed arcs,
and whether the variables are real or discrete.
Gibbs sampling does not apply to graphs with
deterministic  nodes, however, since these put zeroes
in the full joint.
This section describes Gibbs sampling without plates,
as a precursor to discussing
Gibbs sampling with plates in Section~\ref{ssct-gibbsplates}.
On challenging problems, other forms of Markov chain Monte Carlo
sampling can and should be tried.  The literature is extensive.

Gibbs sampling corresponds to a probabilistic 
version of gradient ascent,
although their goals of averaging as opposed to
maximizing are fundamentally different.
Gradient ascent
in real valued problems corresponds to simple methods
from function optimization \cite{Gill.murray.etal.81}
and in discrete problems corresponds to local repair
or local search \cite{johnson.etal:focs85,minton:aaai90,selman:aaai92}.
Gibbs sampling varies gradient ascent by introducing
a random component.  The algorithm usually tries to ascend,
 but will sometimes descend, as a strategy for exploring
further around the search space.  So the algorithm tends to wander
around local maxima with occasional excursions to other regions of the
space.
Gibbs sampling is also the core algorithm of simulated annealing
if temperature is held equal to one \cite{vanLaar}.

To sample a set of variables
$X$ according to some non-zero distribution $p(X)$,
initialize $X$ to some value
and then repeatedly resample each variable $x \in X$ according
to its conditional probability  $p(x|X-\{x\})$.
For the simple medical problem of Figure~\ref{fig-smp},
suppose the value of symptoms is known, and 
the remaining variables are to be sampled, then do as follows:
\begin{enumerate}
\item
Initialize the remaining variables somehow.
\item
Repeat the following for $i=1,2,3,\ldots$,
and record the sample of $Age_i, Occ_i, Clim_i, Dis_i$
at the end of each cycle.
\begin{enumerate}
\item
Reassign $Age$ by sampling it according to the
conditional:
\[
p(Age|Occ, Clim, Dis, Symp)~.
\]
That is, take the values of $Occ, Clim, Dis, Symp$ as given
and compute the resulting conditional
 distribution on $Age$.
Then sample $Age$ according to that distribution.
\item
Reassign $Occ$ by sampling it according to the
conditional:
\[ p(Occ|Age, Clim, Dis, Symp)~.\]
\item
Reassign $Clim$ by sampling it according to the
conditional:
\[ 
p(Clim|Age, Occ, Dis, Symp)~.
\]
\item
Reassign $Dis$ by sampling it according to the
conditional:
\[
p(Dis|Age, Clim, Occ, Symp)~.
\]
\end{enumerate}
\end{enumerate}
This sequence of steps is depicted in Figure~\ref{fig-diseasegibbs}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{diseasegibbs.ps}{6in}
\end{center}
\vskip -0.2in
\caption{Gibbs sampling on the medical example}
\label{fig-diseasegibbs}
\end{figure}
In this figure, the basic graph has been rearranged for
each step to represent the dependencies that arise
during the sampling process.  This uses the arc
reversal and conditioning operators introduced previously.

The effect of sampling is not immediate.
$Age_2, Occ_2, Clim_2, Dis_2$ is
conditionally dependent on $Age_1, Occ_1, Clim_1, Dis_1$,
and in general so is $Age_i, Occ_i, Clim_i, Dis_i$
for any $i$.
However the effect of the sampling scheme is that in
the long run, for large $i$,
$Age_i, Occ_i, Clim_i, Dis_i$ is approximately
generated according to $p(Age, Occ, Clim, Dis | Symp)$
independently of $Age_1, Occ_1, Clim_1, Dis_1$.
In Gibbs sampling, all the conditional
sampling is done in accordance with the original distribution,
and since this is a stationary
process, in the long run the samples converge to the
stationary distribution or fixed-point of the process.
Methods for making subsequent
samples independent are known as regenerative simulation
\cite{Ripley.ssim} and correspond to sending the temperature
back to zero occasionally.

With this sample different quantities such as
the probability a patient will have $Age>20$
and $Clim=tropical$ given $Symp$ can be estimated.
This is done by looking at the frequency of this event
in the generated sample.
The justification for this is 
the subject of Markov process theory \cite[Theorem~2.26]{Cinlar.sp}.
The following result, presented informally,  applies:
\begin{comment}
\label{comm-gibbs}
Let $x_1,x_2,\ldots,x_I$ be a sequence of 
discrete variables from a
Gibbs sampler for the distribution $p(x)>0$.
Then the average of $g(x_i)$ approaches the expected
value with probability 1 as $I$ approaches infinity:
\[
\frac{1}{I}\sum_{i=1}^I g(x_i) ~\longrightarrow ~
   \overline{g(x)}=\expec{}{g(x)}  ~.
\]
Further, for a second function $h(x_i)$,
the ratio of two sample averages for $g$ and $h$  approaches their 
``true'' ratio:
\[
\frac{\sum_{i=1}^I g(x_i)} {\sum_{i=1}^I h(x_i)} ~\longrightarrow ~
   \frac{\overline{g(x)}} {\overline{h(x)}}  ~.
\]
This is used to approximate conditional expected values.
\end{comment}
To complete this procedure it is necessary to know how many Gibbs
samples to take, how large to make
$I$, and how to estimate the error in the estimate.
Both these questions have no easy answer
but heuristic strategies exist \cite{Ripley.ssim,Neal.review}.
%   ?????

For Bayesian networks this scheme is  easy in general since the 
only requirement 
when sampling from $p(x|X-\{x\})$ is
the conditional distribution for nodes connected to $x$,
and the global probabilities do not need to
be calculated.
Notice, for instance, that
 in Figure~\ref{fig-diseasegibbs} some sampling
operations do not require all five variables.
The general form for Bayesian networks given
in Equation~(\ref{eq-dags}) goes as follows:
\[
p(x|X-\{x\}) ~=~ \frac{p(X)}{\sum_x p(X)}
        ~=~ \frac{ p(x|parents(x)) \, \prod_{y~:~x \in parents(y)}
                        p(y|parents(y)) }
         { \sum_x p(x|parents(x)) \, \prod_{y~:~x \in parents(y)}
                        p(y|parents(y)) } ~.
\]
Notice the product is over a subset of variables.
Only include conditional distributions for 
variables that have $x$ as a parent.
Thus, the formula only involves examining the parents, children and children's parents of $x$,
the so-called Markov blanket \cite{Pearl.88}.
Also, notice normalization is only required 
over the single dimension
changed in the current cycle,
done in the denominator.
For $x$ discrete, these conditional probabilities can be
enumerated and direct sampling done for $x$.

The kind of simplification above for
Bayesian networks also applies to undirected graphs and chain graphs,
with or without plates.
Here, modify Equation~(\ref{eq-ugs}):
\begin{equation}
\label{eq-ugsgibbs}
p(x|X-\{x\}) ~=~ \frac{ \exp 
     \left( \sum_{C~:~x\in C \in  cliques(G)} f_C (C) \right)}
        {\sum_x  \exp \left( \sum_{C~:~x\in C \in  cliques(G)} f_C (C) \right)} ~.
\end{equation}
In this formula, ignore all cliques not containing $x$
so, again, Gibbs sampling only computes with information local to the node.
Also, the troublesome normalization constant
does not have to be computed
because the probability is a ratio of functions
and so cancels out.
As before, normalization is only required over the single dimension $x$.

\subsection{Gibbs sampling on plates}
\label{ssct-gibbsplates}

Many learning problems can be represented
as Bayesian networks.
For instance, the simple unsupervised learning problem
represented in Figure~\ref{fig-sbnl}
is a Bayesian network once the plate is expanded out.
It follows that Gibbs sampling is readily applied to
learning as a general inference algorithm
\cite{Gilks.jrssb,Gilks.etal.stat}.

Consider a simplified example of this unsupervised learning problem.
In this model, assume that each variable $var_1$ and $var_2$
belongs to a mixture of Gaussians of known variance equal to 1.0.
This simple model is given in 
Figure~\ref{fig-auto2}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{auto2.ps}{2in}
\end{center}
\vskip -0.2in
\caption{Unsupervised learning in two dimensions}
\label{fig-auto2}
\end{figure}
For a given class, $class=c$, the variables
$var_1$ and $var_2$ are distributed as Gaussian with
means $\mu_{1,c}$ and $\mu_{2,c}$.
In the uniform, unit-variance case the distribution for each sample 
is given by:
\[
p(var_1,var_2|,\phi,\mu,M) ~=~
     \sum_c \phi_c N(var_1-\mu_{1,c}) N(var_2-\mu_{2,c}) ~,
\]
where $N(,)$ is the one-dimensional  
Gaussian probability density function with standard deviation of 1.
This model might seem trivial, but if 
the standard deviation were to vary as well, 
the model corresponds to a Kernel density estimate
so can approximate any other distribution arbitrarily well
using sufficient number of tiny Gaussians.

In this simplified Gaussian mixture model, 
the sequence of steps for Gibbs sampling goes as follows:
\begin{enumerate}
\item
Initialize the variables $\phi_c, \mu_{1,c},\mu_{2,c}$ 
for each class $c$.
\item
Repeat the following 
and record the sample of $\phi_c, \mu_{1,c},\mu_{2,c}$ for each class $c$
at the end of each cycle.
\begin{enumerate}
\item
For $i=1,\ldots, N$, reassign $class_{i}$
according to the conditional:
\[
p(class_i ~|~var_{1,i},var_{2,i},\phi,\mu_1,\mu_2)
        ~.
\]
\item
Reassign the vector $\phi$ by sampling 
according to the
conditional:
\[
p(\phi|class_i~:~i=1,\ldots, N)~.
\]
\item
Reassign the vector $\mu_{1}$ 
(and $\mu_{2}$) by sampling 
according to the
conditional:
\[
p(\mu_{1}|var_{1,i},class_i~:~i=1,\ldots, N)~.
\]
\end{enumerate}
\end{enumerate}
Figure~\ref{fig-auto2gibbs}(a) illustrates Step~2(a)
in the language of graphs.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{auto2gibbs.ps}{4.5in}
\end{center}
\vskip -0.2in
\caption{Gibbs sampling in the unsupervised learning problem}
\label{fig-auto2gibbs}
\end{figure}
Figure~\ref{fig-auto2gibbs}(b) illustrates Steps~2(b) and~2(c).
Step~2(a) represents the standard sampling operation
using inference on Bayesian networks without plates.
Step~2(b) and~2(c) are also easy to perform because
in this case the distributions 
are exponential family, and the graph matches
the conditions for Lemma~\ref{cor-modfacfull}.
Therefore, each of the model parameters
$\phi$, $\mu_1,\mu_2 $ are {\it a posteriori}
independent and their distribution is known in closed form,
with the sufficient statistics calculated in $O(N)$ time.

One important caveat in the use of Gibbs sampling for
learning is the problem of symmetry.
In the above description, there is nothing to distinguish
class~1 from class~2.
Initially, the class centers $\mu$ for the above process
will remain distinct.
Asymptotically, since there is nothing in the problem definition
to distinguish between class~1 and class~2, they will appear
indistinguishable.
This problem is handled by symmetry breaking:
force $\mu_{1,c} < \mu_{2,c}$.

Gibbs sampling applies whenever there are variables
associated with the data that are not given.
Hidden or latent variables are an example.
Incomplete data (or missing values) \cite{Quinlan.ML6},
robust methods and modeling of outliers,
and various density estimation and non-parametric methods
all fall in this family of models
\cite{Titter.smith.85}.
Gibbs sampling generalizes to virtually any graphical model with
plates and unshaded nodes inside the plate;  the sequence
of sampling operations will be much the same as in
Figure~\ref{fig-auto2gibbs}.
If the underlying distribution is exponential family,
for instance, Lemma~\ref{lem-exp}
applies after shading all nodes inside the plate;
each full cycle is guaranteed to be linear time  in the sample size.
The algorithm in the exponential family case is summed up in
Figure~\ref{fig-gibbsgen}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{gibbsgen.ps}{5in}
\end{center}
\vskip -0.2in
\caption{Gibbs sampling with the exponential family}
\label{fig-gibbsgen}
\end{figure}
Figure~\ref{fig-gibbsgen}(I) shows the general learning problem,
extending the mixture model of
Figure~\ref{fig-generic}.  This same structure appears in
Figure~\ref{fig-auto2}. 
However, in Figure~\ref{fig-gibbsgen}(I) the sufficient statistics
$T(x_*,u_*)$
are also shown.
Figure~\ref{fig-gibbsgen}(II) shows more of the algorithm, 
generalizing Figure~\ref{fig-auto2gibbs}(a) and (b).
Again the role of sufficient statistics is shown.
The sampling in Figure~\ref{fig-auto2gibbs}(b) first
computes the sufficient statistics and then sampling 
applies from that. 

Thomas, Spiegelhalter, and Gilks
\citeyear{Spiegelhalter.BayesS4}\cite{Gilks.etal.stat}
have taken advantage of this general applicability of sampling
to create a compiler
that converts a graphical representation of a data analysis problem,
with plates, into a matching Gibbs sampling algorithm.
This scheme applies to a broad variety
of data analysis problems.

\subsection{A closed form approximation}
\label{ssct-cfa}

What would happen to Gibbs sampling if the number of cases in the training sample,
$N$, was large compared to the number of unknown parameters
in the problem?
A good way to think of this is:
first, if $N$ is sufficiently large, then
the samples of the model parameters $\phi$ and
$\mu_1,\mu_2$ will tend to drift around their
mean because their posterior variance would be $O(I/N)$.
That is, after the $i$-th step, the 
sample is $\phi_i,\mu_{1,i},\mu_{2,i}$.
The $i+1$-th sample $\phi_{i+1},\mu_{1,i+1},\mu_{2,i+1}$
would be conditionally dependent on these,
but because $N$ is large, the posterior variance
of the $i+1$-th sample given the $i$-th sample
would be small so that:
\[
\phi_{i+1}
~\approx ~  \expec{\phi_i,\mu_{1,i},\mu_{2,i}}{\phi_{i+1} } ~.
\]
This  approximation is
used with Markov chain
Monte Carlo methods
by the \myindex{mean field method} from statistical physics,
popular in neural networks \cite{Hertz.krogh.palmer}.
Rather than sampling a sequence of parameters
$\theta$,
$\theta_1,\theta_2,\ldots, \theta_i$,
according to some scheme, use the deterministic update:
\begin{equation}
\label{eq-meanfield}
\overline{\theta_{i+1}}
~=~  \expec{\theta_i=\overline{\theta_i}}{\theta_{i+1} } ~,
\end{equation}
where the expected value is according to the sampling distribution.
This instead generates a deterministic
sequence 
$\overline{\theta_1},\overline{\theta_2},\ldots, \overline{\theta_i}$
that under reasonable conditions converges to some maxima.
This kind of approach leads naturally to the EM
algorithm, which will be discussed  in Section~\ref{ssct-EM}.

\subsection{The expectation maximization (EM) algorithm}
\label{ssct-EM}

The expectation maximization algorithm,
widely known as the EM algorithm,
corresponds to a deterministic version of Gibbs sampling
used to search for the MAP estimate
for model parameters \cite{Dempster.etal.EM}.
It is generally considered to be faster than gradient descent.
Convergence is slow near a local maxima so some
implementations switch to conjugate gradient
or other methods \cite{meilijson:jrssb} when near a solution.
The computation used to find the derivative is similar to the
computation used for the EM algorithm, so
this does not require a great deal of additional code.
Also, the determinism means the EM algorithm
 no longer generates
unbiased  posterior estimates of model parameters.
The intended gain is speed, not accuracy.
The EM algorithm can generally be applied to exponential
family models wherever
Gibbs sampling can be applied.  The correspondence between EM and Gibbs
is shown below.

Consider again the simple unsupervised learning problem
represented in Figure~\ref{fig-auto2}
(Section~\ref{ssct-gibbsplates}).
In this case, the sequence of steps for the
EM algorithm is similar to that for the Gibbs sampler.
The EM algorithm works
on the means or modes of unknown variables instead of sampling
them.  
Rather than sampling the set of classes
and thereby computing sufficient statistics so that a distribution for
$\phi$ and $\mu_1,\mu_2$ can be found, 
a sequence of class means are generated and thereby
used to compute \myindex{expected sufficient statistics}.
Likewise, instead of sampling new parameters $\phi$ and
$\mu_1,\mu_2$, modes are then computed
from the expected sufficient statistics.

Consider again the unsupervised learning problem in Figure~\ref{fig-sbn}.
Suppose there are 10 classes and that
the three variables $var_1, var_2, var_3$ are finite valued
and discrete and modeled with a  multinomial with probabilities
conditional on the class value $class_i$.
The sufficient statistics in this case are all counts:
$n_j$ is the number of cases where the class is $j$;
$n_{v,k|j}$ is the number of cases where $class=j$ and
$var_{v}=k$:
\begin{eqnarray*}
n_{j} &=&   \sum_{i=1}^N 1_{class_i=j} ~, \\
n_{v,k|j} &=& \sum_{i=1}^N 1_{class_i=j}1_{var_{v,i}=k}   
        ~.
\end{eqnarray*}
The {\em expected} sufficient statistics
computed from the rules of probability  for a given set of
parameters $\phi$ and
$\theta_1,\theta_2,\theta_3$ are given by:
\begin{eqnarray*}
\overline{n_{j}} &=&   \sum_{i=1}^N p(class_i=j|var_{1,i},
                            var_{2,i},var_{3,i},
                                \phi,\theta_1,\theta_2,\theta_3) ~, \\
\overline{n_{v,k|j}} &=& \sum_{i=1}^N p(class_i=j|var_{1,i},var_{2,i},var_{3,i},
                                \phi,\theta_1,\theta_2,\theta_3) 
                                        \,1_{var_{v,i}=k}  ~.
\end{eqnarray*}
Thanks to Lemma~\ref{lem-efsolve},
these kinds of expected sufficient statistics can be computed 
for most exponential family distributions.
Once sufficient statistics are computed for any of the distributions
posterior means or modes of the model parameters
(in this case $\phi$ and
$\theta_1,\theta_2,\theta_3$)
can be found.
\begin{enumerate}
\item
Initialize the parameters $\phi$ and
$\theta_1,\theta_2,\theta_3$.
\item
Repeat the following 
until some convergence criteria is met:
\begin{enumerate}
\item
Compute the expected sufficient statistics 
$\overline{n_{j}}$  and $\overline{n_{v,k|j}} $.
\item
Recompute $\phi$ and
$\theta_1,\theta_2,\theta_3$ to be equal to their mode 
conditioned on the sufficient statistics.
For many posterior distributions,
these can be found in standard tables,
and in most cases found via Lemma~\ref{lem-efsolve}.
For instance, using the mean for $\phi$ gives:
\[
\phi_j ~=~ \frac{ \overline{n_{j}} + \alpha_j}
                {\sum_j (\overline{n_{j}} + \alpha_j) } ~.
\]
\end{enumerate}
\end{enumerate}

All of the other Gibbs sampling algorithms discussed in  
Section~\ref{ssct-gibbsplates}
can be similarly placed in this EM framework.
When the mode is used in Step~2(b), ignoring numerical
problems,
the EM algorithm converges on a local maxima of the
posterior distribution for the parameters \cite{Dempster.etal.EM}.
The general method is summarized in the following comment
\cite{Dempster.etal.EM}.
\begin{comment}
The conditions of Lemma~\ref{lem-exp} apply with data variables $X$
inside the plate 
and model parameters $\theta$ outside.
In addition, some of the variables $U\subset X$ are latent,
so they are unknown and unshaded.
Some of the remaining variables are sometimes missing,
so, for the data $X_i$, variables $V_i\subset ( X-U)$ are
not given.
This means the data given for the $i$-th datum is
$X-U-V_i$ for $i=1,\ldots N$.
The EM algorithm goes as follows:
\begin{description}
\item{E-step: }
The contribution to the expected sufficient statistics for each datum is:
\[
ET_i ~=~  \expec{U_i,V_i|X-U_i-V_i,\theta}{ t(X_i)}~.
\]
The expected sufficient statistic is then $ET ~=~ \sum_{i=1}^N  ET_i$.
\item{M-step: }
Maximize the conjugate posterior using the expected sufficient
statistics $ET$ in place of the sufficient statistics
using the MAP approach for this distribution.
\end{description}
The fixed point of this algorithm is a local maxima of the
posterior for $\theta$.
\end{comment}
Figure~\ref{fig-gibbsgen}(II)  for Gibbs sampling
in Section~\ref{ssct-gibbsplates}
illustrates the steps in EM nicely.
Figure~\ref{fig-gibbsgen}(II)(a) corresponds to the
E-step.  The expected sufficient statistics are found from
the parameters $\theta$.   Rather than sampling $u$,
and therefore computing the sufficient statistics,
the expected sufficient statistics are computed.
Figure~\ref{fig-gibbsgen}(II)(b) corresponds to the M-step.
Here, given the sufficient statistics, the mean or mode
of the parameters $\theta$ are computed instead of being
sampled.
EM is therefore Gibbs with a mean/mode approximation
done at the two major sampling steps of the algorithm.

In some cases, the expected sufficient statistics can be computed in
closed form.
Assume the exponential family distribution for $p(X|\theta)$
has a known normalization constant $Z(\theta)$ and
the link function $w^{-1}$ exists.
For some $i$, the normalizing constant for the exponential
family distribution  $p(U_i,V_i|X-U_i-V_i,\theta)$
is known in closed form.
Denote it by $Z_i(\theta)$.
Then using the notation of Theorem~\ref{thm-plate}:
\[
ET_i ~=~  \left\{  \begin{array}{lr}
                t(X_i) & \mbox{if } U=V_i=\emptyset ~, \\
                \frac{\diff w(\theta)}{\diff \theta}^{-1}  
                        \frac{\diff \log Z_i(\theta)}{\diff \theta} & \mbox{otherwise} ~.
        \end{array} \right.
\]

\section{Partial exponential models}
\label{ssct-partial}

In some cases, only an initial, inner part of a learning problem can be
handled using the recursive arc reversal 
theorem of Section~\ref{ssct-expon}.
In this case,
simplify what can be simplified, and then solve the
remainder of the problem using a  generic method
like the MAP approximation.
This section presents several examples:
linear regression with heterogeneous variance,
feed-forward networks with a linear output layer,
and Bayesian networks.

This general process was depicted in the graphical
model for the partial exponential family
of Figure~\ref{fig-generic}.
This is an abstraction used to
represent the general process.
Consider the problem of learning a Bayesian network,
both structure and parameters,
where the distribution is exponential family given the
network structure.
The variable $\Tau$ is a discrete variable indicating
which graphical structure is chosen for the Bayesian network.
The variable $X$ represents the full set of variables
given for the problem.
The variable $\Theta_\Tau$ represents the distributional parameters
for the Bayesian network, and is the part of the model that
is conveniently exponential family.
That is, $p(X|\Theta_\Tau,\Tau)$ will be treated as exponential for different
        $\Theta_\Tau$ and hold $\Tau$ fixed.
The sufficient statistics in this case are given by $ss(X*,\Tau)$.
In this case, the subproblem that conveniently
falls in the exponential family, 
$p(\Theta_{\Tau} | X*, \Tau)$ is simplified,
but it is necessary to resort
to the more general learning techniques of previous sections
to solve the remaining part of the problem, $p(\Tau|X*)$.

\subsection{Linear regression with heterogeneous variance}
\label{ssct-linreghet}

Consider the heterogeneous variance problem given in
Figure~\ref{fig-lnhv}.
This shows a graphical model for the linear regression
problem of Section~\ref{ssct-linreg}
modified to the situation
where the standard deviation is heterogeneous,
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{linreghet.ps}{4in}
\end{center}
\vskip -0.2in
\caption{Linear regression with heterogeneous variance}
\label{fig-lnhv}
\end{figure}
so it is a function of the inputs $x$ as well.
In this case, the standard deviation
$s$ is not given but is computed via:
\[
s ~=~ \exp \left( \sum_{i=1}^m weights\mbox{-}\sigma_i \,\,basis_i(x)
        \right) ~.
\]
The
exponential transformation
guarantees that the standard deviation $s$ will also be positive.

The corresponding
learning  model can be simplified to the graph in Figure~\ref{fig-linreghets}.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{linreghets.ps}{4in}
\end{center}
\vskip -0.2in
\caption{The heterogeneous variance problem with the plate simplified}
\label{fig-linreghets}
\end{figure}
Compare this with the model given in Figure~\ref{fig-linregls}.
What is the difference?
In this case, the sufficient statistics
exist, but they are shown to be deterministically
dependent on the sample and ultimately on the unknown
parameters for the standard deviation $weights\dash\sigma$.
If the parameters for the standard deviation were known
then the graph could be reduced to Figure~\ref{fig-linregls}.
Computationally, this is an important gain.
It says that for a given set of values for
$weights\dash\sigma$, calculation can be done that is
linear time in the sample size to arrive at a characterization
of what the parameters from the mean,
$weights\dash\mu$,  should be.
In short, one half of a problem,
$p(weights\dash\mu | weights\dash\sigma, y_i, x_{.,i} ~:~i=1,\ldots , N)$,
is well understood.
To do a search for the MAP  solution,
search the space of parameters for
the standard deviation, $weights\dash\sigma$,
since the remaining ($weights\dash\mu$) is given.

Another variation of linear regression replaces
the Gaussian error function with a more robust error function
such as Student's $t$ distribution, or an $L_q$ norm for
$1<q<2$.
By introducing a convolution, these robust regression models
can be handled by combining the EM algorithm with
standard least squares \cite{Lange.sins.compgraph}.
% ??? check this

\subsection{Feed-forward networks with a linear output layer}
\label{ssct-lffn}

A similar example is the standard feed-forward network
where the final output layer is linear.
This situation is given by Figure~\ref{fig-ffnnl}
if we change the deterministic functions for
$m_1$ and $m_2$ to be linear instead of sigmoidal.
In this case Lemma~\ref{lem-exp}
identifies that when the weight vectors $w_3,w_4$, and $w_5$
are assumed given, the distribution is in the exponential
family.  Thus the simplification to Figure~\ref{fig-ffnnlinsimp}
is possible using the standard sufficient statistics for 
multivariate linear regression.
\begin{figure}[htbp]
\begin{center}
\
\myepsfwd{feedforwllsimp.ps}{4in}
\end{center}
\vskip -0.2in
\caption{Simplified learning of a feed-forward network with linear output}
\label{fig-ffnnlinsimp}
\end{figure}
Algorithmically, this implies that given values for the internal
weight vectors $w_3, w_4$, and $w_5$, and assuming
a conjugate prior holds for the output weight vectors $w_1$ and $w_2$,
the posterior distribution for the output weight vectors $w_1$ and $w_2$,
and their means and variances, can be found in closed form.
The evidence for $w_1$ and $w_2$ given $w_3, w_4$ and $w_5$,
$p(y_*|x_{1,*},\ldots,x_{n,*},w_3, w_4,w_5,M)$,
can also be computed using the exact method of Lemma~\ref{lem-exactbf},
so therefore the posterior for $w_3, w_4$, and $w_5$,
\[
p(w_3, w_4,w_5| y_*,x_{1,*},\ldots,x_{n,*},M) ~\propto~ p(w_3, w_4,w_5|M) \,
        p(y_*|x_{1,*},\ldots,x_{n,*},w_3, w_4,w_5,M) ~,
\]
can be computed in closed form up to a constant.
This effectively cuts the problem into two pieces---$w_3, w_4$ and $w_5$
followed by $w_1$ and $w_2$ given $w_3, w_4$, and $w_5$---and provides a
 clean solution to the second piece.

\subsection{Bayesian networks with missing variables}
\label{ssct-lbn}

Class probability trees and
discrete Bayesian networks can be learned
efficiently by noticing that their basic form is exponential family
\cite{Buntine.AISIII,Buntine.IJCAI91,Buntine.UAI91,Cooper.mlj,Spiegel.dawid.ss}.
Take, for instance, the family of models specified by
the Bayesian network given in Figure~\ref{fig-modsgen}.
In this case, the local evidence corollary,
 Corollary~\ref{coroll-decomp2}, applies.
The evidence for Bayesian networks generated from this graph is
therefore a product over the nodes in the Bayesian network.
If we change a Bayesian network by adding or removing an arc,
the Bayes factor is therefore simply the local Bayes factor
for the node, as mentioned in the
incremental decomposition lemma, Lemma~\ref{lem-localbf}.
Local search is then quite fast,
and Gibbs sampling over the space of Bayesian networks is possible.
A similar situation exists with trees \cite{Buntine.AISIII}.
The same results apply to any Bayesian network 
with exponential family distributions
at each node, such as Gaussian or Poisson.
Results for Gaussians are presented, for instance,
in \cite{Geiger.Heckerman.UAI94}.

This local search approach is a MAP approach because it
searches for the network structure maximizing posterior probability.
More accurate approximation can be done by generating
a Markov chain of Bayesian networks from the search space of Bayesian networks.
Because the Bayes factors are readily computed in this case,
Gibbs sampling or  Markov chain Monte Carlo schemes can be used.
The scheme given below is the 
Metropolis algorithm \cite{Ripley.ssim}.
This only looks at single neighbors
until a successor is found.
This is done be repeating the following steps:
\begin{enumerate}
\item
For the initial Bayesian network $G$, randomly select a neighboring
Bayesian network $G'$ differing only by an arc.
\item
Compute $Bayes\mbox{-}factor(G',G)$ by making the decompositions
described in Theorem~\ref{thm-decomp1},
doing a local computation as described
in Lemma~\ref{lem-localbf}, and using the 
Bayes factors computed with Lemma~\ref{lem-exactbf}.
\item
Accept the new Bayesian network $G'$ with probability given by:
\[
\min \left( 1, Bayes\mbox{-}factor(G',G) \frac{p(G')}{p(G)} \right)~.
\]
If accepted, assign $G'$ to $G$, otherwise $G$ remains unchanged.
\end{enumerate}
A local maxima Bayesian network could be found concurrently,
however, this scheme generates a set of Bayesian networks appropriate for
model averaging and expert evaluation of the space of potential
Bayesian networks.  Of course, initialization might 
 search for local maxima to use as a reference. 
This sampling scheme was illustrated in the context of
averaging in Figure~\ref{fig-multnets}.

This scheme is readily adapted to learn the
structure and parameters of a Bayesian network with missing or latent
variables.  
For the Metropolis algorithm, add Step~4,  which resamples
the missing data and latent variables.
\begin{description}
\item{4.}
For the current complete data and Bayesian network $G$,
compute the predictive distribution for the missing
data or latent variables.
Use this to resample  the missing
data or latent variables to construct a new set
of complete data (for subsequent use in computing Bayes factors).
\end{description} 

\section{Conclusion}
The marriage of learning and graphical models  presented
here provides a framework for understanding learning.
It also provides a framework for
developing a learning or data analysis toolkit,
or more ambitiously, a software generator for learning algorithms.
Such a toolkit combines two
important components:
a language for representing a learning problem
together with techniques
for generating a matching algorithm.
While a working toolbox is not demonstrated,
a blueprint is provided to show how it could be constructed,
and the construction of some
well-known learning algorithms has been demonstrated.
Table~\ref{tab-algs} lists some standard problems,
the derivation of algorithms using the operations
from the previous chapters, and where in the text they are
considered.
\begin{table}
\begin{tabular}{|p{2in}|p{2.5in}|p{1in}|}
\hline\hline
Problem & Method & Sections   \\
\hline
\parbox{2in}{Bayesian networks and exponential
        family conditionals\rule[-3pt]{0cm}{1pt}} &   
         \parbox{2.5in}{\rule[10pt]{0cm}{1pt}Decomposition
        of exact Bayes factors, with local search or Gibbs
        sampling to generate alternative models\rule[-3pt]{0cm}{1pt}} 
        & \parbox{1.5in}{\ref{ssct-decomp}, \ref{ssct-select},
                        \ref{ssct-lbn} }  \\
\hline
\parbox{2in}{\rule[10pt]{0cm}{1pt}Bayesian networks with missing/latent variables, and other
        unsupervised learning models\rule[-3pt]{0cm}{1pt}}
        & \parbox{2.5in}{\rule[10pt]{0cm}{1pt}Gibbs sampling or EM  making use of  above
         techniques where possible\rule[-3pt]{0cm}{1pt}}
                     & \parbox{1.5in}{\ref{ssct-gibbsplates}, 
         \ref{ssct-EM}, \ref{ssct-lbn}} \\
\hline
\parbox{2in}{Feed-forward networks} &  \parbox{2.5in}{\rule[10pt]{0cm}{1pt}MAP method
                by exact computation of derivatives\rule[-3pt]{0cm}{1pt}}
        &  \ref{ssct-deriv} \\
\hline
\parbox{2in}{\rule[10pt]{0cm}{1pt}Feed-forward networks with linear output} &
        \parbox{2.5in}{\rule[10pt]{0cm}{1pt}As above with an initial
                removal of the linear component\rule[-3pt]{0cm}{1pt}}  &
        \ref{ssct-lffn} \\
\hline
\parbox{2in}{\rule[10pt]{0cm}{1pt}Linear regression and extensions\rule[-3pt]{0cm}{1pt}} &
        \parbox{2.5in}{\rule[10pt]{0cm}{1pt}Least squares, EM, and MAP}  & 
      \ref{ssct-linreg}, \ref{ssct-linreghet} \\
\hline
\parbox{2in}{Generalized linear models} &
        \parbox{2.5in}{\rule[10pt]{0cm}{1pt}MAP method
                by exact computation of derivatives\rule[-3pt]{0cm}{1pt}} &
          \ref{ssct-cond}, \ref{ssct-deriv} \\
\hline\hline
\end{tabular}
\caption{Derivation of learning algorithms}
\label{tab-algs}
\end{table}
The notion of a learning toolkit is not new, and
can be seen in the BUGS system by Thomas, Spiegelhalter, and Gilks
\citeyear{Spiegelhalter.BayesS4}\cite{Gilks.etal.stat},
in the work of Cohen \citeyear{cohen:92} for 
inductive logic programming,
and emerging in software for handling generalized linear models
  \cite{Mccullagh.nelder.glm,NewS}.

There is an important role for a data analysis toolkit.
Every problem has its own quirks and requirements.
Knowledge discovery, for instance,
can vary in many ways depending
on the user-defined notion of interestingness.
Learning is often an embedded task in a larger system.
So while there are some easy applications of learning,
generally learning applications require
special purpose development of learning systems or related
support software.
Sometimes, this can be achieved
by patching together some existing techniques
or by decomposing a problem into subproblems.
Nevertheless, the decomposition and patching of 
learning algorithms with
inference and decision making
can be formalized and understood
within graphical models.
In some ways the S system
plays the role of a toolkit \cite{StatsinS}.
It provides a system for prototyping learning algorithms,
includes the ability to handle generalized linear
models, does automatic differentiation of expressions, 
and includes  many statistical and mathematical 
functions useful as primitives.
The language of graphical models is best viewed
as an additional layer on top of this kind of system.
Note, also, that it is impractical to assume that
a software generator could create algorithms competitive
with current finely tuned algorithms, for instance,
for hidden Markov models.  However, a software toolkit for
learning could be used to prototype an algorithm
that could later be refined by hand.

The combination of learning and graphical models
shares some of the superior aspects of each of the different
learning fields.
Consider the philosophy of neural networks.
These nonparametric systems are composed of
simple computational components, usually readily parallelizable,
and often nonlinear.
The components can be pieced together to tailor systems for
specific applications.
Graphical models for learning have these same features.
Graphical models also
have the expressibility of probabilistic knowledge representations
that were developed in artificial intelligence
to be used in knowledge acquisition contexts.
They therefore form an important basis of knowledge refinement.
Finally, graphical models for learning allow the powerful
tools of statistics
to be applied to the problem.

Once learning problems are specified in the common language
of graphical models,
their associated learning algorithms, their derivation,
and their interrelationships can be explored.
This allows commonalities between seemingly diverse pairs of
algorithms---such
as k-means clustering versus approximate methods for learning
hidden Markov models, learning decision trees versus 
learning Bayesian networks
\cite{Buntine.IJCAI91},
and Gibbs sampling versus the expectation maximization algorithm
in Section~\ref{ssct-EM}---to be
understood as variations of one another.
The framework is important as an educational tool.


\appendix

\section{Proofs of Lemmas and Theorems}
\label{app-proofs}

\subsection{Proof of Theorem \protect{\ref{th-ugs}}}
A useful property of independence is that   %???? check with ripley ref.
$A$ is independent of $B$ given $C$ if and only if
$p(A,B,C) =f(A,C)g(B,C)$ for some functions $f$ and $g$.
The only if result follows directly from this property.
The proof of the if result makes use of the following simple lemma.
If $A$ is independent of
$B$ given $X-A-B$, and $p(X) = \prod_i f_i(X_i)$ for
some functions $f_i>0$ and variable sets $X_i\subseteq X$, then:
\begin{equation}
\label{eq-ugprod}
p(X) = \prod_i g_i(X_i -B) \, h_i(X_i-A) 
\end{equation}
for some functions $g_i,h_i>0$.
Notice that it is known $p(X)=g(X-B) \, h(X-A)$
for some functions $g$ and $h$ by independence.
Instantiate the variables in $B$
to some value $b$.
Then:
\[
g(X-B) \, h(X-A,B=b)
   ~=~  \prod_i f_i(X_i,B=b)  ~.
\]
Similarly, instantiate $A$ to $a$, then:
\[
g(X-B,A=a) \, h(X-A) ~=~ \prod_i f_i(X_i,A=a) ~.
\]
Multiplying both sides of the two equalities together, and substitute in 
\[
g(X-B,A=a) \, h(X-A,B=b) ~=~p(X,A=a,B=b)~=~\prod_i f_i(X_i,A=a,B=b)~.
\]
Get:
\[
p(X) \, \prod_i f_i(X_i,A=a,B=b)  ~=~
     \prod_i f_i(X_i,B=b)  \, \prod_i f_i(X_i,A=a)  ~.
\]
This is defined for all $X$ since the domain is a cross product.
The lemma holds because all functions are strictly positive if 
\[
g_i(X_i-B) ~=~ \frac{f_i(X_i,B=b)  }{ f_i(X_i,B=b,A=a) } 
\]
and
\[
h_i(X_i-A) ~=~ f_i(X_i,A=a)~.
\]

The final proof of the if result follows by applying 
Equation~(\ref{eq-ugprod}) repeatedly.
Suppose the variables in $X$ are $x_1,\ldots,x_v$.
Now $p(X)=f_0(X)$ for some strictly positive function $f_0$.
Therefore:
\[
p(X)~=~g_0(X-\{x_1\})h_0(\{x_1\}\cup neighbors(x))~.
\]
Denote $A_{i,0}=X-\{x_i\}$ and $A_{i,1}=\{x_i\}\cup neighbors(x_i)$.
Repeating the application of Equation~(\ref{eq-ugprod}) for each variable
yields:
\[
p(X)~=~\prod_{i_1=0,1}\ldots \prod_{i_v=0,1}
                    g_{i_1,\ldots,i_v} ( X - \bigcup_{j=1,\ldots,v} A_{j,i_j} ) ~.
\]
for strictly positive functions $g_{i_1,\ldots,i_v} $.
Now,  consider these functions.
It is only necessary to keep the function $g_{i_1,\ldots,i_v} $ if the
set $X - \bigcup_{j=1,\ldots,v} A_{j,i_j} $ is maximal:
it is not contained in any such set.
Equivalently, keep the function $g_{i_1,\ldots,i_v} $ if the
set $\bigcup_{j=1,\ldots,v} A_{j,i_j} $ is minimal.
The minimal such sets are the set of cliques on the undirected graph.
The result follows.

\subsection{Proof of Theorem \protect{\ref{thm-decomp1}}}
It takes $O(|X|^2)$ operations to remove 
the deterministic nodes from a graph using Lemma~\ref{lem-nodet}.
These nodes can be removed from the graph
and then reinserted at the end.
Hereafter, assume the graph contains no
deterministic nodes.    
Also, denote the unknown variables in a set $Y$ by 
$unknown(Y)=Y-known(Y)$.
Then, without loss of generality, assume
that $X_i$ contains all known variables
in the Markov blanket of $unknown(X_i)$.

Showing the independence model is equivalent
amounts to showing for $i=1,\ldots,P$ that 
$unknown(X_i)$
is independent of $\bigcup_{j\neq i} unknown(X_j)$
given $known(X)$.
To test independence using the method of Frydenberg
\citeyear{Frydenberg.scand}, 
 each plate must be expanded (that is, duplicate it
the right number of times),
moralize the graph, removing the given nodes, 
and then test for separability.
The Markov blanket for each node in this expanded
graph corresponds to those nodes
directly connected in the moralized expanded graph.
Suppose we have the finest unique
partition $unknown(X_i)$ of the 
unknown nodes.  $X_i$'s are then reconstructed by adding
known variables in the Markov blankets for variables in $unknown(X_i)$.
Suppose $V$ is an unknown variable in a plate,
and $V_j$ are its instances once the plate is expanded.
Now, by symmetry, every $V_j$ is either
in the same element of the finest partition, or they are all
in separate elements.  If $V_j$ has a certain unknown
variable in its Markov boundary outside the plate,
then so must $V_k$ for $k\neq j$ by symmetry.
Therefore $V_j$ and $V_k$ are in the same element of the partition.
Hence by contradiction,
if $V_j$ is in a separate element,
that element occurs wholly within the plate
boundaries.
Therefore, this finest partition
can be represented using plates,
and the finest partition identified
from the graph ignoring plates.
The operation of finding the finest separated sets
in a graph is quadratic in the size of the graph,
hence the $O(|X|^2)$ complexity.

Assume the condition holds and consider Equation~(\ref{eq-modfac}).
Let $cliques(\tau)$ denote the subsets of variables in
$\tau\cup parents(\tau)$ that form cliques in the graph
formed by restricting $G$ to $\tau\cup parents(\tau)$
and placing an undirected arc between all parents.
Let $\tau(X)$ be the set of chain components in $X$.
From Frydenberg \citeyear{Frydenberg.scand}, we have:
\[
p(X| M) ~=~  \prod_{\tau\in \tau(X)}  \prod_{i\in ind(\tau)}
        \prod_{C\in cliques(\tau)} g_C(C_i) ~.
\]
Furthermore, if $u\in X_i$ is not known, then
the variables in $u$'s Markov blanket will occur in 
$X_i$, and therefore, if $u\in C$ for some clique $C$,
then $C\subseteq X_i$.
Therefore cliques containing an unknown
variable can be partitioned according to which subgraph they belong in.
Let:
\[
cliques_j' ~=~ \{ C ~:~C\in cliques(\tau(X)),~
        unknown(X_j)\cap C \neq \emptyset \}~,
\]
and add to this any remaining cliques wholly
contained in the set so far:
\[
cliques_j ~=~ cliques_j' \cup \left\{ C ~:~C\in cliques(\tau(X)),~
        C\subseteq \bigcup_{C'\in cliques_j'}C' \right\}~.
\]
Call any remaining cliques $cliques_0$.
Therefore:
\[
p(X | M) ~=~  
                \prod_{j=0}^P 
        \prod_{C\in cliques_j } \prod_{i\in ind(C)}  
          g_C(C_i) ~.
\]
Results in:
\[
f_j( known(X_{j,*}) ) ~=~ 
        \int_{unknown(X_j,*)} \prod_{C\in cliques_j }  
           \prod_{i\in ind(C)}  g_C(C_i) \diff
                                unknown(X_j,*)   ~.
\]
Furthermore, the potential functions on the cliques 
in $G_i$ are well defined as described.

\subsection{Proof of Corollary \protect{\ref{coroll-decomp2}}}
If $X_j ~=~ \tau_j \cup ndparents(\tau_j)$, then
every clique in a chain component in $\tau_j $
will occur in $cliques_j$.
Therefore:
\begin{eqnarray*}
\prod_{C\in cliques_j } \prod_{i\in ind(C)}  
          g_C(C_i) &=&
 \prod_{\tau\in\tau_j} p(\tau |  ndparents(\tau) )~, \\
  &=& p(\tau_j |  ndparents(\tau_j) ) ~.
\end{eqnarray*}

\subsection{Proof of Lemma \protect{\ref{lem-localbf}}}
Consider the definition of the Markov blanket.
If a directed arc is added between the nodes,
then the Markov blanket will only change for
an unknown node $X$ if $U$ now enters the 
set  of non-deterministic parents
of the chain-components containing 
non-deterministic children of $X$.   
This will not  effect the subsequent graph separability,
however, because it will only subsequently add arcs
between $U$, a given node, and other nodes. 

\section{The exponential family}
\label{app-ef}

The exponential family of distributions 
was described in Definition~\ref{defn-expfam}.
The common use of the exponential family exists because of 
Theorem~\ref{thm-plate}.
Table~\ref{tab-ef} gives a few
exponential family distributions and their functional form.
Further details and more extensive tables
can be found in most Bayesian textbooks on probability
distributions \cite{Degroot.70,Bernardo.Smith.94}.
\begin{table}
\begin{tabular}{|l|p{3.4in}|}
\hline\hline
Distribution & Functional form  \\
\hline\rule[12pt]{0cm}{1pt}
$j~\sim$ $C$-dim multinomial($\theta_{1},\ldots,\theta_{C}$) &
        $\theta_{j}$ for $j\in [1,\ldots,C]$ \newlinewithspace
$y|x~\sim$ Gaussian($x^\dagger\theta,\sigma$) 
        & $\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2\sigma}
        (y-x^\dagger\theta)^2 \right)$ for $y\in\Re, x\in\Re^d$  \newlinewithspace
$x~\sim$ Gamma($\alpha>0,\beta>0$) &
        $\frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}$
                for $x\in\Re^+$ \newlinewithspace
$\theta~\sim$ $C$-dim Dirichlet($\alpha_1,\ldots,\alpha_C$) &
        $\frac{1}{Beta(\alpha_1,\ldots,\alpha_C)}\prod_{i=1}^C 
                \theta_i^{\alpha_i-1}  $ \newlinewithbigspace
$x~\sim$ 
        $d$-dim Gaussian($\mu\in \Re^d,\Sigma\in \Re^{d\times d}$) 
        &   $\frac{\det^{1/2}\Sigma}{(2\pi)^{d/2}}\exp\left(-\frac{1}{2}
        (x-\mu)^\dagger \Sigma (x-\mu) \right)$
                for $x\in\Re^d$ \newlinewithbigspace
$S~\sim$ 
        $d$-dim Wishart($\alpha\geq d,\Sigma\in \Re^{d\times d}$)  
        & \parbox[t]{4in}{\rule[12pt]{0cm}{1pt}
        $\frac{\det^{-\alpha/2} \Sigma \, \det^{\alpha-d-1/2} S}
                { 2^{d\alpha/2}\, \pi^{d(d-1)/4}\prod_{i=1}^d 
                        \Gamma((\alpha+1-i)/2) }
                \exp\left( -\frac{1}{2} \mbox{trace} \Sigma^{-1} S \right)$ 
                \newlinewithspace
        ~~~~~for $S,\Sigma$ symmetric positive definite 
        \rule[-5pt]{0cm}{1pt}} \\
\hline\hline
\end{tabular}
\caption{Distributions and their functional form}
\label{tab-ef}
\end{table}
Table~\ref{tab-conj} 
gives some standard conjugate prior distributions
for those in Table~\ref{tab-ef},
and Table~\ref{tab-conjpost} gives their matching posteriors
\cite{Degroot.70,Bernardo.Smith.94}.
\begin{table}
\begin{tabular}{|l|p{4.1in}|}
\hline\hline
Distribution & Conjugate prior  \\
\hline\rule[12pt]{0cm}{1pt}
$j~\sim$ $C$-dim multinomial &
        $\theta~\sim$ Dirichlet($\alpha_1, \ldots, \alpha_C$) \newlinewithspace
$y|x~\sim$ Gaussian
        &  $\theta|\sigma \sim$ $d$-dim Gaussian($\theta_0,
                        \frac{1}{\sigma^2}\Sigma_0$), 
             $\sigma^{-2}\sim$ Gamma($\alpha_0/2,2/\beta_0$)\newlinewithspace
$x~\sim$ Gamma &
        $ \beta|\alpha~\sim$ Gamma($\alpha_0,\beta_0$) \newlinewithspace
$x~\sim$ $d$-dim Gaussian  &
        $\mu|\Sigma~ \sim$ Gaussian($\mu_0,N_0\Sigma$), 
                $\Sigma~\sim$ Wishart($\delta_0,S_0$) \\
\hline\hline
\end{tabular}
\caption{Distributions and their conjugate priors}
\label{tab-conj}
\end{table}
\begin{table}
\begin{tabular}{|l|p{4.3in}|}
\hline\hline
Distribution & Conjugate posterior   \\
\hline\rule[12pt]{0cm}{1pt}
$j~\sim$ $C$-dim multinomial &
        \parbox{4.3in}{\rule[10pt]{0cm}{1pt}$\theta~\sim$ 
        Dirichlet($n_{1}+\alpha_1, \ldots, n_{C}+\alpha_C$)
        \newlinewithspace
                     ~~~~~~for $n_c=\sum_{i=1}^N 1_{j_i=c}=\#<j's=c>$ 
                     \rule[-5pt]{0cm}{1pt}}\\
\hline
$y|x~\sim$ Gaussian
        &  \parbox{4.5in}{ \rule[12pt]{0cm}{1pt}$\theta|\sigma~ \sim$ $d$-dim 
   Gaussian($\overline{\theta},\frac{1}{\sigma^2}\Sigma$),
        $\sigma^{-2}~\sim$ Gamma($(\alpha_0+N)2,2/\beta$), \newlinewithspace
        ~~~~~~for $\Sigma=\Sigma_0+\sum_{i=1}^Ny_iy_i^\dagger$,
            $\overline{\theta} =
                \Sigma^{-1}\left(\Sigma_0\theta_0+\sum_{i=1}^Nx_iy_i\right)$, \newlinewithspace
        ~~~~~~for $\beta = 
           \sum_{i=1}^N (x_i-\overline{\theta}^\dagger y_i)^2
                + (\overline{\theta}-\theta_0)^\dagger \Sigma_0
                                (\overline{\theta}-\theta_0)
                + \beta_0$               \rule[-5pt]{0cm}{1pt}} \\\hline
              \rule[12pt]{0cm}{1pt}
$x~\sim$ Gamma &
     $ \beta|\alpha~\sim$ Gamma($N\alpha+\alpha_0,\sum_{i=1}^Nx_i+\beta_0$) 
\rule[-5pt]{0cm}{1pt}\\\hline
$x~\sim$ $d$-dim Gaussian  &
          \parbox{4.5in}{  \rule[12pt]{0cm}{1pt}     
        $\mu|\Sigma~ \sim$ Gaussian($\overline{\mu},(N+N_0)\Sigma$), 
        $\Sigma~\sim$ Wishart($N+\delta_0,S+S_0$)  \newlinewithspace
       ~~~~~~for  $\overline{\mu} = \overline{x} + 
                        \frac{N_0}{N+N_0} (\mu_0-\overline{x})$, \newlinewithspace
        ~~~~~~for $S = \sum_{i=1}^N (x_i-\overline{x})(x_i-\overline{x})^\dagger$ 
        \rule[-5pt]{0cm}{1pt} } \\\hline
\hline
\end{tabular}
\caption{Distributions and matching conjugate posteriors}
\label{tab-conjpost}
\end{table}
For the distributions in Table~\ref{tab-ef} with
priors in Table~\ref{tab-conj} ,
Table~\ref{tab-bf} gives their matching evidence
derived using Lemma~\ref{lem-exactbf} and cancelling
a few common terms.
\begin{table}
\begin{tabular}{|l|p{3.75in}|}
\hline\hline
Distribution & Evidence  \\
\hline\rule[12pt]{0cm}{1pt}
$j~\sim$ $C$-dim multinomial &
        $Beta(n_1+\alpha_1, \ldots, n_C + \alpha_C) /
                Beta(\alpha_1, \ldots, \alpha_C)$ \newlinewithspace
$y|x~\sim$ Gaussian
        &  $\frac{\det^{1/2}\Sigma_0 }{\pi^{N/2}\det^{1/2}\Sigma}
                \frac{\Gamma((\alpha_0+N)/2) \beta^{(\alpha_0+N)/2}}
                {\Gamma(\alpha_0/2) \beta_0^{\alpha_0/2} }  $
         \newlinewithbigspace
$x~\sim$ Gamma &
         $\frac{\beta_0^{\alpha_0}\Gamma(N\alpha+\alpha_0)}
           {\Gamma(\alpha_0)
                \left(\sum_{i=1}^Nx_i+\beta_0\right)^{N\alpha+\alpha_0} }$ 
        for $\alpha$ fixed \newlinewithbigspace
$x~\sim$ $d$-dim Gaussian  & 
        $\frac{\det^{\delta_0/2}S_0}
                {(\pi)^{dN/2}\det^{(\delta_0+N)/2}(S+S_0)}
        \frac{N_0^d}{(N+N_0)^d}
        \prod_{i=1}^d\frac{\Gamma((\delta_0+N-1-i)/2)}
                        {\Gamma((\delta_0-1-i)/2)}  $
           \rule[-5pt]{0cm}{1pt}\\
\hline\hline
\end{tabular}
\caption{Distributions and their evidence}
\label{tab-bf}
\end{table}

In the case where the functions $w_i$ are {\em full rank} in
$\theta$
(dimension of $\theta$ is $k$, same as $w$, and the
Jacobian of $w$ with respect to $\theta$ is invertible, 
$\det\left(\frac{\diff w(\theta)}{\diff \theta}\right)\neq 0$),
then various moments of the distribution can be easily found:
\begin{equation}
\label{eq-expeqs1}
        \expec{x|y,\theta}{t(x,y)}
 ~=~\left(\frac{\diff w(\theta)}{\diff \theta}\right)^{-1}
  \frac{\diff Z(\theta)}{\diff \theta}   ~.
\end{equation}
The vector function $w(\theta)$ now has an inverse
and it is referred to as the {\em link function} 
\cite{Mccullagh.nelder.glm}.
This yields:
\begin{equation}
\label{eq-expeqs2}
\expec{x|y,\theta}{ \exp \left( \sum_{i=1}^k
                \phi_i \,t_i(x,y) \right)  } ~=~
        \frac{Z(w^{-1}(\phi+w(\theta)))}{Z(\theta)} ~.
\end{equation}
These are important because if the normalization
constant $Z$ can be found in closed form, then
it can be differentiated and divided,
for instance, symbolically,
to construct formula for
various moments of the distribution
such as 
$\expec{x|\theta}{t_i(x)}$ and $\expec{x|\theta}{t_i(x)t_j(x)}$.
Furthermore,
Equation~(\ref{eq-expeqs1}) implies that
derivatives of the normalization constant,
${\diff Z(\theta)}/{\diff \theta} $, can be found
by estimating moments of the sufficient statistics
(for instance, by Markov chain Monte Carlo methods).

\acks{The general program presented here is shared
by many, including Peter Cheeseman,
who encouraged this development from its inception.
These ideas were presented in formative stages
at Snowbird 1993 (Neural Networks for Computing),
April 1993, and
to the Bayesian Analysis in Expert Systems (BAIES) group
in Pavia, Italy, June 1993.
Feedback from that group
helped further develop these ideas.
Graduate students at Stanford
and Berkeley have also received
various incarnations of this ideas.
Thanks also to 
George John, Ronny Kohavi, Scott Schmidler, Scott Roy, 
Padhraic Smyth,
and Peter Cheeseman for their
feedback on drafts, and to the JAIR reviewers.
Brian Williams pointed out the extension of
the decomposition theorems to the deterministic case.
Brian Ripley reminded me of the extensive features of S.   }

\bibliographystyle{theapa}
\bibliography{../../handbook/Version/Bibliography/master,../../bib}

%TR \newpage
%TR \input{lwgmja.ind}

\end{document}


