Bayesian estimation of a decreasing density
Geurt Jongbloed, Frank van der Meulen and Lixue Pang

arXiv:1801.02539v1 [math.ST] 8 Jan 2018

Delft University of Technology

Abstract: Suppose X1 , . . . , Xn is a random sample from a bounded and decreasing
density f0 on [0, ∞). We are interested in estimating such f0 , with special interest
in f0 (0). This problem is encountered in various statistical applications and has
gained quite some attention in the statistical literature. It is well known that the
maximum likelihood estimator is inconsistent at zero. This has led several authors
to propose alternative estimators which are consistent. As any decreasing density
can be represented as a scale mixture of uniform densities, a Bayesian estimator
is obtained by endowing the mixture distribution with the Dirichlet process prior.
Assuming this prior, we derive contraction rates of the posterior density at zero
by carefully revising arguments presented in Salomond (2014). Various methods
for estimating the density are compared using a simulation study. We apply the
Bayesian procedure to the current durations data described in Keiding et al. (2012).
Key words and phrases: Bayesian nonparametric approach, posterior consistency,
concentrate rate.

1

Introduction

Consider an independent and identically distributed (i.i.d.) sample X1 , . . . , Xn from a bounded
decreasing density f0 on [0, ∞). The problem of estimating f0 based on the sample, only using
the information that it is decreasing, has attracted quite some attention in the literature. One
of the reasons for this is that the estimation problem arises naturally in several applications.
See for instance the introductory section of Kulikov & Lopuhaä (2006), Vardi (1989), Watson
(1971) and Keiding et al. (2012). In all these examples, the sampling density f0 can be
expressed in terms of an underlying distribution function of interest H0 :
1 − H0 (x)
f0 (x) = R ∞
, x ≥ 0.
0 y dH0 (y)

(1)

In words: the sampling density is proportional to a survival function of interest, which is by
definition decreasing. Having an estimate of the bounded decreasing density f0 , an estimate
of H0 is obtained using the ‘inverse relation’
H0 (x) = 1 −

f0 (x)
, x ≥ 0,
f0 (0)
1

(2)

expressing H0 in terms of f0 .
The most commonly used estimator for f0 is the maximum likelihood estimator derived
in Grenander (1956). This estimator is defined as the maximizer of the log likelihood `(f ) =
P
n
ˆ
i=1 log f (Xi ) over all decreasing density functions on (0, ∞). The solution fn of this maximization problem can be graphically constructed. Starting from the empirical distribution Fn
based on X1 , . . . , Xn , the least concave majorant of Fn can be constructed. This is a concave
distribution function. The left-continuous derivative of this piecewise linear concave function
yields the maximum likelihood (or Grenander) estimator for f0 . For more details on the derivation of this estimate, see Section 2.2 in Groeneboom & Jongbloed (2014). As can immediately
be inferred from the characterization of the Grenander estimator,
Fn (X(1) )
Fn (Xi )
1
fˆn (0) := lim fˆn (x) = max
≥
=
,
1≤i≤n
x↓0
Xi
X(1)
nX(1)
where X(i) denotes the i-th order statistic of the sample. Denoting convergence in distribution
d

by →,
d

nf0 (0)X(1) → Y

as

n→∞

where Y has the standard exponential distribution, it is clear that fˆn (0) does not converge in
probability to f0 (0). This inconsistency of fˆn (0) was first studied in Woodroofe & Sun (1993).
There it is also shown that
N (t) d 1
fˆn (0) d
→ sup
=
f0 (0)
t
U
t>0

as n → ∞,

where N is a standard Poisson process on [0, ∞) and U is a standard uniform random variable.
It is clear from (2) that this inconsistency is undesirable, as estimating the distribution function of interest, H0 , at any point x > 0, requires estimation of f0 (0). Various approaches have
been taken to obtain a consistent estimator of f0 (0). The idea in Kulikov & Lopuhaä (2006) is
to estimate f0 (0) by fˆn evaluated at a small positive (but vanishing) number: fˆn (cn−1/3 ) for
some c > 0. There it is shown that the estimator is n1/3 -consistent, assuming f0 (0) < ∞ and
|f00 (0)| < ∞.
A likelihood related approach was taken in Woodroofe & Sun (1993). There a penalized
log likelihood function is introduced, where the estimator is defined as maximizer of
`α (f ) =

n
X

log f (Xi ) − αnf (0).

i=1

For fixed α ≥ 0, this estimator can be computed explicitly by first transforming the data
using a data dependent affine transformation and then applying the basic concave majorant
algorithm to the empirical distribution function based these transformed data. It is shown
(again, assuming f0 (0) < ∞ and |f00 (0)| < ∞) that the optimal rate to choose α is n−2/3 .
P
1/3 -consistent.
Then, the maximum penalized estimator fˆn,
α̂n (0) is n
Groeneboom & Jongbloed (2014) proposed to estimate f0 (0) by the histogram estimator
b−1
F
n (bn ), where {bn } is a sequence of positive numbers with bn → 0 if n → ∞. The bin
n
widths bn can e.g. be chosen by estimating the asymptotically Mean Squared Error-optimal
choice.
2

In this paper we take a Bayesian nonparametric approach to the problem. An advantage
of the Bayesian setup is the ease of constructing credible regions. To construct frequentist
analogues of these, confidence regions, can be quite cumbersome, relying on either bootstrap
simulations or asymptotic arguments.
To formulate a Bayesian approach for estimating a decreasing density, note that any decreasing density on [0, ∞) can be represented as a scale mixture of uniform densities:
Z ∞
ψx (θ)dGf (θ), where ψx (θ) = θ−1 1[0,θ] (x),
(3)
fG (x) =
0

where G is a distribution function concentrated on the positive half line. Therefore, by endowing the mixing measure with a prior distribution we obtain the posterior distribution of
the decreasing density, and in particular of f0 (0). A convenient and well studied prior for
distribution functions on the real line is the Dirichlet process (DP) prior (see for instance Ferguson (1973) and Van der Vaart and Ghosal (2017)). This prior contains two parameters: the
concentration parameter, usually denoted by α, and the base measure, which we will denote
by G0 . The approach where a prior is obtained by putting a Dirichlet process prior on G in
(3) was previously considered in Salomond (2014). In that paper, the asymptotic properties
of the posterior in a frequentist setup were studied. More specifically, contraction rates are
derived to quantify the performance of the Bayesian procedure. This is a rate for which we
can shrink balls around the true parameter value, while maintaining most of the posterior
mass. More formally, if L is a semimetric on the space of density functions, a contraction
rate εn is a sequence of positive numbers εn ↓ 0 for which the posterior mass of the balls
{f : L(f, f0 ) ≤ εn } converges in probability to 1 as n → ∞, when assuming X1 , X2 , . . . are
independent and identically distributed with density f0 . A general discussion on contraction
rates is given in Van der Vaart and Ghosal (2017).
In this paper, we show that the posterior of f (0) contracts at rate (log n/n)2/9 . The same
result is also claimed in Salomond (2014), though under stronger assumptions on the base
measure than ours. The reason for considering these different conditions is that we believe
there is a difficulty in the argument presented in his paper (see Appendix A.1). Nevertheless,
many results from Salomond (2014) are important ingredients to the proof we present. One
of the key contributions of this paper is to show how these can be used, together with new
arguments to derive the claimed contraction rate. We also address computational aspects
of the problem and show how draws from the posterior can be obtained using the algorithms
presented in Neal (2000). By a simulation study, we compare the behaviour of various proposed
methods for estimating f0 (0).
The outline of this paper is as follows. In section 2 we derive pointwise contraction rates
for the density evaluated at x, for any x ≥ 0. In section 3 a Markov Chain Monte Carlo
method for obtaining draws from the posterior is given, based on the results of Neal (2000).
This is followed by a review of some existing methods to consistently estimate f0 at zero. In
section 5 these methods are compared with the Bayesian approach. The Bayesian method is
also applied to the current durations data described in Keiding et al. (2012). In the appendix
we point out a difficulty in the proof of (Salomond (2014)) to establish posterior consistency
at zero.

3

2

Pointwise posterior contraction rates for f0

Let F denote the collection of all bounded decreasing densities on [0, ∞). Denote the distribution of X n = (X1 , . . . , Xn ) under f by Pf and expectation under Pf by Ef . In this section
we are interested in the asymptotic behaviour of the posterior in a frequentist setup. This
entails that we study the behaviour of the posterior distribution on F while assuming a true
underlying density f0 . Set P0 = Pf0 and E0 = Ef0 . Denote the prior measure on F by Π and
the posterior measure by Π(· | X n ).
Given a loss function L on F, we say that the posterior is consistent with respect to L
if for any ε > 0, E0 Π(L(f, f0 ) > ε | X n ) → 0 when n → ∞. If {εn } is a sequence that
tends to zero, then we say that the posterior contracts at rate εn (with respect to L) if
E0 Π(L(f, f0 ) > εn | X n ) → 0 when n → ∞. The rate {εn } is called a contraction rate.
Salomond (2014) derived contraction rates based on the Dirichlet process prior for the
L1 −, Hellinger- and pointwise loss function. As we are specifically interested in estimating
the density at zero, we focus on Theorem 4 in that paper, which gives sufficient conditions for
posterior consistency under pointwise loss at strictly positive values of x. As we will explain
more in Appendix A.1, there is a problem with the argument in Salomond (2014) for proving
posterior contraction for f0 (0).
In the following theorem we show how the proof of Salomond (2014) can be adapted to
yield sufficient conditions (different from Salomond’s) for posterior contraction in terms of the
behaviour of the density of the base measure near zero. Although the argument in Salomond
(2014) for proving posterior concentrate rate n for f0 (x) with x > 0 is correct, we prove the
theorem below for x ≥ 0 rather than only for x = 0 case.
0

0

2.1 Theorem Fix f0 ∈ F and x ∈ [0, ∞) such that f0 (x) < 0 and supx∈[0,x0 ] | f0 (x) |< ∞ for
some x0 > 0. Let Xi , i = 1, 2, . . . be independent random variables, each with density f0 . Take
Π the distribution on F that is obtained via (3), putting the DP (G0 , α) on the distribution
functions G. We assume that G0 has a strictly positive Lebesgue density g0 on (0, ∞) such
that there exist constants θ0 > 0, a > 0 and t1 > 0 with
e−a/θ . g0 (θ) . θt1

for all

θ ∈ (0, θ0 ).

(4)

Then there exists a constant C > 0 such that, when n → ∞,
n
o

E0 Π f ∈ F : |f (x) − f0 (x)| > C(log n/n)2/9
X n → 0.
In the proof we will use the following three lemmas. In each of these lemmas, Π is the prior
as in Theorem 2.1 and f0 satisfies the assumptions stated within this theorem. Throughout,
we take
n = (log n/n)1/3
and
ηn = 2/3
(5)
n .
Proofs are given in section 2.1.
2.2 Lemma Define
Dn =

Z Y
n
f (Xi )
dΠ(f ).
f0 (Xi )
i=1

4

(6)

There exists a constant c > 0 such that


2
P0 Dn < e−cnn = o(1)

as n → ∞,

(7)

where εn as just defined in (5).
2.3 Lemma Let c0 > 0 and {ξn } be any sequence with ξn → ∞ as n → ∞. Define Mn =
exp (t1 + 1)−1 c0 ξn with t1 as it occurs in condition (4) and the sieves Fn ⊂ F by
Fn = {f ∈ F : f (0) ≤ Mn }.
Then
Π(Fnc ) . e−c0 ξn .
2.4 Lemma Let {ξn } be a sequence of positive numbers tending to zero. Then for any  > 0,
there exists a set Ω with probability Π(Ω) ≥ 1 −  and N0 (depending on G0 and ) such that
when n ≥ N0 , for any G ∈ Ω, we have fG (0) − fG (ξn ) ≤ ξn .
We now give the proof of Theorem 2.1.
Proof The posterior measure of a measurable set E ⊂ F is given by
n

Π(E | X ) =

Dn−1

Z Y
n
f (Xi )
dΠ(f ),
f0 (Xi )
E
i=1

2

where Dn is as defined in (6). Recalling (5), let Dn = {Dn < e−cnn } (the constant c being as
in lemma 2.2). For a constant C > 0 and x ≥ 0, define the following subsets of F:
Bn (x) = {f ∈ F : |f (x) − f0 (x) |> Cηn },
and
Bn1 = Fnc ,

Bn2 (x) = {f ∈ Fn :| f (x) − f0 (x) |> Cηn }.

Clearly, as Bn (x) ⊂ Bn1 ∪ Bn2 (x), Π(Bn (x) | X n ) ≤ Π(Bn1 | X n ) + Π(Bn2 (x) | X n ). We have
E0 Π(Bn1 | X n ) = E0 Π(Bn1 | X n )1Dn + E0 Π(Bn1 | X n )1Dnc
Z Y
n
f (Xi )
cn2n
≤ E0 (1Dn ) + e
E0
dΠ(f )
f0 (Xi )
Fnc
i=1

= P0 (Dn ) + e

cn2n

Π(Fnc )

= o(1)

as

n → ∞,

where the final equality follows from lemma 2.2 and lemma 2.3 (applied with ξn = nε2n and
c0 > c).

5

Consider (test-) functions Φn : R → [0, 1]. Then we can bound E0 Π(Bn2 (x) | X n ) as follows
E0 Π(Bn2 (x) | X n )
= E0 Π(Bn2 (x) | X n )1Dn + E0 Π(Bn2 (x) | X n )1Dnc Φn (x)
+ E0 Π(Bn2 (x) | X n )1Dnc (1 − Φn (x))
"
Z

n
Y
f (Xi )
(1 − Φn (x))dΠ(f )1Dnc
f0 (Xi )
Bn2 (x)

Dn−1

≤ E0 [1Dn ] + E0 (Φn (x)) + E0
2

≤ P0 (Dn ) + E0 (Φn (x)) + ecnn E0

i=1
n
Y

Z

Bn2 (x) i=1
2

= E0 (Φn (x)) + ecnn

#

f (Xi )
(1 − Φn (x))dΠ(f )
f0 (Xi )

Z
Ef (1 − Φn (x))dΠ(f ) + o(1).

(8)

Bn2 (x)

Here we use lemma 2.2 again to obtain the final equality.
To construct the specific test functions Φn (x), we distinguish between x > 0 and x = 0.
√
First consider x > 0 and fix a sequence {cn } ↓ 0 such that ncn → ∞ and cn ≤ C 2 ηn2 /16 (one
could e.g. choose cn ∼ n−4/9 ). Fix a constant c0 > 0 and define
)
(
Z x
n
X
+
−1
f0 (t)dt > cn
ϕn (x) = 1 n
1[x−c0 Cηn ,x] (Xi ) −
(
−1
ϕ−
n (x) = 1 n

i=1
n
X

x−c0 Cηn

1[x,x+c0 Cηn ] (Xi ) −

i=1

)

x+c0 Cηn

Z

f0 (t)dt < −cn

.

(9)

x

(recall that the constant C is as it appears in the statement of theorem 2.1).
−
Define the tests by Φn (x) = max(ϕ+
n (x), ϕn (x)). As derived in the proof of theorems 3 and
5 in Salomond (2014), there exists a constant C 0 > 0 such that the tests Φn (x) satisfy
E0 Φn (x) = o(1)
sup

0

3

0

Ef (1 − Φn (x)) ≤ e−C n(Cηn ) = e−C C

3 n2
n

.

f ∈Bn2 (x)

Substituting these bounds into (8) and choosing C > (c/C 0 )1/3 shows that E0 Π(Bn2 (x) |
X n ) → 0 as n → ∞. This finishes the proof for x > 0.
We now consider the case x = 0. Define subsets
+
Bn2
(0) = {f ∈ Fn : f (0) − f0 (0) > Cηn }
−
Bn2
(0) = {f ∈ Fn : f (0) − f0 (0) < −Cηn }.
+
−
+
−
As Bn2 (0) = Bn2
(0) ∪ Bn2
(0), Π(Bn2 (0) | X n ) ≤ Π(Bn2
(0) | X n ) + Π(Bn2
(0) | X n ). For
−
−
−
−
bounding E0 Π(Bn2
(0) | X n ), we define the test function Φ−
n (0) by Φn (0) := ϕn (0), with ϕn (0)
−
as defined in (9). Then it follows from the inequalities in (8), applied with Bn2 (0) instead of
−
Bn2 (x), that E0 Π(Bn2
(0) | X n ) = o(1) as n → ∞.
+
For bounding E0 Π(Bn2
(0) | X n ), we use a similar but slightly different test, given by
(
)
Z ξn
n
X
−1
Φ+
1[0,ξn ] (Xi ) −
f0 (t)dt > cn ,
n (0) = 1 n
0

i=1

6

where ξn = Cηn /4. By Markov’s inequality E0 Φ+
n (0) = o(1).
+
We aim to derive an upper bound on Ef (1 − Φ+
n (0)) for f ∈ Bn2 (0) for n sufficiently large.
Suppose that f (0) − f0 (0) > Cηn . Then it follows that
Z ξn
(f (t) − f0 (t)) dt ≥ (f (ξn ) − f0 (0))ξn ≥ (f (ξn ) − f (0))ξn + Cξn ηn .
0

Here, the first inequality follows from the monotonicity of both f and f0 . This implies
!
Z ξn
n
X
−1
f0 (t)dt ≤ cn
Ef (1 − Φ+
1[0,ξn ] (Xi ) −
n (0)) = Pf n
0

i=1


Z
n
X
−1
≤ Pf n
1[0,ξn ] (Xi ) −

ξn

f (t)dt ≤ cn

0

i=1


+ (f (0) − f (ξn ))ξn − Cξn ηn .
By lemma 2.4, for any  > 0, there exists a set Ω with probability at least 1 −  and for any
f ∈ Ω we have f (0) − f (ξn ) ≤ ξn when n ≥ N0 . Hence
Z
Z
+
Ef (1 − Φn (0)) =
Ef (1 − Φ+
n (0))dΠ(f )
+
Bn2
(0)

+
Bn2
(0)∩Ω

Z
+
+
Bn2
(0)∩Ωc

Z
≤
+
Bn2
(0)∩Ω

Ef (1 − Φ+
n (0))dΠ(f )

Ef (1 − Φ+
n (0))dΠ(f ) + .

+
(0) ∩ Ω, because cn ≤ ξn2 , we find that
For f ∈ Bn2

Ef (1 −

Φ+
n (0))

≤ Pf

n

−1

n
X

Z
1[0,ξn ] (Xi ) −

i=1

ξn

!
2

f (t)dt ≤ −(Cηn ) /8 .
0

It is easily seen that the random variable appearing on the right hand side of this display
Rξ
is bounded, as n−1 (1[0,ξn ] (X1 ) − 0 n f (t)dt) ≤ n1 . The next step is to obtain the bound on
R
2
Rξ
ξ
variance Vf := 0 n f (t)dt − 0 n f (t)dt . Note that E0 Π(k f0 − f kL1 > Cn | X n ) → 0
(Theorem 1 in Salomond (2014)). As any set A can be written as union of two disjoint subsets
A1 = A ∩ {k f0 − f kL1 > Cn } and A2 = A ∩ {k f0 − f kL1 ≤ Cn }, we can see that the
expectation under P0 of Π(A1 | X n ) tends to zero. Hence we can just consider f such that
k f0 − f kL1 ≤ Cn , then
Z ξn
f (t)dt ≤k f0 − f kL1 +ξn f0 (0) ≤ Cn + M ξn . ξn .
0

Hence, Vf . ξn . By Bernstein’s inequality (Van der Vaart (1998), lemma 19.32), we have
!
1
2 /8)2
(c
(Cη
)
1
n
2
≤ exp(−c2 n(Cηn )3 )
Ef (1 − Φ+
n (0)) ≤ exp −
2
n · n12 Vf + 3n
(c1 (Cηn )2 /8)
2

for some c2 > 0. Substituting these bounds in (8) and choosing C > (c/c2 )1/3 ,  = o(e−cnn )
+
shows that E0 Π(Bn2
(0) | X n ) → 0 as n → ∞.

7

2.1

Proof of Lemmas

Proof of Lemma 2.2
Let F and F0 denote the cumulative distribution functions corresponding to the probability densities
respectively. Define the Kullback-Leiber divergence from f to f0 by
R f and f0f (x)
KL(f, f0 ) = f (x) log f0 (x) dx. For n as defined in (5), let
θn = F0−1 (1 − n /(2n)),

fn (·) =

f (·)I[0,θn ] (·)
,
F (θn )

f0,n (·) =

f0 (·)I[0,θn ] (·)
.
F0 (θn )

In appendix B of Salomond (2014) it is shown that if suffices to prove that
Π(Sn ) ≥ exp(−C1 n2n )

(10)

for some positive constant C1 , where
)
(


Z θn
Z
f (x) 2
2
2
2
f (x)dx ≥ 1 − n .
dx ≤ n ,
Sn = f : KL(f0,n , fn ) ≤ n , f0,n (x) log
f0 (x)
0
The main idea for proving (10) is to construct a subset Nn ⊂ Sn such that the prior mass of
2
Nn is bounded below by e−C1 nn (Cf. Salomond (2014) lemma 8). Here
Nn = {fP 0 : P 0 ([0, ∞)) = 1, | P 0 (Ui ) − p̃i |≤ 2n m−1 , i = 1, . . . , m0 + m}
satisfies

Pm+m0
i=1

p̃i = 1 and the sets U0 , U1 , . . . , Um0 +m are defined as
Ui = (in , (i + 1)n ]

for i = 1, . . . , m0

Um0 +i = [θ0 ∨ (yi − 2n ), yi + 2n ]
U0 = [0, ∞) ∩

for i = 1, . . . , m

0
c
(∪m+m
i=1 Ui ) .

−1
The numbers satisfy m0 ∼ ε−1
n , m . εn and yi ∈ [θ0 , ∞) for all i = 1, . . . , m. In Supplementary A.2 we show that U0 , U1 , . . . , Um0 +m is a partition of [0, ∞). By the definition of Dirichlet
Process, if P 0 ∼ DP (α, G0 ), then P 0 (U0 ), . . . , P 0 (Um ) ∼ Dir(αG0 (U0 ), . . . , αG0 (Um )). Without loss of generality, assume αG0 (Ui ) < 1, for i = 0, 1, . . . , m0 + m. Similar to Lemma 6.1 in
Ghosal et al. (2000), we have

Π(Nn ) ≥ Dir(P 0 (Ui ) ∈ [p̃i ± 2n m−1 ], 1 ≤ i ≤ m + m0 )
0 +m
Z p̃i +2n m−1
mY
1
αG (U )−1
≥ Γ(α)
x 0 i dxi .
Γ(αG0 (Ui )) 0∧(p̃i −2n m−1 ) i
i=1

αG0 (Ui )−1

Here we use (P 0 (U0 ))αG0 (U0 )−1 ≥ 1. As xi
Z

pi +2n m−1

0∧(pi −2n m−1 )

≥ 1 we have

αG0 (Ui )−1

xi

8

dxi ≥ 22n m−1 .

Substituting this bound into the lower bound on Π(Nn ), combined with the inequalities
βΓ(β) = Γ(β + 1) ≤ 1 for 0 < β ≤ 1 and m . −1
n , we obtain
!
0 +m
0 +m
mY
mX
3(m0 +m)
0
Π(Nn ) & n
G0 (Ui ) = exp 3(m + m) log n +
log G0 (Ui ) .
i=1

i=1

When m0 < i ≤ m0 + m it is trivial that G0 (Ui ) & 2n and therefore
0 +m
mX

log G0 (Ui ) & m log n .

i=m0 +1

For bounding G0 (Ui ) when i ≤ m0 , we use the lower bound of (4): g0 (θ) & e−a/θ . Without
loss of generality, we assume the associated constant of this lower bound is bigger than 1. In
this case we have
Z
e−a/θ dθ & n exp(−a/(in )).
G0 (Ui ) &
Ui

Implying
0

m
X

0

0

log G0 (Ui ) ≥ m log n −

a−1
n

i=1

Since

Pm0

−1
i=1 i

m
X

i−1 .

i=1

∼ log(m0 ) ∼ log −1
n , we therefore have (up to an additive constant)
0

m
X

log G0 (Ui ) & −1
n log n .

i=1

Therefore, we obtain
−1

Π(Sn ) ≥ Π(Nn ) & eC1 n

log n

2

& e−C1 nn

for some C1 > 0. This is exactly as is required.
Proof of Lemma 2.3
any decreasing density f , the mixture representation given in (3) gives fG (0) =
R ∞ For
−1 dG (θ) and
θ
f
0
Z ∞
Z ∞
θ−1 dGf (θ) ≤ 2−1 Mn
dGf (θ) ≤ Mn /2.
2Mn−1

2Mn−1

Hence, we have
Π (Fnc )

Z

∞

=Π

θ

−1


≤Π

dGf (θ) > Mn

≤

!
θ

0

2Mn−1

2Mn−1

Z

−1

dGf (θ) > Mn /2

0

Z

2Mn−1

θ

−1

g0 (θ)dθ .

0

Mn−1

Z

2Mn−1

θt1 −1 dθ . e−c0 ξn .

0

Here, the second inequality follows from Markov’s inequality and the third inequality follows
from g0 (θ) . θt1 , which is in the prior mass condition (4).
9

Proof of Lemma 2.4
Let Π be as described in the statement of Theorem 2.1 and let h be a real-valued function
on [0, 1] which is strictly increasing and convex on (0, ζ) for some ζ > 0. Let γt be a standard
Gamma process on [0, ∞), then
E exp(−λγt ) = exp(−tΦ(λ)) = (λ + 1)−t ,
where Φ is the Laplace transform of γ1 . By the proof of proposition 10 in Bertoin (1998), we
have for any δ ∈ (0, 1),
P (γt ≥ δh(t)) ≤ (1 − 1/e)−1 (1 − exp(−tΦ(1/δh(t))))
= (1 − 1/e)−1 [1 − (1/(δh(t)) + 1)−t ]
= (1 − 1/e)−1 [1 − exp(−t log(1 + 1/δh(t)))].
Take h(t) = exp(−(t | log t |r )−1 ), r > 1. When t ↓ 0, we have
t log(1 + 1/δh(t)) ∼ −t log δh(t) = −t log δ+ | log t |−r → 0.
Hence, for any  > 0, exists t0 (, δ) > 0 such that for any t ≤ t0 ,
P (γt ≥ δh(t)) ≤ .
0 (θ))
. The above statement gives that for any  > 0, there exists
We can represent G(θ) as γ(αG
γ(α)
a set Ω with probability P (Ω) ≥ 1 −  and θ1 (, G0 ) such that for any G ∈ Ω, we have

G(θ) ≤ h(αG0 (θ))

for θ ≤ θ1 .

(11)
R ξn

From the mixture representation in (3), it follows that fG (0) − fG (ξn ) = 0 θ−1 dG(θ). Under
condition (4), G0 (θ) . θt1 +1 , hence we have limθ↓0 θ−1 G(θ) ≤ limθ→0 θ−1 h(αG0 (θ)) = 0. Then
it follows from integration by parts followed by (11) that exists N0 (depends on ), when
n ≥ N0 ,
Z ξn
Z ξn
Z ξn
θ−2 h(αG0 (θ))dθ.
θ−2 G(θ)dθ ≤ ξn−1 G(ξn ) +
θ−1 dG(θ) = ξn−1 G(ξn ) +
0

0

0

Integrating by parts and using (11) once again, we get that
Z ξn
Z ξn
θ−1 dG(θ) ≤ ξn−1 G(ξn ) − ξn−1 h(αG0 (ξn )) +
θ−1 dh(αG0 (θ))
0
0
Z ξn
≤
θ−1 dh(αG0 (θ)).
0

Combining these inequalities and using that for t close to zero
h0 (t) = exp(−(t | log t |r )−1 )

| log t | −r
≤ t.
| log t |r+1

t2

As well as (4), we obtain for n ≥ N0
Z ξn
Z
−1
0
fG (0) − fG (ξn ) ≤
αθ g0 (θ)h (αG0 (θ))dθ ≤
0

0

10

ξn

θ−1 θt1 θt1 +1 dθ . ξn2t1 +1 . ξn .

3

Gibbs Sampling in the DMP model

Since a decreasing density can be represented as a scale mixture of uniform densities (see
(3)) and the mixing measure is chosen according to a Dirichlet process, the model is called a
Dirichlet Mixture Model. Algorithms for drawing from the posterior in such models have been
studied in many papers over the past two decades, a key reference being Neal (2000). Here we
shortly discuss the algorithm coined ”algorithm 2” in that paper. We assume G0 has a density
g0 with respect to Lebesgue measure.
Let #(x) denote the number of distinct values in the vector x and let x−i denote the vector
obtained by removing the i-th element of x. Denote by ∨(x) and ∧(x) the maximum and
minimum of all elements in the vector x respectively.
The starting point for deriving the algorithm is a construction to sample from the DPM
model:
Z := (Z1 , . . . , Zn ) ∼ CRP (α)
i.i.d

Θ1 , . . . , Θ#(Z) ∼ G0

(12)

ind

X1 , . . . , Xn | Θ1 , . . . , Θ#(Z) , Z1 , . . . , Zn ∼ U nif (0, ΘZi ).
Here CRP(α) denotes the “Chinese Restaurant Process” prior, which is a distribution on the
space of partitions of the integers {1, 2, . . . , n}. This distribution is most easily described in a
recursive way. Initialize by setting Z1 = 1. Next, given Z1 , . . . , Zi , let Li = #(Z1 , . . . , Zi ) and
set
(
Li + 1 with probability α/(i + α)
Zi+1 =
k
with probability Nk /(i + α).
P
where k varies over {1, . . . , Li } and Nk = ij=1 1{Zj = k}. In principle this process can be
continued indefinitely, but for our purposes it ends after n steps. One can interpret the vector Z
as a partitioning of the index set {1, . . . , n} (and hence the data X = (X1 , . . . , Xn )) into #(Z)
disjoint sets (sometimes called “clusters”). For ease of notation, write Θ = (Θ1 , . . . , Θ#(Z) ).
An algorithm for drawing from the posterior of (Z, Θ) is obtained by successive substitution
sampling (also known as Gibbs sampling), where the following two steps are iterated:
1. sample Θ | (X, Z);
2. sample Z | (X, Θ).
The first step entails sampling from the posterior within each cluster. For the k−th component
of Θ, Θk , this means sampling from
Y
Y
fΘk |X,Z (θk | x, z) ∝ fΘk (θk )
fXj |Θk (xj | θk ) = g0 (θk )
ψ(xj | θk ).
j:zj =k

j:zj =k

Sampling Z | (X, Θ) is done by cycling over all Zi (1 ≤ i ≤ n) iteratively. For i ∈ {1, . . . , n}
and k ∈ {1, . . . , 1 + ∨(Z)} we have
fZi |Z−i ,X,Θ (k | z−i , x, θ) ∝ fXi |Zi ,Z−i ,Θ (xi | k, z−i , θ)fZi |Z−i ,Θ (k | z−i , θ)
= fXi |ΘZ (xi | θk )fZi |Z−i (k | z−i )
i

11

The right-hand-side of this display equals
Nk,−i
ψ(xi | θk )
nZ− 1 + α

α
n−1+α

1 ≤ k ≤ ∨(Z),

if

ψ(xi | θ)dG0 (θ)

if k = 1 + ∨(Z),

P
where Nk,−i = j∈{1,...,n}\{i} 1{Zj = k}. The expression for k = 1 + ∨(Z) follows since in that
case sampling from Xi | Θk boils down to sampling from the marginal distribution of Xi .
It may happen that over subsequent iterations of the Gibbs sampler certain clusters disappear. Hence #(Z) and ∨(Z) will typically not be the same. If this happens, the Θj corresponding to the disappearing cluster is understood to be removed from the vector Θ (because
the cluster becomes “empty”, the
R prior and posterior distribution of such a Θj are equal). In
this step we need to evaluate ψ(xi | θ)dG0 (θ). One option is to numerically evaluate this
quantity for i = 1, . . . , n (it only needs to be evaluated once). Alternatively, the “no-gaps”
algorithm of MacEachern and Müller (1998) or “algorithm 8” of Neal (2000) can be used and
refer for further details to these papers.
P
It is well known (see e.g. Ferguson (1973)) that G | Θ1 , . . . , Θn ∼ DP (αG0 + ni=1 δΘi , α +
n), where δx denotes the Dirac measure on {x}. Moreover, given Θ1 , . . . , Θn , the random
variables X1 , . . . , Xn and random distribution function G are independent. Hence, for any
measurable function h
!
Z

Z
n
X
1
E
hdG | Θ1 , . . . , Θn , X1 , . . . , Xn =
α hdG0 +
h(Θi ) .
α+n
i=1

Now suppose in the j-th iteration
of the Gibbs
sampler (possibly after neglecting “burn in”


(j)
(j)
samples) we have obtained ΘZ1 , . . . , ΘZn . At iteration j, the mean of the posterior density
is given by
fˆ(j) (x) := E

Z


ψx dG | Θ =

α
α+n

Z

n

1 X
(j) 
ψx dG0 +
ψx ΘZi .
α+n

(13)

i=1

Assuming J iterations,
a Rao-Blackwellized estimator for the posterior mean is obtained by
P
computing J −1 Jj=1 fˆ(j) (x).
If g0 is chosen from a class conjugate to the uniform distribution, then sampling of both
Θ | (X, Z) and Z | (X, Θ) is simplified. It is well known that the class of Pareto-distributions
is conjugate to the uniform distribution. Write Y ∼ P ar(s, ξ) with s > 0 and ξ > 0 if
fY (y) = ξsξ y −ξ−1 1{y ≥ s}. As choosing the base measure to be the P ar(s, ξ) distribution
would exclude support points in (0, s), we propose to further endow S with the Gamma(λ, 1)distribution as prior. In this way, the base measure is a mixture of P ar(s, ξ)-distributions
where ξ is fixed and the mixing measure over s is the Gamma(λ, 1)-distribution:
Z
g0 (θ) =

ξ −ξ−1

ξs θ

1{θ ≥ s}s

λ−1 −s

e

−1 −ξ−1

ds = ξ(Γ(λ))

Z

θ

0

12

θ

sξ+λ−1 e−s ds.

(14)

By l’Hospital’s rule, we have g0 (θ) ∼ θλ−1 as θ ↓ 0. Therefore, if λ > 1 then this choice of base
measure satisfies condition (4) for posterior contraction.
For the Gibbs sampling algorithm we then have to iterate the steps
1. sample Θ | (X, Z, S);
2. sample Z | (X, Θ, S);
3. sample S | (X, Z, Θ).
Steps (1) and (2) remain unchanged then. Clearly, in step (3) we have
fS|X,Z,Θ (s | x, z, θ) = fS|Θ,Z (s | θ, z) ∝ fΘ|S,Z (θ | s, z)fS (s)
#(z)

= fS (s)

Y

fΘi |S (θi | s) ∝ fS (s)sξ#(z) 1{s ≤ ∧(θ)},

i=1

the product being taken over all (distinct) values in the vector θ. This implies that updating
S in step (3) requires sampling from a compactly supported density.

4

Review of existing methods for estimating the decreasing
density at zero

In this section we review some consistent estimators for a decreasing density f0 at zero that
have appeared in the literature. These will be compared with the Bayesian method of this
paper using a simulation study in section 5.

4.1

Maximum penalised likelihood

In Woodroofe & Sun (1993), the maximum penalised likelihood estimator is defined as the
maximiser of the following penalised log likelihood function:
`α (f ) =

n
X

log f (Xi ) − αnf (0).

i=1

Here α ≥ 0 is the penalty parameter. This estimator has the same form as the maximum
likelihood estimator (MLE), being piecewise constant with at most n discontinuities. For fixed
α ≥ 0, let x1 < · · · < xn < ∞ denote the ordered observed values and
w0 = 0

and wk = α + γxk ,

k = 1, . . . , n

where γ is the unique solution of the equation


αs/n
.
γ = min 1 −
1≤s≤n
α + γxs
Denote by f P (α, ·) the penalized estimator with penalty parameter α. When xn > α,
is a step function with

f P (α, ·)

f P (α, x) = f P (α, xk ),

∀xk−1 < x ≤ xk ,
13

∀k = 1, . . . , n.

At zero it is defined by right continuity and for x 6∈ [0, xn ] as f P (α, x) = 0. Here
(j − i)/n
.
0≤i<k k≤j≤n wj − wi

f P (α, xk ) = min max

Geometrically, for k = 1, 2, . . . , n, f P (α, xk ) is the left derivative of the least concave majorant
of the empirical distribution function of the transformed data wi , i = 1, . . . , n evaluated at wk .
Note that an alternative expression for f P (α, 0) is (1 − γ)/α which can be easily calculated.
Theorem 4 in Woodroofe & Sun (1993) states that
n1/3 {f P (αn , 0) − f0 (0)} ⇒d sup
t>0

W (t) − (c + βt2 )
t

where αn = cn−2/3 , β = −f0 (0)f00 (0)/2 and W (t) denotes the standard Brownian motion. In
Woodroofe & Sun (1993), the theoretically optimal constant c is determined by minimizing
the expected absolute value of the limiting distribution f P , resulting in c = 0.649 · β −1/3 .

4.2

Simple and ‘adaptive’ estimators

In Kulikov & Lopuhaä (2006), f0 (0) is estimated by the maximum likelihood estimator fˆn
evaluated at a small positive (but vanishing) number: fˆn (cn−1/3 ) for some c > 0. Of course,
the estimator depends on the choice of the parameter c.
In Kulikov & Lopuhaä (2006), Theorem 3.1, it is shown that
n
o
A21 n1/3 (fˆn (cB21 n−1/3 ) − f0 (cB21 n−1/3 )) + cB21 f00 (0)
converges in distribution to DR [W (t) − t2 ](c) when n → ∞. Here DR [Z(t)](c) is the right
derivative of the least concave majorant on [0, ∞)
p of the process Z(t), evaluated at c. Furthermore, B21 = 41/3 f0 (0)1/3 |f00 (0)|−2/3 and A21 = B21 /f0 (0).
Based on this asymptotic result, two estimators are proposed, denoted as f S and f A (‘S’ for
simple, ‘A’ for adaptive). The first is a simple one with cB21 = 1, then f S (0) = fˆn (n−1/3 ). The
second is f A (0) = fˆn (c∗ B21 n−1/3 ), where c∗ ≈ 0.345 is taken such that the the second moment
of the limiting distribution is minimized. Of course, to really turn this into an estimator, B21
has to be estimated. More on this in section 5.2.

4.3

Histogram estimator

In chapter 2 of Groeneboom & Jongbloed (2014) a natural and simple histogram-type estimator
for f0 (0) is proposed. Let {bn } be a vanishing sequence of positive numbers and consider the
estimator f H (0) = b−1
n Fn (bn ), where Fn is the empirical distribution of X1 , . . . , Xn . It can
be shown that Ef H (0) − f0 (0) behaves like bn f00 (0)/2 and the variance of f H (0) behaves like
f0 (0)/(nbn ) as n → ∞. Then the asymptotic mean square error (aMSE) optimal choice for bn
is (2f0 (0)/f00 (0)2 )1/3 n−1/3 = 2−1/3 B21 n−1/3 , where B21 is as defined in the Section 4.2.

14

5

Simulation study

In this section, we first study the influence of the base measure on the posterior mean by
considering two choices for this measure. Next, we compare our Bayesian method with various
consistent methods for estimating a decreasing density at zero as described in section 4. To
approximate the posterior mean estimator, we use the algorithm described in Section 3. These
simulation results are shown in section ?? of Supplementary Materials. Finally, in Section 5.3
we apply the Bayesian estimator to the fertility dataset from Keiding et al. (2012).

5.1

Influence of the base measure on the posterior mean

To assess the influence of the base measure in the Dirichlet-process prior, we ran simulations
using the following two choices for the base density g0 .
1. The mixture of Pareto densities:
Z
g0 (θ) =

ξsξ θ−s−1 1[s,∞) (θ)fS (s) ds,

where ξ > 0 is fixed and fS is the density of the Gamma(λ, 1)-distribution. As shown
in section 3 this is a convenient choice of base measure from a computational point of
view. Moreover, by the results in section 2 the resulting posterior is consistent at zero if
λ > 1. Here we take ξ = 0.1 and λ = 5.
2. The density g0 (θ) ∝ e−2θ−θ
bound of assumption (4).

−1

that vanishes exponentially fast near zero, as the lower

Figure 1 shows estimates of the posterior mean under the two proposed base measures and
various sample sizes (n = 50, 200, 1000 and 2000) from the standard exponential distribution.
The total number of MCMC iterations was chosen to be 2000, with 500 burn-in iterations.
The concentration parameter was α = 1. From these figures we can see that for increasing
sample sizes the posterior mean approximates f0 better. The estimates of the posterior mean
under the two choices of the base measure are similar for large sample sizes.

5.2

Comparing various methods for estimating f0 at 0

In this section we present a simulation study comparing our Bayesian estimator with the
various estimators available for f0 (0) discussed in section 4. We simulated 50 samples of sizes
n = 50, 200, 1000 from the standard exponential distribution and half-normal distribution. For
each sample, the following estimators are calculated: the posterior mean estimator f B , the
penalized NPMLE f P , the two estimators f S and f A and the histogram type estimator f H .
All these estimators require choosing some input parameters.
1. The posterior mean estimator f B (0) is computed using the DPM prior with concentration
parameter α = 1 and base measure the mixture of Pareto distributions as described in
section 5.1. For the hyperparameters ξ and λ, the values ξ = 0.1 and λ = 5 were taken
for samples from the exponential distribution and ξ = 0.1 and λ = 10 for the samples
from the half-normal distribution. The posterior mean was computed as described in
Section 3, in particular equation (13).
15

1.0
0.0

0.2

0.4

0.6

0.8

1.0
0.8
0.6
0.4
0.2
0.0

0

1

2

3

4

0

1

2

4

5

6

0.8
0.6
0.4
0.2
0.0

0.0

0.2

0.4

0.6

0.8

1.0

(b) n=200

1.0

(a) n=50

3

0

1

2

3

4

5

6

7

0

(c) n=1000

1

2

3

4

5

6

7

(d) n=2000

Figure 1: Estimated posterior mean in case the data are sampled from the standard exponential
distribution. The red and blue piecewise constant functions represent the posterior means based
on two choices of the base measure in section 5.1 (red: mixture of pareto; blue: exponential).
The black lines depicts the true density. The concentration parameter is α = 1 in all cases.

16

−1/3

2. For the penalized estimator f P (αn , 0) the parameter αn = 0.649β̂n n−2/3 was taken
with


f P (α0 , 0) − f P (α0 , xm ) −1/3
P
.
β̂n = max f (α0 , 0)
,n
2xm
Here xm is the second point of jump of f P (α0 , ·) and α0 = 0.0516, 0.0205 for n = 50, 200
(listed in Woodroofe & Sun (1993)).
3. For f S (0) = fˆn (n−1/3 ) no tuning is needed. For the other estimator we take f A (0) =
fˆn (0.345B̂21 n−1/3 ), where
B̂21 = 41/3 f S (0)1/3 |fˆn0 (0)|−2/3 ,

(15)

a consistent estimator of B21 where
fˆn0 (0) = min{n1/6 (fˆn (n−1/6 ) − fˆn (n−1/3 )), −n−1/3 }.
4. For the histogram estimator f H (0) = Fn (b̂n )/b̂n , b̂n = 2−1/3 B̂21 n−1/3 was chosen with
B̂21 as in (15).
Figure 2 shows, for each combination of sample size and estimation method described, the
boxplots of the 50 realized values based on samples from the standard exponential distribution.
Figure 3 shows these boxplots for the samples from the half normal distribution.
In table 1 we compare the bias, variance and mean squared error of these consistent estimators based on data from the standard exponential distribution. For the standard exponential
data, the Bayesian estimator f B performs best in the MSE sense for sample sizes n = 200 and
1000. When n = 50, the penalized estimator f P (0) has smallest MSE.

50

Bias
Var
MSE

fP
-0.059
0.033
0.036

200

Bias
Var
MSE

-0.033
0.023
0.024

-0.165
0.015
0.042

-0.140
0.026
0.045

-0.068
0.017
0.044

-0.077
0.015
0.021

1000

Bias
Var
MSE

-0.034
0.008
0.009

-0.101
0.005
0.016

-0.084
0.009
0.016

-0.102
0.008
0.018

-0.009
0.007
0.007

n

fS
-0.216
0.030
0.077

fA
-0.175
0.054
0.085

fH
-0.189
0.028
0.063

fB
-0.217
0.028
0.071

Table 1: Simulated bias, variance and mean squared error for the five estimators from standard
exponential distribution.
Table 2 lists the bias, variance and MSE values of the estimators with observations sampled
from the half-normal distribution. For the half-normal data, the histogram estimator f H
behaves best in the MSE sense. The posterior mean has smaller variance than the others,
except for the big sample size n = 1000. The penalized MLE is the most biased.
17

1.6

n=1000

1.2

0.8

0.4
1.6

1.2
n=200

0.8

0.4
1.6

1.2
n=50

0.8

0.4
P

S

A

H

B

Estimator
Estimator

P

S

A

H

B

Figure 2: Boxplots based on 50 replications, where a sample of size n is drawn from the standard
exponential distribution. The rows correspond to the sample sizes n = 50, 200 and 1000; the
columns to the different estimators. Here P ,S,A,H,B correspond to the penalized maximum
likelihood-, simple-, adaptive-, histogram- and posterior mean- estimator respectively. The
horizontal lines indicate the true value of f0 (0) = 1.

18

1.2

n=1000

1.0

0.8

0.6

0.4

1.2

1.0
n=200

0.8

0.6

0.4

1.2

1.0
n=50

0.8

0.6

0.4
P

S

A

H

B

Estimator
Estimator

P

S

A

H

B

Figure 3: Boxplots based on 50 replications, where a sample of size n is drawn from the half
normal distribution. The rows correspond to the sample sizes n = 50, 200 and 1000; the
columns to the different estimators. Here P ,S,A,H,B correspond to the penalized maximum
likelihood-, simple-, adaptive-, histogram- and posterior
mean- estimator respectively. The
p
horizontal lines indicate the true value of f0 (0) = 2/π.

19

50

Bias
Var
MSE

fP
0.073
0.020
0.025

200

Bias
Var
MSE

0.075
0.009
0.015

0.039
0.012
0.013

0.035
0.012
0.013

0.001
0.010
0.010

0.015
0.005
0.005

1000

Bias
Var
MSE

0.041
0.006
0.008

0.015
0.003
0.003

0.009
0.002
0.003

-0.014
0.002
0.002

0.023
0.006
0.006

n

fS
-0.019
0.014
0.015

fA
-0.010
0.020
0.021

fH
-0.030
0.011
0.012

fB
-0.078
0.006
0.012

Table 2: Simulated bias, variance and mean squared error for the five estimators based on
samples from the standard half-normal distribution.

5.3

Application to fertility data

In Keiding et al. (2012) data concerning the fertility of a population are analyzed. The aim
is to estimate the distribution of the duration for women to become pregnant from when they
start attempting, based on data from so-called current durations. For more information on the
design of this study see Keiding et al. (2012).
What is important, is that the sampling density of current durations, f0 , is decreasing on
[0, ∞) and that the distribution function of interest H0 can be expressed in terms of f0 via
relation (2). In this section we estimate the distribution function H0 based on the n = 618
measured current durations using our Bayesian approach. In Groeneboom & Jongbloed (2015)
chapter 9, pointwise confidence bands for f0 and H0 are constructed based on the smoothed
maximum likelihood estimator. Having derived the estimators, producing such confidence
bands needs quite some fine tuning. In this section, we construct the Bayesian counterpart of
the confidence bands, credible regions for H0 . Contrary to the frequentist approach, having
the machinery available for computing the posterior mean, the pointwise credible sets can be
obtained directly from the MCMC output. The settings for computing the posterior mean
are as described in Section 5.1, with α = 1 and the mixture of Pareto distributions as base
measure, using as hyperparameters λ = 5 and ξ = 0.1. The result fo ther fertility data is
shown in figure ??.

A
A.1

Appendix
A difficulty in Salomond’s proof for posterior consistency at zero

In this section we point out a problem in the proof of theorem 4 in Salomond (2014). There,
the following condition on g0 is used
e−a1 /θ

t2

. g0 (θ) . e−a2 /θ

t2

for

0 < θ ≤ θ0

(16)

for some constants 1 < a1 ≤ a2 and t2 = 1. Note that this is different from what we assume in
(4) in Theorem 2.1. These conditions is used to show inequality (10). However, (16) does not
20

1.0

0.10

0.8

0.08

0.6

0.06

0.4

0.04

0.2

0.02

0.0

0.00

0

5

10

15

20

25

30

35

0

5

10

(a)

15

20

25

30

35

(b)

Figure 4: 95% pointwise credible sets for probability density function f0 (a) and distribution
function H0 (x) = 1 − f0 (x)/f0 (0) (b) with fertility data.
2

imply ecnn Π(f (0) − f0 (0) ≥ Cηn ) → 0. As shown in the following, changing (16) requiring
that t2 > 3, which seems to be natural, does not solve the problem because then inequality (10)
will not hold anymore. Our proof uses the weaker assumption (4) and posterior consistency is
supported by the new part of the proof.
In Salomond (2014) a different way is used to show that E0 Π(f (0) − f0 (0) > Cηn | X n )
tends to zero. We now review the arguments given there. By monotonicity of both f and f0 ,
for any sequence of positive numbers {hn },
f (0) − f0 (0) = f (0) − f (hn ) + f (hn ) − f0 (0)
Z hn
−1
(f (x) − f0 (x))dx
≤ f (0) − f (hn ) + hn
0
Z hn
≤
θ−1 dGf (θ) + h−1
n L1 (f, f0 ),
0

where L1 (f, f0 ) =

R∞
0

|f (x) − f0 (x)|dx. Therefore,

Π(f (0) − f0 (0) ≥ Cηn | X n ) ≤ Π(An1 | X n ) + Π(An2 | X n ),
where,

An1 =

Z
f:

hn

θ

−1


dGf (θ) ≥ Cηn /2

and An2 = {f : h−1
n L1 (f, f0 ) ≥ Cηn /2}.

0

Upon taking hn = εn /ηn , it follows from theorem 1 in Salomond (2014) that E0 Π[An2 | X n ] → 0
as n → ∞. Using (7), we get for n → ∞,
2

2

E0 Π(An1 | X n ) ≤ P0 (Dn ) + ecnn Π(An1 ) = ecnn Π(An1 ) + o(1).
21

(17)

By Markov’s inequality,
Π(An1 ) ≤ 2(Cηn )

−1

Z

hn

θ

E0

−1

−1

Z

Z

θ−1 g0 (θ)dθ.

0

0

Using (16), we have
Z hn

hn

dGf (θ) ≤ 2(Cηn )

hn

t2

θ−1 e−a2 /θ dθ
0
Z ∞
t
−1 t2 −a2 /hn2
= t hn e
(uhtn2 + 1)−1 e−a2 u du

θ−1 g0 (θ)dθ .

0

0
t

.

2
htn2 e−a2 /hn .

Here we use the substitution u = 1/θt2 − 1/htn2 at the equality sign and the simple bound
t

1/hn2
t
u+1/hn2

≤ 1 at the second inequality. Hence, upon substituting ηn and εn as given in (5) and

using hn = εn /ηn we obtain
Z hn
t2
2
2
θ−1 g0 (θ)dθ . ecnn · ηn−1 · htn2 e−a2 /hn
ecnn · ηn−1
0

. (n/ log n)−(t2 −2)/9 exp(cn1/3 (log n)2/3
− a2 nt2 /9 (log n)−t2 /9 ).
Provided that t2 > 3 (this is essential), the right hand side tends to zero.
The prior mass condition (4) used in theorem 2.1 is different from the one given in (16).
Moreover, as we have shown above, the latter condition should be satisfied with t2 > 3. This
raises the question at which places in his proof the prior mass condition is used. A careful
inspection reveals that this is the case in lemma 2.3 and 2.2. We verify whether these two
lemmas are valid under (16) with t2 > 3. It turns out that lemma 2.3 also holds. However,
lemma 2.2 is not valid anymore.
Assume the lower bound in (16) with t2 > 3. Following the proof of lemma 2.2, we see that
in this case we have
Z
t2
G0 (Ui ) &
e−a1 /θ dθ & n exp(−a1 /(in )t2 )
Ui

for i ≤

m0

and hence, apart from an additive constant,
0

m
X

0

0

log G0 (Ui ) ≥ m log n −

i=1

As

P∞

i=1

i−t2

2
a1 −t
n

m
X

i−t2 .

i=1

converges, we have, up to an additive constant
0

m
X

−t2
log G0 (Ui ) & −1
n log n − n .

i=1

This cannot be lower bounded by −C1 n2n for any t2 > 1 and henceforth in this case the “right”
lower bound on Π(Sn ) cannot be obtained. The method of proof is sharp, in the sense that
t
assuming g0 (θ) & e−a1 /θ 2 will only give the right bound if t2 ∈ (0, 1].
22

A.2

Supplementary proof for lemma 2.2

We review some results needed to define set Nn . By lemma 11 in Salomond (2014), there exists
a discrete measure
P =

N
X

zi ∈ [θ0 , ∞),

pi δzi ,

pi ∈ [0, 1]

and

N
X

i=1

Z

∞

pi =

i=1

f0 (x)dx
θ0

R∞ √
√
2
2
where the number of support points satisfies N . −1
n , such that θ0 ( f0 − fP ) . n . Based
on this P , we can take a maximal 2-separated subset of {z1 , z2 , . . . , zN }. Let P̃ be the discrete
measure on this 2-net with weights p̃i obtained by moving the masses pi to the closest point
in thePsupport of P̃ . Denote yi , i = 1, . . . , m are the support points of P̃ , m ≤ N . −1
n . Then
2 for i 6= j. We have
P̃ = m
p̃
δ
with
|
y
−
y
|>
2
i
j
n
i=1 i yi
Z
|

kfP − fP̃ kL1 =

N
X

pi ψx (zi ) −

i=1

Since for any θ0 < θ1 < θ2 ,
Z
Z
| ψx (θ1 ) − ψx (θ2 ) | dx =
x≤θ1

m
X

p̃j ψx (yj ) | dx.

j=1

Z
+

Z
| ψx (θ1 ) − ψx (θ2 ) | dx

+
θ1 <x≤θ2

x>θ2

= 2(θ2 − θ1 )/θ2 . θ2 − θ1 .
By the construction of P̃ this implies that kfP − fP̃ kL1 .  . Upon taking  = 2n we get
Z ∞ p
q 2
fP − fP̃ ≤ 2kfP − fP̃ kL1 . ε2n .
θ0

We now give the construction of the set Nn . On [0, θ0 ], consider the regular partition
with bθ0 −1
n c points. Define the above discrete measure P̃ support on these points such that
fP̃ (x) = f0 (in ) in interval (in , (i + 1)n ], i = 1, . . . , m0 , m0 = bθ0 −1
n c − 1. Hence we construct
0
a discrete probability P̃ on [0, ∞) by support on in , i = 1, . . . , m and yi , i = 1, . . . , m.

References
Balabdaoui, F., Jankowski, H., Pavlides, M., Seregin, A. and Wellner, J.A. (2011). On the
Grenander estimator at zero. Statist. Sinica 21, p. 873–899.
Bertoin J. (1998). Lévy Processes. Cambridge University Press.
Blackwell, D. and MacQueen J.B. (1973). Ferguson distributions via Polya urn schemes.
Ann. Statist. 1, p. 353-355.
Doss, H. and Sellke T. (1982). The Tails of Probabilities Chosen From A Dirichlet Prior.
Ann. Statist. 10, p. 1302–1355.
Ferguson, T.S. (1973). A Bayesian Analysis of Some Nonparametric Problem. Ann. Statist.
1, p. 209–230.
23

Ghosal, S., Ghosh, J.K. and Van der Vaart, A.W. (2000). Convergence rates of posterior
distributions. Ann. Statist. 28, p. 500-531.
Grenander, U. (1956). On the theory of mortality measurement. II. Skand. Aktuarietidskr.,
39, p. 125-153.
Groeneboom, P. and Jongbloed, G. (2014). Nonparametric estimation under shape constraints. Cambridge University Press.
Groeneboom, P. and Jongbloed, G. (2015). Nonparametric confidence intervals for monotone functions. Ann. Statist. 43, p.2019-2054.
Groeneboom, P., Jongbloed, G. and Wellner, J.A. (2001). Estimation of a convex function:
characterizations and asymptotic theory. Ann. Statist. 29, p. 1653-1698.
Hoffmann, M., Rousseau, J., and Schmidt-Hieber, J. (2015). On adaptive posterior concentration rates. Ann. Statist. 43, p. 2259-2295.
Kulikov, V.N. and Lopuhaä, H.P. (2006). The behavior of the NPMLE of a decreasing
density near the boundaries of the support. Ann. Statist. 34, p. 742-768.
MacEachern, S.N. and Müller, P. (1998). Estimating Mixture of Dirichlet Process Models.
J. Comput. Graph. Statist. 7(2), 223–238.
Meyer, M.C. and Woodroofe, M. (2004). Consistent maximum likelihood estimation of a
unimodal density using shape restrictions. Can. J. Statist. 32, p. 55-100.
Neal, R.M. (2000). Markov Chain Sampling Methods for Dirichlet Process Mixture Models. Journal of Comp. and Graph. Statist. 9, p. 249-265.
Orbanz, P. (2014). Lecture Notes on Bayesian Nonparametrics. Version: May 16, 2014.
Salomond, JB. (2014). Concentration rate and Consistency of the posterior distribution
for selected priors under monotonicity constraints. Electron. J. Statist. 8, p. 1380-1404.
Sethuraman, J. (1994) A constructive definition of Dirichlet priors. Statist. Sinica, 4(2),
639–650.
Slama, R. and Højbjerg Hansen, O.K. and Ducot, B. and Bohet, A. and Sorensen, D.
and Allemand, L. and Eijkemans, M.J. and Rosetta, L. and Thalabard, J.C. and
Keiding, N. and others. (2012). Estimation of the frequency of involuntary infertility
on a nation-wide basis. Human Reproduction 27, p. 1489–1498.
Van der Vaart, A.W. and Ghosal, S. (2017) Fundamentals of Nonparametric Bayesian
Inference. Cambridge Series in Statistical and Probabilistic Mathematics.
Van der Vaart, A.W. (1998). Asymptotic Statistics. Cambridge University Press.
Vardi, Y. (1989). Multiplicative censoring, renewal processes, deconvolution and decreasing density: nonparametric estimation. Biometrika 76, p.751–761.
Watson, G.S. (1971). Estimating functionals of particle size distributions. Biometrika 58,
p. 483-490.
Woodroofe, M. and Sun, J. (1993). A penalized maximum likelihood estimate of f (0+)
when f is non-increasing. Statist. Sinica 3, p. 501–515.

24

