arXiv:1703.03353v1 [math.ST] 9 Mar 2017

A Note on Bayesian Model Selection for Discrete Data Using
Proper Scoring Rules
A. Philip Dawid∗

Monica Musio†

Silvia Columbu‡

March 10, 2017

Abstract
We consider the problem of choosing between parametric models for a discrete observable,
taking a Bayesian approach in which the within-model prior distributions are allowed to be
improper. In order to avoid the ambiguity in the marginal likelihood function in such a
case, we apply a homogeneous scoring rule. For the particular case of distinguishing between
Poisson and Negative Binomial models, we conduct simulations that indicate that, applied
prequentially, the method will consistently select the true model.
Keywords: Consistent model selection; homogeneous score; discrete data; prequential

1

Introduction

It is well known that Bayesian model selection with improper within-model prior distributions
is not well-defined, owing to the presence of an arbitrary multiplicative constant in each term
of the marginal likelihood function. Recently (Dawid and Musio, 2015) it has been shown how
this problem can be overcome if one replaces negative log-likelihood (the log score) by another,
homogeneous, proper scoring rule (Parry et al., 2012)—since then the arbitrary constants do not
enter into the formulae. That paper considered the case of continuous variables and, in particular,
the Hyvärinen scoring rule (Hyvärinen, 2005), and showed that this approach will generally lead
to consistent selection of the correct model.
The above approach can not be applied directly when the data are discrete, since then we
need to use scoring rules specifically adapted to the discrete case, as characterised in Dawid et al.
(2012). Here we investigate, by example, such a discrete data problem. In particular we consider the problem of distinguishing between the Poisson and the Negative Binomial distributions.
Simulations indicate that the method will again deliver consistent selection of the true model.

2

Local scoring rules

Let X be a discrete sample space endowed with a structure whereby with each x ∈ X is associated
a neighbourhood Nx ⊆ X , containing x. In Dawid et al. (2012) it was shown how to define a proper
local scoring rule S(x, P ) on X , where x ∈ X , and P is a distribution over X . The rule is proper if,
for all P , S(P, Q) := EX∼P S(X, Q) is minimised for Q = P , and local if S(x, P ) depend on P only
through the probabilities it assigns to points in Nx . Under a condition on the neighbourhoods,
we can define an undirected graph G on X such that we can take y ∈ Nx just when x and y
are identical or are adjacent in G. Then all proper local scorings can be characterised, and (on
excluding the log score, yielding what are termed key local proper scoring rules) any of these will
∗ Department

of Pure Mathematics and Mathematical Statistics, University of Cambridge, U.K.
of Mathematics, University of Cagliari, Italy
‡ Dipartiment of Mathematics, University of Cagliari, Italy
† Dipartiment

1

be homogeneous in the sense that its value is unchanged when all probabilities in Nx are scaled
by the same positive constant.
In particular, suppose the sample space X is the set Z≥0 of non-negative integers, and we
regard x and y as neighbours if and only if they differ by at most 1. It is shown in Dawid et al.
(2012) that a key local scoring rule adapted to this structure has the form






p(x)
p(x + 1)
p(x + 1) 0 p(x + 1)
0
S(x, P ) = Gx−1
+ Gx
−
Gx
(x = 0, 1, . . .) (1)
p(x − 1)
p(x)
p(x)
p(x)
where, for each x ∈ Z≥0 , p(x) = P (X = x), Gx is a concave function on R+ , and the first term in
(1) is absent if x = 0. It is clear from the way in which ratios enter (1) that such a scoring rule is
homogeneous.
The cumulative score (1) based on an independent and identically distributed sample (x1 , . . . , xn )
in which the frequency of y is fy (y = 0, 1, . . .) is
∞
X

fy Gy (vy ) + (fy+1 − fy vy ) G0y (vy )

(2)

y=0

with vy := p(y + 1)/p(y). If for example we wished to fit the Poisson model p(x) ∝ θx /x!, we
might estimate θ by minimising the total empirical score
 




∞
X
fy
θ
θ
0
+ fy+1 −
θ Gy
.
(3)
fy Gy
y+1
y+1
y+1
y=0
In the sequel we shall use the special case of (1) with
Gx (v) = −(x + 1)a v m /m(m − 1)

(m > 0, m 6= 1).

This gives the scoring rule
 −1
m
m {p(1)/p(0)}




m
S(x, P ) =
{m(m − 1)}−1 [(m − 1)(x + 1)a {p(x + 1)/p(x)}

i


m−1

−mxa {p(x)/p(x − 1)})

3

(4)

(x = 0)
(5)
(x > 0).

Bayesian Model Selection

Let M be a finite or countable class of statistical models for the same observable X ∈ X . Each
M ∈ M is a parametric family, with parameter θM ∈ ΘM , a dM -dimensional Euclidean space;
when M obtains, with parameter value θM , then X has distribution PθM , with density function
(probability mass function) pM (x | θM ). Having observed data X = x, we wish to make inference
about which model M ∈ M (and possibly which parameter-value θM ) actually generated the data.
The Bayesian approach assigns, within each model M , a prior distribution ΠM , with density
πM (·) say, for its parameter θM . The associated predictive distribution PM of X (given only the
validity of model M , but no information on its parameter) has density function
Z
pM (x) =
pM (x | θM ) πM (θM ) dθM .
(6)
ΘM

Any function over M proportional to pM (x) (considered as a function of M , for fixed x) supplies
the marginal likelihood function, L(M ), based on data X = x. In typical asymptotic scenarios,
selection of the model maximising L(M ), or, equivalently, minimising the log score SL (x, PM ) :=
− log pM (x), will consistently select the true model (Dawid, 2011).
“Objective” Bayesian inference attempts to use standardised within-model priors ΠM intended
to represent “prior ignorance”. In many applications, such an “ignorance prior” for θM is not a
2

genuine distribution, but rather an “improper” σ-finite but not finite measure, with a “density”
πM (·) that does not have a finite integral and so can not be normalised to be a proper probability
density. Typically one writes πM (θM ) ∝ fM (θM ), where fM is a given non-integrable function
and the constant of proportionality is not specified. Even without that specification, this allows
mechanical computation of a formal within-model-M posterior density πM (θM | x), by application
of Bayes’s formula: πM (θM | x) ∝ pM (x | θM ) πM (θM ) ∝ pM (x | θM ) fM (θM ). This will often yield
an integrable function and hence the possibility of normalisation to supply a genuine probability
density.
However things do not work out so well when we turn to model selection. We have, for each
model M ,
πM (θM ) = cM fM (θM ),
where cM is the unspecified proportionality constant. This formally leads to the marginal likelihood function
Z
L M ∝ cM
pM (x | θM ) fM (θM ) dθM .
ΘM

But since this involves the unspecified constants cM , which could vary arbitrarily with M , it is no
longer meaningful to compare models by means of their marginal likelihoods.
A way round this problem was proposed in Dawid and Musio (2015): instead of attempting
to minimise the log score SL (x, PM ) := − log pM (x), we replace that with another proper scoring
rule S(x, PM ). And if that scoring rule is homogeneous, it will simply not involve the unspecified
constant cM . In Dawid and Musio (2015) a detailed analysis of this approach was conducted for
the case of continuous data and the Hyvärinen scoring rule, and it was shown that it will typically
deliver consistent selection of the true model.

4

Discrete model selection

We shall investigate empirically, for a simple example, the validity of the above results when
generalised to the case of discrete data. We shall use the scoring rule (5), and apply this to
the choice between a Poisson and a Negative Binomial model. For this purpose we first need to
compute, for each of these models separately, the appropriate score.

5

Poisson model

Consider the Poisson model X ∼ P(kΛ):
p(x | λ) = e−kλ (kλ)x /x!

(x = 0, 1, . . .),

(7)

with conjugate prior Λ ∼ Γ(α, β):
π(λ) =

β α α−1 −βλ
λ
e
.
Γ(α)

(8)

For propriety we require α > 0, β > 0.
The predictive distribution P has density function
p(x) =

Γ(α + x)
(1 − φ)α φx
Γ(α)x!

(9)

with φ := k/(β + k).
Then p(x + 1)/p(x) = φ(x + α)/(x + 1), and so
S(0, P )
S(x, P )

= m−1 αm φm

(10)
−1

= {m(m − 1)}


(m − 1)φm (x + 1)a−m (x + α)m

− mφm−1 xa−m+1 (x + α − 1)m−1
3

(x > 0).

(11)

5.1

Multiple observations

Suppose now we have N independent and identically distributed observations X N = (X1 , . . . , XN )
from the above Poisson distribution. We can apply the above score in two different ways:
(a). Apply direct to the sufficient statistic.
(b). Apply prequentially to all observations.
5.1.1

Sufficient statistic

PN
The sufficient statistic is TN = i=1 Xi , with distribution P(N kΛ). So the score computed this
way is simply obtained from (10) and (11) on replacing x by tN and k by N k. This gives
= m−1 αm φm
N

SN (0, P )

(12)
−1

= {m(m − 1)}

SN (x, P )


a−m
(m − 1)φm
(tN + α)m
N (tN + 1)

− mφm−1
ta−m+1
(tN + α − 1)a−m+1
N
N

(tN > 0)

(13)

where φN := N k/(β + N k).
5.1.2

Prequential

Now suppose we have already observed X n−1 = xn−1 . The posterior distribution of Λ is
Λ | X n−1 = xn−1 ∼ Γ {α + tn−1 , β + (n − 1)k} .
So the predictive distribution of Xn , given the previous observations X n−1 = xn−1 , is obtained
from (10) and (11) on replacing x with xn , α with α + tn−1 , and β with β + (n − 1)k. The
incremental contribution to the prequential score is thus given by:
Sn∗ (0, P )

=

Sn∗ (xn , P )

=

m−1 (φ∗n )m (α + tn−1 )m

{m(m − 1)}−1 (m − 1)(φ∗n )m (xn + 1)a−m (tn + α)m

(14)

− m(φ∗n )m−1 xa−m+1
(tn + α − 1)a−m+1
n

(15)

(xn > 0)

with φ∗n := k/(β + nk).
The total prequential score is obtained by summing this from n = 1 to N .

5.2

Improper prior

The usual improper prior is the formal limit with α, β ↓ 0. In this case (12) and (13) become:
SN (0, P )
SN (x, P )

=

0

(16)
−1

= {m(m − 1)}


(m − 1)(tN + 1)a−m tm
N

− mta−m+1
(tN − 1)a−m+1
N

(tN > 0).

(17)

Note that the score is well-defined even when all observations are 0, in which case the posterior is
improper.
For the prequential version, we obtain, from (14) and (15):
Sn∗ (0, P )
Sn∗ (xn , P )

m
= tm
n /mn

=

(18)

m
(xn + 1)a−m tm
n /mn
− xa−m+1
(tn − 1)a−m+1 /(m
n

− 1)nm−1

(xn > 0).

(19)

An alternative improper prior is the Jeffreys prior, having α = 1/2, β ↓ 0, which is easily handled
similarly.
4

6

Negative Binomial model

Now we consider an alternative model, the Negative Binomial X ∼ N B(s; Θ), having
p(x | θ) =

(s + x − 1)!
(1 − θ)s θx
x!(s − 1)!

(x = 0, 1, . . .),

(20)

with conjugate prior Θ ∼ β(p, q):
Γ(p + q) p−1
θ
(1 − θ)q−1 .
Γ(p)Γ(q)

(21)

Γ(p + q) (s + x − 1)! Γ(p + x)Γ(q + s)
.
Γ(p)Γ(q) x!(s − 1)! Γ(p + q + s + x)

(22)

π(θ) =
For propriety we require p > 0, q > 0.
The predictive density is
p(x) =
Then

(x + s)(x + p)
p(x + 1)
=
,
p(x)
(x + 1)(x + p + q + s)

and so we have:

6.1

S(0, P )

=

S(x, P )

=

m−1 (sp)m (p + q + s)−m

{m(m − 1)}−1 (m − 1)(x + 1)a−m {(x + s)(x + p)}m (x + p + q + s)−m

− mxa−m+1 {(x + s − 1)(x + p − 1)}m−1 (x + p + q + s − 1)−m+1 .

(23)
(24)

Multiple observations

Again, we can handle multiple observations either by restricting to the sufficient statistic, or by
cumulating the prequential score.
6.1.1

Sufficient statistic

PN
The sufficient statistic is TN = i=1 Xi , with distribution N B(N s, Θ). So the score computed
this way is simply obtained from (23) and (24) on replacing x by tN and s by N s. This gives
SN (0, P )

=

SN (x, P )

=

6.1.2

m−1 (N sp)m (p + q + N s)−m
(25)

−1
a−m
m
−m
(m − 1)(tN + 1)
{(tN + N s)(tN + p)} (tN + p + q + N s)
{m(m − 1)}

a−m+1
− mtN
{(tN + N s − 1)(tN + p − 1)}m−1 (tN + p + q + N s − 1)−m+1 . (26)

Prequential

Now suppose we have already observed X n−1 = xn−1 . The posterior distribution of Θ is
Θ | X n−1 = xn−1 ∼ β {p + tn−1 , q + (n − 1)s} .
So the predictive density of Xn , given the previous observations X n−1 = xn−1 , is obtained
from (23) and (24) on replacing x with xn , p with p+tn−1 , and q with q +(n−1)s. The incremental
contribution to the prequential score is thus given by:
Sn∗ (0, P )

=

Sn∗ (xn , P )

=

m−1 sm (p + tn−1 )m (p + q + tn−1 + ns)−m
(27)

m
{m(m − 1)}−1 (m − 1)(xn + 1)a−m {(xn + s)(p + tn )} (p + q + tn + ns)−m
i
m−1
−m+1
− mxa−m+1
{(x
+
s
−
1)(p
+
t
−
1)}
(p
+
q
+
t
+
ns
−
1)
. (28)
n
n
n
n

The total prequential score is obtained by summing this from n = 1 to N .
5

6.2

Improper prior

The usual improper prior is the formal limit with p, q ↓ 0. In this case (25) and (26) become:

SN (0, P )
S(x, P )

=

0

(29)
−1

= {m(m − 1)}


(m − 1)(tN + 1)a−m tm
N

− mta−m+1
(tN − 1)m−1 .
N

(30)

The score is well-defined even when all observations are 0, in which case the posterior is
improper.
For the prequential version, we obtain, from (27) and (28):
Sn∗ (0, P )

−m
= m−1 sm tm
n−1 (tn−1 + ns)

−m
Sn∗ (xn , P ) = {m(m − 1)}−1 (m − 1)(xn + 1)a−m (xn + s)m tm
n (tn + ns)

− mxa−m+1
(xn + s − 1)m−1 (tn − 1)m−1 (tn + ns − 1)−m+1 .
n

(31)

(32)
The total prequential score is obtained by summing this from n = 1 to N .
Again, similar expressions can be found using the improper Jeffreys prior, which has p ↓ 0,
q = 1/2.

7

Simulations

We generated observations from either the Poisson distribution (7) with k = 1, λ = 10, or the
Negative Binomial distribution (20) with s = 81, θ = 0.1. These both have variance 10, the former
having mean 10, and the latter mean 9. We used, as the scoring rule, the special case of (5) having
a = m = 2, namely
1
S(x, P ) = (x + 1)2
2



p(x + 1)
p(x)

2
−x

2



p(x)
p(x − 1)


δ(x > 0).

For each generating distribution we computed the excess of the cumulative prequential score
for the wrong model over that for the correct model. These differences are shown, as a function
of increasing data, in Figures 1 and 2 respectively. Each figure displays 10 sample sequences
generated from the indicated distribution, as well as the average taken over a sample Areof 100
sequences.
In each case we see a clear linear upward trend, supporting the expectation of consistent model
selection, although even with 1000 observations there is a non-negligible probability of a negative
value, giving a misleading preference for the wrong model.

8

Conclusions

We have extended the Bayesian model selection methodology of Dawid and Musio (2015) to apply
to problems with discrete data. We have conducted a simulation study to compare Poisson and
Negative Binomial distributions. The results suggest that the method will consistently select the
correct model as the number of data points increases.

Acknowledgements
Philip Dawid’s research was supported through an Emeritus Fellowship from the Leverhulme Trust.
6

(Negbinomial − Poisson), under Poisson
150
10 Single differences in score
Average difference in score
over 100 replicates

Differences in score

100

50

0

−50

0

200

400

600

800

1000

n

Figure 1: Data from Poisson distribution P(10)

(Poisson− Negbinomial), under Negbinomial

Differences in score

150

10 Single differences in score
Average difference in score
over 100 replicates

100

50

0

−50

0

200

400

600

800

n

Figure 2: Data from Negative Binomial distribution N B(81; 0.1)

7

1000

References
Dawid, A. P. (2011). Posterior model probabilities. In Philosophy of Statistics, (ed. P. S. Bandyopadhyay and M. Forster), pp. 607–30. Elsevier, New York.
Dawid, A. P., Lauritzen, S. L. and Parry, M. (2012). Proper local scoring rules on discrete
sample spaces. Ann. Statist. 40 593–608.
Dawid, A. P. and Musio, M. (2015). Bayesian model selection based on proper scoring rules
(with Discussion). Bayesian Analysis 10 479–521.
Hyvärinen, A. (2005). Estimation of non-normalized statistical models by score matching. Journal of Machine Learning 6 695–709.
Parry, M. F., Dawid, A. P., and Lauritzen, S. L. (2012). Proper local scoring rules. Annals
of Statistics 40 561–92.

8

