Kullback-Leibler Divergence
for the Normal-Gamma Distribution
Joram Soch1,3,• & Carsten Allefeld1,2
1
2

arXiv:1611.01437v1 [] 4 Nov 2016

3

Bernstein Center for Computational Neuroscience, Berlin, Germany
Berlin Center for Advanced Neuroimaging, Berlin, Germany
Department of Psychology, Humboldt-Universität zu Berlin, Germany

• Corresponding author: joram.soch@bccn-berlin.de.

Abstract
We derive the Kullback-Leibler divergence for the normal-gamma distribution and
show that it is identical to the Bayesian complexity penalty for the univariate
general linear model with conjugate priors. Based on this finding, we provide two
applications of the KL divergence, one in simulated and one in empirical data.

Contents
1 Introduction

1

2 Theory
2.1 Multivariate normal KL divergence
2.2 Univariate gamma KL divergence .
2.3 Normal-gamma KL divergence . . .
2.4 The Bayesian model evidence . . .
2.5 The general linear model . . . . . .

.
.
.
.
.

2
2
2
3
4
5

3 Application
3.1 Polynomial basis functions . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Neuroimaging model selection . . . . . . . . . . . . . . . . . . . . . . . .

6
6
8

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

4 Conclusion

10

5 References

10

0

1

Introduction

Let x be a k × 1 random vector and y > 0 be a random variable. Then, x and y are
said to follow a normal-gamma distribution (NG distribution), if their joint probability
density function is given by
p(x, y) = N(x; µ, (yΛ)−1 ) · Gam(y; a, b)

(1)

where N(x; µ, Σ) denotes a multivariate normal density with mean µ and covariance Σ
and Gam(x; a, b) denotes a gamma density with shape a and rate b. In full, the density
function is given by (Koch, 2007, p. 55)
s
i ba
h y
|yΛ|
T
(x
−
µ)
Λ(x
−
µ)
·
y a−1 exp[−by] .
(2)
p(x, y) =
exp
−
(2π)k
2
Γ(a)
The Kullback-Leibler divergence (KL divergence) is a non-symmetric distance measure
for two probability distributions P and Q and is defined as
KL[P ||Q] =

X
i∈Ω

P (i) ln

P (i)
.
Q(i)

(3)

For continuous probability distribtions P and Q with probability density functions p(x)
and q(x) on the same domain X, it is given by (Bishop, 2006, p. 55)
Z
p(x)
dx .
(4)
KL[P ||Q] =
p(x) ln
q(x)
X
The KL divergence becomes important in information theory and statistical inference.
Here, we derive the KL divergence for two NG distributions with vector-valued x and
real-positive y and provide two examples of its application.

1

2
2.1

Theory
Multivariate normal KL divergence

First, consider two multivariate normal distributions over the k × 1 vector x specified by
p(x) = N(x; µ1 , Σ1 )
q(x) = N(x; µ2 , Σ2 )
According to equation (4), the KL divergence of P from Q is defined as
Z
N(x; µ1 , Σ1 )
dx .
KL[P ||Q] =
N(x; µ1 , Σ1 ) ln
N(x; µ2 , Σ2 )
Rk
Using the multivariate normal density function


1
1
T −1
N(x; µ, Σ) = p
exp − (x − µ) Σ (x − µ) ,
2
(2π)n |Σ|
it evaluates to (Duchi, 2014)


1
|Σ1 |
T −1
−1
(µ2 − µ1 ) Σ2 (µ2 − µ1 ) + tr(Σ2 Σ1 ) − ln
−k .
KL[P ||Q] =
2
|Σ2 |

2.2

(5)

(6)

(7)

(8)

Univariate gamma KL divergence

Next, consider two univariate gamma distributions over the real-positive y specified by
p(y) = Gam(y; a1 , b1 )
q(y) = Gam(y; a2 , b2 )
According to equation (4), the KL divergence of P from Q is defined as
Z ∞
Gam(y; a1 , b1 )
dy .
KL[P ||Q] =
Gam(y; a1 , b1 ) ln
Gam(y; a2 , b2 )
0

(9)

(10)

Using the univariate gamma density function
Gam(y; a, b) =

ba a−1
y
exp[−by] for y > 0 ,
Γ(a)

(11)

it evaluates to (Penny, 2001)
KL[P ||Q] = a2 ln

b1
Γ(a1 )
a1
− ln
+ (a1 − a2 ) ψ(a1 ) − (b1 − b2 )
b2
Γ(a2 )
b1

where ψ(x) is the digamma function.

2

(12)

2.3

Normal-gamma KL divergence

Now, consider two normal-gamma distributions over x and y specified by
p(x, y) = N(x; µ1 , (yΛ1 )−1 ) Gam(y; a1 , b1 )
q(x, y) = N(x; µ2 , (yΛ2 )−1 ) Gam(y; a2 , b2 )
According to equation (4), the KL divergence of P from Q is defined as
Z ∞Z
p(x, y)
p(x, y) ln
KL[P ||Q] =
dx dy .
q(x, y)
0
Rk
Using the law of conditional probability, it can be evaluated as follows:
Z

∞

(13)

(14)

Z

p(x|y)p(y)
p(x|y)p(y) ln
dx dy
q(x|y)q(y)
0
Rk
Z
Z ∞
p(x|y)
p(y)
p(x|y) ln
=
dx dy
q(x|y)
0
Rk
Z ∞
Z
p(y)
+
p(y) ln
p(x|y) dx dy
q(y) Rk
0
= hKL[p(x|y)||q(x|y)]ip(y) + KL[p(y)||q(y)]

KL[P ||Q] =

(15)

In other words, the KL divergence for two normal-gamma distributions over x and y
is equal to the sum of a multivariate normal KL divergence regarding x conditional on
y, expected over y, and a univariate gamma KL divergence regarding y. Together with
equation (8), the first term becomes
hKL[p(x|y)||q(x|y)]ip(y)

 

|(yΛ1 )−1 |
1
T
−1
(µ2 − µ1 ) (yΛ2 )(µ2 − µ1 ) + tr (yΛ2 )(yΛ1 )
−k
− ln
=
2
|(yΛ2 )−1 |
p(y)


y
1
1 |Λ2 | k
=
(µ2 − µ1 )T Λ2 (µ2 − µ1 ) + tr(Λ2 Λ−1
ln
−
.
1 )−
2
2
2 |Λ1 | 2 p(y)

(16)

Using the relation y ∼ Gam(a, b) ⇒ hyi = a/b, we have
hKL[p(x|y)||q(x|y)]ip(y) =

1
1 |Λ2 | k
1 a1
(µ2 − µ1 )T Λ2 (µ2 − µ1 ) + tr(Λ2 Λ−1
ln
− .
1 )−
2 b1
2
2 |Λ1 | 2
(17)

Thus, from (17) and (12), the KL divergence in (15) becomes
 1
1 a1 
1 |Λ2 | k
(µ2 − µ1 )T Λ2 (µ2 − µ1 ) + tr(Λ2 Λ−1
ln
−
1 )−
2 b1
2
2 |Λ1 | 2
b1
a1
Γ(a1 )
+ a2 ln − ln
+ (a1 − a2 ) ψ(a1 ) − (b1 − b2 )
.
b2
Γ(a2 )
b1

KL[P ||Q] =

3

(18)

2.4

The Bayesian model evidence

Consider Bayesian inference on data y using model m with parameters θ. In this case,
Bayes’ theorem is a statement about the posterior density:
p(θ|y, m) =

p(y|θ, m) p(θ|m)
.
p(y|m)

(19)

The denominator p(y|m) acts as a normalization constant on the posterior density p(θ|y, m)
and according to the law of marginal probability is given by
Z
p(y|m) = p(y|θ, m) p(θ|m) dθ .
(20)
This is the probability of the data given only the model, regardless of any particular
parameter values. It is also called “marginal likelihood” or “model evidence” and can act
as a model quality criterion in Bayesian inference, because parameters are integrated out
of the likelihood.
For computational reasons, only the logarithmized or log model evidence (LME) L(m) =
ln p(y|m) is of interest in most cases. By rearranging equation (19), the model evidence
can be represented as
p(y|m) =

p(y|θ, m) p(θ|m)
.
p(θ|y, m)

(21)

Logarithmizing both sides of the equation and taking the expectation with respect to the
posterior density over model parameters θ gives the LME
Z
Z
p(θ|y, m)
dθ .
(22)
L(m) = p(θ|y, m) ln p(y|θ, m) dθ − p(θ|y, m) ln
p(θ|m)
Using this reformulation, the LME as a model quality measure can be naturally decomposed into an accuracy term, the posterior expected likelihood, and a complexity term,
the KL divergence between the posterior and the prior distribution:
L(m) = Acc(m) − Com(m)
Acc(m) = hlog p(y|θ, m)ip(θ|y,m)

(23)

Com(m) = KL [p(θ|y, m)||p(θ|m)]
Intuitively, the accuracy acts increasing and the complexity acts decreasing on the log
model evidence. This reflects the capability of the LME to select models that achieve the
best balance between accuracy and complexity, i.e. models that explain the observations
sufficiently well (high accuracy) without employing too many principles (low complexity).
The fact that the complexity term is a KL divergence between posterior and prior means
that models with prior assumptions that are close to the posterior evidence receive a low
complexity penalty, because one is not surprised very much when accepting such a model
which renders the Bayesian complexity a measure of surprise.

4

2.5

The general linear model

Consider multiple linear regression using the univariate general linear model (GLM)
y = Xβ + ε, ε ∼ N (0, σ 2 V )

(24)

where y is an n × 1 vector of measured data, X is an n × p matrix called the design
matrix, β is a p × 1 vector of weight parameters called regression coefficients and ε is
an n × 1 vector of errors or noise. These residuals are assumed to follow a multivariate
normal distribution whose covariance matrix is the product of a variance factor σ 2 and
an n × n correlation matrix V . Usually, X and V are known while β and τ are unknown
parameters to be inferred via model estimation.
For mathematical convenience, we rewrite σ 2 = 1/τ and V = P −1 so that equation (24)
implies the following likelihood function:
p(y|β, τ ) = N(y; Xβ, (τ P )−1 ) .

(25)

The conjugate prior relative to this likelihood function is a normal-gamma distribution
on the model parameters β and τ (Koch, 2007, ch. 2.6.3):
p(β|τ ) = N(β; µ0 , (τ Λ0 )−1 )
p(τ ) = Gam(τ ; a0 , b0 )

(26)

Due to the conjugacy of (26) to (25), the posterior is also a normal-gamma distribution
p(β|τ, y) = N(β; µn , (τ Λn )−1 )
p(τ |y) = Gam(τ ; an , bn )

(27)

where the posterior parameters in (27) are given by (Koch, 2007, ch. 4.3.2)
T
µn = Λ−1
n (X P y + Λ0 µ0 )

Λn = X T P X + Λ 0
n
an = a0 +
2
1 T
bn = b0 + (y P y + µT0 Λ0 µ0 − µTn Λn µn )
2
From (23), the complexity for the model defined by (25) and (26) is given by
Com(m) = KL [p(β, τ |y)||p(β, τ )] .

(28)

(29)

In other words, the complexity penalty for a general linear model with normal-gamma
priors (GLM-NG) is identical to a KL divergence between two NG distributions and using
(18) can be written in terms of the prior and posterior parameters as
 1
1 an 
1 |Λ0 | p
(µn − µ0 )T Λ0 (µn − µ0 ) + tr(Λ0 Λ−1
ln
−
n )−
2 bn
2
2 |Λn | 2
Γ(an )
an
bn
+ a0 ln − ln
+ (an − a0 ) ψ(an ) − (bn − b0 )
.
b0
Γ(a0 )
bn

Com(m) =

5

(30)

3
3.1

Application
Polynomial basis functions

Consider a linear model with polynomial basis functions (Bishop, 2006, p. 5) given by
y=

p
X

ci xi + ε with

− 1 ≤ x ≤ +1 .

(31)

i=0

Essentially, this model assumes that y is an additive mixture of polynomial terms xi
weighted with the coefficients ci with i = 1, . . . , p where the natural number p is called
the model order. This means that p = 0 corresponds to a constant value (plus noise ε);
p = 1 corresponds to a linear function; p = 2 corresponds to a quadratic pattern; p = 3
corresponds to a 3rd degree polynomial etc.
Given that x is an n × 1 vector of real numbers between −1 and +1, this model can be
rewritten as a GLM given in equation (24) with

 

c0
x01 x11 · · · xp1
 c1 
 x0 x1 · · · xp 
2
2
 
 2
(32)
X =  ..
.. . .
..  and β =  ..  .
.
.
. .
.
cp
x0n x1n · · · xpn
Based on this reformulation, we simulate polynomial data. We perform N = 100 simulations with n = 100 data points in each simulation. We generate simulated data based on
a true model order ptrue = 5 and analyze these data using a set of models ranging from
pmin = 0 to pmax = 20. The predictor x is equally spaced between −1 and +1 and design
matrices Xp are created according to equation (32).
In each simulation, six regression coefficients βtrue are drawn independently from the standard normal distribution N(0, 1). Then, Gaussian noise ε is sampled from the multivariate
normal distribution N(0, σε2 In ) with a residual variance of σε2 = 1. Finally, simulated data
is generated as as y = X5 βtrue + ε.
Then, for each p ∈ {0, . . . , 20}, Bayesian model estimation is performed using the design
matrix Xp , a correlation matrix V = In and the prior distributions (26) with the prior
parameters µ0 = 0p , Λ0 = Ip invoking a standard multivariate normal distribution and
a0 = 1, b0 = 1 invoking a relatively flat gamma prior. Posterior parameters are calculated
using equation (28) and give rise to the model complexity via (30) as well as model
accuracy and the log model evidence via (23).
Average LMEs, accuracies and complexities are shown in Figure 1. One can see that
the true model order is correctly identified by the maximal log model evidence. This is
achieved by an increasing complexity penalty which outweighs the saturating accuracy
gain for models with p > 5. This demonstrates that the KL divergence for the NG
distribution can be used to select polynomial basis functions when basis sets cannot be
separated based on model accuracy alone.

6

Log Model Evidence
−150

L(m)

−155
−160
−165
−170
0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

14

15

16

17

18

19

20

14

15

16

17

18

19

20

p

Model Accuracy
−140

Acc(m)

−145
−150
−155
−160
−165
0

1

2

3

4

5

6

7

8

9

10

11

12

13

p

Model Complexity
10

Com(m)

9
8
7
6
5
4
0

1

2

3

4

5

6

7

8

9

10

11

12

13

p

Figure 1. Bayesian model selection for polynomial basis functions. All displays have
model order on the x-axis and average model quality measures (across simulations) on
the y-axis. Intuitively, the model accuracy (middle panel) increases with model order,
but saturates at around −140 with no major increase after p = 5. Moreover, the model
complexity (lower panel) – which is the KL divergence between prior and posterior distribution – also grows with model order, but switches to a linear increase at around
p = 5 reaching a value of 10 at p = 20. Together, this has the consequence that the log
model evidence (upper panel) is maximal for p = 5 (black cross) where exact values are:
Acc(m) = −140.77, Com(m) = 8.34, L(m) = −149.11.

7

3.2

Neuroimaging model selection

In neuroimaging, especially functional magnetic resonance imaging (fMRI), GLMs as
given by equation (24) are applied to time series of neural data y (Friston et al., 1995). The
design matrix X is specified by the temporal occurrence of experimental conditions and
the covariance structure V is estimated from residual auto-correlations. Model estimation
and statistical inference are performed “voxel-wise”, i.e. separately for each measurement
location in the brain, usually referred to as the “mass-univariate GLM”.
Here, we analyze data from a study on orientation pop-out processing (Bogler et al., 2013).
During the experimental paradigm, the screen showed a 3 × 7 array of homogeneous bars
oriented either 0°, 45°, 90° or 135° relative to the vertical axis. This background stimulation
changed every second and was interrupted by trials in which one target bar on the left and
one target bar on the right were independently rotated either 0°, 30°, 60° or 90° relative to
the rest of the stimulus display. Those trials of orientation contrast (OC) lasted 4 seconds
and were alternated with inter-trial intervals of 7, 10 or 13 seconds. Each combination of
OC on the left side and OC on the right side was presented three times resulting in 48
trials in each of the 5 sessions lasting 672 seconds.
After fMRI data preprocessing (slice-timing, realignment, normalization, smoothing), two
different models of hemodynamic activation were applied to the fMRI data. The first
model (GLM I) considers the experiment a factorial design with two factors (left vs. right
OC) having four levels (0°, 30°, 60°, 90°). This results in 4 × 4 = 16 possible combinations
or experimental conditions modelled by 16 onset regressors convolved with the canonical
hemodynamic response function (HRF). The second model (GLM II) puts all trials from
all conditions into one HRF-convolved regressor and encodes orientation contrast using a
parametric modulator (PM) that is given as PM = deg/90° with deg = (0°, 30°, 60°, 90°),
resulting in PM = (0, 1/3, 2/3, 1), such that the parametric modulator is proportional to
orientation contrast. There was one PM for each factor of the design, i.e. one PM for left
OC and one PM for right OC. Note that both models encode the same information and
that every signal that can be identified using GLM II can also be characterized using
GLM I, but not vice versa, because the first model allows for a greater flexibility of
activation patterns across experimental conditions than the second.
For these two models, we performed Bayesian model estimation. To overcome the challenge of having to specify prior distributions on the model parameters, we apply crossvalidation across fMRI sessions. This gives rise to a cross-validated log model evidence
(cvLME) as well as cross-validated accuracies and complexities for each model in each
subject. We then performed a paired t-test to find voxels where GLM II has a significantly
higher cvLME than GLM I. Due to the specific assumptions in GLM II and the higher
flexibility of GLM I, we assumed that these differences might be primarily based on a
complexity advantage of GLM II over GLM I.
We focus on visual area 4 (V4) that is known to be sensitive to orientation contrast.
Within left V4, specified by a mask from a separate localizer paradigm (Bogler et al.,
2013), we identified the peak voxel ([x y z] = [−15, −73, −5] mm) defined by the maximal
t-value (t = 3.95) and extracted log model evidence as well as model accuracy and model
complexity from this voxel for each subject. Differences in LME, accuracy and complexity
are shown in Figure 2. Again, the model complexity enables a model selection that would
not be possible based on the model accuracy alone.

8

Log Model Evidence
L(m2) − L(m1)

10

5

0

−5

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 avg

subject

Model Accuracy
Acc(m2) − Acc(m1)

4
2
0
−2
−4
−6
−8
−10
1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 avg

subject

Model Complexity
Com(m2) − Com(m1)

0
−2
−4
−6
−8
−10
1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 avg

subject

Figure 2. Bayesian model selection for orientation pop-out processing. All displays have
subject on the x-axis and difference in model qualities (m1 = GLM I, m2 = GLM II)
on the y-axis. Interestingly, there is a slight disadvantage for GLM II regarding only the
model accuracy (middle panel), its mean difference across subjects being smaller than
zero. However, model complexity (lower panel) – measured as the KL divergence between
prior and posterior distribution – is consistently higher for GLM I. Together, this has
the consequence that the log model evidence (upper panel) most often favors GLM II.
Average values are: ∆Acc = −1.78, ∆Com = −6.02, ∆L = 4.24.

9

4

Conclusion

We have derived the Kullback-Leibler divergence of two normal-gamma distributions
using earlier results on the KL divergence for multivariate normal and univariate gamma
distributions. Moreover, we have shown that the KL divergence for the NG distribution
occurs as the complexity term in the univariate general linear model when using conjugate
priors.
Analysis of simulated and empirical data demonstrates that the complexity penalty has
the desired theoretical features, namely to quantify the relative informational content of
two generative models and to detect model differences that cannot be detected by just
relying on model accuracy, e.g. given by the maximum log-likelihood (as in information
criteria like AIC or BIC) or the posterior log-likelihood (as in the Bayesian log model
evidence).

5

References

[1] Bishop CM (2006): Pattern Recognition and Machine Learning. 1st Edition, Springer,
New York.
[2] Bogler C, Bode S, Haynes JD (2013): “Orientation pop-out processing in human visual
cortex”. NeuroImage, vol. 81, pp. 73-80.
[3] Duchi J (2014): “Derivations for Linear Algebra and Optimization”. University of
California, Berkeley; URL: http://www.eecs.berkeley.edu/∼jduchi/projects/general
notes.pdf.
[4] Friston KJ, Holmes AP, Worsley KJ, Poline JP, Frith CD, Frackowiak RSJ (1995):
“Statistical parametric maps in functional imaging: A general linear approach”. Human Brain Mapping, vol. 2, iss. 4, pp. 189-210.
[5] Koch KR (2007): Introduction to Bayesian Statistics. 2nd Edition, Springer, Berlin/
Heidelberg, 2007.
[6] Penny WD (2001): “KL-Divergences of Normal, Gamma, Dirichlet and Wishart
densities”. University College, London; URL: www.fil.ion.ucl.ac.uk/∼wpenny/
publications/densities.ps.

10

