Bayesian parameter estimation via variational methods.
We consider a logistic regression model with a Gaussian prior distribution over the parameters. We show that an accurate variational transformation can be used to obtain a closed form approximation to the posterior distribution of the parameters thereby yielding an approximate posterior predictive model. This approach is readily extended to binary graphical model with complete observations. For graphical models with incomplete observations we utilize an additional variational transformation and again obtain a closed form approximation to the posterior. Finally, we show that the dual of the regression problem gives a latent variable density model, the variational formulation of which leads to exactly solvable EM updates.
Introduction
Bayesian methods have a number of virtues, particularly their uniform treatment
of uncertainty at all levels of the modeling process. The formalism also allows
ready incorporation of prior knowledge and the seamless combination of such
knowledge with observed data (Bernardo & Smith 1994, Gelman 1995, Heckerman et
al. 1995). The elegant semantics, however, often comes at a sizable computational
cost|posterior distributions resulting from the incorporation of observed data must
be represented and updated, and this generally involves high-dimensional integra-
tion. The computational cost involved in carrying out these operations can call into
question the viability of Bayesian methods even in relatively simple settings, such
as generalized linear models (McCullagh & Nelder 1983). We concern ourselves in
this paper with a particular generalized linear model|logistic regression|and we
focus on Bayesian calculations that are computationally tractable. In particular we
describe a
exible deterministic approximation procedure that allows the posterior
distribution in logistic regression to be represented and updated e-ciently. We also
show how our methods permit a Bayesian treatment of a more complex model|a directed
graphical model (a \belief network") in which each node is a logistic regression
model.
The deterministic approximation methods that we develop in this paper are known
generically as variational methods. Variational techniques have been used extensively
in the physics literature (see, e.g., Parisi 1988, Sakurai 1985) and have also found
applications in statistics (Rustagi 1976). Roughly speaking, the objective of these
methods is to transform the problem of interest into an optimization problem via the
introduction of extra degrees of freedom known as variational parameters. For xed
values of the variational parameters the transformed problem often has a closed form
solution, providing an approximate solution to the original problem. The variational
parameters are adjusted via an optimization algorithm to yield an improving sequence
of approximations. For an introduction to variational methods in the context of
graphical models see Jordan et al. (1999).
Let us brie
y sketch the variational method that we develop in this paper. We
study a logistic regression model with a Gaussian prior on the parameter vector. Our
variational transformation replaces the logistic function with an adjustable lower
bound that has a Gaussian form; that is, an exponential of a quadratic function of
the parameters. The product of the prior and the variationally transformed likelihood
thus yields a Gaussian expression for the posterior (conjugacy), which we optimize
variationally. This procedure is iterated for each successive data point.
Our methods can be compared to the Laplace approximation for logistic regression
(cf. Spiegelhalter & Lauritzen 1990), a closely related method which also utilizes a
Gaussian approximation to the posterior. To anticipate the discussion in following
sections, we will see that the variational approach has an advantage over the Laplace
approximation; in particular, the use of variational parameters gives the variational
approach greater
exibility. We will show that this
exibility translates into improved
accuracy of the approximation.
Variational methods can also be contrasted with sampling techniques, which have
become the method of choice in Bayesian statistics (Thomas et al. 1992, Neal 1993,
Gilks et al. 1996). Sampling techniques enjoy wide applicability and can be powerful
in evaluating multi-dimensional integrals and representing posterior distributions.
They do not, however, yield closed form solutions nor do they guarantee monotonically
improving approximations. It is precisely these features that characterize
variational methods.
The paper is organized as follows. First we describe in some detail a variational
approximation method for Bayesian logistic regression. This is followed by an evaluation
of the accuracy of the method and a comparison to Laplace approximation. We
then extend the framework to belief networks, considering both complete data and
incomplete data. Finally, we consider the dual of the regression problem and show
that our techniques lead to exactly solvable EM updates.
Bayesian logistic regression
We begin with a logistic regression model given by:
is the logistic function, S the binary response variable, and
the set of explanatory variables. We represent the uncertainty in
the parameter values  via a prior distribution P () which we assume to be a Gaussian
with possibly full covariance structure. Our predictive distribution is therefore:
Z
In order to utilize this distribution we need to be able to compute the posterior parameter
distribution assume that each D
r g
is a complete observation. This calculation is intractable for large n or T , thus we
consider a variational approximation.
Our approach involves nding a variational transformation of the logistic function
and using this transformed function as an approximate likelihood. In particular
we wish to consider transformations that combine readily with a Gaussian prior in
the sense that the Gaussian prior becomes the conjugate prior to the transformed
likelihood. We begin by introducing the type of variational transformations we will
use for this purpose.
2.1 A brief introduction to variational methods
Consider any continuously dierentiable convex function f(z). Figure 1 provides
an example of a convex function that we will make use of later on. Convexity of
this function guarantees by denition that any tangent line always remains below
the function itself. We may thus interpret the collection of all the tangent lines as
a parameterized family of lower bounds for this convex function (cf. convex duality,
Rockafellar 1976). The tangents in this family are naturally parameterized by their
locations. From the point of view of approximating the convex non-linear function
f , it seems natural to use one of the simpler tangent lines as a lower bound. To
formulate this a little more precisely, let L(z; z 0 ) be the tangent line at
@z
then it follows that f(z)  L(z; z 0 ) for all z; z 0 and f(z 0 In the terminology
of variational methods, L(z; z 0 ) is a variational lower bound of f(z) where the
parameter z 0 is known as the variational parameter. Since the lower bound L(z; z 0 )
is considerably simpler (linear in this case) than the non-linear function f(z), it may
be attractive to substitute the lower bound for f . Note that we are free to adjust
the variational parameter z 0 , the location of the tangent, so as to make L(z; z 0 ) as
accurate an approximation of f(z) as possible around the point of interest, i.e., when
z  z 0 . The quality of this approximation degrades as z receeds from z 0 ; the rate at
which this happens depends on the curvature of f(z). Whenever the function f has
relatively low curvature as is the case in Figure 1, the adjustable linear approximation
seems quite attractive.
-5

Figure

1: A convex function f and its two tangent lines. The locations of the tangents
are indicated with short vertical line segments.
2.2 Variational methods in Bayesian logistic regression
Here we illustrate how variational methods, of the type described above, can be
used to transform the logistic likelihood function into a form that readily combines
with the Gaussian prior (conjugacy). More precisely, the transformed logistic function
should depend on the parameters  at most quadratically in the exponent. We begin
by symmetrizing the log logistic function:
log
and noting that function in the variable
x 2 . (This is readily veried by taking second derivatives; the behavior of f(x) as
a function of x 2 is shown in Figure 1). As discussed above, a tangent surface to a
convex function is a global lower bound for the function and thus we can bound f(x)
globally with a rst order Taylor expansion in the variable x
Note that this lower bound is exact whenever  . Combining this result with Eq.
(4) and exponentiating yields the desired variational transformation of the logistic
where tanh(=2)=(4). We also introduce the following
that is, P denotes the variational lower bound on the logistic function
As a lower bound it is no longer normalized. We refer to eq. (8) as a -
transformation of the conditional probability.
For each xed value of H S we can in fact recover the exact value of the logistic
function via a particular choice of the variational parameter. Indeed, maximizing the
lower bound with respect to  yields substituting this value back into the
lower bound recovers the original conditional probability. For all other values of  we
obtain a lower bound.
The true posterior P (jD) can be computed by normalizing P (SjX; )P (). Given
that this calculation is not feasible in general, we instead form the bound:
and normalize the variational approximation P (SjX; ; )P (). Given that P () is
Gaussian and given our choice of a Gaussian variational form for P (SjX; ; ), the
normalized variational distribution is a Gaussian. Note that although P
a lower bound on the true conditional probability, our variational posterior approximation
is a proper density and thus no longer a bound. This approximate Bayesian
update amounts to updating the prior mean  and the prior covariance matrix
into the posterior mean and the posterior covariance matrix. Omitting the algebra
we nd that the updates take the following form:
for a single observation (S; X), where Successive observations can
be incorporated into the posterior by applying these updates recursively.
Our work is not nished, however, because the posterior covariance matrix depends
on the variational parameter  through () and we have yet to specify . We
choose  via an optimization procedure; in particular, we nd a value of  that yields
a tight lower bound in eq. (9). The fact that the variational expression in eq. (9) is
a lower bound is important|it allows us to use the EM algorithm to perform the
optimization. We derive such an EM algorithm in Appendix A; the result is the
following (closed form) update equation for :
post
where the expectation is taken with respect to P (jD;  old ), the variational posterior
distribution based on the previous value of . Owing to the EM formulation, each
update for  corresponds to a monotone improvement to the posterior approximation.
Empirically we nd that this procedure converges rapidly; only a few iterations are
needed. The accuracy of the approximation is considered in the following two sections.
To summarize, the variational approach allows us to obtain a closed form expression
for the posterior predictive distribution in logistic regression:
Z
where the posterior distribution P (jD) comes from making a single pass through
the data set applying the updates in eq. (10) and eq. (11) after
optimizing the associated variational parameters at each step. The predictive lower
bound P (S t jX t ; D) takes the form:
log
for any complete observation D t , where  and  signify the parameters in P (jD)
and the subscript t refers to the posterior P (jD; D t ) found by augmenting the data
set to include the point D t .
We note nally that the variational Bayesian calculations presented above need
not be carried out sequentially. We could compute a variational approximation to
the posterior probability P (jD) by introducing (separate) transformations for each
of the logistic functions in
Y
Y
The resulting variational parameters would have to be optimized jointly rather than
one at a time. We believe the sequential approach provides a cleaner solution.
3 Accuracy of the variational method
The logistic function is shown in Figure 2(a), along with a variational approximation
for 2. As we have noted, for each value of the variational parameter , there
is a particular point x where the approximation is exact; for the remaining values of
x the approximation is a lower bound.
a) -4 -2 0 2 40.20.61

Figure

2: a) The logistic function (solid line) and its variational form (dashed line)
for 2. b) The dierence between the predictive likelihood and its variational
approximation as a function of g( 0 ), as described in the text.
Integrating eq. (9) over the parameters we obtain a lower bound on the predictive
probability of an observation. The tightness of this lower bound is a measure of
accuracy of the approximation. To assess the variational approximation according
to this measure, we compared the lower bound to the true predictive likelihood that
was evaluated numerically. Note that for a single observation, the evaluation of the
predictive likelihood can be reduced to a one-dimensional integration problem:
Z
Z
where the eective prior P Gaussian with mean
where the actual prior distribution P () has mean  and covariance
. This reduction has no eect on the accuracy of the variational transformation and
thus it can be used in evaluating the overall accuracy. Figure 2(b) shows the dierence
between the true predictive probability and the variational lower bound for various
settings of the eective mean  0 and variance  2 , with  optimized separately for each
dierent values of  0 and  2 . The fact that the variational approximation is a lower
bound means that the dierence in the predictive likelihood is always positive.
We emphasize that the tightness of the lower bound is not the only relevant
measure of accuracy. Indeed, while a tight lower bound on the predictive probability
assures us that the associated posterior distribution is highly accurate, the converse is
not true in general. In other words, a poor lower bound does not necessarily imply a
poor approximation to the posterior distribution at the point of interest, only that we
no longer have any guarantees of good accuracy. In practice, we expect the accuracy of
the posterior to be more important than that of the predictive probability since errors
in the posterior run the risk of accumulating in the course of the sequential estimation
procedure. We defer the evaluation of the posterior accuracy to the following section
where comparisons are made to related methods.
4 Comparison to other methods
There are other sequential approximation methods that yield closed form posterior
parameter distributions in logistic regression models. The method most closely
related to ours is that of Spiegelhalter and Lauritzen (1990), which we refer to as
the S-L approximation in this paper. Their method is based on the Laplace ap-
proximation; that is, they utilize a local quadratic approximation to the complete
log-likelihood centered at the prior mean . The parameter updates that implement
this approximation are similar in spirit to the variational updates of eq. (10) and
eq.
post
post X (18)
X). Since there are no additional adjustable parameters in this
approximation, it is simpler than the variational method; however, we would expect
this lack of
exibility to translate into less accurate posterior estimates.
We compared the accuracy of the posterior estimates for the two methods in the
context of a single observation. To simplify the comparison we utilized the reduction
described in the previous section. Since the accuracy of neither method is aected
by this reduction, it su-ces for our purposes here to carry out the comparison in
this simpler setting. 1 The posterior probability of interest was therefore P ( 0 jD) /
computed for various choices of values for the prior mean  0 and the prior
standard deviation . The correct posterior mean and standard deviations were
obtained numerically. Figures 3 and 4 present the results. We plot signed dierences
in comparing the obtained posterior means to the correct ones; relative errors were
used for the posterior standard deviations. The error measures were left signed to
reveal any systematic biases. Note that the posterior mean from the variational
method is not guaranteed to be a lower bound on the true mean. Such guarantees
can be given only for the predictive likelihood. As can be seen in Figures 3(a) and 4(a)
the variational method yields signicantly more accurate estimates of the posterior
means, for both values of the prior variance. For the posterior variance, the S-L
estimate and the variational estimate appear to yield roughly comparable accuracy
for the small value of the prior variance (Figure 3(b)); however, for the larger prior
variance, the variational approximation is superior (Figure 4(b)). We note that the
variational method consistently underestimates the true posterior variance; a fact that
could be used to rene the approximation. Finally, in terms of the KL-divergences
between the approximate and true posteriors, the variational method and the S-L
approximation are roughly equivalent for the small prior variance; and again the
1 Note that the true posterior distribution over  can be always recovered from the posterior
computed for the one-dimensional reduced parameter
variational method is superior for the larger value of the prior variance. This is
shown in Figure 5.
a)
error
in
mean
S-L approximation
Variational
-0.020.02relative
error
in
stdv S-L approximation
Variational

Figure

3: a) The errors in the posterior means as a function of g( 0 ), where  0 is the
prior mean. Here for the prior. b) The relative errors in the posterior standard
deviations as a function of g( 0 ). Again for the prior distribution.
a)
-0.4
error
in
mean
S-L approximation
Variational
-0.2
-0.10.1relative
error
in
stdv S-L approximation
Variational

Figure

4: The plots are the same as in Figure 3, but now for the prior
distribution.
5 Extension to belief networks
A belief network is a probabilistic model over a set of variables fS i g that are
identied with the nodes in an acyclic directed graph. Letting (i) denote the set of
parents of node S i in the graph, we dene the joint distribution associated with the
belief network as the following product:
Y
KL-divergence S-L approximation
Variational
KL-divergence S-L approximation
Variational

Figure

5: KL-divergences between the approximate and the true posterior distribution
as a function of g( 0 ). a) for the prior. b) 3. The two approximation
methods have (visually) identical curves for
We refer to the conditional probabilities P (S i jS (i) ) as the \local probabilities" associated
with the belief network.
In this section we extend our earlier work in this paper and consider belief networks
in which logistic regression is used to dene the local probabilities (such models
have been studied in a non-Bayesian setting by Neal 1992 and by Saul, Jaakkola, &
Jordan 1994). Thus we introduce parameter vectors  i , one for each binary variable
consider models in which each local probability P (S i jS (i) ;  i ) is a logistic
regression of node S i on its parents S (i) .
To simplify the arguments in the following sections, we will consider augmented
belief networks in which the parameters themselves are treated as nodes in the belief
network (see Figure 6). This is a standard device in the belief network literature and
is of course natural within the Bayesian formalism.
5.1 Complete cases
A \complete case" refers to a data point in which all of the variables fS i g are
observed. If all of the data points are complete cases, then the methods that we
developed in the previous section apply immediately to belief networks. This can be
seen as follows. Consider the Markov blankets associated with each of the parameters

Figure

6(a)). For complete cases each of the nodes within the Markov blanket for
each of the parameters is observed (shaded in the diagram). By the independence
semantics of belief networks, this implies that the posterior distributions for the
parameters are independent of one another (conditioned on the observed data). Thus
the problem of estimating the posterior distributions for the parameters reduces to
a set of n independent subproblems, each of which is a Bayesian logistic regression
problem. We apply the methods developed in the previous sections directly.
SS
SS

Figure

a) A complete observation (shaded variables) and the Markov blanket
(dashed line) associated with the parameters  4 . b) An observation where the value
of S 4 is missing (unshaded in the gure).
5.2 Incomplete cases
The situation is substantially more complex when there are incomplete cases in
the data set. Incomplete cases imply that we no longer have all the Markov blankets
for the parameters in the network. Thus dependencies can arise between the parameter
distributions in dierent conditional models. Let us consider this situation in
some detail. A missing value implies that the observations arise from a marginal distribution
obtained by summing over the missing values of the unobserved variables.
The marginal distribution is thus a mixture distribution, where each mixture component
corresponds to a particular conguration of the missing variables. The weight
assigned to that component is essentially the posterior probability of the associated
conguration (Spiegelhalter & Lauritzen 1990). Note that the dependencies arising
from the missing values in the observations can make the network quite densely connected
(a missing value for a node eectively connects all of the neighboring nodes in
the graph). The dense connectivity leaves little structure to be exploited in the exact
probabilistic computations in these networks and tends to make exact probabilistic
calculations intractable.
Our approach to developing Bayesian methods for belief networks with missing
variables combines two variational techniques. In particular, we augment the -
transformation introduced earlier with a second variational transformation that we
refer to as a q-transformation. While the purpose of the -transformation is to convert
a local conditional probability into a form that can be integrated analytically, the
purpose of the q-transformation is to approximate the eect of marginalizing across
missing values associated with one or more parents. 2 Intuitively, the q-transformation
2 Treating the parameter as a parent node helps to emphasize the similarity between these two
variational transformations. The principal dierence is that a parameter node has only a single
\lls in" the missing values, allowing the variational transformation for complete data
to be invoked. The overall result is a closed-form approximation to the marginal
posterior.
The correct marginalization across missing variables is a global operation that
aects all of the conditional models that depend on the variables being marginalized
over. Under the variational approximation that we describe below, marginalization
is a local operation that acts individually on the relevant conditional models.
5.2.1 Approximate marginalization
Consider the problem of marginalizing over a set of variables S 0 under a joint
distribution:
Y
If we performed the marginalization exactly, then the resulting distribution would
not retain the same factorization as the original joint (assuming S 0 is involved in
more than one of the conditionals); this can be seen from:
Y
Y
where we have partitioned the product into the set of factors that depend on S 0
(indexed by i 0 ) and those that do not (indexed by i 00 ). Marginalization is not generally
a local operation on the individual node probabilities P
such locality, a desirable goal for computational reasons, can be achieved if we forgo
exact marginalization and instead consider approximations. In particular, we describe
a variational approximation that preserves locality at the expense of providing a lower
bound on the marginal probability instead of an exact result.
To obtain the desired variational transformation, we again exploit a convexity
property. In particular, for a given sequence consider the geometric
average
probability distribution. It is well known that the
geometric average is less than or equal to the arithmetic average
(This can be
easily established via an invocation of Jensen's inequality). We can exploit this fact
as follows. Consider an arbitrary distribution q(S 0 ), and rewrite the marginalization
operation in the following way:
child, while in general parents have multiple children.
Y
Y
where the inequality comes from transforming the average over the bracketed term
(with respect to the distribution q) into a geometric average. The third line follows
from plugging in the form of the joint distribution and exchanging the order of the
products. The logarithm of the multiplicative constant C(q) is the entropy of the
variational distribution q:
Y
and therefore log
Let us now make a few observations about the result in Eq. (24). First, note that
the lower bound in this equation has the same factored form as the original joint
probability. In particular, we dene the q-transformation of the ith local conditional
probability as follows:
~
Y
the lower bound in Eq. (24) is then a product of these q-transformations. Second,
note that all the conditionals are transformed by the same distribution q. A change in
q can thus aect all the transformed conditionals. This means that the dependencies
between variables S that would have resulted from exact marginalization over S 0 have
been replaced with \eective dependencies" through a shared variational distribution
q.
While the bound in Eq. (24) holds for an arbitrary variational distribution q(S 0 ),
to obtain a tight bound we need to optimize across q(S 0 ). In practice this involves
choosing a constrained class of distributions and optimizing across the class. The
simplest form of variational distribution is the completely factorized distribution:
Y
which yields a variational bound which is traditionally referred to as the \mean eld
approximation." This simplied approximation is appropriate in dense models with a
relatively large number of missing values. More generally, one can consider structured
variational distributions involving partial factorizations that correspond to tractable
substructures in the graphical model (cf. Saul & Jordan, 1996). We consider this
topic further in the following two sections.
Although the main constraint on the choice of q(S 0 ) is the computational one
associated with evaluation and optimization, there is one additional constraint that
must be borne in mind. In particular, the q-transformed conditional probabilities
must be in a form such that a subsequent -transformation can be invoked, yielding
as a result a tractable Bayesian integral. A simple way to meet this constraint is to
require that the variational distribution q(S 0 ) should not depend on the parameters .
As we discuss in the following section, in this case all of the q-transformations simply
involve products of logistic functions, which behave well under the -transformation.
5.2.2 Bayesian parameter updates
The derivation presented in the previous section shows that approximate variational
marginalization across a set of variables S 0 can be viewed as a geometric
average of the local conditional probabilities:
Y
where q(S 0 ) is the variational distribution over the missing values. Note that while
the -transformations are carried out separately for each relevant conditional model,
the variational distribution q associated with the missing values is the same across
all the q-transformations.
Given the transformation in eq. (28), the approximate Bayesian updates are obtained
readily. In particular, when conditioning on a data point that has missing
components we rst apply the q-transformation. This eectively lls in the missing
values, resulting in a transformed joint distribution that factorizes as in the case
of complete observations. The posterior parameter distributions therefore can be
obtained independently for the parameters associated with the transformed local
probabilities.
Two issues need to be considered. First, the transformed conditional probabilities
(cf. Eq. (28)) are products of logistic functions and therefore more complicated than
before. The -transformation method, however, transforms each logistic function
into an exponential with quadratic dependence on the parameters. Products of such
transforms are also exponential with quadratic dependence on the parameters. Thus
the approximate likelihood will again be Gaussian and if the prior is a multivariate
Gaussian the approximate posterior will also be Gaussian.
The second issue is the dependence of the posterior parameter distributions on the
variational distribution q. Once again we have to optimize the variational parameters
(a distribution in this case) to make our bounds as tight as possible; in particular, we
set q to the distribution that maximizes our lower bound. This optimization is carried
out in conjunction with the optimization of the  parameters for the transformations
of the logistic functions, which are also lower bounds. As we show in Appendix B.1,
the fact that all of our approximations are lower bounds implies that we can again
devise an EM algorithm to perform the maximization. The updates that are derived
in the Appendix are as follows:
pos i
where S (i) is the vector of parents of S i , and the expectations are taken with respect
to the variational distribution q.
5.2.3 Numerical evaluation
In this section, we provide a numerical evaluation of our proposed combination
of q-transformation and -transformation. We study a simple graph that consists
of a single node S and its parents S  . In contrast to the simple logistic regression
case analyzed earlier, the parents S  are not observed but instead are distributed
according to a distribution P (S  ). This distribution, which we manipulate directly
in our experiments, essentially provides a surrogate for the eects of a pattern of
evidence in the ancestral graph associated with node S (cf. Spiegelhalter & Lauritzen
1990).
Our interest is in the posterior probability over the parameters  associated with
the conditional probability P (SjS  ; ).
Suppose now that we observe 1. The exact posterior probability over the
parameters  in this case is given by
Our variational method focuses on lower bounding the evidence term in brackets. It
is natural to evaluate the overall accuracy of the approximation by evaluating the
accuracy of the marginal data likelihood:
We consider two dierent variational approximations. In the rst approximation
the variational distribution q is left unconstrained; in the second we use an approximation
that factorizes across the parents S  (the \mean eld" approximation). We
emphasize that in both cases the variational posterior approximation over the parameters
is a single Gaussian.
The results of our experiment are shown in Figures 7 and 8. Each gure displays
three curves, corresponding to the exact evaluation of the data likelihood P (D) and
the two variational lower bounds. The number of parents in S  was 5 and the prior
distribution P () was taken to be a zero mean Gaussian with a variable covariance
matrix. By the symmetry of the Gaussian distribution and the sigmoid function, the
exact value of P (D) was 0:5 in all cases. We considered several choices of P (S  )
and P (). In the rst case, the P (S  ) were assumed to factorize across the parents
and for each leaving a single parameter p that species the
stochasticity of P (S  ). A similar setting would arise when applying the mean eld
approximation in the context of a more general graph. Figure 7 shows the accuracy of
the variational lower bounds as a function of p where in Figure 7(a) P
i.e., the covariance matrix is diagonal with diagonal components set to 1=5, and in

Figure

is a sample covariance matrix of 5 Gaussian
random vectors distributed according to N(0; I=5). The results of Figure 7(b) are
averaged over 5 independent runs. The choice of scaling in N(0; I=5) is made to
insure that j
gures indicate that the variational approximations
are reasonably accurate and that there is little dierence between the two methods.
In

Figure

8 we see how the mean eld approximation (which is unimodal) deteriorates
as the distribution P (S  ) changes from a factorized distribution toward
a mixture distribution. More specically, let P f (S  jp) be the (uniform) factorized
distribution discussed above with parameter p and let Pm (S  ) be a pure mixture
distribution that assigns a probability mass 1=3 to three dierent (randomly chosen)
congurations of the parents S  . We let P (S
where the parameter p m controls the extent to which P (S  ) resembles a (pure) mixture
distribution. Figure 8 illustrates the accuracy of the two variational methods
as a function of p m where in Figure 8(a) As expected,
the mean eld approximation deteriorates with an increasing p m whereas our rst
variational approximation remains accurate.
6 The dual problem
In the logistic regression formulation (eq. (1)), the parameters  and the explanatory
variables X play a dual or symmetric role (cf. Nadal and Parga 1994). In
the Bayesian logistic regression setting, the symmetry is broken by associating the
same parameter vector  with multiple occurences of the explanatory variables X
as shown in Figure 9. Alternatively, we may break the symmetry by associating a
single instance of the explanatory variable X with multiple realizations of . In this
sense the explanatory variables X play the role of parameters while  functions as a
continuous latent variable. The dual of the Bayesian regression model is thus a latent
a)
likelihood
likelihood

Figure

7: Exact data likelihood (solid line), variational lower bound 1 (dashed line),
and variational lower 2 (dotted line) as a function of the stochasticity parameter p
of P (S  ). In (a) P and in (b) P is a sample
covariance of 5 random vectors distributed according to N(0; I=5).
a)
likelihood
likelihood

Figure

8: Exact data likelihood (solid line) and the two variational lower bounds
(dashed and dotted lines respectively) as a function of the mixture parameter p m . In
0:1 and in (b)
variable density model over a binary response variable S. Graphically, in the dual
interpretation we have a single \parameter" node for X whereas separate nodes are
required for dierent realizations of  (illustrated as  (i) in the gure) to explain successive
observations S (i) . While a latent variable density model over a single binary
variable is not particularly interesting, we can generalize the response variable S to a
vector of binary variables each component S i has a distinct
set of \parameters" associated with it. The latent variables ,
however, remain in this dual interpretation the same for all S i . We note that strictly
speaking the dual interpretation would require us to assign a prior distribution over
the new \parameters" vectors X i . For simplicity, however, we omit this consideration
and treat X i simply as adjustable parameters. The resulting latent variable
(1) (2) (3)
(1) (2) (3)
(1) (2) (3)
x
(1) (2) (3)
a) b)

Figure

9: a) Bayesian regression problem. b) The dual problem.
density model over binary vectors is akin to the standard factor analysis model (see
e.g. Everitt 1984). This model has already been used to facilitate visualization of
high dimensional binary vectors (Tipping 1999).
We now turn to a more technical treatment of this latent variable model. The
joint distribution is given by
where the conditional probabilities for the binary observables are logistic regression
models
We would like to use the EM- algorithm for parameter estimation. To achieve
this we again exploit the variational transformations. The transformations can be
introduced for each of the conditional probability in the joint distribution and optimized
separately for every observation D
n g in the database consisting
only of the values of the binary output variables. As in the logistic regression case,
the transformations change the unwieldy conditional models into simpler ones that
depend on the parameters only quadratically in the exponent. The variational ev-
idence, which is a product of the transformed conditional probabilities, retains the
same property. Consequently, under the variational approximation, we can compute
the posterior distribution over the latent variables  in closed form. The mean and
the covariance of this posterior can be obtained analogously to the regression case
giving
The variational parameters  t
associated with each observation and the conditional
model can be updated using eq. (12) where X is replaced with X i , now the vector of
parameters associated with the i th conditional model.
We can solve the M-step of the EM-algorithm by accumulating su-cient statistics
for the parameters X based on the closed form posterior distributions corresponding
to the observations in the data set. Omitting the algebra, we obtain the
following explicit updates for the parameters:
where
and the subscript t denotes the quantities pertaining to the observation D t . Note that
since the variational transformations that we expoited to arrive at these updates
are all lower bounds, the M-step necessarily results in a monotonically increasing
lower bound on the log-probability of the observations. This desirable monotonicity
property is unlikely to arise with other types of approximation methods, such as the
Laplace approximation.
We have exemplied the use of variational techniques in the setting of Bayesian
parameter estimation. We found that variational methods can be exploited to yield
closed form expressions that approximate the posterior distributions for the parameters
in logistic regression. The methods apply immediately to a Bayesian treatment
of logistic belief networks with complete data. We also showed how to combine mean
eld theory with our variational transformation and thereby treat belief networks
with missing data. Finally, our variational techniques lead to an exactly solvable
EM algorithm for a latent variable density model|the dual of the logistic regression
problem.
It is also of interest to note that our variational method provides an alternative
to the standard iterative Newton-Raphson method for maximum likelihood estimation
in logistic regression (an algorithm known as \iterative reweighted least
squares" or \IRLS"). The advantage of the variational approach is that it guarantees
monotone improvement in likelihood. We present the derivation of this algorithm in


Appendix

C.
Finally, for an alternative perspective on the application of variational methods
to Bayesian inference, see Hinton and van Camp (1993) and MacKay (1997). These
authors have developed a variational method known as \ensemble learning," which
can be viewed as a mean eld approximation to the marginal likelihood.



--R

Bayesian Theory.
Learning Bayesian networks: the combination of knowledge and statistical data.
An Introduction to Latent Variable Models.
Bayesian Data Analysis.
Markov Chain Monte Carlo in Practice.
Keeping neural networks simple by minimizing the description length of the weights.
An introduction to variational methods in graphical models.
Generalized Linear Models.
Ensemble learning for hidden Markov models.
Duality between learning machines: A bridge between supervised and unsupervised learning.

Technical report CRG-TR-93-1

Convex Analysis.
Variational Methods in Statistics.
Exploiting tractable substructures in intractable net- works

Sequential updating of conditional probabilities on directed graphical structures.
Probabilistic visualisation of high-dimensional binary data
BUGS: A program to perform Bayesian inference using Gibbs sampling.
A Optimization of the variational parameters To optimize the variational approximation of eq.
The form of this posterior, however, remains at least as unwieldy as the Bayesian logistic regression problem considered earlier in the paper. Proceeding analogously, we transform the logistic functions as in Eq. (7) corresponding to each of the conditional probabilities in the product and obtain P

The optimization of these parameters is shown in Appendix B.

The metric for optimizing the parameters comes from the fact that the transformations associated with these parameters introduce a lower bound on the probability of the observations.
Similarly to the case of the simple Bayesian logistic regression considered previously (see
--TR

--CTR
Edward Snelson , Zoubin Ghahramani, Compact approximations to Bayesian predictive distributions, Proceedings of the 22nd international conference on Machine learning, p.840-847, August 07-11, 2005, Bonn, Germany
R. Ksantini , D. Ziou , B. Colin , F. Dubeau, Weighted pseudo-metric for a fast CBIR method, Machine Graphics & Vision International Journal, v.15 n.3, p.471-480, January 2006
Aaron D'Souza , Sethu Vijayakumar , Stefan Schaal, The Bayesian backfitting relevance vector machine, Proceedings of the twenty-first international conference on Machine learning, p.31, July 04-08, 2004, Banff, Alberta, Canada
Dmitry Kropotov , Dmitry Vetrov, On one method of non-diagonal regularization in sparse Bayesian learning, Proceedings of the 24th international conference on Machine learning, p.457-464, June 20-24, 2007, Corvalis, Oregon
Wang , D. M. Titterington, Lack of Consistency of Mean Field and Variational break Bayes Approximations for State Space Models, Neural Processing Letters, v.20 n.3, p.151-170, November  2004
Stephen Roberts , Evangelos Roussos , Rizwan Choudrey, Hierarchy, priors and wavelets: structure and signal modelling using ICA, Signal Processing, v.84 n.2, p.283-297, February 2004
Shinichi Nakajima , Sumio Watanabe, Variational Bayes Solution of Linear Neural Networks and Its Generalization Performance, Neural Computation, v.19 n.4, p.1112-1153, April 2007
R. A. Choudrey , S. J. Roberts, Variational mixture of Bayesian independent component analyzers, Neural Computation, v.15 n.1, p.213-252, January
Perlich , Foster Provost , Jeffrey S. Simonoff, Tree induction vs. logistic regression: a learning-curve analysis, The Journal of Machine Learning Research, 4, p.211-255, 12/1/2003
variational methods;logistic regression;graphical models;incomplete data;belief networks;bayesian estimation