arXiv:1701.07572v1 [] 26 Jan 2017

Adaptive posterior convergence rates in non-linear latent
variable models
Shuang Zhou∗, Debdeep Pati
Department of Statistics, Florida State University
Anirban Bhattacharya
Department of Statistics, Texas A&M University
David Dunson
Department of Statistical Science, Duke University
January 27, 2017

Abstract
Non-linear latent variable models have become increasingly popular in a variety of applications. However, there has been little study on theoretical properties of these models.
In this article, we study rates of posterior contraction in univariate density estimation for
a class of non-linear latent variable models where unobserved U(0, 1) latent variables are
related to the response variables via a random non-linear regression with an additive error.
Our approach relies on characterizing the space of densities induced by the above model
as kernel convolutions with a general class of continuous mixing measures. The literature
on posterior rates of contraction in density estimation almost entirely focuses on finite or
countably infinite mixture models. We develop approximation results for our class of continuous mixing measures. Using an appropriate Gaussian process prior on the unknown
regression function, we obtain the optimal frequentist rate up to a logarithmic factor under
standard regularity conditions on the true density.

Keywords: Bayesian nonparametrics; Density estimation; Gaussian process; One factor
∗

Corresponding author: shuang.zhou@stat.fsu.edu

1

model; Rate of convergence

1

Introduction

Latent variable models are popular in statistics and machine learning for dimension reduction, parsimonious dependence modeling and data visualization. Linear latent variable
models, such as factor models or probabilistic principal components, assume a linear relationship between the observed and latent variables. A number of non-linear latent variable
models have been proposed in the literature for structured dimension reduction and manifold learning; example include the generative tomographic mapping (GTM; [2, 1]) and
the Gaussian process latent variable model (GP-LVM; [12, 11, 13, 5]). These models flexibly model the relationship between observed and latent variables, notably using Gaussian
process (GP) priors. In spite of their empirical success, a general theoretical framework
studying the properties of the induced density of the data after marginalizing out the latent
variables seems lacking.
[10] proposed an NL-LVM (non linear latent variable model) approach for univariate
density estimation in which the response variables are modeled as unknown functions (referred to as the transfer function) of uniformly distributed latent variables with an additive
Gaussian error. Operationally similar to a univariate GP-LVM model, the latent variable
specification allows straightforward posterior computation via conjugate posterior updates.
Since inverse c.d.f. transforms of uniform random variables can generate draws from any
distribution, a prior with large support on the space of transfer functions can approximate draws from any continuous distribution function arbitrarily closely. One can also
conveniently center the non-parametric model on a parametric family by centering the
prior on the transfer function on a parametric class of quantile (or inverse c.d.f.) functions
{Fθ−1 : θ ∈ Θ}. While such centering on parametric guesses can be achieved in Dirichlet

process (DP; [3, 4]) mixture models by appropriate choice of the base measure G0 , posterior
computation becomes complicated unless the base measure is conjugate to the kernel K.

Although there is an increasingly rich literature on asymptotic properties of Bayesian

density estimation, this literature mainly focuses on discrete mixture models that have a
fundamentally different form from the NL-LVM models. Hence, it is unclear what types of
asymptotic properties NL-LVMs have for density estimation, and technical tools developed
in the existing literature cannot be fully utilized to study this problem. Our focus is
2

on closing this gap focusing in particular on studying how the posterior measure for the
unknown density concentrates around the true density f0 as the sample size n increases.
Assuming f0 belongs to a Hölder space of univariate functions with smoothness β, it is
known that the minimax optimal rate of convergence for an estimate of the density is
n−β/(2β+1) . Assuming the prior for the unknown density is induced through a discrete
mixture of exponential power distributions, [9] showed that the posterior measure for the
density concentrates at the optimal rate up to a logarithmic factor. Their result shows
rate adaptivity to any degree of smoothness of the true density f0 , generalizing previous
results, such as posterior consistency [7] or optimal rates for a particular smoothness level
[6]. We seek to obtain an adaptive rate result for a class of NL-LVM models, and in the
process significantly advance technical understanding of this relatively new class of models.
The main contributions of this article are as follows. We provide an accurate characterization of the prior support in terms of kernel convolution with a class of continuous
mixing measures. We provide conditions for the mixing measure to admit a density with
respect to Lebesgue measure and show that the prior support of the NL-LVM is at least as
large as that of DP mixture models. We then develop approximation results for the above
class of continuous mixing measures, and show adaptive convergence rates. This involves
some novel issues and technical details relative to the existing literature.
The rest of the article is organized as follows. We introduce relevant notations and
terminologies in Section 2. To make the article self-contained, we also provide a brief
background on Gaussian process priors. In Section 3, we formulate our assumptions on the
true density f0 and in the following section, we describe the NL-LVM model and relate it
to convolutions. We state our main theorem on convergence rates in Section 5. Section 6
provides auxiliary results and Section 7 proves the main theorem of posterior concentration
rate. We discuss some implications of our results and outlines possible future directions
Section 8.

2

Notations

Throughout the article, Y1 , . . . , Yn are independent and identically distributed with density
f0 ∈ F, the set of all densities on R absolutely continuous with respect to the Lebesgue
measure λ. The supremum and L1 -norm are denoted by
let

·

p,ν

·

∞

and

·

1

, respectively. We

denote the norm of Lp (ν), the space of measurable functions with ν-integrable
3

pth absolute power. For two density functions f, g ∈ F, let h denote the Hellinger distance
R
√
√
defined as h2 (f, g) =
f − g 2,λ = (f 1/2 − g1/2 )2 dλ, K(f, g) the Kullback-Leibler
R
R
divergence given by K(f, g) = log(f /g)f dλ and V (f, g) = log(f /g)2 f dλ. The notation

C[0, 1] is used for the space of continuous functions f : [0, 1] → R endowed with the

supremum norm. For β > 0, we let C β [0, 1] denote the Hölder space of order β, consisting
of the functions f ∈ C[0, 1] that have ⌊β⌋ continuous derivatives with the ⌊β⌋th derivative

f ⌊β⌋ being Lipschitz continuous of order β − ⌊β⌋. The ǫ-covering number N (ǫ, S, d) of a

semi-metric space S relative to the semi-metric d is the minimal number of balls of radius
ǫ needed to cover S. The logarithm of the covering number is referred to as the entropy.
By near-optimal rate of convergence we mean optimal rate of convergence slowed down by
a logarithmic factor.
We
φ(x) =

write
(2π)−1/2

“-”

for

inequality

up

to

a

constant

multiple.

Let

exp(−x2 /2)

denote the standard normal density, and let φσ (x) = (1/σ)φ(x/σ).
R
(i)
Let an asterisk denote a convolution e.g., (φσ ∗ f )(y) = φσ (y − x)f (x)dx, and let φσ ∗ f

denote the i-fold convolution. The support of a density f is denoted by supp(f ).

We briefly recall the definition of the RKHS of a Gaussian process prior; a detailed
review can be found in [14]. A Borel measurable random element W with values in a
separable Banach space (B, · ) (e.g., C[0, 1]) is called Gaussian if the random variable

b∗ W is normally distributed for any element b∗ ∈ B∗ , the dual space of B. The reproducing

kernel Hilbert space (RKHS) H attached to a zero-mean Gaussian process W is defined as

the completion of the linear space of functions t 7→ EW (t)H relative to the inner product
hEW (·)H1 ; EW (·)H2 iH = EH1 H2 ,
where H, H1 and H2 are finite linear combinations of the form
and si in the index set of W .

P

i ai W (si )

with ai ∈ R

Let W = (Wt : t ∈ R) be a Gaussian process with squared exponential covariance kernel.

The spectral measure mw of W is absolutely continuous with respect to the Lebesgue
measure λ on R with the Radon-Nikodym derivative given by
dmw
2
1
(x) = 1/2 e−x /4 .
dλ
2π
Define a scaled Gaussian process W a = (Wat : t ∈ [0, 1]), viewed as a map in C[0, 1].

Let Ha denote the RKHS of W a , with the corresponding norm
the RKHS is denoted

Ha1 .
4

·

Ha

. The unit ball in

Throughout the paper, C, C1 , C2 , . . . denote global constants whose value may change
one line to another.

3

Assumptions on the true density

It has been widely recognized that one needs certain smoothness assumptions and tail
conditions on the true density f0 to derive posterior convergence rates. We make the
following assumptions in our case:
Assumption 3.1. Assume log f0 ∈ C β [0, 1]. Let lj (x) =

dj
dxj

log f0 (x) be the derivatives

for j = 1, . . . , r with r = ⌊β⌋. For any β > 0, assume that there exists a constant L > 0
such that

|lr (x) − lr (y)| ≤ L|x − y|β−r ,

(3.1)

for all x 6= y.
Assumption 3.2. Assume f0 is compactly supported on [a0 , b0 ], for −∞ < a0 < b0 < ∞,

and that there exists some interval [a, b] ⊂ [a0 , b0 ] such that f0 is nondecreasing on [a0 , a],
bounded away from 0 on [a, b] and non-increasing on [b, b0 ].

Assumption 3.1 is useful in simplifying expressions for f0 as convolutions with a given
density, providing a key piece in our theoretical developments. Similar assumption on the
local smoothness appeared in [9], while in our case the global smoothness assumption is
sufficient since f0 is assumed to be compactly supported. Assumption 3.2 guarantees that
for every δ > 0, there exists a constant C > 0 such that f0 ∗ φσ ≥ Cf0 for every σ < δ.

4

The NL-LVM model

Consider the nonlinear latent variable model,
yi = µ(ηi ) + ǫi ,
µ ∼ Πµ ,

ǫi ∼ N(0, σ 2 ), (i = 1, . . . , n)

σ ∼ Πσ ,

ηi ∼ U(0, 1),

(4.1)
(4.2)

where ηi ’s are latent variables, µ ∈ C[0, 1] is a transfer function relating the latent variables

to the observed variables and ǫi is an idiosyncratic error. The density of y conditional on
5

the transfer function µ and scale σ is obtained on marginalizing out the latent variable as
Z 1
def
φσ (y − µ(x))dx.
(4.3)
f (y; µ, σ) = fµ,σ (y) =
0

Define a map g : C[0, 1] × [0, ∞) → F with g(µ, σ) = fµ,σ . One can induce a prior Π

on F via the mapping g by placing independent priors Πµ and Πσ on C[0, 1] and [0, ∞)

respectively, with Π = (Πµ ⊗ Πσ )◦g−1 . [10] assumed a Gaussian process prior with squared
exponential covariance kernel on µ and an inverse-gamma prior on σ 2 .

It is not immediately clear whether the class of densities fµ,σ in the range of g encompass
a large subset of the density space. The following intuition relates the above class with
continuous convolutions which plays a key role in our proofs. Let f0 be a continuous
Rt
density with cumulative distribution function F0 (t) = −∞ f0 (x)dx. Assume f0 to be
non-zero almost everywhere within its support, so that F0 : supp(f0 ) → [0, 1] is strictly

monotone and hence has an inverse F0−1 : [0, 1] → supp(f0 ) satisfying F0 {F0−1 (t)} = t for
all t ∈ supp(f0 ). Letting µ0 (x) = F0−1 (x), one obtains
fµ0 ,σ (y) =

Z

0

1

φσ (y

− F0−1 (x))dx

=

Z

∞
−∞

φσ (y − t)f0 (t)dt,

(4.4)

where the second equality follows from the change of variable theorem. Thus, fµ0 ,σ (y) =
φσ ∗ f0 , i.e., fµ0 ,σ is the convolution of f0 with a normal density having mean 0 and
standard deviation σ. It is well known that the convolution φσ ∗ f0 can approximate f0

arbitrary closely as the bandwidth σ → 0 in the sense that for f0 ∈ Lp (λ) with p ≥ 1,
φσ ∗ f0 − f0

p,λ

→ 0 as σ → 0. For Holder-smooth functions, the order of approximation

can be characterized in terms of the smoothness. If f0 ∈ C β [0, 1] with β ≤ 2, it follows

from standard Taylor series expansion that ||φσ ∗ f0 − f0 ||∞ = O(σ β ). For β > 2, a similar

Taylor series expansion yields a sub-optimal error ||φσ ∗ f0 − f0 ||∞ = O(σ 2 ). In this case,

we can refine the approximation by convoluting with a sequence of functions fj constructed

by the procedure,
fj+1 = f0 − △σ fj ,

△σ fj = φσ ∗ fj − fj ,

For f0 ∈ C β [0, 1] with β ∈ (2j, 2j +2] we have φσ ∗fj −f0

∞

j ≥ 1.

(4.5)

= O(σ β ) [9]. Although the fj s

need not be non-negative in general, we show that they are non-negative on supp(f0 ) when

f0 is compactly supported. It can be additionally shown that the normalizing constant is
6

1 + O(σ β ); let hj denote the density obtained by normalizing fj . We then approximate f0
by φσ ∗ hβ , where hβ = hj for β ∈ (2j, 2j + 2].

Let λ̃ denote the Lebesgue measure on [0, 1], or equivalently, the U[0, 1] distribution.

For any measurable function µ : [0, 1] → R, let νµ denote the induced measure on (R, B),

with B denoting the Borel sigma-field on R. Then, for any Borel measurable set B, νµ (B) =

λ̃(µ−1 (B)), where µ−1 (B) = {x ∈ [0, 1] : µ(x) ∈ B}. By the change of variable theorem
for induced measures,

Z

0

1

φσ (y − µ(x))dx =

Z

φσ (y − t)dνµ (t),

(4.6)

so that fµ,σ can be expressed as a kernel mixture form with mixing distribution νµ . It turns
out that this mechanism of creating random distributions is very general. Depending on the
choice of µ, one can create a large variety of mixing distributions based on this specification.
For example, if µ is a strictly monotone function, then νµ is absolutely continuous with
respect to the Lebesgue measure, while choosing µ to be a step function, one obtains a
discrete mixing distribution. However, it is easier to place a prior on µ supported on the
space of continuous functions C[0, 1] without further shape restrictions and Theorem 4.1
assures us that this specification leads to large L1 support on the space of densities.
Suppose the prior Πµ on µ has full sup-norm support on C[0, 1] so that Pr( µ−µ∗

∞

<

ǫ) > 0 for any ǫ > 0 and µ∗ ∈ C[0, 1], and the prior Πσ on σ has full support on [0, ∞). If

f0 is compactly supported so that the quantile function µ0 ∈ C[0, 1], then it can be shown
that under mild conditions, the induced prior Π assigns positive mass to arbitrarily small
L1 neighborhoods of any density f0 . We summarize the above discussion in the following
theorem, with a proof provided in the appendix.
Theorem 4.1. If Πµ has full sup-norm support on C[0, 1] and Πσ has full support on
[0, ∞), then the L1 support of the induced prior Π on F contains all densities f0 which
have a finite first moment and are non-zero almost everywhere on their support.

Remark 4.1. The conditions of Theorem 4.1 are satisfied for a wide range of Gaussian
process priors on µ (for example, a GP with a squared exponential or Matérn covariance
kernel).
Remark 4.2. When f0 has full support on R, the quantile function µ0 is unbounded near
R1
R
0 and 1, so that µ0 ∞ = ∞. However, 0 |µ0 (t)| dt = R |x| f0 (x)dx, which implies that
7

µ0 can be identified as an element of L1 [0, 1] if f0 has finite first moment. Since C[0, 1] is
dense in L1 [0, 1], the previous conclusion regarding L1 support can be shown to hold in
the non-compact case too.

5

The main theorem

We consider the case where f0 satisfies Assumption 3.1 and Assumption 3.2. For β-Hölder
density f0 , we consider density hβ as defined after expression (4.5). Denote µ0 the quantile
: [0, 1] → [a0 , b0 ], a continuous monotone function inheriting the smoothness
function Fh−1
β

of hβ . Note that hβ has the same smoothness of f0 based on the construction of fj , therefore
with the fundamental theorem of calculus it is easy to see that µ0 ∈ C β+1 [0, 1].
We now mention our choices for the prior distributions Πµ and Πσ .

Assumption 5.1. We assume µ follows a centered and rescaled Gaussian process denoted
by GP(0, cA ), where A denotes the rescaled parameter, and assume A is a density satisfying
for a > 0,
C1 ap exp (−D1 a logq a) ≤ g(a) ≤ C2 ap exp (−D2 a logq a)

(5.1)

for positive constants C1 , C2 , D1 , D2 , nonnegative constant p and q, and every sufficiently
large a > 0.
Assumption 5.2. We assume σ ∼ IG(aσ , bσ ).
Note that contrary to the usual conjugate choice of an inverse-Gamma prior for σ 2 , we
have assumed an inverse-Gamma prior for σ. This enables one to have slightly more prior
mass near zero compared to an inverse-Gamma prior for σ 2 , leading to the optimal rate of
posterior convergence. Refer also to [9] for a similar prior choice for the bandwidth of the
kernel in discrete location-scale mixture priors for densities.
We state below the main theorem of posterior convergence rates.
Theorem 5.1. If f0 satisfies Assumption 3.1 and the priors Πµ and Πσ are as in Assumptions 5.1 and 5.2 respectively, the best obtainable rate of posterior convergence relative to
Hellinger metric h is
ǫn = max(ǫ̃n , ǭn ),
8

(5.2)

where ǫ̃n = n

β
− 2β+1

(log n)t1 , ǭn = n

β
− 2β+1

q)/(2β + 1), t2 = t1 + 1.

(log n)t2 , with nonnegative constants t1 = β(2 ∨

Unlike the treatment in discrete mixture models [6] where a compactly supported density is approximated with a discrete mixture of normals, the main trick here is to approximate the true density f0 by the convolution φσ ∗ f0 and allow the prior on the transfer

function to appropriately concentrate around the true quantile function µ0 ∈ C[0, 1].

6

Auxiliary results

To guarantee that the above scheme leads to the optimal rate of convergence, we first derive
sharp bounds for the Hellinger distance between fµ1 ,σ1 and fµ2 ,σ2 for µ1 , µ2 ∈ C[0, 1] and

σ1 , σ2 > 0. We summarize the result in the following Lemma 6.1.

Lemma 6.1. For µ1 , µ2 ∈ C[0, 1] and σ1 , σ2 > 0,
s

2 
µ1 − µ2 ∞
2σ1 σ2
2
h (fµ1 ,σ1 , fµ2 ,σ2 ) ≤ 1 −
exp −
.
σ12 + σ22
4(σ12 + σ22 )

(6.1)

Proof. Note that by Hölder’s inequality,
2
Z 1
p
p
φσ1 (y − µ1 (x)) φσ2 (y − µ2 (x))dx .
fµ1 ,σ1 (y)fµ2 ,σ2 (y) ≥
0

Hence,

2

Z Z

1

Z

1

φσ2 (y − µ2 (x))dx

Z 1p
p
φσ1 (y − µ1 (x)) φσ2 (y − µ2 (x))dx dy.
− 2

h (fµ1 ,σ1 , fµ2 ,σ2 ) ≤

0

φσ1 (y − µ1 (x))dx +

0

0

By changing the order of integration (applying Fubini’s theorem since the function within
the integral is jointly integrable) we get,
Z 1
2
h2 (fµ1 (x),σ1 , fµ2 (x),σ2 )dx
h (fµ1 ,σ1 , fµ2 ,σ2 ) ≤
0
s


Z 1
2σ1 σ2
(µ1 (x) − µ2 (x))2
=
1−
dx
exp
−
σ12 + σ22
4(σ12 + σ22 )
0
s

2 
µ1 − µ2 ∞
2σ1 σ2
exp −
.
≤ 1−
σ12 + σ22
4(σ12 + σ22 )
9

Remark 6.1. When σ1 = σ2 = σ, h2 (fµ1 ,σ , fµ2 ,σ ) ≤ 1 − exp

implies that h2 (fµ1 ,σ , fµ2 ,σ ) - µ1 −

2
µ2 ∞ /σ 2 .

Remark 6.2. The standard inequality h2 (fµ1 ,σ1 , fµ2 ,σ2 ) ≤



µ1 − µ2

2
/8σ 2
∞

fµ1 ,σ1 − fµ2 ,σ2

1

, which

relating the

Hellinger distance to the total variation distance leads to the cruder bound
h2 (fµ1 ,σ1 , fµ2 ,σ2 ) ≤ C1
which is linear in µ1 − µ2

∞

µ1 − µ2 ∞
|σ2 − σ1 |
+ C2
,
(σ1 ∧ σ2 )
(σ1 ∧ σ2 )

. This bound is less sharp than what is obtained in Lemma

6.1 and does not suffice for obtaining the optimal rate of convergence.
To control the Kullback-Leibler divergence between the true density f0 and the model
fµ,σ , we derive an upper bound for log

f0
fµ,σ ∞

in Lemma 6.2.

Lemma 6.2. If f0 satisfies Assumption 3.2,
f0
log
fµ,σ

∞

µ − µ0
≤C+
σ2

2
∞

(6.2)

for some constant C > 0.
Proof. Note that
fµ,σ (y) =
≥
≥
≥

1


(y − µ(x))2
dx
exp −
2σ 2
0



Z 1
µ − µ0
(y − µ(x))2
1
√
dx exp −
exp −
2
σ
σ2
2πσ 0

2 
µ − µ0 ∞
√
Cφσ/ 2 ∗ f0 (y) exp −
σ2


2
µ − µ0 ∞
,
Cf0 (y) exp −
σ2
1
√
2πσ

Z



2 
∞

where the last inequality follows from Lemma 6 of [6] since f0 is compactly supported by
Assumption 3.2. This provides the desired inequality.
Lemma 6.3. For β ∈ (2j, 2j + 2], j ≥ 0 and fj constructed by 4.5, we have the expression
 (i)
P
fj = ji=0 (−1)i j+1
i+1 φσ f0 .
10

The proof can be found in Appendix A.2. The expression of fj as a linear combination
of true density and the folded convolutions indicates that fj is as smooth as f0 . One can
get infinitely differentiable function by convoluting with the kernel, so for the true density
with higher regularity degree, we need add the ”smoother” function into the approximation
fj to ensure the approximation error remains of order O(σ β ).
Lemma 6.4. For any β > 0, let f0 satisfy Assumption 3.1 and 3.2, integer j be such that
for β ∈ (2j, 2j + 2], fj constructed by (4.5). For any constant L and all x ∈ [a0 , b0 ], we
have

φσ ∗ fβ (x) = f0 (x)(1 + O(D(x)σ β )),
where
D(x) =

r
X
i=1

β

ci |lj (x)| i + cr+1 ,

for nonnegative constants ci , i = 1, . . . , r, and cr+1 a multiple of L.
Proof. Following the proof of Lemma 1 in [9], for any x, y ∈ [a0 , b0 ],
log f0 (y) ≤ log f0 (x) +

r
X
lj (x)

(y − x)j + L|y − x|β ,

log f0 (y) ≥ log f0 (x) +

r
X
lj (x)

(y − x)j − L|y − x|β .

i=1

i=1

j!

j!

Define
Bfu0 ,r (x, y) =

r
X
lj (x)

(y − x)j + L|y − x|β ,

r
X
lj (x)

(y − x)j − L|y − x|β .

i=1

Bfl 0 ,r (x, y) =

i=1

j!

j!

Then we have
Bfu

e

0 ,r

≤ 1 + Bfu0 ,r +

1 u 2
(B ) + · · · + M |Bfu0 ,r |r+1 ,
2! f0 ,r
11

(6.3)

Bfl

e

0 ,r

≥ 1 + Bfl 0 ,r +

1 l 2
(B ) + · · · − M |Bfl 0 ,r |r+1 .
2! f0 ,r

where
M=

1
exp
(r + 1)!



sup
x,y∈[a0 ,b0 ],x6=y

(|

r
X
lj (x)
j=1

j!


(y − x)j | + L|y − x|β ) .

Note that f0 is bounded on [a0 , b0 ], we consider the convolution on the whole real line by
extending f0 analytically outside [a0 , b0 ]. For β ∈ (1, 2], r = 1 and x ∈ (a0 , b0 ),
Z
B u (x,y)
φσ (y − x)dy
φσ ∗ f0 (x) ≤ f0 (x) e f0 ,r
Z
≤ f0 (x) φσ (y − x)[1 + L|y − x|β + M (l12 (x)(y − x)2 + L2 |y − x|2β )]dy.
R

Since lj (x)’s are all continuous on [a0 , b0 ], there exist finite constants Mj such that |lj | ≤ Mj

and |y − x| ≤ |b0 − a0 |. The integral in the last inequality can be bounded by
Z
φσ (y − x)[1 + L|y − x|β + M ((M1 |b0 − a0 |)2−β |l1 (x)(y − x)|β + (L2 |b0 − a0 |β )|y − x|β )]dy
R

Therefore,
φσ ∗ f0 (x) ≤ f0 (x)(1 + (r1 |l1 (x)|β + r2 )σ β ),
where r1 = M (M1 |b0 − a0 |)2−β µβ , r2 = (L + M L2 )µβ .

In the other direction,
Z
φσ ∗ f0 (x) ≥ f0 (x) φσ (y − x)[1 − L|y − x|β − M (l12 (x)(y − x)2 + L2 |y − x|2β )]dy.
Thus we achieve (6.3).
(2)

For any β > 2, suppose β ∈ (2j, 2j + 2], j > 1. First we calculate φσ ∗ f0 , φσ ∗ f0 , . . . ,

(j)

(i)

φσ ∗ f0 (x), by Lemma 6.3 to get φσ ∗ fβ (x). The calculation of φσ ∗ f0 (x) is the same

as φσ ∗ f0 (x) except taking the convolution with φ√iσ . The terms σ 2 , σ 4 , . . . , σ 2j caused
(i)

by the factors containing |y − x|k , k < β in φσ f0 can be canceled out by Lemma 6.3. For
terms containing |y − x|k , k ≥ β, we take out |y − x|β and bound the rest by a certain
power of |lj (x)| or some constant.

12

Lemma 6.5. Let f0 satisfy Assumption 3.1 and 3.2. With Aσ = {x : f0 (x) ≥ σ H } , we

have

Z

2β

Acσ

f0 (x)dx = O(σ ),

Z

Acσ

φσ ∗ fj (x)dx = O(σ 2β ),

(6.4)

for all non-negative integer j, sufficiently small σ and sufficiently Large H.
Proof. Under Assumption 3.2 there exists (a, b) ⊂ [a0 , b0 ] such that Acσ ⊂ [a0 , a) ∪ (b, b0 ]
R
if we choose σ sufficiently small, so that f0 (x) ≤ σ H for x ∈ Acσ . Therefore, Ac f0 (x) ≤

σ H |b0 − a0 | ≤ O(σ 2β ) if we choose H ≥ 2β. Using Lemma 6.4,
Z
Z
f0 (x)(1 + O(D(x)σ β )) ≤ O(σ H ),
φσ ∗ fj (x)dx =

σ

Acσ

Acσ

with bounded D(x) and H ≥ 2β it is easy to bound the second integral by O(σ 2β ).
Lemma 6.6. Suppose f0 (x) satisfies Assumption 3.1 and 3.2. For β > 2 and j such that
β ∈ (2j, 2j + 2], we can construct the density hβ from (4.5) and show that hβ satisfies
Lemma 6.4 and Lemma 6.5.

Proof. To get the density function we first show that fj is nonnegative and compute the
R
normalizing constant fj (x) = 1 + O(σ β ). Following the proof of Lemma 2 in [9], we treat
log f0 as a function in C 2 [0, 1] and obtain the same form of φσ ∗ f0 as (6.3). For small

enough σ we can find ρ1 ∈ (0, 1) very close to 0 such that

φσ ∗ f0 (x) = f0 (x)(1 + O(D (2) (x)σ 2 )) < f0 (x)(1 + ρ1 ),
where D (2) contains |l1 (x)| and |l2 (x)| with certain power, so D (2) is bounded. Then we

have

f1 (x) = 2f0 (x) − Kσ f0 (x) > 2f0 (x) − f0 (x)(1 + ρ1 ) = f0 (x)(1 − ρ1 ).
Then we treat log f0 as a function with β = 4, j = 1. Similarly, we can get
φσ ∗ f1 (x) = f0 (x)(1 + O(D (4) (x)σ 4 )),
where D (4) contains |l1 (x)|, . . . .|l4 (x)|. We can find 0 < ρ2 < ρ1 such that φσ ∗ f1 (x) <

f0 (x)(1 + ρ2 ), then can get

f2 (x) = f0 (x) − (φσ ∗ f1 (x) − f1 (x)) > f0 (x)(1 − ρ1 − ρ2 ) > f0 (x)(1 − 2ρ1 ).
13

Continuing this procedure, we can get fj (x) > f0 (x)(1 − jρ1 ), with sufficiently small σ,
and 1 − jρ1 ∈ (0, 1) but very close to 1. Obviously fj is nonnegative.

Now we calculate the normalizing constant for fj . When β < 2, with Lemma 6.3,
Z
Z
Z
Z
f1 (x) = f0 − (φσ ∗ f0 − f0 ) ≤ f0 + | (φσ ∗ f0 − f0 )| ≤ 1 + O(σ β ).
For β ∈ (2, 4],
Z
Z
Z
Z
f2 (x) = f0 − (φσ ∗ f1 − f1 ) ≤ f1 + | (φσ ∗ f1 − f0 )| ≤ 1 + O(σ β ).
Then by induction, we have

R

fj = 1 + O(σ β ), so that we have the density

hβ =

fj
,
1 + O(σ β )

β ∈ (2j, 2j + 2].

(6.5)

Now to show hβ satisfying (6.3), note that
φσ ∗ hβ =

φσ ∗ fj
f0 (x)(1 + O(D(x)σ β ))
R
.
=
1 + O(σ β )
fj

Since D(x) is bounded and for sufficiently small σ we consider term
(1 + O(D(x)σ β ))
(1 + O((D(x) + 1)σ β ) + O(D(x)σ 2β ))
= 1 + O(D(x)σ β ),
1 + O(σ β )
1 + O(σ β )
which directly leads to the same form as (6.3), and obviously Lemma 6.5 is satisfied.
Lemma 6.7. Let f0 satisfy Assumption 3.1 and 3.2, and integer j be such that β ∈
(2j, 2j + 2]. Then we can show that the density hβ defined by (6.5) satisfies,
Z
f0 (x)
= O(σ 2β ),
f0 (x) log
φσ ∗ hβ (x)

(6.6)

for sufficiently small σ and all x ∈ [a0 , b0 ].
Proof. Again consider the set Aσ = {x : f0 (x) ≥ σ H } with arbitrarily large H. We separate
the Kullback-Leibler divergence into
Z
Z
Z
Z
(f0 − φσ ∗ hβ )2
f0
f0
f0 log
(φσ ∗ hβ − f0 ).
≤
+
+
f0 log
φσ ∗ hβ
φσ ∗ hβ
φσ ∗ hβ
Aσ
Acσ
Acσ
[a0 ,b0 ]

Under Assumption 3.2 and by Remark 3 in [7], for small enough σ there exists constant C
such that for all x ∈ [a0 , b0 ], φσ ∗ f0 ≥ Cf0 , especially on set Acσ f0 satisfies φσ ∗ f0 ≥ f0 /3.
14

Also in the proof of Lemma 6.6 we can find ρ ∈ (0, 1) such that fβ > ρf0 . Then we have

on set Aσ with sufficiently small σ
φσ ∗ hβ =

φσ ∗ fβ
ρφσ ∗ f0
≥
≥ Kf0 ,
β
1 + O(σ )
1 + O(σ β )

for some positive and finite constant K. Applying Lemma 6.4, the first integral on the
right side of 6.6 can be bounded by
Z
Z
(f0 − φσ ∗ hβ )2
[f0 (x) − f0 (x)(1 + O(D(x)σ β ))]2
≤
φσ ∗ hβ
K1 f0 (x)
Aσ
ZAσ
f0 (x)O(D 2 (x)σ 2β ) = O(σ 2β ).
Aσ

To bound the second integral of r.h.s again by Remark 3 in [7] we get φσ ∗hβ ≥

ρ
f ,
3(1+O(σβ )) 0

so easily we can find a constant C < 1 such that φσ ∗ hβ ≥ Cf0 . With Lemma 6.5 clearly

the second and third term can be bounded by O(σ 2β ).

Lemma 6.8. Let Ha1 denote the unit ball of RKHS of the Gaussian process with rescaled
parameter a and B1 be the unit ball of C[0, 1]. For r > 1, there exists a constant K, such
that for ǫ < 1/2,
log N (ǫ, ∪a∈[0,r] Ha1 ,



1
· ∞ ) ≤ Kr log
ǫ

2

.

(6.7)

Proof. Since we can write any element of Ha1 as a function of Re(z) by [15] which can
be analytically extended to some interval containing Ωa = {z ∈ C : | Re(z)| ≤ R} with

δ
R = δ/(6 max(a, 1)), so for any h ∈ Ωa , |2a Re(z)| ≤ | 6 max(a,1)
· 2a| = δ/3. Consider any b

with |b − a| ≤ 1, we can show that any element of Hb1 can be extended analytically to Ωa

noting that for z ∈ Ωa related to the maximum norm,

2b Re(z) ≤ 2a Re(z) + 2(b − a) Re(z) ≤

δ
2
+ 2 Re(z) ≤ δ.
3
3

Therefore, F a forms one ǫ-net over Hb1 . We find one set Γ = {ai , i = 1, . . . , k} with

k = ⌊r⌋ + 1 and ak = r, such that for any b ∈ [0, r] there exists some ai satisfying

|b − ai | ≤ 1, so that ∪i≤k F ai forms an ǫ-net over ∪a≤r Ha1 . Since the covering number of
∪i≤k F ai is bounded by summation of covering number of F ai , we obtain
log N (ǫ, ∪a∈[0,r] Ha1 ,

·

∞

) ≤ log

X
k
i=1




1 2
r
#(F ) ≤ log(k · #(F )) ≤ Kr log
.
ǫ
ai

15

To complete the proof, note that the piecewise polynomials are constructed on the partition of [0, 1], ∪i≤m Bi , where Bi ’s are disjoint interval with length shorter than R =

δ/(6 max(a, 1)), so the total number of polynomials is a non-decreasing in a. Also we find
that when building the mesh grid of the coefficients of polynomials in each Bi , both the
approximation error and tail estimate are invariant to interval length R, therefore we have
#(F a ) ≤ #(F b ) if a ≤ b, for a, b ∈ [0, r].

Remark 6.3. With larger a we need a finer partition on [0, 1] while the grid of coefficients of
piecewises polynomial remains the same except the range and the meshwidth will change
together along with a. Since we can see the element h of RKHS ball as a function of it and
with Cauchy formula we can bound the derivatives of h by C/Rn , where |h|2 ≤ C 2 .

7

Proof of the main theorem

Proof of Theorem 5.1: Following [8], we need to find sequences ǭn , ǫ̃n → 0 with

n min{ǭ2n , ǫ̃2n } → ∞ such that there exist constants C1 , C2 , C3 , C4 > 0 and sets Fn ⊂ F so

that,

log N (ǫn , Fn , d) ≤ C1 nǭ2n

(7.1)

Π(Fnc ) ≤ C3 exp{−nǫ̃2n (C2 + 4)}




Z
Z
f0
f0 2
2
2
Π fµ,σ : f0 log
≤ ǫ̃n , f0 log
≤ ǫ̃n ≥ C4 exp{−C2 nǫ̃2n }.
fµ,σ
fµ,σ

(7.2)
(7.3)

Then we can conclude that for ǫn = max{ǭn , ǫ̃n } and sufficiently large M > 0, the posterior
probability

Πn (fµ,σ : d(fµ,σ , f0 ) > M ǫn |Y1 , . . . , Yn ) → 0 a.s. Pf0 .
We consider the Gaussian process µ ∼ W A given A, with A satisfying Assumption 5.1.
We will first verify (7.3) along the lines of [6]. Note that

h2 (f0 , fµ,σ ) - h2 (f0 , fµ0 ,σ ) + h2 (fµ0 ,σ , fµ,σ ).
Since fµ0 ,σ = φσ ∗ hβ , by Lemma 6.7, one obtains under Assumptions 3.1 and 3.2,


Z
f0
2
- O(σ 2β ).
h (f0 , fµ0 ,σ ) ≤ f0 log
fµ0 ,σ
16

(7.4)

(7.5)

From Lemma 6.1 and the following remark, we obtain
h2 (fµ0 ,σ , fµ,σ ) -

µ − µ0
σ2

2
∞

.

From Lemma 8 of [6], one has



Z
f0 i
f0
2
≤ h (f0 , fµ,σ ) 1 + log
f0 log
fµ,σ
fµ,σ

(7.6)

∞

i

(7.7)

for i = 1, 2.
From (7.4)-(7.7), for any b ≥ 1 and ǫ̃2n = σn2β ,
σ ∈ [σn , σn + σnb ], µ − µ0 ∞ - σnβ+1 ⊂



Z
Z
f0 2
f0
2β
2β
- σn .
- σn , f0 log
f0 log
fµ,σ
fµ,σ


Since µ0 ∈ C β+1 [0, 1], from Section 5.1 of [15],
P( µ − µ0

1

∞

≤ 2δn ) ≥ C4 exp{−C5 (1/δn ) β+1 log(

1 2∨q
) }(C6 /δn )p+1/β+1 ,
δn

for δn → 0 and constants C4 , C5 , C6 > 0. Letting δn = σn3 , we obtain

 


1
1
2∨q
log
P( µ − µ0 ∞ ≤ 2δn ) ≥ exp − C7
,
σn
σnβ+1
for some constant C7 > 0. Since σ ∼ IG(aσ , bσ ), we have
Z 2σn
baσσ
P(σ ∈ [σn , 2σn ]) =
x−(aσ +1) e−bσ /x dx
Γ(aσ ) σn
Z 2σn
baσσ
e−2bσ /x dx
≥
Γ(aσ ) σn
baσσ
σn exp{−bσ /σn }
≥
Γ(aσ )
≥ exp{−C8 /σn },
for some constant C8 > 0. Hence
P{σ ∈ [σn , 2σn ], µ − µ0

∞

-

σnβ+1 }




1
1
2∨q
exp{−C8 /σn }
≥ exp − C7
log
σn
σnβ+1

 


1
1
≥ exp − 2C7
.
log2∨q
σn
σnβ+1


17



Then (7.3) will be satisfied with ǫ̃n = n

β
− 2β+1

β(2∨q)
2β+1

logt1 (n), where t1 =

and some C9 > 0.

Next we construct a sequence of subsets Fn such that 7.1 and 7.2 are satisfied with
β

ǭn = n− 2β+1 logt2 n and ǫ̃n for some global constant t2 > 0.

Letting B1 denote the unit ball of C[0, 1] and given positive sequences Mn , rn , define
Bn = ∪a<rn (Mn Ha1 ) + δ̄n B1
as in [15], with δ̄n = ǭn ln /K1 , K1 = 2(2/π)1/2 and let
Fn = {fµ,σ : µ ∈ Bn , ln < σ < hn }.
First we need to calculate N (ǭn , Fn , ·
fµ1 ,σ1 − fµ2 ,σ2

1

1

). Observe that for σ2 > σ1 >

 1/2
µ1 − µ2
2
≤
π
σ1

∞

+

σ2
2 ,

3(σ2 − σ1 )
.
σ1

n = l (1 + κ )m , m ≥ 0, we obtain a partition of [l , h ]
Taking κn = min{ ǭ6n , 1} and σm
n
n
n n

n
n with
as ln = σ0n < σ1n < · · · < σm
< hn ≤ σm
n
n −1


1
hn
+ 1.
mn = log
ln log(1 + κn )

One can show that

n −σ n
3(σm
m−1 )
n
σm−1

(7.8)

= 3κn ≤ ǭn /2. Let {µ̃nk , k = 1, . . . , N (δ̄n , Bn , ·

δ̄n -net of Bn . Now consider the set

n
) : k = 1, . . . , N (δ̄n , Bn , ·
{(µ̃nk , σm

∞

), 0 ≤ m ≤ mn }.

n ) such that µ − µ̃n
Then for any f = fµ,σ ∈ Fn , we can find (µ̃nk , σm
k

if one has σ ∈

n ],
n
, σm
(σm−1

then

n
fµ,σ − fµnk ,σm

1

∞

∞

)} be a

(7.9)
< δ̄n . In addition,

≤ ǭn .

Hence the set in (7.9) is an ǭn -net of Fn and its covering number is given by
mn N (δ̄n , Bn , ·

∞

).

From the proof of Theorem 3.1 in [15], for any Mn , rn with rn > a0 , we obtain
log N (2δ̄n , Bn , ·

∞



) ≤ K2 rn log
18



Mn
δ̄n

2

.

(7.10)

Again from the proof of Theorem 3.1 in [15], for rn > 1 and for Mn2 > 16K3 rn (log(rn /δ̄n ))2 ,
we have
P(W A ∈
/ Bn ) ≤

K4 rnp e−K5 rn log
K5 logq rn

q

rn

+ exp{−Mn2 /8}

(7.11)

for constants K3 , K4 , K5 > 0.
Next we calculate P (σ ∈
/ [ln , hn ]). Observe that
−1
P(σ ∈
/ [ln , hn ]) = P (σ −1 < h−1
> ln−1 )
n ) + P (σ
−1
Z ∞
∞
k
X
e−bσ hn (bσ h−1
baσ
n )
≤
e−bσ x/2 dx
+ σ
k!
Γ(aσ ) ln−1
k=ασ

≤ e−aσ log(hn ) +

baσσ −bσ ln−1 /2
e
.
Γ(aσ )

(7.12)

Thus with hn = O(exp{n1/2β+1 (log n)2t1 }), ln = O(n−1/2β+1 (log n)−2t1 ), rn = O(n1/2β+1 (log n)2t1 ), Mn =

O(n1/2β+1 (log n)t1 +1 ), (7.11) and (7.12) implies

Π(Fnc ) = exp{−K6 nǫ̃2n }
for some constant K6 > 0 guaranteeing that (7.2) is satisfied with ǫ̃n = n−β/2β+1 (log n)t1 .
Also with ǭn = n−β/2β+1 (log n)t1 +1 , it follows from (7.8) and (7.10) that
log N (ǭn , Fn , ·

1

) ≤ K7 n1/2β+1 (log n)2t1 +2

for some constant K7 > 0.
Hence max{ǭn , ǫ̃n } = n−β/2β+1 (log n)t1 +1 .

8

Discussion

Non-linear latent variable models offer a flexible modeling framework in a broad variety of
problems and improved practical performance has been demonstrated by [12, 11, 13, 5, 10]
among others. The univariate density estimation model studied here can be extended to
multivariate density estimation, latent factor models and density regression problems.
When the density is compactly supported, the quantile function is a continuous function on [0, 1]. Hence one can use standard results on concentration bounds for Gaussian
processes [14]. However, for densities supported on R, the results fail as the corresponding
19

quantile function is unbounded near zero and one. In this case, assumptions on the tail
behavior of the true density as well as a careful analysis on the behavior of the corresponding quantile function near boundary are required. We propose to address this problem
elsewhere.

A

Appendix

A.1

Proof of Theorem 4.1

Proof. Let f0 be a density with quantile function µ0 that satisfies the conditions of Theorem
R1
R∞
4.1. Observe that µ0 1 = t=0 |µ0 (t)| dt = −∞ |z| f0 (z)dz < ∞ since f0 has a finite first
moment, and thus µ0 ∈ L1 [0, 1]. Fix ǫ > 0. We want to show that Π{Bǫ (f0 )} > 0, where
Bǫ (f0 ) = {f :

f − f0

1

< ǫ}.

Note that µ0 ∈
/ C[0, 1], so that pr( µ − µ0

∞

< ǫ) can be zero for small enough ǫ. The

main idea is to find a continuous function µ˜0 close to µ0 in L1 norm and exploit the fact
that the prior on µ places positive mass to arbitrary sup-norm neighborhoods of µ˜0 . The
details are provided below.
Since φσ ∗ f0 − f0

1

→ 0 as σ → 0, find σ1 such that φσ ∗ f0 − f0

1

< ǫ/2 for σ < σ1 .

Pick any σ0 < σ1 . Since C[0, 1] is dense in L1 [0, 1], for any δ > 0, we can find a continuous
function µ˜0 such that µ0 − µ˜0

1

constant C. Thus, for δ = ǫσ0 /4,


< δ. Now, fµ,σ − fµ˜0 ,σ

fµ,σ : σ0 < σ < σ1 , µ − µ˜0

since f0 − fµ,σ

1

< f0 − fµ0 ,σ

Thus, Π{Bǫ (f0 )} > pr( µ − µ˜0

1
∞

∞

1


< δ ⊂ fµ,σ :

+ fµ0 ,σ − fµ˜0 ,σ

1

≤ C µ − µ˜0
f0 − fµ,σ

+ fµ˜0 ,σ − fµ,σ

1

1

1

/σ for a global

<ǫ ,

and fµ0 ,σ = φσ ∗ f0 .

< δ) pr(σ0 < σ < σ1 ) > 0, since Πµ has full sup-norm

support and Πσ has full support on [0, ∞).

20

A.2

Proof of Lemma 6.3

Proof. Consider fj constructed by (4.5). When i = 1, f1 = 2f0 − φσ ∗ f0 , so the form holds.
By induction, suppose this form holds for j > 1, then

fj+1 = f0 − (φσ ∗ fj − fj )




j
j
X
X
i j+1
i+1 j + 1
(i+1)
(−1)
(−1)
φ(i) ∗ f0
φ
∗ f0 +
= f0 +
i+1 σ
i+1 σ
i=0
i=0




j+1
j
X
X
i j+1
i j+1
(i)
(−1)
(−1)
= (j + 2)f0 +
φ ∗ f0 +
φ(i)
σ ∗ f0
i+1 σ
i
i=1
i=1

 

j
X
j+1
j +1
i
j+1 (i+1)
(−1)
= (j + 2)f0 +
+
φ(i)
φσ
∗ f0
σ ∗ f0 + (−1)
i+1
i
i=1


j
X
i j+2
(−1)
= (j + 2)f0 +
φ(i) ∗ f0 + (−1)j+1 φ(i+1)
∗ f0
σ
i+1 σ
i=1


j+1
X
i j+2
(−1)
φ(i) ∗ f0 .
=
i+1 σ
i=0

It holds for j + 1, which completes the proof.

References
[1] Christopher M Bishop, Markus Svensén, and Christopher KI Williams. Developments
of the generative topographic mapping. Neurocomputing, 21(1):203–224, 1998.
[2] Christopher M Bishop, Markus Svensén, and Christopher KI Williams. Gtm: The
generative topographic mapping. Neural computation, 10(1):215–234, 1998.
[3] T.S. Ferguson. A Bayesian analysis of some nonparametric problems. The Annals of
Statistics, 1(2):209–230, 1973.
[4] T.S. Ferguson. Prior distributions on spaces of probability measures. The Annals of
Statistics, 2(4):615–629, 1974.
[5] B. Ferris, D. Fox, and N. Lawrence. Wifi-slam using gaussian process latent variable models. In Proceedings of the 20th International Joint Conference on Artificial
Intelligence, pages 2480–2485, 2007.
21

[6] S. Ghosal and A. van der Vaart. Posterior convergence rates of Dirichlet mixtures at
smooth densities. The Annals of Statistics, 35(2):697–723, 2007.
[7] S. Ghosal, JK Ghosh, and RV Ramamoorthi. Posterior consistency of Dirichlet mixtures in density estimation. The Annals of Statistics, 27(1):143–158, 1999.
[8] S. Ghosal, J.K. Ghosh, and A.W. van der Vaart. Convergence rates of posterior
distributions. Annals of Statistics, 28(2):500–531, 2000.
[9] W. Kruijer, J. Rousseau, and A. van der Vaart. Adaptive bayesian density estimation
with location-scale mixtures. Electronic Journal of Statistics, 4:1225–1257, 2010.
[10] Suprateek Kundu and David B Dunson. Latent factor models for density estimation.
Biometrika, 101(3):641–654, 2014.
[11] N. Lawrence. Probabilistic non-linear principal component analysis with gaussian
process latent variable models. The Journal of Machine Learning Research, 6:1783–
1816, 2005.
[12] N.D. Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. In Advances in neural information processing systems 16: proceedings
of the 2003 conference, volume 16, page 329. The MIT Press, 2004.
[13] N.D. Lawrence and A.J. Moore. Hierarchical gaussian process latent variable models.
In Proceedings of the 24th international conference on Machine learning, pages 481–
488. ACM, 2007.
[14] AW van der Vaart and JH van Zanten. Reproducing kernel Hilbert spaces of Gaussian
priors. IMS Collections, 3:200–222, 2008.
[15] A.W. van der Vaart and J.H. van Zanten. Adaptive Bayesian estimation using a
Gaussian random field with inverse Gamma bandwidth. The Annals of Statistics, 37
(5B):2655–2675, 2009.

22

