# Parameter Priors for Directed Acyclic Graphical Models and the
Characterization of Several Probability Distributions

Dan Geiger
<EMAIL_ADDRESS>

David Heckerman
<EMAIL_ADDRESS>
(Annals October 2002 version with corrections and updates made June 2021)

###### Abstract

We develop simple methods for constructing parameter priors for model choice
among Directed Acyclic Graphical (DAG) models. In particular, we introduce
several assumptions that permit the construction of parameter priors for a
large number of DAG models from a small set of assessments. We then present a
method for directly computing the marginal likelihood of every DAG model given
a random sample with no missing observations. We apply this methodology to
Gaussian DAG models which consist of a recursive set of linear regression
models. We show that the only parameter prior for complete Gaussian DAG models
that satisfies our assumptions is the normal-Wishart distribution. Our
analysis is based on the following new characterization of the Wishart
distribution: let $W$ be an $n\times n$, $n\geq 3$, positive-definite
symmetric matrix of random variables and $f(W)$ be a pdf of $W$. Then, f$(W)$
is a Wishart distribution if and only if
$W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}$ is independent of
$\\{W_{12},W_{22}\\}$ for every block partitioning
$W_{11},W_{12},W^{\prime}_{12},W_{22}$ of $W$. Similar characterizations of
the normal and normal-Wishart distributions are provided as well.

Keywords: Bayesian network, Directed Acyclic Graphical Model, Dirichlet
distribution, Gaussian DAG model, learning, linear regression model, normal
distribution, Wishart distribution.

Corrections to the original text in red are taken from J. Kuipers, G. Moffa,
and D. Heckerman, Addendum on the scoring of Gaussian directed acyclic
graphical models. Annals of Statistics 42, 1689-1691, Aug 2014. Other updates
to the original are in blue.

##  1 Introduction

Directed Acyclic Graphical (DAG) models have increasing number of applications
in Statistics (Cowell, Dawid, Lauritzen, and Spiegelhalter, 1999) as well as
in Decision Analysis and Artificial Intelligence (Howard and Matheson, 1981;
Heckerman, Mamdani, and Wellman, 1995b; Pearl, 1988). A DAG model $m=(s,{\cal
F}_{s})$ for a set of variables ${\bf X}=\\{X_{1},\ldots,X_{n}\\}$ each
associated with a set of possible values $\mbox{D}_{i}$, respectively, is a
set of joint probability distributions for
$\mbox{D}_{1}\times\cdots\times\mbox{D}_{n}$ specified via two components: a
structure $s$ and a set of local distribution families ${\cal F}_{s}$. The
structure $s$ for ${\bf X}$ is a directed graph with no directed cycles (i.e.,
a Directed Acyclic Graph) having for every variable $X_{i}$ in ${\bf X}$ a
node labeled $X_{i}$ with parents labeled by ${\bf Pa}^{m}_{i}$. The structure
$s$ represents the set of conditional independence assertions, and only these
conditional independence assertions, which are implied by a factorization of a
joint distribution for ${\bf X}$ given by $p({\bf
x})=\prod_{i=1}^{n}p(x_{i}|{\bf pa}^{m}_{i})$, where ${\bf
x}=(x_{1},\ldots,x_{n})$ is a value for ${\bf X}$ (an $n$-tuple) and $x_{i}$
is a value for $X_{i}$ and where ${\bf pa}^{m}_{i}$ is the value for ${\bf
Pa}^{m}_{i}$ as in ${\bf x}$. When $x_{i}$ has no incoming arcs in $m$ (no
parents), $p(x_{i}|{\bf pa}^{m}_{i})$ stands for $p(x_{i})$. The local
distributions are the $n$ conditional and marginal probability distributions
that constitute the factorization of $p({\bf x})$. Each such distribution
belongs to the specified family of allowable probability distributions ${\cal
F}_{s}$. A DAG model is often called a Bayesian network, although the later
name sometimes refers to a specific joint probability distribution that
factorizes according to a DAG, and not, as we mean herein, a set of joint
distributions each factorizing according to the same DAG. A DAG model is
complete if it has no missing arcs. Note that any two complete DAG models for
${\bf X}$ encode the same assertions of conditional independence—namely, none.
Also note that a complete DAG determines a unique ordering of the variables in
which $X_{i}$ precedes $X_{j}$ if and only if $X_{i}\rightarrow X_{j}$ is an
arc in this DAG.

In this paper, we assume that each local distribution is selected from a
family ${\cal F}_{s}$ which depends on a finite set of parameters
$\mbox{$\theta$}_{m}\in\Theta_{m}$ (a parametric family). The parameters for a
local distribution is a set of real numbers that completely determine the
functional form of $p(x_{i}|{\bf pa}^{m}_{i})$ when $x_{i}$ has parents and of
$p(x_{i})$ when $x_{i}$ has no parents. We denote by $m^{h}$ the model
hypothesis that the true joint probability distribution of ${\bf X}$ is
perfectly represented by a structure $s$ of a DAG model $m$ with local
distributions from ${\cal F}_{s}$—namely, that the joint probability
distribution satisfies only the conditional independence assertions implied by
this factorization and none other. Consequently, the true joint distribution
for a DAG model $m$ is given by,

$p({\bf x}|\mbox{$\theta$}_{m},m^{h})=\prod_{i=1}^{n}p(x_{i}|{\bf
pa}^{m}_{i},\mbox{$\theta$}_{i},m^{h})$ (1)

where $\theta_{1},\ldots\theta_{n}$ are subsets of $\mbox{$\theta$}_{m}$.
Whereas in a general formulation of DAG models, the subsets
$\\{\theta_{i}\\}_{i=1}^{n}$ could possibly overlap allowing several local
distributions to have common parameters, in this paper, we shall shortly
exclude this possibility (Assumption 5). Note that $\theta_{m}$ denotes the
union of $\theta_{1},\ldots,\theta_{n}$ for a DAG model $m$.

We consider the Bayesian approach when the parameters $\mbox{$\theta$}_{m}$
and the model hypothesis $m^{h}$ are uncertain but the parametric families are
known. Given data $d=\\{{\bf x}_{1},\ldots,{\bf x}_{N}\\}$, a random sample
from $p({\bf x}|\mbox{$\theta$}_{m},m^{h})$ where $\mbox{$\theta$}_{m}$ and
$m^{h}$ are the true parameters and model hypothesis, respectively, we can
compute the posterior probability of a model hypothesis $m^{h}$ using

$p(m^{h}|d)=c\ p(m^{h})\ p(d|m^{h})=c\ p(m^{h})\int
p(d|\mbox{$\theta$}_{m},m^{h})\ p(\mbox{$\theta$}_{m}|m^{h})\
d\mbox{$\theta$}_{m}$ (2)

where $c$ is a normalization constant. We can then select a DAG model that has
a high posterior probability or average several good models for prediction.

The problem of selecting an appropriate DAG model, or sets of DAG models,
given data, posses a serious computational challenge, because the number of
DAG models grows faster than exponential in $n$. Methods for searching through
the space of model structures are discussed (e.g.) by Cooper and Herskovits
(1992), Heckerman, Geiger, and Chickering (1995a), and Friedman and Goldszmidt
(1997).

From a statistical viewpoint, an important question which needs to be
addressed is how to specify the quantities $p(m^{h})$,
$p(d|\mbox{$\theta$}_{m},m^{h})$, $p(\mbox{$\theta$}_{m}|m^{h})$, needed for
evaluating $p(m^{h}|d)$ for every DAG model $m$ that could conceivably be
considered by a search algorithm. Buntine (1991) and Heckerman et al. (1995a)
discuss methods for specifying the priors $p(m^{h})$ via a small number of
direct assessments.

Herein, we develop practical methods for assigning parameter priors,
$p(\mbox{$\theta$}_{m}|m^{h})$, to every candidate DAG model $m$ via a small
number of direct assessments. Our method is based on a set of assumptions, the
most notable of which is the assumption that complete DAG models represent the
same set of distributions, which implies that data cannot distinguish between
two complete DAG models. Multivariate Gaussian, multinomial, and multivariate
$t$ distributions satisfy this assumption. Another assumption is likelihood
and prior modularity, which says that the local distribution for $x_{i}$ and
its parameter priors depend only on the parents of $x_{i}$ but not on the
entire description of the structure. These assumptions, together with global
parameter independence, introduced by Spiegelhalter and Lauritzen (1990), are
the heart of the proposed methodology.

The methodology described herein for setting priors to DAG models and
consequently calculating their marginal likelihoods is an extension of the
results by Dawid and Lauritzen (1993) for decomposable graphical models. For
decomposable graphical models, which is a set of models that can be regarded
both as DAG models as well as undirected graphical models, the two
methodologies are identical. Our specification of a formal set of assumptions
followed by a technical derivation of this methodology provides an easy access
to examine the validity of the approach and devise alternatives when needed.

The contributions of this paper are as follows: A methodology for specifying
parameter priors for many DAG structures using a few direct assessments
(Section 2). A formula that computes the marginal likelihood for every DAG
model (Section 3). A specialization of this formula to an efficient
computation for Gaussian DAG models (Section 4). An analysis of complete
Gaussian DAG models which shows that the only parameter prior that satisfies
our assumptions is the normal-Wishart distribution (Section 5). The analysis
is based on the following new characterization of the Wishart, normal, and
normal-Wishart distributions.

Theorem Let $W$ be an $n\times n$, $n\geq 3$, positive-definite symmetric
matrix of real random variables such that no entry in $W$ is zero, $\mu$ be an
$n$-dimensional vector of random variables, $f_{W}(W)$ be a pdf of $W$,
$f_{\mu}(\mu)$ be a pdf of $\mu$, and $f_{\mu,W}(\mu,W)$ be a pdf of
$\\{\mu,W\\}$. Then, $f_{W}(W)$ is a Wishart distribution, $f_{\mu}(\mu)$ is a
normal distribution, and $f_{\mu,W}(\mu,W)$ is a normal-Wishart distribution
if and only if global parameter independence holds for unknown $W$, unknown
$\mu$, or unknown $\\{\mu,W\\}$, respectively.

The assumption of global parameter independence is expressed differently for
each of the three cases treated by this theorem and the proof follows from
Theorems 7, 9 and 10, respectively, proven in Section 5. It should be noted
that a single principle, global parameter independence, is used to
characterize three different distributions.

A similar characterization for the bivariate Wishart, bivariate normal, and
bivariate normal-Wishart distributions has recently been obtained under the
assumption that the pdf is strictly positive, and assuming also some
additional independence constraints—termed standard local parameter
independence (Geiger and Heckerman, 1998). Another related result is the
characterization of the Dirichlet distribution via global and local parameter
independence (Geiger and Heckerman, 1997; Járai , 1998).

##  2 The Construction of Parameter Priors

In this section, we present assumptions that simplify the assessment of
parameter priors and a method of assessing these priors. The assumptions are
as follows:

######  Assumption 1 (Complete model equivalence)

Let $m_{1}=(s_{1},{\cal F}_{s1})$ be a complete DAG model for ${\bf X}$. The
family ${\cal F}_{s2}$ of every complete DAG model $m_{2}=(s_{2},{\cal
F}_{s2})$ for ${\bf X}$ is such that $m_{1}$ and $m_{2}$ represent the same
set of joint probability distributions, namely, that for every
$\mbox{$\theta$}_{m1}$ there exists $\mbox{$\theta$}_{m2}$ such that $p({\bf
x}|\mbox{$\theta$}_{m1},m^{h}_{1})=p({\bf x}|\mbox{$\theta$}_{m2},m^{h}_{2})$
and vice versa.

Two examples where this assumption holds are quite common. One happens when
$p({\bf x}|\mbox{$\theta$}_{m1},m^{h}_{1})$ and $p({\bf
x}|\mbox{$\theta$}_{m2},m^{h}_{2})$ are multivariate normal distributions and
the other happens when ${\bf X}$ consists of variables with finite domains and
$p({\bf x}|\mbox{$\theta$}_{m1},m^{h}_{1})$ and $p({\bf
x}|\mbox{$\theta$}_{m2},m^{h}_{2})$ are unrestricted discrete distributions.
In these two cases, all the local distributions have the same functional form
in every ordering of the variables. If the joint distribution for ${\bf X}$ is
a multivariate $t$ distribution, then too, all local conditional distributions
have the same functional form (e.g., DeGroot, 1970), however, unlike the
unrestricted discrete and multivariate normal distributions, for $t$
distributions, the parameters of the local distributions are dependent which
violates assumption 5 discussed below.

We now provide an example where this assumption fails. Suppose the set of
variables ${\bf X}=\\{X_{1},X_{2},X_{3}\\}$ consists of three variables each
with possible values $\\{x_{i},\overline{x}_{i}\\}$, respectively, and $s_{1}$
is the complete structure with arcs $X_{1}\rightarrow X_{2}$,
$X_{1}\rightarrow X_{3}$, and $X_{2}\rightarrow X_{3}$. Suppose further, that
the local distributions ${\cal F}_{s1}$ of model $m_{1}$ are restricted to the
logit

$p(x_{i}|{\bf pa}^{m}_{i},\mbox{$\theta$}_{i},m^{h})=\frac{1}{1+{\rm
exp}\left\\{a_{i}+\sum_{x_{j}\in{\bf pa}^{m}_{i}}b_{ji}x_{j}\right\\}}$

where $\theta_{1}=\\{a_{1}\\}$, $\theta_{2}=\\{a_{2},b_{12}\\}$, and
$\theta_{3}=\\{a_{3},b_{13},b_{23}\\}$.

Consider now a second complete model $m_{2}$ for ${\bf
X}=\\{X_{1},X_{2},X_{3}\\}$ whose structure consists of the arcs
$X_{1}\rightarrow X_{2}$, $X_{1}\rightarrow X_{3}$, and $X_{3}\rightarrow
X_{2}$. Assumption 1 asserts that the families of local distributions for
$m_{1}$ and $m_{2}$ are such that the set of joint distributions for ${\bf X}$
represented by these two complete models is the same. In this example,
however, if we specify the local families for $m_{2}$ by also restricting them
to be logit distributions, then the two models will represent different sets
of joint distributions over $\\{X_{1},X_{2},X_{3}\\}$. Hence, Assumption 1
will be violated. Using Bayes rule one can always determine a set of local
distribution families that will satisfy Assumption 1, however, their
functional form will usually involve an integral (and will often violate
Assumption 5 below).

Note that whenever two DAG models represent the same set of probability
distributions for ${\bf X}$, they must also specify the same set of
independence assumptions. The example with the logit distributions highlights
that the converse does not hold because every complete DAG represents the same
independence assumptions, namely none, and yet complete DAG models can
represent different sets of probability distributions.

Our definition of $m^{h}$, that the true joint pdf of a set of variables ${\bf
X}$ is perfectly represented by $m$, and Assumption 1, which says that two
complete models represent the same set of joint pdfs for ${\bf X}$, imply that
for two complete models $m^{h}_{1}=m^{h}_{2}$. This is a strong assumption. It
implies that
$p(\mbox{$\theta$}_{m2}|m^{h}_{2})=p(\mbox{$\theta$}_{m2}|m^{h}_{1})$ because
two complete models represent the same set of distributions. It also implies
$p(d|m^{h}_{1})=p(d|m^{h}_{2})$ which says that the marginal likelihood for
two complete DAG models is the same for every data set, or equivalently, that
complete DAG models cannot be distinguished by data. Obviously, in the example
with the logit distributions, the two models can be distinguished by data
because they do not represent the same set of joint distributions.

######  Assumption 2 (Regularity)

For every two complete DAG models $m_{1}$ and $m_{2}$ for ${\bf X}$ there
exists a one-to-one mapping $k_{1,2}$ between the parameters
$\mbox{$\theta$}_{m1}$ of $m_{1}$ and the parameters $\mbox{$\theta$}_{m2}$ of
$m_{2}$ such that the likelihoods satisfy $p({\bf
x}|\mbox{$\theta$}_{m1},m^{h}_{1})=p({\bf x}|\mbox{$\theta$}_{m2},m^{h}_{2})$
where $\mbox{$\theta$}_{m2}=k_{1,2}(\mbox{$\theta$}_{m1}$). The Jacobian
$|\partial{\mbox{$\theta$}_{m1}}/\partial{\mbox{$\theta$}_{m2}}|$ exists and
is non-zero for all values of $\Theta_{m1}$.

Assumption 2 implies
$p(\mbox{$\theta$}_{m2}|m^{h}_{1})=\left|\frac{\partial\mbox{$\theta$}_{m1}}{\partial\mbox{$\theta$}_{m2}}\right|\
p(\mbox{$\theta$}_{m1}|m^{h}_{1})$ where
$\mbox{$\theta$}_{m2}=k_{1,2}(\mbox{$\theta$}_{m1}$). Furthermore, due to
Assumption 1,
$p(\mbox{$\theta$}_{m2}|m^{h}_{2})=p(\mbox{$\theta$}_{m2}|m^{h}_{1})$, and
thus

$p(\mbox{$\theta$}_{m2}|m^{h}_{2})=\left|\frac{\partial\mbox{$\theta$}_{m1}}{\partial\mbox{$\theta$}_{m2}}\right|\
p(\mbox{$\theta$}_{m1}|m^{h}_{1}).$ (3)

For example, suppose $x=\\{x_{1},x_{2}\\}$ have a non-singular Bivariate
normal pdf $f({\bf x})=N({\bf x}|\mu,W)$ where $\mu$ is the vector of means
and $W=(w_{ij})$ is the inverse of a positive definite covariance matrix. If
we write $f({\bf x})=f_{x_{1}}(x_{1})f_{x_{2}|x_{1}}(x_{2}|x_{1})$ where
$f_{x_{1}}(x_{1})=N(x_{1}|e_{1},1/v_{1})$ and
$f_{x_{2}|x_{1}}(x_{2}|x_{1})=N(x_{2}|e_{2|1}+b_{12}x_{1},1/v_{2|1})$, then
the following well known relationships are satisfied:

$w_{11}=\frac{1}{v_{1}}+\frac{b_{12}^{2}}{v_{2|1}}\ \ \ \
w_{12}=-\frac{b_{12}}{v_{2|1}}\ \ \ \ w_{22}=\frac{1}{v_{2|1}}\ \ \ \
e_{1}=\mu_{1}\ \ \ \ e_{2|1}=\mu_{2}-b_{12}\mu_{1}\ \ \ \ $ (4)

Note that the transformation between $\\{\mu,W\\}$ and
$\\{e_{1},v_{1},e_{2|1},v_{2|1},b_{12}\\}$ is one to one and onto as long as
$W$ is the inverse of a covariance matrix and the conditional variances
$v_{1},v_{2|1}$ are positive. The Jacobian of this transformation is given by,

$|\frac{\partial{w_{11},w_{12},w_{22},\mu_{1},\mu_{2}}}{\partial{v_{1},v_{1|2},b_{12},e_{1},e_{2|1}}}|=v_{1}^{-2}v_{2|1}^{-3}$
(5)

Symmetric equations hold when $f({\bf x})$ is written as
$f_{x_{2}}(x_{2})f_{x_{1}|x_{2}}(x_{1}|x_{2})$ and so there is a one-to-one
and onto mapping between $\\{e_{1},v_{1},e_{2|1},v_{2|1},b_{12}\\}$ and
$\\{e_{2},v_{2},e_{1|2},v_{1|2},b_{21}\\}$. Note that the parameters $\mu,W$
for the joint space are instrumental for decomposing the needed mapping into a
composition of two mappings.

######  Assumption 3 (Likelihood Modularity)

For every two DAG models $m_{1}$ and $m_{2}$ for ${\bf X}$ such that $X_{i}$
has the same parents in $m_{1}$ and $m_{2}$, the local distributions for
$x_{i}$ in both models are the same, namely, $p(x_{i}|{\bf
pa}^{m}_{i},\mbox{$\theta$}_{i},m^{h}_{1})=p(x_{i}|{\bf
pa}^{m}_{i},\mbox{$\theta$}_{i},m^{h}_{2})$ for all $X_{i}\in{\bf X}$.

######  Assumption 4 (Prior Modularity)

For every two DAG models $m_{1}$ and $m_{2}$ for ${\bf X}$ such that $X_{i}$
has the same parents in $m_{1}$ and $m_{2}$,
$p(\mbox{$\theta$}_{i}|m^{h}_{1})=p(\mbox{$\theta$}_{i}|m^{h}_{2})$.

######  Assumption 5 (Global Parameter Independence)

For every DAG model $m$ for ${\bf X}$,
$p(\mbox{$\theta$}_{m}|m^{h})=\prod_{i=1}^{n}p(\mbox{$\theta$}_{i}|m^{h})$.

The likelihood and prior modularity assumptions have been used implicitly in
the work of (e.g.) Cooper and Herskovits (1992), Spiegelhalter, Dawid,
Lauritzen, and Cowell (1993), and Buntine (1994). Heckerman et al. (1995a)
made Assumption 4 explicit in the context of discrete variables under the name
parameter modularity. Spiegelhalter and Lauritzen (1990) introduced Assumption
5 in the context of DAG models under the name global independence.

Note that the first three assumptions concern the distribution of ${\bf X}$
whereas the last two assumptions concern the distribution of the parameters.
Obviously, when the parameters $\theta_{1},\ldots\theta_{n}$ are not variation
independent for every complete DAG model for $\bf X$, the assumption of global
parameter independence is inconsistent with the model and can not be true.
Hence, assumption 5 excludes, for example, the possibility that two local
distributions share a common parameter. On the other hand, even when the
parameters are variation independent, it is possible to specify a prior
distribution for $\theta$ that violates global parameter independence. Cowell
et al. (1999, pp. 191-2) highlight this point.

The assumptions we have made lead to the following significant consequence:
When we specify a parameter prior $p(\mbox{$\theta$}_{mc}|m^{h}_{c})$ for one
complete DAG model $m_{c}$, we also implicitly specify a prior
$p(\mbox{$\theta$}_{m}|m^{h})$ for any DAG model $m$ among the super
exponentially many possible DAG models. Consequently, we have a framework in
which a manageable number of direct assessments leads to all the priors needed
to search the model space. In the rest of this section, we explicate how all
parameter priors are determined by the one elicited prior. In Section 4, we
show how to elicit the one needed prior $p(\mbox{$\theta$}_{mc}|m^{h}_{c})$
under specific distributional assumptions.

Due to the complete model equivalence and regularity assumptions, we can
compute $p(\mbox{$\theta$}_{mc}|m^{h}_{c})$ for one complete model for ${\bf
X}$ from the prior of another complete model for ${\bf X}$. In so doing, we
are merely performing coordinate transformations between parameters for
different variable orderings in the factorization of the joint likelihood (Eq.
3). Thus by specifying parameter prior for one complete model, we have
implicitly specified a prior for every complete model.

It remains to examine how the prior $p(\mbox{$\theta$}_{m}|m^{h})$ is computed
for an incomplete DAG model $m$ for ${\bf X}$ from the prior
$p(\mbox{$\theta$}_{mc}|m^{h}_{c})$ for some complete model $m_{c}$. Due to
global parameter independence we have
$p(\mbox{$\theta$}_{m}|m^{h})=\prod_{j=1}^{n}p(\mbox{$\theta$}_{j}|m^{h})$ and
therefore it suffices to examine each of the $n$ terms separately. To compute
$p(\mbox{$\theta$}_{i}|m^{h})$, we identify a complete DAG model $m_{c(i)}$
such that ${\bf Pa}^{m}_{i}={\bf Pa}^{m_{c(i)}}_{i}$. The prior
$p(\mbox{$\theta$}_{mc(i)}|m^{h}_{c(i)})$ is obtained from
$p(\mbox{$\theta$}_{mc}|m^{h}_{c})$, as we have shown for every pair of
complete DAG models. Due to global parameter independence
$p(\mbox{$\theta$}_{mc(i)}|m^{h}_{c(i)})$ is a product one term of which is
$p(\mbox{$\theta$}_{i}|m^{h}_{c(i)})$. Finally, due to prior modularity
$p(\mbox{$\theta$}_{i}|m^{h})$ is equal to
$p(\mbox{$\theta$}_{i}|m^{h}_{c(i)})$.

This methodology of constructing priors is described by Heckerman et al.
(1995a) for discrete DAG models and in Section 4 for Gaussian DAG models. Our
method is equivalent to the method of compatible priors devised for
decomposable graphical models (Dawid and Lauritzen, 1993). Our arguments, via
a set of assumptions, can be regarded as an axiomatic justification for
compatible priors, and as an extension of this method to general DAG models
and to any probability distributions that satisfy Assumptions 1 – 5. We are
currently unaware, however, of additional probability distributions that
satisfy these five assumptions.

The following theorem summarizes the general construction which was formulated
to cover both cases—the discrete and the Gaussian.

###### Theorem 1

Given Assumptions 1 through 5, the parameter prior
$p(\mbox{$\theta$}_{m}|m^{h})$ for every DAG model $m$ is determined by a
specified parameter prior $p(\mbox{$\theta$}_{mc}|m^{h}_{c})$ for an arbitrary
complete DAG model $m_{c}$.

Theorem 1 shows that once we specify the parameter prior for one complete DAG
model all other priors can be generated automatically and need not be
specified manually. Consequently, together with Eq. 2 and due to the fact that
also likelihoods can be generated automatically in a similar fashion, we have
a manageable methodology to automate the computation of $p(d|m^{h})$ for any
DAG model of ${\bf X}$ which is being considered by a search algorithm as a
candidate model. Next we show how this computation can be done efficiently.

##  3 Computation of the Marginal Likelihood for Complete Data

For a given ${\bf X}$, consider a DAG model $m$ and a complete random sample
$d$. Assuming global parameter independence, the parameters remain independent
given complete data. That is,

$p(\mbox{$\theta$}_{m}|d,m^{h})=\prod_{i=1}^{n}p(\mbox{$\theta$}_{i}|d,m^{h})$
(6)

In addition, assuming global parameter independence, likelihood modularity,
and prior modularity, the parameters remain modular given complete data. In
particular, if $X_{i}$ has the same parents in $s_{1}$ and $s_{2}$, then

$p(\mbox{$\theta$}_{i}|d,m^{h}_{1})=p(\mbox{$\theta$}_{i}|d,m^{h}_{2})$ (7)

Also, for any ${\bf Y}\subseteq{\bf X}$, define $d^{{\bf Y}}$ to be the random
sample $d$ restricted to observations of ${\bf Y}$. For example, if ${\bf
X}=\\{X_{1},X_{2},X_{3}\\}$, ${\bf Y}=\\{X_{1},X_{2}\\}$, and $d=\\{{\bf
x}_{1}=\\{x_{11},x_{21},x_{31}\\},{\bf x}_{2}=\\{x_{12},x_{22},x_{32}\\}\\}$,
then we have $d^{{\bf Y}}=\\{\\{x_{11},x_{21}\\},\\{x_{12},x_{22}\\}\\}$. Let
${\bf Y}$ be a subset of ${\bf X}$, and $s_{c}$ be a complete structure for
any ordering where the variables in ${\bf Y}$ come first. Then, assuming
global parameter independence and likelihood modularity, it is not difficult
to show that

$p({\bf y}|d,m^{h}_{c})=p({\bf y}|d^{{\bf Y}},m^{h}_{c})$ (8)

Given these observations, we can compute the marginal likelihood as follows,
yielding an important component for searching DAG models via a Bayesian
methodology.

###### Theorem 2

Given any complete DAG model $m_{c}$ for ${\bf X}$, any DAG model $m$ for
${\bf X}$, and any complete random sample $d$, Assumptions 1 through 5 imply

$p(d|m^{h})=\prod_{i=1}^{n}\frac{p(d^{{\bf
Pa}_{i}\cup\\{X_{i}\\}}|m^{h}_{c})}{p(d^{{\bf Pa}_{i}}|m^{h}_{c})}$ (9)

Proof: From the rules of probability, we have

$p(d|m^{h})=\prod_{l=1}^{{\color[rgb]{0,0,1}N}}\int p({\bf
x}_{l}|\mbox{$\theta$}_{m},m^{h})\ p(\mbox{$\theta$}_{m}|d_{l},m^{h})\
d\mbox{$\theta$}_{m}$ (10)

where $d_{l}=\\{{\bf x}_{1},\ldots,{\bf x}_{l-1}\\}$. Using Equations 1 and 6
to rewrite the first and second terms in the integral, respectively, we obtain

$p(d|m^{h})=\prod_{l=1}^{{\color[rgb]{0,0,1}N}}\int\prod_{i=1}^{n}p(x_{il}|{\bf
pa}_{il},\mbox{$\theta$}_{i},m^{h})\ p(\mbox{$\theta$}_{i}|d_{l},m^{h})\
d\mbox{$\theta$}_{m}$

where $x_{il}$ is the value of $X_{i}$ in the $l$-th data point.

Using likelihood modularity and Equation 7, we get

$p(d|m^{h})=\prod_{l=1}^{{\color[rgb]{0,0,1}N}}\int\prod_{i=1}^{n}p(x_{il}|{\bf
pa}_{il},\mbox{$\theta$}_{i},m^{h}_{c(i)})\
p(\mbox{$\theta$}_{i}|d_{l},m^{h}_{c(i)})\ d\mbox{$\theta$}_{m}$ (11)

where $s_{c(i)}$ is a complete structure with variable ordering ${\bf
Pa}_{i}$, $X_{i}$ followed by the remaining variables. Decomposing the
integral over $\mbox{$\theta$}_{m}$ into integrals over the individual
parameter sets $\mbox{$\theta$}_{i}$, and performing the integrations, we have

$p(d|m^{h})=\prod_{l=1}^{{\color[rgb]{0,0,1}N}}\prod_{i=1}^{n}p(x_{il}|{\bf
pa}_{il},d_{l},m^{h}_{c(i)})$

Using Equation 8, we obtain

$\displaystyle p(d|m^{h})$ $\displaystyle=$
$\displaystyle\prod_{l=1}^{{\color[rgb]{0,0,1}N}}\prod_{i=1}^{n}\frac{p(x_{il},{\bf
pa}_{il}|d_{l},m^{h}_{c(i)})}{p({\bf pa}_{il}|d_{l},m^{h}_{c(i)})}$ (12)
$\displaystyle=$
$\displaystyle\prod_{l=1}^{{\color[rgb]{0,0,1}N}}\prod_{i=1}^{n}\frac{p(x_{il},{\bf
pa}_{il}|d^{{\bf Pa}_{i}\cup\\{X_{i}\\}}_{l},m^{h}_{c(i)})}{p({\bf
pa}_{il}|d^{{\bf Pa}_{i}}_{l},m^{h}_{c(i)})}$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{n}\frac{p(d^{{\bf
Pa}_{i}\cup\\{X_{i}\\}}|m^{h}_{c(i)})}{p(d^{{\bf Pa}_{i}}|m^{h}_{c(i)})}$

By the likelihood modularity, complete model equivalence, and regularity
assumptions, we have that $p(d|m^{h}_{c(i)})=p(d|m^{h}_{c}),\ i=1,\ldots,n$.
Consequently, for any subset ${\bf Y}$ of ${\bf X}$, we obtain $p(d^{{\bf
Y}}|m^{h}_{c(i)})=p(d^{{\bf Y}}|m^{h}_{c})$ by summing over the variables in
${\bf X}\setminus{\bf Y}$. Consequently, using Equation 12, we get Equation 9.
$\Box$

An equivalent approach for computing the marginal likelihood (Equation 9) for
decomposable discrete and Gaussian DAG models has been developed by Dawid and
Lauritzen (1993) using compatible priors.

An important feature of Equation 9, which we now demonstrate, is that two DAG
models that represent the same assertions of conditional independence have the
same marginal likelihood. We say that two structures for ${\bf X}$ are
independence equivalent if they represent the same assertions of conditional
independence. Independence equivalence is an equivalence relation, and induces
a set of equivalence classes over the possible structures for ${\bf X}$.

Verma and Pearl (1990) provide a simple characterization of independence
equivalent structures using the concept of a v-structure. Given a structure
$s$, a v-structure in $s$ is an ordered node triple $(X_{i},X_{j},X_{k})$
where $s$ contains the arcs $X_{i}\rightarrow X_{j}$ and $X_{j}\leftarrow
X_{k}$, and there is no arc between $X_{i}$ and $X_{k}$ in either direction.
Verma and Pearl show that two structures for ${\bf X}$ are independence
equivalent if and only if they have identical edges and identical
v-structures. This characterization makes it easy to identify independence
equivalent structures.

An alternative characterization developed by Chickering (1995) and
independently by Andersson, Madigan, and Perlman (1997, Lemma 3.2), is useful
for proving our claim that independence equivalent structures have the same
marginal likelihood. An arc reversal is a transformation from one structure to
another, in which a single arc between two nodes is reversed. An arc between
two nodes is said to be covered if those two nodes would have the same parents
if the arc were removed.

######  Theorem 3 (Chickering, 1995; Andersson, Madigan, and Perlman, 1997)

Two structures for ${\bf X}$ are independence equivalent if and only if there
exists a set of covered arc reversals that transform one structure into the
other.

A proof of this theorem can also be found in Heckerman et al. (1995a).

Theorem 3 implies that if every pair of DAGs that differ by a single covered
arc represent the same set of distributions, then every two independence
equivalent DAGs represent the same set of distributions. Furthermore, a
consequence of the next theorem is that Assumptions 1 through 5 imply that
indeed every two independence equivalent DAGs represent the same set of
distributions. Without these assumptions, two independence equivalent DAGs can
represent different sets of distributions.

###### Theorem 4

Given Assumptions 1 through 5, every two independence equivalent DAG models
have the same marginal likelihood.

Proof: Theorem 3 implies that we can restrict the proof to two DAG models that
differ by a single covered arc. Say the arc is between $X_{i}$ and $X_{j}$ and
that the joint parents of $X_{i}$ and $X_{j}$ are denoted by $\pi$. For these
two models, Equation 9 differs only in terms $i$ and $j$. For both models the
product of these terms is
$p(d^{\pi\cup\\{X_{i},X_{j}\\}}|m^{h}_{c})/p(d^{\pi}|m^{h}_{c})$. $\Box$

The conclusions of Theorem 2, and, consequently, of Theorem 4 are not
justified when our assumptions are violated. In the example of the logit
distributions, discussed in the previous subsection, which violates assumption
1, the structures $s_{1}$ and $s_{2}$ differ by the reversal of a covered arc
between $X_{2}$ and $X_{3}$, but, given that all local distribution families
are logit, there are certain joint distributions that can be represented by
one structure, but not the other, and so their marginal likelihood will be
different.

The implication of Theorem 4 is quite strong: all models in the same
independence equivalence class are scored equivalently. This severely
constrains possible parameter priors as shown in the next two sections. One
possible approach to bypass our assumptions is to select one representative
DAG model from each class of independence equivalent DAG models, assume global
parameter independence only for these representatives, and evaluate the
marginal likelihood only for these representatives. The search can then be
conducted in the space of representative models as suggested in Spirtes and
Meek (1995), Chickering (1996), and Madigan, Andersson, Perlman, and Volinsky
(1996). The difficulty with this approach is that when projecting a prior from
a complete DAG model to a DAG model with missing edges, one needs to perform
additional high dimensional integrations before using the parameter modularity
property (see Section 2). Another approach is to modify the definition of
$m^{h}$ to allow independence equivalent DAG models to have different
parameter priors. This alternative is needed when arcs have a causal
interpretation. However, when choosing this alternative, the parameter prior
for each model examined by a search procedure must be provided by a user as
the search is being conducted, or a new mechanism to produce acceptable priors
on-the-fly must be devised.

##  4 Gaussian Directed Acyclic Graphical Models

We now apply the methodology of previous sections to Gaussian DAG models. A
Gaussian DAG model is a DAG model as defined by Eq 1, where each variable
$X_{i}\in{\bf X}$ is continuous, and each local likelihood is the linear
regression model

$p(x_{i}|{\bf
pa}^{m}_{i},\mbox{$\theta$}_{i},m^{h})=N(x_{i}|m_{i}+\sum_{x_{j}\in{\bf
pa}_{i}}b_{ji}x_{j},1/v_{i})$ (13)

where $N(x_{i}|\mu,\tau)$ is a normal distribution with mean $\mu$ and
precision $\tau>0$. Given this form, a missing arc from $X_{j}$ to $X_{i}$ is
equivalent to $b_{ji}=0$ in the DAG model. The local parameters are given by
$\mbox{$\theta$}_{i}=(m_{i},{b}_{i},v_{i})$, where ${b}_{i}$ is the column
vector $(b_{1i},\ldots,b_{i-1,i})$ of regression coefficients. Furthermore,
$m_{i}$ is the conditional mean of $X_{i}$ and $v_{i}$ is the conditional
variance of $X_{i}$.

For Gaussian DAG models, the joint likelihood $p({\bf
x}|\mbox{$\theta$}_{m},m^{h})$ obtained from Eqs 1 and 13 is an
$n$-dimensional multivariate normal distribution with a mean vector $\mu$ and
a symmetric positive definite precision matrix $W$,

$p({\bf x}|\mbox{$\theta$}_{m},m^{h})=\prod_{i=1}^{n}p(x_{i}|{\bf
pa}^{m}_{i},\mbox{$\theta$}_{i},m^{h})=N({\bf x}|\mu,W).$

For a complete model $m_{c}$ with ordering $(X_{1},\ldots,X_{n})$ there is a
one-to-one mapping between
$\mbox{$\theta$}_{mc}=\bigcup_{i=1}^{n}\mbox{$\theta$}_{i}$ where
$\mbox{$\theta$}_{i}=(m_{i},{b}_{i},v_{i})$ and $\\{\mu,W\\}$ which has a
nowhere singular Jacobian matrix. Consequently, assigning a prior for the
parameters of one complete model induces a parameter prior, via the change of
variables formula, for $\\{\mu,W\\}$ and in turn, induces a parameter prior
for every complete model. Any such induced parameter prior must satisfy,
according to our assumptions, global parameter independence. Not many prior
distributions satisfy such a requirement. In fact, in the next section we show
that the parameter prior $p(\mu,W|m^{h}_{c})$ must be a normal-Wishart
distribution.

For now we proceed by simply choosing $p(\mu,W|m^{h}_{c})$ to be a normal-
Wishart distribution. In particular, $p(\mu|W,m^{h}_{c})$ is a multivariate
normal distribution with mean $\nu$ and precision matrix $\alpha_{\mu}W$
($\alpha_{\mu}>0$); and $p(W|m^{h}_{c})$ is a Wishart distribution, given by,

$p(W|m^{h}_{c})=c(n,\alpha_{w})|T|^{\alpha_{w}/2}|W|^{(\alpha_{w}-n-1)/2}e^{-1/2\mbox{\it
tr}\\{TW\\}}\equiv\mbox{\cal Wishart}(W|\alpha_{w},T)$ (14)

with $\alpha_{w}$ degrees of freedom $(\alpha_{w}>n-1)$ and a positive-
definite precision matrix $T$ and where $c(n,\alpha_{w})$ is a normalization
constant given by

$c(n,\alpha_{w})=\left[2^{\alpha_{w}n/2}\pi^{n(n-1)/4}\prod_{i=1}^{n}\Gamma\left(\frac{\alpha_{w}+1-i}{2}\right)\right]^{-1}$
(15)

(DeGroot, 1970, pp. 57). We provide interpretations for $\nu$, $\alpha_{\mu}$,
$T$, and $\alpha_{w}$ later in this section. Note that in some expositions of
the Wishart distribution, the inverse of $T$ is used for the parameterization;
$T^{-1}$ is called the scale matrix (e.g., Press 1971, pp. 101).

This choice of a prior satisfies global parameter independence due to the
following well known theorem.

Let $M^{\prime}$ denote the transpose of matrix $M$. Define a block
partitioning $\\{W_{11},W_{12},W^{\prime}_{12},W_{22}\\}$ of an $n$ by $n$
matrix $W$ to be compatible with a partitioning $\mu_{1},\mu_{2}$ of an $n$
dimensional vector $\mu$, if the indices of the rows that correspond to block
$W_{11}$ are the same as the indices of the terms that constitute $\mu_{1}$
and similarly for $W_{22}$ and $\mu_{2}$. Also define
$W_{11.2}=((W^{-1})_{11})^{-1}=W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}$.

###### Theorem 5

If $f(\mu,W)$ is an $n$ dimensional normal-Wishart distribution, $n\geq 2$,
with parameters $\nu,\alpha_{\mu}$, $\alpha_{w}$, and $T$, then
$\\{\mu_{1},W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}\\}$ is independent of
$\\{\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1},$ $W_{12},W_{22}\\}$ for every
partitioning $\mu_{1},\mu_{2}$ of $\mu$ where $W_{11}$,$W_{12}$,
$W^{\prime}_{12}$, $W_{22}$ is a block partitioning of $W$ compatible with the
partitioning $\mu_{1},\mu_{2}$. Furthermore, the pdf of
$\\{\mu_{1},W_{11.2}\\}$ is normal-Wishart with parameters $\nu_{1}$,
$\alpha_{\mu}$, $T_{11}$, and $\alpha_{w}-n+l$ where $T_{11}$,$T_{12}$,
$T^{\prime}_{12}$, $T_{22}$ is a compatible block partitioning of $T$,
$\nu_{1},\nu_{2}$ is a compatible partitioning of $\nu$, and $l$ is the size
of the vector $\nu_{1}$.

The proof of Theorem 5 requires a change of variables from $(\mu,W)$ to
$(\mu_{1}$, $\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1})$ and
$(W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12},$ $W_{12},W_{22})$. Press (1971, p.
117-119) carries out these computations for the Wishart distribution. Standard
changes are needed to obtain the claim for the normal-Wishart distribution. A
consequence of Theorem 5 is the following.

###### Corollary 6

Let $W$ be a $n\times n$ positive-definite matrix of random variables. Let
$a$, $b$, and $c$ be three sets of indices of $W$. If $f(W_{ab.c})=\mbox{\cal
Wishart}(W_{ab.c}|\alpha_{1},T_{1})$ and $f(W_{bc.a})=\mbox{\cal
Wishart}(W_{bc.a}|\alpha_{2},T_{2})$, then
$\alpha_{1}-l_{ab}=\alpha_{2}-l_{bc}$ where $l_{ab}$ is the number of indices
in the block $a,b$ and $l_{bc}$ is the number of indices in the block $b,c$.

Proof: The pdf for $W_{b.ac}=(W_{ab.c})_{b.a}=(W_{cb.a})_{b.c}$ is a Wishart
distribution, and from the two alternative ways by which this pdf can be
formed, using Theorem 5, it follows that
$\alpha_{1}-l_{ab}=\alpha_{2}-l_{bc}$. $\Box$

To see why the independence conditions in Theorem 5 imply global parameter
independence, consider the partitioning in which the first block contains the
first $n-1$ coordinates which correspond to $X_{1},\ldots,X_{n-1}$ while the
second block contains the last coordinate which corresponds to $X_{n}$. For
this partitioning, $b_{n}=-W_{22}^{-1}W^{\prime}_{12}$, $v_{n}=W_{22}^{-1}$,
and $m_{n}=\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1}$. Furthermore,
$((W^{-1})_{11})^{-1}=W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}=W_{11.2}$ is the
precision matrix associated with $X_{1},\ldots,X_{n-1}$. Consequently,
$\\{m_{n},b_{n},v_{n}\\}$ is independent of $\\{\mu_{1},W_{11.2}\\}$. We now
recursively repeat this argument with $\\{\mu_{1},W_{11.2}\\}$ instead of
$\\{\mu,W\\}$, to obtain global parameter independence. The converse, namely
that global parameter independence implies the independence conditions in
Theorem 5, is established similarly.

Our choice of prior implies that the posterior $p(\mu,W|d,m^{h}_{c})$ is also
a normal-Wishart distribution (DeGroot, 1970, p. 178). In particular,
$p(\mu|W,d,m^{h}_{c})$, where $d$ is a sample of $N$ complete cases, is
multivariate normal with mean vector given by

$\frac{\alpha_{\mu}\nu+N\,\overline{{\bf x}}_{N}}{\alpha_{\mu}+N}$ (16)

and precision matrix $(\alpha_{\mu}+N)W$, where $\overline{{\bf x}}_{N}$ is
the sample mean of $d$, and $p(W|d,m^{h}_{c})$ is a Wishart distribution with
$\alpha_{w}+N$ degrees of freedom and precision matrix $R$ given by

$R=T+S_{N}+\frac{\alpha_{\mu}N}{\alpha_{\mu}+N}(\nu-\overline{{\bf
x}}_{N})(\nu-\overline{{\bf x}}_{N})^{\prime}$ (17)

where $S_{N}=\sum_{l=1}^{N}({\bf x}_{l}-\overline{{\bf x}}_{N})({\bf
x}_{l}-\overline{{\bf x}}_{N})^{\prime}$. From these equations, we see that
$\alpha_{\mu}$ and $\alpha_{w}$ can be thought of as effective sample sizes
for the normal and Wishart components of the prior, respectively.

In order to calculate the marginal likelihood of a Gaussian DAG model, we can
work in the parametric space $(\mu,W)$. According to Theorem 5, if
$p(\mu,W|m^{h}_{c})$ is a normal-Wishart distribution with the parameters
given by the theorem, then $p(\mu_{{\bf Y}},((W^{-1})_{{\bf Y}{\bf
Y}})^{-1}|m^{h}_{c})$ is also a normal–Wishart distribution with effective
sample sizes $\alpha_{\mu}$ and $\alpha_{w}-n+l$, and parameters $\nu_{{\bf
Y}}$ and $T_{{\bf Y}{\bf Y}}$, where ${\bf Y}$ is a subset of $l$ coordinates
and $M_{{\bf Y}{\bf Y}}$ is the matrix $M$ with elements restricted to the
corresponding variables in ${\bf Y}$. Thus, we obtain the terms in Equation 9:

$p(d^{{\bf
Y}}|m^{h}_{c})=(2\pi)^{-lN/2}\left(\frac{\alpha_{\mu}}{\alpha_{\mu}+N}\right)^{l/2}\frac{c(l,\alpha_{w}-n+l)}{c(l,\alpha_{w}-n+l+N)}\
|{\color[rgb]{1,0,0}T_{{\bf Y}{\bf Y}}}|^{\frac{\alpha_{w}-n+l}{2}}\
|{\color[rgb]{1,0,0}R_{{\bf Y}{\bf Y}}}|^{-\frac{\alpha_{w}-n+l+N}{2}}$

(See Geiger and Heckerman, 1994, for a derivation when $l=n$.) Equations 9 and
4 together provide a way to compute the marginal likelihood for Gaussian DAG
models given the assessment of the parameter prior $p(\mu,W|m^{h}_{c})$.

Rather than a direct assessment of the parmeter prior, we consider a partially
indirect approach. We start with the observation that when
$p(\mu,W|m^{h}_{c})$ is normal–Wishart as we have described, then then $p({\bf
x}|m^{h}_{c})$ is a multivariate $t$ distribution with $\alpha_{w}-n+1$
degrees of freedom, location vector $\nu$, and precision matrix
$\alpha_{\mu}(\alpha_{w}-n+1)/(\alpha_{\mu}+1)T^{-1}$. This result can be
derived by first integrating over $\mu$ using Equation 6 on p. 178 of DeGroot
with sample size equal to one, and then integrating over $W$ following an
approach similar to that on pp. 179–180. Next, when $\alpha_{w}>n+1$, it
follows that

${\rm E}({\bf x}|m^{h}_{c})=\nu\ \ \ \ \ \ {\rm Cov}({\bf
x}|m^{h}_{c})=\frac{\alpha_{\mu}+1}{\alpha_{\mu}}\ \frac{1}{\alpha_{w}-n-1}\
T$ (18)

(e.g., DeGroot, 1970, pp. 61). Thus, a person can assess the parameter prior
by assessing $\alpha_{\mu}$ and $\alpha_{w}$, driectly, and by assessing a DAG
model for E$({\bf x}|m^{h}_{c})$ and Cov$({\bf x}|m^{h}_{c})$ and then
computing $\nu$ and $T$ using Equations 18. We call this model a prior DAG
model. The unusual aspect of this assessment is the conditioning hypothesis
$m^{h}_{c}$ (see Heckerman et al. [1995b] for a discussion). This indirect
approach provides a suitable Bayesian alternative for many of the examples
discussed in Spirtes, Glymour, and Scheines (2001).

##  5 Characterization of Several Probability Distributions

We now characterize the Wishart distribution as the only pdf that satisfies
global parameter independence for an unknown precision matrix $W$ with $n\geq
3$ coordinates (Theorem 7). This theorem is phrased and proven in a
terminology that relates to known facts about the Wishart distribution. We
proceed with similar characterizations of the normal and normal-Wishart
distributions (Theorems 9 and 10).

We will use $\mbox{\it tr}\\{A+B\\}$ to denote the sum of traces $\mbox{\it
tr}\\{A\\}+\mbox{\it tr}\\{B\\}$ even when the dimensions of the square
matrices $A$ and $B$ are different.

###### Theorem 7

Let $W$ be an $n\times n$, $n\geq 3$, positive-definite symmetric matrix of
random variables and $f(W)$ be a pdf of $W$. Then, f$(W)$ is a Wishart
distribution if and only if $W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}$ is
independent of $\\{W_{12},W_{22}\\}$ for every block partitioning
$W_{11},W_{12},W^{\prime}_{12},W_{22}$ of $W$.

Proof: That $W_{11.2}=W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}$ is independent
of $\\{W_{12},W_{22}\\}$ whenever $f(W)$ is a Wishart distribution is a well
known fact (Press 1971, p. 117-119). It is also expressed by Theorem 5. The
other direction is proven by induction on $n$. The base case $n=3$ is treated
at the end.

The pdf of $W$ can be written in $n!$ orderings. In particular, due to the
assumed independence conditions, and since the transformations from
$\\{W_{11},W_{12},W_{22}\\}$ to $\\{W_{11.2},W_{12},W_{22}\\}$ and to
$\\{W_{22.1},W_{11},W_{12}\\}$ both have a Jacobian determinant of 1, we
obtain the following functional equation:

$f(W)=f_{1}(W_{11}-W_{12}W_{22}^{-1}W_{12}^{\prime})f_{2|\\!|1}(W_{22},W_{12})=f_{2}(W_{22}-W_{12}^{\prime}W_{11}^{-1}W_{12})f_{1|\\!|2}(W_{11},W_{12})$
(19)

where a subscripted $f$ denotes a pdf.

We divide the indices of $W$ into two blocks, the first block (say, block 1)
contains $n-1$ indices and the second block (say, block 2) consists of one
index. By the induction hypothesis, and since the independence conditions on
$W$ also hold for $W_{11.2}$ we conclude that $W_{11.2}$ is distributed
according to $\mbox{\cal Wishart}(W_{11.2}|\;\alpha_{1},T_{1})$. Since this
argument holds for every block of size $n-1$ of $W$, and since if a matrix $V$
is distributed Wishart so does $V_{11.2}$ for any block of indices (Theorem
5), it follows that $W_{11.2}$ is distributed according to $\mbox{\cal
Wishart}(W_{11.2}|\;\alpha_{1},T_{1})$ also for blocks of size smaller than
$n-1$.

Thus,

$\displaystyle
c_{1}|W_{11}-W_{12}W_{22}^{-1}W_{12}^{\prime}|^{\beta_{1}}e^{\mbox{\it
tr}\\{T_{1}(W_{11}-W_{12}W_{22}^{-1}W_{12}^{\prime})\\}}f_{2|\\!|1}(W_{22},W_{12})=$
$\displaystyle
c_{2}|W_{22}-W_{12}^{\prime}W_{11}^{-1}W_{12}|^{\beta_{2}}e^{\mbox{\it
tr}\\{T_{2}(W_{22}-W_{12}^{\prime}W_{11}^{-1}W_{12})\\}}f_{1|\\!|2}(W_{11},W_{12})$
(20)

where $c_{1}$ and $c_{2}$ are normalizing constants.

We now argue that $\beta_{1}=\beta_{2}$. Divide the indices of $W$ to three
non empty sets $a,b,c$ such that block 1 consists of the indices in $a,b$ and
block 2 consists of the indices in $c$. The matrices $W_{ab.c}$ and $W_{bc.a}$
have a Wishart distribution, with, say, degrees of freedom $\alpha_{1}$ and
$\alpha_{2}$ respectively, and so according to Corollary 6,
$\alpha_{1}-l_{ab}=\alpha_{2}-l_{bc}$. Furthermore,
$W_{c.ab}=(W_{bc.a})_{c.b}$ has a Wishart distribution with
$\alpha_{2}-l_{bc}+l_{c}$ degrees of freedom. Consequently,
$\beta_{1}=(\alpha_{1}-l_{ab}-1)/2$ is equal to
$\beta_{2}=(\alpha_{2}-l_{bc}+l_{c}-l_{c}-1)/2$. Let
$\beta=\beta_{1}=\beta_{2}$.

Define

$\displaystyle F_{2|\\!|1}(W_{22},W_{12})$ $\displaystyle=$ $\displaystyle
c_{1}f_{2|\\!|1}(W_{22},W_{12})/|W_{22}|^{\beta}e^{\mbox{\it
tr}\\{T_{2}W_{22}+T_{1}(W_{12}W_{22}^{-1}W_{12}^{\prime})\\}}$ (21)
$\displaystyle F_{1|\\!|2}(W_{11},W_{12})$ $\displaystyle=$ $\displaystyle
c_{2}f_{1|\\!|2}(W_{11},W_{12})/|W_{11}|^{\beta}e^{\mbox{\it
tr}\\{T_{1}W_{11}+T_{2}(W_{12}^{\prime}W_{11}^{-1}W_{12})\\}},$ (22)

substitute into Equation 5, and obtain, using
$|W_{11}-W_{12}W_{22}^{-1}W_{12}^{\prime}||W_{22}|=|W|$, that
$F_{2|\\!|1}(W_{22},W_{12})=F_{1|\\!|2}(W_{11},W_{12})$. Consequently,
$F_{2|\\!|1}$ and $F_{1|\\!|2}$ are functions only of $W_{12}$ and thus, using
Equation 19, we obtain

$f(W)=|W|^{\beta}e^{\mbox{\it tr}\\{T_{1}W_{11}+T_{2}W_{22}\\}}H(W_{12})$ (23)

for some function $H$.

To show that $f(W)$ is Wishart we must find the form of $H$ and show that it
is proportional to $e^{2\mbox{\it tr}\\{T_{12}W_{12}\\}}$ for some matrix
$T_{12}$.

Considering the three possible pairs of blocks formed with the sets of indices
$a$, $b$, and $c$, Equation 23 can be rewritten as follows.

$\displaystyle f(W)=|W|^{\beta_{1}}e^{\mbox{\it
tr}\\{T_{aa}W_{aa}+T_{bb}W_{bb}+T_{cc}W_{cc}\\}}e^{2\mbox{\it
tr}\\{T^{\prime}_{ab}W_{ab}+T^{\prime}_{ac}W_{ac}+T^{\prime}_{bc}W_{bc}\\}}H_{1}(W_{ac},W_{bc})$
(24) $\displaystyle f(W)=|W|^{\beta_{2}}e^{\mbox{\it
tr}\\{S_{aa}W_{aa}+S_{bb}W_{bb}+S_{cc}W_{cc}\\}}e^{2\mbox{\it
tr}\\{S^{\prime}_{ab}W_{ab}+S^{\prime}_{ac}W_{ac}+S^{\prime}_{bc}W_{bc}\\}}H_{2}(W_{ab},W_{bc})$
(25) $\displaystyle f(W)=|W|^{\beta_{3}}e^{\mbox{\it
tr}\\{R_{aa}W_{aa}+R_{bb}W_{bb}+R_{cc}W_{cc}\\}}e^{2\mbox{\it
tr}\\{R^{\prime}_{ab}W_{ab}+R^{\prime}_{ac}W_{ac}+R^{\prime}_{bc}W_{bc}\\}}H_{3}(W_{ab},W_{ac})$
(26)

By setting $W_{ab}=W_{ac}=W_{bc}=0$, we get $\beta_{1}=\beta_{2}=\beta_{3}$
and $T_{ii}=S_{ii}=R_{ii}$, for $i=a,b,c$. By comparing Equations 24 and 25 we
obtain

$e^{2\mbox{\it
tr}\\{(T^{\prime}_{ac}-S^{\prime}_{ac})W_{ac}\\}}H_{1}(W_{ac},W_{bc})=e^{2\mbox{\it
tr}\\{(S^{\prime}_{ab}-T^{\prime}_{ab})W_{ab}+(S^{\prime}_{bc}-T^{\prime}_{bc})W_{bc}\\}}H_{2}(W_{ab},W_{bc})$
(27)

Each side of this equation must be a function only of $W_{bc}$. We denote this
function by $H_{12}$. Hence,

$H_{1}(W_{ac},W_{bc})=H_{12}(W_{bc})e^{2\mbox{\it
tr}\\{(S^{\prime}_{ac}-T^{\prime}_{ac})W_{ac}\\}}$

and by symmetric arguments, comparing Equations 24 and 26,

$H_{1}(W_{ac},W_{bc})=H_{13}(W_{ac})e^{2\mbox{\it
tr}\\{(R^{\prime}_{bc}-T^{\prime}_{bc})W_{bc}\\}}$

Consequently, $H_{12}(W_{bc})$ is proportional to $e^{2\mbox{\it
tr}\\{(R^{\prime}_{bc}-T^{\prime}_{bc})W_{bc}\\}}$ and so, substituting into
Equation 23, $f(W)$ is found to be a Wishart distribution, as claimed.

It remains to examine the case $n=3$. We first assume $n=2$ in which case
$f(W)$ is not necessarily a Wishart distribution. In the appendix we show that
given the independence conditions for two coordinates, $f$ must have the form

$f(W)=c|W|^{\beta}e^{\mbox{\it tr}\\{TW\\}}H(W_{12})$ (28)

where $H$ is an arbitrary function, and that the marginal distributions of
$W_{11.2}$ and $W_{22.1}$ are one dimensional Wishart distributions.

We now treat the case $n=3$ using these assertions about the case $n=2$.
Starting with Equation 19, and proceeding with blocks $a,b,c$ each containing
exactly one coordinate, we get, due to the given independence conditions for
two coordinates, that $f_{1}$ has the form given by Equation 28, and that
$f_{2}$ is a one dimensional Wishart distribution. Proceeding parallel to
Equations 5 through 22, we obtain,

$H(a_{12}-b_{1}b_{2}/W_{22})F_{2|\\!|1}(W_{22},W_{12})=F_{1|\\!|2}(W_{11},W_{12})$
(29)

where $(b_{1},b_{2})$ is the matrix $W_{12}$, $a_{12}$ is the off-diagonal
element of $W_{11}$, $a_{12}-b_{1}b_{2}/W_{22}$ is the off diagonal element of
$W_{11}-W_{12}W^{-1}_{22}W^{\prime}_{12}$, and $W_{22}$ is a $1\times 1$
matrix. Note that the left hand side depends on $W_{11}$ only through
$a_{12}$. Thus also the right hand side depends on $W_{11}$ only through
$a_{12}$. Let $b_{1}$ and $b_{2}$ be fixed, $y=b_{1}b_{2}/W_{22}$, and
$x=a_{12}$. Also let $F(t)=F_{2|\\!|1}(b_{1}b_{2}/t,(b_{1},b_{2}))$ and
$G(a_{12})=F_{1|\\!|2}(W_{11},(b_{1},b_{2}))$. We can now rewrite Equation 29
as $H(x-y)F(y)=G(x)$. Now set $z=x-y$, and obtain for every $y$ and $z$

$H(z)F(y)=G(y+z)$ (30)

the only measurable solution of which for $H$ is $H(z)=ce^{bz}$ (e.g., Aczél ,
1966).

Substituting this form of $H$ into Equation 28, we see that $W_{11.2}$ has a
two dimensional Wishart distribution. Recall that $W_{22.1}$ has a one
dimensional Wishart distribution. Consequently, we can apply the induction
step starting form Equation 5 and prove the Theorem for $n=3$. $\Box$

We now treat the situation when only the means are unknown, characterizing the
normal distribution. The two-dimensional case turns out to be covered by the
Skitovich-Darmois theorem (e.g., Kagan, Linnik, and Rao, 1973).

######  Theorem 8 (Skitovich-Darmois)

Let $z_{1},\ldots,z_{k}$ be independent random variables and
$\alpha_{i},\beta_{i}$, $1<i<k$, be constant coefficients. If
$L_{1}=\sum\alpha_{i}z_{i}$ is independent of $L_{2}=\sum\beta_{i}z_{i}$, then
each $z_{i}$ for which $\alpha_{i}\beta_{i}\neq 0$ is normal.

The Skitovich-Darmois theorem is used in the proof of the base case of our
next characterization. Several generalizations of the Skitovich-Darmois
theorem are described by Kagan et al. (1973).

###### Theorem 9

Let $W$ be an $n\times n$, $n\geq 2$, positive-definite symmetric matrix such
that no entry in $W$ is zero, $\mu$ be an $n$-dimensional vector of random
variables, and $f(\mu)$ be a pdf of $\mu$. Then, f$(\mu)$ is an $n$
dimensional normal distribution $N(\mu|\eta,\gamma W)$ where $\gamma>0$ if and
only if $\mu_{1}$ is independent of
$\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1}$ for every partitioning
$\mu_{1},\mu_{2}$ of $\mu$ where $W_{11}$,$W_{12}$, $W^{\prime}_{12}$,
$W_{22}$ is a block partitioning of $W$ compatible with the partitioning
$\mu_{1},\mu_{2}$.

Proof: The two independence conditions, $\mu_{1}$ independent of
$\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1}$ and $\mu_{2}$ independent of
$\mu_{1}+W_{11}^{-1}W_{12}\mu_{2}$, are equivalent to the following functional
equation

$f(\mu)=f_{1}(\mu_{1})f_{2|\\!|1}(\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1})=f_{2}(\mu_{2})f_{1|\\!|2}(\mu_{1}+W_{11}^{-1}W_{12}\mu_{2})$
(31)

where a subscripted $f$ denotes a pdf. We show that the only solution for $f$
that satisfies this equation is the normal distribution. Consequently, both
the if and only if portions of the theorem will be established.

For $n\geq 3$, we can divide the indices of $W$ into three non-empty sets
$a,b$ and $c$. We group $a$ and $b$ to form a block and $b$ and $c$ to form a
block. For each of the two cases, let $W_{11}$ be the block consisting of the
indices in $\\{a,b\\}$ or $\\{b,c\\}$, respectively, and $W_{22}$ be the block
consisting of the indices of $c$ or $a$, respectively. By the induction
hypothesis applied to both cases and marginalization we can assume that
$f_{1}(\mu_{1})$ is a normal distribution
$N(\mu_{1}|\eta_{1},\gamma_{1}((W^{-1})_{11})^{-1})$ and that
$f_{2}(\mu_{2})=N(\mu_{2}|\eta_{2},\gamma_{2}((W^{-1})_{22})^{-1})$.
Consequently, the pdf of the block corresponding to the indices in $b$ is a
normal distribution, and from the two alternative ways by which this pdf can
be formed, it follows that $\gamma_{1}=\gamma_{2}$.

Let $\gamma=\gamma_{i}$, $i=1,2$, and define

$\displaystyle F_{2|\\!|1}(x)$ $\displaystyle=$ $\displaystyle
f_{2|\\!|1}(x)/N(x|\eta_{2}+W_{22}^{-1}W^{\prime}_{12}\eta_{1},\gamma W_{22})$
$\displaystyle F_{1|\\!|2}(x)$ $\displaystyle=$ $\displaystyle
f_{1|\\!|2}(x)/N(x|\eta_{1}+W_{11}^{-1}W_{12}\eta_{2},\gamma W_{11}).$

By substituting these definitions into Equation 31, substituting the normal
form for $f_{1}(\mu_{1})$ and $f_{2}(\mu_{2})$, and canceling on both sides of
the equation the term $N(\mu|\eta,\gamma W)$ (which is formed by standard
algebra pertaining to quadratic forms (E.g., DeGroot, p. 55)), we obtain a new
functional equation,

$F_{2|\\!|1}(\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1})=F_{1|\\!|2}(\mu_{1}+W_{11}^{-1}W_{12}\mu_{2}).$

By setting $\mu_{2}=-W_{22}^{-1}W^{\prime}_{12}\mu_{1}$, we obtain
$F_{1|\\!|2}((I-(W_{11}^{-1}W_{12})(W_{22}^{-1}W^{\prime}_{12}))\mu_{1})=F_{2|\\!|1}(0)$
for every $\mu_{1}$. Hence, the only solution to this functional equation is
$F_{1|\\!|2}=F_{2|\\!|1}\equiv\mbox{constant}$. Consequently,
$f(\mu)=N(\mu|\eta,\gamma W)$.

It remains to prove the theorem for $n=2$. Let $z_{1}=\mu_{1}$,
$z_{2}=\mu_{2}+w_{22}^{-1}w_{12}\mu_{1}$,
$L_{1}=\mu_{1}+w_{11}^{-1}w_{12}\mu_{2}$, and $L_{2}=\mu_{2}$. By our
assumptions $z_{1}$ and $z_{2}$ are independent and $L_{1}$ and $L_{2}$ are
independent. Furthermore, rewriting $L_{1}$ and $L_{2}$ in terms of $z_{1}$
and $z_{2}$, we get,
$L_{1}=w^{-1}_{11}w^{-1}_{22}(w_{11}w_{22}-w^{2}_{12})z_{1}+w_{11}^{-1}w_{12}z_{2}$
and $L_{2}=z_{2}-w_{22}^{-1}w_{12}z_{1}$. All linear coefficients in this
transformation are non zero due to the fact that $W$ is positive definite and
that $w_{12}$ is not zero. Consequently, due to the Skitovich-Darmois theorem,
$z_{1}$ is normal and $z_{2}$ is normal. Furthermore, since $z_{1}$ and
$z_{2}$ are independent, their joint pdf is normal as well. Finally,
$\\{\mu_{1},\mu_{2}\\}$ and $\\{z_{1},z_{2}\\}$ are related through a non-
singular linear transformation and so $\\{\mu_{1},\mu_{2}\\}$ also have a
joint normal distribution $f(\mu)=N(\mu|\eta,A)$ where $A=(a_{ij})$ is a
$2\times 2$ precision matrix. Substituting this solution into Equation 31 and
comparing the coefficients of $\mu_{1}^{2}$, $\mu_{2}^{2}$, and
$\mu_{1}\mu_{2}$, we obtain $a_{12}/a_{11}=w_{12}/w_{11}$ and
$a_{12}/a_{22}=w_{12}/w_{22}$. Thus $A=\gamma W$ where $\gamma>0$. $\Box$

The proofs of Theorems 7 and 9 can be combined to form the following
characterization of the normal-Wishart distribution.

###### Theorem 10

Let $W$ be an $n\times n$, $n\geq 3$, positive-definite symmetric matrix of
real random variables such that no entry in $W$ is zero, $\mu$ be an
$n$-dimensional vector of random variables, and $f(\mu,W)$ be a joint pdf of
$\\{\mu,W\\}$. Then, $f(\mu,W)$ is an $n$ dimensional normal-Wishart
distribution if and only if
$\\{\mu_{1},W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}\\}$ is independent of
$\\{\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1},W_{12},W_{22}\\}$ for every
partitioning $\mu_{1},\mu_{2}$ of $\mu$ where $W_{11}$,$W_{12}$,
$W^{\prime}_{12}$, $W_{22}$ is a block partitioning of $W$ compatible the
partitioning $\mu_{1},\mu_{2}$.

Proof sketch: The two independence conditions,
$\\{\mu_{1},W_{11}-W_{12}W_{22}^{-1}W^{\prime}_{12}\\}$ independent of
$\\{\mu_{2}+W_{22}^{-1}W^{\prime}_{12}\mu_{1},$ $W_{12},W_{22}\\}$ and
$\\{\mu_{2},W_{22}-W^{\prime}_{12}W_{11}^{-1}W_{12}\\}$ independent of
$\\{\mu_{1}+W_{11}^{-1}W_{12}\mu_{2},$ $W^{\prime}_{12},W_{11}\\}$, are
equivalent to the following functional equation

$\displaystyle f(\mu,W)$ $\displaystyle=$ $\displaystyle
f_{1}(\mu_{1},W_{11}-W_{12}W_{22}^{-1}W_{12}^{\prime})f_{2|\\!|1}(\mu_{2}+W_{22}^{-1}W_{12}^{\prime}\mu_{1},W_{22},W_{12})$
(32) $\displaystyle=$ $\displaystyle
f_{2}(\mu_{2},W_{22}-W_{12}^{\prime}W_{11}^{-1}W_{12})f_{1|\\!|2}(\mu_{1}+W_{11}^{-1}W_{12}\mu_{2},W_{11},W_{12})$

where a subscripted $f$ denotes a pdf. We show that the only solution for $f$
that satisfies this functional equation is the normal-Wishart distribution.
Setting $W$ to a fixed value yields Equation 31 the solution of which is

$\displaystyle f(\mu,W)$ $\displaystyle\propto$ $\displaystyle
N(\mu|\eta(W),\gamma(W)\cdot W)$ $\displaystyle=$ $\displaystyle
N(\mu_{2}|\eta_{2}(W),\gamma(W)\cdot[W_{22}-W^{\prime}_{12}W_{11}^{-1}W_{12}])$
$\displaystyle\ \ \ \ \ \ \ \ \ \cdot N(\mu_{1}|\eta_{1}(W)+\eta_{2}(W)\cdot
W^{-1}_{11}W_{12}-W_{11}^{-1}W_{12}\mu_{2},\gamma(W)\cdot W_{11})$

where both $\gamma$ and $\eta=(\eta_{1},\eta_{2})$ potentially can be
functions of $W$. To see that these quantities in fact do not depend on $W$,
first note that the normal distributions for $\mu_{2}$ and $\mu_{1}$ in Eq. 5
must be proportional to the functions $f_{2}$ and $f_{1|\\!|2}$ in Eq. 32,
respectively. Comparing the form of $f_{2}$ with the normal distribution for
$\mu_{2}$, we see that $\gamma(W)$ and $\eta_{2}(W)$ can only depend on
$W_{22}-W^{\prime}_{12}W_{11}^{-1}W_{12}$. Comparing the form of $f_{1|\\!|2}$
with the normal distribution for $\mu_{1}$, we see that $\gamma(W)$ and
$\eta_{2}(W)$ can only depend on $\\{W_{11},W_{12}\\}$. Consequently,
$\gamma(W)$ and $\eta_{2}(W)$ must be constant. Similarly, $\eta_{1}(W)$ must
be a constant. Substituting these solutions into Equation 32 and dividing by
the common terms which are equal to $f(\mu|W)$ yields Equation 19, the
solution of which for $f$ is a Wishart pdf. $\Box$

Note that the conditions set on $W$ in Theorem 10, namely, a positive-definite
symmetric matrix of real random variables such that no entry in $W$ is zero,
are necessary and sufficient in order for $W$ to be a precision matrix of a
complete Gaussian DAG model.

##  6 Local versus Global Parameter Independence

We have shown that the only pdf for $\\{\mu,W\\}$ which satisfies global
parameter independence, when the number of coordinates is greater than two, is
the normal-Wishart distribution. We now discuss additional independence
assertions implied by the assumption of global parameter independence.

Consider the parameter prior for $\\{m_{n},b_{n},v_{n}\\}$ when the prior for
$\\{\mu,W\\}$ is a normal-Wishart as specified by Equations 14 and 15. By a
change of variables, we get

$\displaystyle f_{n}(m_{n},b_{n},v_{n})=$ $\displaystyle\mbox{\cal
Wishart}(1/v_{n}\;|\;\alpha+n-1,T_{22}-T^{\prime}_{12}T_{11}^{-1}T_{12})\cdot
N(b_{n}\;|\;T_{11}^{-1}T_{12},T_{22}/v_{n})\cdot
N(m_{n}\;|\;\nu_{n},\alpha_{\mu}/v_{n})$

where the first block ($T_{11}$) corresponds to $X_{1},\ldots,X_{n-1}$ and the
second one-dimensional block ($T_{22}$) corresponds to $X_{n}$. We note that
the only independence assumption expressed by this product is that $m_{n}$ and
$b_{n}$ are independent given $v_{n}$. However, by standardizing $m_{n}$ and
$b_{n}$, namely defining,
$m^{*}_{n}=(m_{n}-\nu_{n})/(\alpha_{\mu}/v_{n})^{1/2}$ and
$b^{*}_{n}=(T_{22}/v_{n})^{1/2}(b_{n}-T_{11}^{-1}T_{12})$, which is well
defined because $T_{22}$ is positive definite and $v_{n}>0$, we obtain a set
of parameters $(m^{*}_{n},b^{*}_{n},v_{n})$ which are mutually independent.
Furthermore, this mutual independence property holds for every local family
and for every Gaussian DAG model over $X_{1},\ldots,X_{n}$. We call this
property the standard local independence for Gaussian DAG models.

This observation leads to the following corollary of our characterization
theorems.

###### Corollary 11

If global parameter independence holds for every complete Gaussian DAG model
over $X_{1},\ldots,X_{n}$ ($n\geq 3)$, then standard local parameter
independence also holds for every complete Gaussian DAG model over
$X_{1},\ldots,X_{n}$.

This corollary follows from the fact that global parameter independence
implies that, due to Theorem 10, the parameter prior is a normal-Wishart, and
for this prior, we have shown that standard local parameter independence must
hold.

It is interesting to note that when $n=2$, there are distributions that
satisfy global parameter independence but do not satisfy standard local
parameter independence. In particular, a prior for a $2\times 2$ positive
definite matrix $W$ which has the form $\mbox{\cal
Wishart}(W|\alpha,T)H(w_{12})$, where $H$ is some real function and $w_{12}$
is the off-diagonal element of $W$, satisfies global parameter independence
(as shown in the appendix) but need not satisfy standard local parameter
independence. Furthermore, if standard local parameter independence is
assumed, then $H(w_{12})$ must be proportional to $e^{aw_{12}}$, which means
that, for $n=2$, the only pdf for $W$ that satisfies global and standard local
parameter independence is the bivariate Wishart distribution. In contrast, for
$n>2$, global parameter independence alone implies a Wishart prior.

##  7 Discussion

The formula for the marginal likelihood applies whenever Assumptions 1 through
5 are satisfied, not only for Gaussian DAG models. Another important special
case is when all variables in ${\bf X}$ are discrete and all local
distributions are multinomial. This case has been treated in Heckerman et al.
(1995) and Geiger and Heckerman (1997) under the additional assumption of
local parameter independence which was introduced by Spiegelhalter and
Lauritzen (1990). Our generalized derivation herein dispenses this assumption
and unifies the derivation in the discrete case with the derivation needed for
Gaussian DAG models.

Our characterization means that the assumption of global parameter
independence when combined with the definition of $m^{h}$, the assumption of
complete model equivalence, and the regularity assumption, may be too
restrictive. One common remedy for this problem is to use a hierarchical prior
$p(\theta|\eta)p(\eta)$ with hyperparameters $\eta$. When such a prior is used
for Gaussian DAG models, our results show that for every value of $\eta$ for
which global parameter independence holds, $p(\theta|\eta)$ must be a normal-
Wishart distribution. The difficulty with this approach is that the marginal
likelihood no longer has closed form and therefore approximate methods such as
MCMC are usually employed to compute the marginal likelihood. Also the
elicitation of hierarchical priors is often difficult. Other alternative
approaches have been discussed at the end of Section 3.

We conclude with a technical comment. Equation 19, which encodes global
parameter independence for an unknown covariance matrix, is an interesting
example of a matrix functional equation. The domain of each unknown function
is a non-singular matrix and the range is $R$. A well known functional
equation of this sort is the equation

$f(XY)=f(X)f(Y)$ (34)

where $X$ and $Y$ are non-singular matrices. The general solution of this
equation is $f(X)=|X|^{\alpha}$ or $f(X)=|X|^{\alpha}sgn(|X|)$ (e.g., Aczél ,
1966). When the domain of $f$ is the set of positive definite matrices, the
solution is simply $f(X)=|X|^{\alpha}$.

We note that the solution of Equation 34 is obtained for matrices over
arbitrary fields. Only algebraic manipulations are used in its proof. It seems
reasonable to believe and interesting to investigate, whether a solution to
Equation 19 can be obtained via purely algebraic manipulations. The proof
technique that we have employed, however, especially for the base case of the
induction, uses the fact that the matrices are over the real numbers.

## Acknowledgments

We thank Chris Meek for helping us shape the definition of DAG models and
correcting earlier versions of this manuscript, Bo Thiesson for implementing
the proposed scheme, and Jim Kajiya for his help in regard to the
characterization theorems. We also thank János Aczél , Enrique Castillo, Clark
Glymour, Antal Járai , Gerard Letac, Helene Masson, Peter Spirtes, and the
reviewers, for their useful suggestions. The first four sections of this work
with different emphasizes (including, for example, details of the derivation
of Equation 4) have been reported in [Geiger and Heckerman, 1994, Heckerman
and Geiger, 1995]. A short version of this work appeared in [Geiger and
Heckerman, 1999].

## References

  * [Aczél, 1966] Aczél , J. (1966). Lectures on Functional Equations and Their Applications. Academic Press, New York.
  * [Andersson, Madigan, and Perlman, 1997] Andersson, S. A., Madigan, D,. and Perlman, M.D. (1997) A characterization of Markov equivalence classes for acyclic digraphs.  Annals of Statistics, 25(2):1272–1317.
  * [Bernardo and Smith, 1994] Bernardo, J. and Smith, A. (1994).  Bayesian Theory.  John Wiley and Sons, New York.
  * [Buntine, 1994] Buntine, W. (1994).  Operations for learning with graphical models.  Journal of Artificial Intelligence Research, 2:159–225.
  * [Chickering, 1995] Chickering, D. (1995).  A transformational characterization of equivalent Bayesian networks structures.  In Proceedings of Eleventh Conference on Uncertainty in Artificial Intelligence, Montreal, QU, pages 87–98. Morgan Kaufmann.
  * [Chickering, 1996] Chickering, D. (1996).  Learning Bayesian networks from data,  Ph.D dissertation, University of California, Los Angeles, CA.
  * [Cooper and Herskovits, 1992] Cooper, G. and Herskovits, E. (1992).  A Bayesian method for the induction of probabilistic networks from data.  Machine Learning, 9:309–347.
  * [Cowell, Dawid, Lauritzen, and Spiegelhalter, 1999] Cowell, R., Dawid, P., Lauritzen, S. and Spiegelhalter D. (1999).  Probabilistic networks and expert systems.  Springer-Verlag, New York.
  * [Dawid and Lauritzen, 1993] Dawid, P. and Lauritzen, S. (1993).  Hyper Markov laws in statistical analysis of decomposable graphical models.  Annals of Statistics, 21:1272–1317.
  * [DeGroot, 1970] DeGroot, M. (1970).  Optimal Statistical Decisions.  McGraw-Hill, New York.
  * [Friedman and Goldszmidt, 1997] Friedman, N., and Goldszmidt, M. (1997).  Sequential update of Bayesian networks structure.  In Proceedings of Thirteenth Conference on Uncertainty in Artificial Intelligence, Providence, RI, pages 165–174. Morgan Kaufmann.
  * [Geiger and Heckerman, 1994] Geiger, D. and Heckerman, D. (1994).  Learning Gaussian networks.  In Proceedings of Tenth Conference on Uncertainty in Artificial Intelligence, Seattle, WA, pages 235–243. Morgan Kaufmann, arXiv:1302.6808.
  * [Geiger and Heckerman, 1997] Geiger, D. and Heckerman, D. (1997).  A characterization of the Dirichlet distribution through global and local parameter independence.  Annals of Statistics, 25(3):1344–1369.
  * [Geiger and Heckerman, 1998] Geiger, D. and Heckerman, D. (1998).  A characterization of the bivariate Wishart distribution.  Probability and Mathematical Statistics, 18(1):119–131.
  * [Geiger and Heckerman, 1999] Geiger, D. and Heckerman, D. (1999).  Parameter priors for directed graphical models and the characterization of several probability distributions.  In Proceedings of Fifteenth Conference on Uncertainty in Artificial Intelligence, Stockholm, Sweden, pages 216–225. Morgan Kaufmann.
  * [Heckerman and Geiger, 1995] Heckerman, D. and Geiger, D. (1995).  Learning Bayesian networks: A unification for discrete and Gaussian domains.  In Proceedings of Eleventh Conference on Uncertainty in Artificial Intelligence, Montreal, QU, pages 274–284. Morgan Kaufmann.
  * [Heckerman, Geiger, and Chickering, 1995] Heckerman, D., Geiger, D., and Chickering, D. (1995a).  Learning Bayesian networks: The combination of knowledge and statistical data.  Machine Learning, 20, pages 197–243.
  * [Heckerman, Mamdani, and Wellman, 1995] Heckerman, D., Mamdani, A., and Wellman, M. (1995b).  Real-world applications of Bayesian networks.  Communications of the ACM, 38.
  * [Howard and Matheson, 1981] Howard, R. and Matheson, J. (1981).  Influence diagrams.  In Howard, R. and Matheson, J., editors, Readings on the Principles and Applications of Decision Analysis, volume II, pages 721–762. Strategic Decisions Group, Menlo Park, CA.
  * [Járai, 1986] Járai , A. (1986). On regular solutions of functional equations. Aequationes Mathematicae, 30:21–54.
  * [Járai, 1998] Járai , A. (1998). Regularity property of the functional equation of the Dirichlet distribution. Aequationes Mathematicae, 56:37–46.
  * [Kadane,Dickey, Winkler, Smith, and Peters, 1980] Kadane, J.B., Dickey, J.M., Winkler, R.L., Smith, W.S., and Peters, S.C. (1980). Interactive elicitation of opinion for a normal linear model. Journal American Statistical Association, 75:845–854.
  * [Kagan, Linnik, and Rao, 1973] Kagan, A.M., Linnik, Y.V., and Rao, C.R. (1973). Characterization problems in Mathematical Statistics. Wiley & Sons.
  * [Madigan, Andersson, Perlman, and Volinksy, 1996] Madigan, D., Andersson, S. A., Perlman, M.D., and Volinsky, C.T. (1996). Model averaging and model selection for Markov equivalence classes of acyclic digraphs.  Communications in Statistics - Theory and Methods, 25:2493-2520.
  * [Pearl, 1988] Pearl, J. (1988).  Probabilistic reasoning in intelligent systems: Networks of plausible inference.  Morgan Kaufmann, San Mateo, CA.
  * [Press, 1971] Press S. (1971). Applied Multivariate Analysis, Holt, Rinehart, and Winston, New York.
  * [Shachter and Kenley, 1989] Shachter, R. and Kenley, C. (1989).  Gaussian influence diagrams.  Management Science, 35:527–550.
  * [Spiegelhalter, Dawid, Lauritzen, and Cowell, 1993] Spiegelhalter, D., Dawid, A., Lauritzen, S., and Cowell, R. (1993).  Bayesian analysis in expert systems.  Statistical Science, 8:219–282.
  * [Spiegelhalter and Lauritzen, 1990] Spiegelhalter, D. and Lauritzen, S. (1990).  Sequential updating of conditional probabilities on directed graphical structures.  Networks, 20:579–605.
  * [Spirtes, Glymour, and Scheines, 2001] Spirtes, P., Glymour, C., and Scheines, R. (2001).  Causation, Prediction, and Search.  MIT Press, Cambridge, MA.
  * [Spirtes and Meek, 1995] Spirtes, P. and Meek, C. (1995).  Learning Bayesian networks with discrete variables from data.  In Proceedings of First International Conference on Knowledge Discovery and Data Mining, Montreal, QU, pages 294–299. Morgan Kaufmann.
  * [Thiesson, Meek, Chickering, and Heckerman, 1998] Thiesson B., Meek C., Chickering D., and D. Heckerman (1998).  Computationally efficient methods for selecting among mixtures of graphical models.  In Proceedings of the Sixth Valencia International Meeting on Bayesian Statistics, Alcossebre, Spain, pages 223-240,
  * [Verma and Pearl, 1990] Verma, T., and Pearl, J. (1990).  Equivalence and synthesis of causal models.  In Proceedings of Sixth Conference on Uncertainty in Artificial Intelligence, Boston, MA, pages 220–227. Morgan Kaufmann.

## Appendix

We now characterize the pdfs of an unknown $2\times 2$ precision matrix that
satisfy global parameter independence. This result has been obtained in
(Geiger and Heckerman, 1998) under additional regularity conditions.

###### Theorem 12

Let $W$ be a $2\times 2$ positive-definite symmetric matrix with random
entries $w_{11},w_{12}$, and $w_{22}$ and let $f(W)$ be a pdf of $W$. Then,
$f(W)=|W|^{\beta}e^{\mbox{\it tr}\\{TW\\}}H(w_{12})$ where $H$ is a real
function if and only if $w_{11}-w^{2}_{12}/w_{22}$ is independent of
$\\{w_{12},w_{22}\\}$ and $w_{22}-w^{2}_{12}/w_{11}$ is independent of
$\\{w_{12},w_{11}\\}$.

Proof: That $w_{11}-w^{2}_{12}/w_{22}$ is independent of $\\{w_{12},w_{22}\\}$
whenever $f(W)$ is a Wishart distribution (e.g., when $H(x)=\mbox{constant}$)
is a well known fact (Press 1971, p. 117-119). Consequently, this claim holds
for any real function $H$. We prove the other direction by solving the
functional equation, which is implied by the given independence assumptions,

$f(W)=f_{1}(w_{11}-w^{2}_{12}/w_{22})f_{2|\\!|1}(w_{22},w_{12})=f_{2}(w_{22}-w^{2}_{12}/w_{11})f_{1|\\!|2}(w_{11},w_{12})$
(35)

where a subscripted $f$ denotes a pdf. To solve this functional equation,
namely to find all pdfs that satisfy it, we use techniques described in
[Aczél, 1966] and results from [Járai, 1986, Járai, 1998].

Let $w_{12}$ be a value such that the integral of $f_{2|\\!|1}(x,w_{12})$ over
the domain of $x$ is not identically zero. Such a value for $w_{12}$ exists
because $f_{2|\\!|1}(x,w_{12})$ integrates to 1 over its domain. Without loss
of generality, suppose this value of $w_{12}$ is 1, lest we can modify the
scale using the transformations $w_{11}\leftarrow w_{12}w_{11}$ and
$w_{22}\leftarrow w_{12}w_{22}$. We rewrite Equation 35 as

$f_{1}(w_{11}-1/w_{22})f_{2|\\!|1}(w_{22},1)=f_{2}(w_{22}-1/w_{11})f_{1|\\!|2}(w_{11},1).$
(36)

We claim that all density functions satisfying Equation 36 must be positive
everywhere and smooth. This is shown in Lemmas 14 and 16 at the end of the
proof. Consequently, we can take the logarithm of Equation 35 and then take
derivatives. First we take the logarithm and rename the functions. We get

$g_{1}(w_{11}-1/w_{22})+g_{2|\\!|1}(w_{22})=g_{2}(w_{22}-1/w_{11})+g_{1|\\!|2}(w_{11})$
(37)

where $g_{1}(x)=\ln f_{1}(x)$, $g_{2|\\!|1}(x)=\ln f_{2|\\!|1}(x,1)$, and
where $g_{2}$ and $g_{1|\\!|2}$ are defined analogously.

We take a mixed second derivative with respect to $w_{11}$ and $w_{22}$ of
Equation 37. We get

$g^{\prime\prime}_{1}(w_{11}-1/w_{22})/w^{2}_{22}=g^{\prime\prime}_{2}(w_{22}-1/w_{11})/w^{2}_{11}$
(38)

By substituting $w_{11}=w_{22}$ we obtain
$g^{\prime\prime}_{1}=g^{\prime\prime}_{2}$. We denote this function by $h$
and so,

$w^{2}_{11}h(w_{11}-1/w_{22})=w^{2}_{22}h(w_{22}-1/w_{11})$ (39)

It is easy to show, using this functional equation for $h$, that if $h$ were
zero at some point then $h$ must be identically zero, if $h$ is positive at
one point then $h$ is positive everywhere, and if $h$ is negative at one point
then $h$ is negative everywhere. We now take a derivative wrt $w_{11}$ and a
derivative with respect to $w_{22}$

$2w_{11}h(w_{11}-1/w_{22})+w_{11}^{2}h^{\prime}(w_{11}-1/w_{22})=\\{w_{22}/w_{11}\\}^{2}h^{\prime}(w_{22}-1/w_{11})$
$2w_{22}h(w_{22}-1/w_{11})+w_{22}^{2}h^{\prime}(w_{22}-1/w_{11})=\\{w_{11}/w_{22}\\}^{2}h^{\prime}(w_{11}-1/w_{22}).$

From these equations, and using Equation 39, we get

$2(w_{22}+1/w_{11})h(w_{22}-1/w_{11})=-(w^{2}_{22}-1/w^{2}_{11})h^{\prime}(w_{22}-1/w_{11})$

Consequently,

$h^{\prime}(x)/h(x)=-2/x$

where $x=w_{22}-1/w_{11}$. This equation holds for every $x\in R^{+}$.
Assuming $h$ is positive everywhere, we have $(\ln h(x))^{\prime}=-2/x$ and so
$\ln h(x)=\ln x^{-2}+c^{\prime}$ where $c^{\prime}$ is a constant. If $h$ is
negative everywhere, we have $(\ln-h(x))^{\prime}=-2/x$ and so $\ln(-h(x))=\ln
x^{-2}+c^{\prime}$. Consequently, whether $h$ is positive everywhere, negative
everywhere, or identically zero, it has the form $h(x)=c/x^{2}$ where $c$ is a
constant. Recall that $h=\left(\ln f_{1}\right)^{\prime\prime}$. Hence,
$f_{1}(x)=c_{1}x^{-c}e^{c_{2}x}$ and similarly
$f_{2}(x)=c^{\prime}_{1}x^{-c}e^{c^{\prime}_{2}x}$ (i.e., one-dimensional
Wishart distributions with the same degrees of freedom). We conclude by
substituting $f_{1}$ and $f_{2}$ into Equation 35 and proceeding as in
Equations 5 through 23. $\Box$

The next lemma shows that every positive everywhere pdf that satisfies
Equation 36 must be smooth. Our lemma is an immediate consequence of Járai’s
Theorem which we now state.

######  Theorem 13 (Járai , 1986,1998)

Let $X_{i}$ be an open subset of $R^{r_{i}}$ $(i=1,2,\ldots,n)$, $T$ be an
open subset of $R^{s}$, $Y$ be an open subset of $R^{k}$, $Z_{i}$ be an open
subset of $R^{m_{i}}$ $(i=1,2,\ldots,n)$, $D$ be an open subset of $T\times Y$
and let $Z$ be an Euclidean space. Consider the functions $f:T\rightarrow Z$,
$g_{i}:D\rightarrow X_{i}$, $f_{i}:X_{i}\rightarrow Z_{i}$, $h_{i}:D\times
Z_{i}\rightarrow Z$, (i=1,2,…,n). Suppose that $0\leq p\leq\infty$ and

  * (i)

for each $(t,y)\in D$,

$f(t)=\sum_{i=1}^{n}h_{i}(t,y,f_{i}(g_{i}(t,y)));$

  * (ii)

$h_{i}$ is $p+1$ times continuously differentiable $(1\leq i\leq n)$;

  * (iii)

$g_{i}$ is $p+2$ times continuously differentiable and for each $t\in T$ there
exists a $y\in Y$ such that $(t,y)\in D$ and
$\frac{\partial{g_{i}}}{\partial{y}}(t,y)$ has rank $r_{i}$, $1\leq i\leq n$.

Then

  * (iv)

if $f_{i}$ $(i=1,2,\ldots,n)$ is Lebesgue measurable and (ii),(iii) are
satisfied with $p=0$ then $f$ is continuous on $T$;

  * (v)

if $f_{i}$ $(i=1,2,\ldots,n)$ is continuous and (ii),(iii) are satisfied with
$p=0$ then $f$ is continuously differentiable on $T$;

  * (vi)

if $f_{i}$ $(i=1,2,\ldots,n)$ is $p$ times continuously differentiable and
(ii),(iii) are satisfied then $f$ is $p+1$ continuously differentiable on $T$.

This theorem is stated in [Járai, 1998] and its proof is based on Theorems
3.3, 5.2, and 7.2 of [Járai, 1986]. A simple corollary of Járai ’s theorem is
the following.

###### Lemma 14

Every Lebesgue measurable real functions $l_{1},l_{2},l_{1|\\!|2}$ and
$l_{2|\\!|1}$ defined on $R^{+}$ which satisfy

$l_{1}(y-1/t)+l_{2|\\!|1}(t)=l_{2}(t-1/y)+l_{1|\\!|2}(y)$ (40)

for every $y,t>0$ such that $yt>1$, are $p$ times continuously differentiable
where $p$ is arbitrary large.

Proof: The proof follows closely the lines of reasoning that Járai (1998)
applied to another functional equation.

Using statement $(iv)$ of Theorem 13 we show that $l_{2|\\!|1}$ is continuous.
To match Járai ’s theorem notation we define $f=l_{2|\\!|1}$, $f_{1}=-l_{1}$,
$f_{2}=l_{2}$, $f_{3}=l_{1|\\!|2}$, $h_{i}(t,y,w)=w$ for $i=1,2,3$,
$g_{1}(t,y)=(y-1/t)$, $g_{2}(t,y)=(t-1/y)$, and $g_{3}(t,y)=y$. The only non
obvious condition to check is that for each $t\in R^{+}$ there exists a $y\in
R^{+}$ such that $ty>1$ and $\frac{\partial{g_{i}}}{\partial{y}}(t,y)$ has
rank $r_{i}$, $1\leq i\leq n$. But here the rank is 1 and so we just need to
observe that there exists a $y$ such that
$\frac{\partial{g_{i}}}{\partial{y}}(t,y)$ is not zero.

To show that $l_{1}$ is continuous rewrite Equation 40 as

$l_{1}(t)+l_{2|\\!|1}(y)=l_{2}(\frac{ty^{2}}{ty+1})+l_{1|\\!|2}(t+1/y)$ (41)

where $t,y>0$. Now define $f=l_{1}$, $f_{1}=-l_{2|\\!|1}$, $f_{2}=l_{2}$,
$f_{3}=l_{1|\\!|2}$, $h_{i}(t,y,w)=w$ for $i=1,2,3$, $g_{1}(t,y)=y$,
$g_{2}(t,y)=\frac{ty^{2}}{ty+1}$, and $g_{3}(t,y)=t+1/y$. Observe that the
conditions of Járai ’s theorem hold and so $f=l_{1}$ is continuous. By the
symmetry of the equation, $l_{2}$ and $l_{1|\\!|2}$ are also continuous on
$R^{+}$.

Now we can apply statement (v) of Járai ’s Theorem. We obtain, in the same way
as above, that all four functions are continuously differentiable. Finally,
applying statement (vi) of Járai ’s Theorem in the same way, we get that all
four functions are twice continuously differentiable. Repeating this process
shows that all four functions are $p$ times continuously differentiable for
every $p>0$. $\Box$

The next Theorem and lemma show that every pdf that satisfies Equation 36 must
be positive everywhere and so taking the logarithm of this equation, as we
have done, is legitimate. We denote by $\lambda^{s}$ the s-dimensional
Lebesgue measure and by $\lambda$ the one dimensional Lebesgue measure.

######  Theorem 15 (Járai , 1995, 1998)

Let $X_{1},\ldots,X_{n}$ be orthogonal subspaces of $R^{r}$ of dimensions
$r_{1},\ldots r_{n}$, respectively. Suppose that $r_{i}\geq 1$ $(1\leq i\leq
n)$ and $\sum_{i=1}^{n}r_{i}=r$. Let $U$ be an open subset of $R^{r}$ and
$F:U\rightarrow R^{m}$ be a continuously differentiable function. For each
$x\in U$, let $N_{x}$ denote the nullspace of $F^{\prime}(x)$. Let $p_{i}$
denote the orthogonal projection of $X$ onto $X_{i}$. Suppose that
$dimN_{x}=r-m$ and $p_{i}(N_{x})=X_{i}$ $(i=1,\ldots,n)$ for all $x\in U$. Let
$A_{i}$ be a subset of $X_{i}$ $(i=1,\ldots,n)$. If $A_{1}\times
A_{2}\times\ldots\times A_{n}\subset U$, and $A_{i}$ is $\lambda^{r_{i}}$
measurable with $\lambda^{r_{i}}(A_{i})>0$ $(1\leq i\leq n)$, then
$F(A_{1}\times A_{2}\times\ldots\times A_{n})$ contains a non-empty open set.

Recall that if $X_{1},\ldots,X_{n}$ are the standard orthogonal axis of
$R^{n}$, then $p_{i}(X_{1},\ldots,X_{n})=X_{i}$, and
$P_{i}(N_{x})=\\{x|(X_{1},\ldots,X_{i-1},x,X_{i+1},\ldots,X_{n})\in N_{x}\\}$.

###### Lemma 16

Let $f,g,h,k$ be non-negative real functions that are Lebesgue integrable with
integral $c>0$. If these functions satisfy

$f(s-1/t)g(t)=h(t-1/s)k(s)$ (42)

for every $s,t>0$ such that $st>1$, then they are everywhere positive.

Proof: The proof follows closely the lines of reasoning that Járai (1998)
applied to another functional equation.

Let $\\{f=0\\}$ denote the set of points in the domain of $f$ for which $f$ is
zero and let $\\{f\neq 0\\}$ denote the complementary set of all points in the
domain for which $f$ is not zero—namely, the set of points for which $f$ is
positive. Similar notation is used for the functions $g,h$ and $k$. The idea
of the proof is to show that the set $\\{f=0\\}$ and the set $\\{f\neq 0\\}$
are both open and therefore, since the domain of $f$ is connected, one of
these sets must be empty. The set $\\{f\neq 0\\}$ cannot be empty because $f$
is non-negative and integrates to a positive constant and so $\\{f=0\\}$ must
be empty as claimed by the Theorem. Similar arguments show that $g,h$ and $k$
are also positive everywhere.

The proof proceeds in three steps. First we use Theorem 15 to establish that
the set $\\{g\neq 0\\}$ contains a non-empty open set (i.e., it contains an
inner point). Then we show that every point in $\\{f\neq 0\\}$ is an inner
point and so $\\{f\neq 0\\}$ is open. Finally we show that every point in
$\\{f=0\\}$ is an inner point and so $\\{f=0\\}$ is open as well. Similar
arguments work for $g,h$ and $k$.

We start by rewriting Equation 42 in two symmetric ways. First as

$f(y)g(z)=h(x(y,z))k(w(y,z))$ (43)

for all $y>0$, and $z>0$ where $x(y,z)=yz^{2}/(yz+1)$ and $w(y,z)=y+1/z$.
Second as

$f(y(x,w))g(z(x,w))=h(x)k(w)$ (44)

for all $x>0$, and $w>0$ where $y(x,w)=xw^{2}/(xw+1)$ and $z(x,w)=x+1/w$.

Step I. We show that $\\{g\neq 0\\}$ contains an inner point. Since both $h$
and $k$ integrate to a positive constant, there must exist two $\lambda$
measurable sets $A_{h}$ in $\\{h\neq 0\\}$ and $A_{k}$ in $\\{k\neq 0\\}$ such
that $\lambda(A_{h})>0$ and $\lambda(A_{k})>0$. The image of these sets under
$z(x,w)=x+1/w$ contains an inner point $z$ according to Theorem 15. This
theorem is applicable because the nullspace of $z^{\prime}$ is
$\\{a(1/w^{2},1))|a>0\\}$ and its projection on either of the two coordinates
is $R^{+}$. Due to Equation 44, and because the right hand side is not zero
for any $x\in A_{h}$ and $w\in A_{k}$, each term on the left hand side is also
not zero. Consequently, their image under $z(x,w)$, which includes an inner
point, belongs to $\\{g\neq 0\\}$.

Step II. Let $y$ be an arbitrary point in $\\{f\neq 0\\}$. We now show that
$y$ is an inner point and so $\\{f\neq 0\\}$ is open. Let $z$ be an inner
point in $\\{g\neq 0\\}$. It follows that the image of a sufficiently small
open set containing $z$ under $x(y,z)=yz^{2}/(yz+1)$ and under $w(y,z)=y+1/z$
are open sets. These images belong to $\\{h\neq 0\\}$ and $\\{k\neq 0\\}$,
respectively, because the left hand side of Equation 43 is positive. Now we
fix $x$ in the image and vary $w$ in a small open neighborhood. Then $y$ is
varied in a small open neighborhood. Since the right hand side of Equation 43
is positive, the neighborhood of $y$ belongs to $\\{f\neq 0\\}$ and so $y$ is
an inner point. Similar arguments show that $\\{g\neq 0\\}$ is open as well.
By the symmetry of Equation 42 the same claim holds for $h$ and $k$.

Step III. Let $y$ be an arbitrary point in $\\{f=0\\}$. We now show that $y$
is an inner point and so $\\{f=0\\}$ is open. Let $z$ be an inner point in
$\\{g\neq 0\\}$. It follows that the image of a sufficiently small open set
containing $z$ under $x(y,z)=yz^{2}/(yz+1)$ and under $w(y,z)=y+1/z$ are open
sets. Since the left hand side of Equation 43 is zero, at least one term in
the right hand side must be zero. If $x$ is in $\\{h=0\\}$, then fix $x$. As
we vary $w$ in a small open neighborhood in the image, $g$ remains positive
due to continuity. Also $y$ is varied in a small open neighborhood. Since the
right hand side of Equation 43 is zero, the neighborhood of $y$ belongs to
$\\{f=0\\}$ and so $y$ is an inner point. The other case is when $w$ is in
$\\{k=0\\}$, in which case we fix $w$ and vary $x$ in a small neighborhood.
Similar arguments show that $\\{g=0\\}$ is open as well. By the symmetry of
Equation 42, the same claim holds for $h$ and $k$. $\Box$

Note that Theorems 14 and 16 together imply that every pdf that solves
Equation 42 must be positive everywhere and smooth.