Adaptive posterior convergence rates in non-linear latent variable models
Shuang Zhou Debdeep Pati Anirban Bhattacharya David Dunson





















Adaptive posterior convergence rates in non-linear latent variable models
    
Shuang Zhou[Corresponding author: shuang.zhou@stat.fsu.edu],    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

    December 30, 2023
==================================================================================================================================================================================================================================================================================





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 (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 model; 
Rate of convergence



§ 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; <cit.>) and the Gaussian process latent variable model (GP-LVM; <cit.>). 
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.  

<cit.> 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; <cit.>) mixture models by appropriate choice of the base measure G_0, posterior computation becomes complicated unless the base measure is conjugate to the kernel 𝒦. 

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 on closing this gap focusing in particular on studying how the posterior measure for the unknown density concentrates around the true density f_0 as the sample size n increases.  Assuming f_0 belongs to a Hö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, <cit.> 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 f_0, generalizing previous results, such as posterior consistency <cit.> or optimal rates for a particular smoothness level <cit.>. 
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 <ref>. To make the article self-contained, we also provide a brief background on Gaussian process priors. In Section <ref>, we formulate our assumptions on the true density f_0 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 <ref>. Section <ref> provides auxiliary results and Section <ref> proves the main theorem of posterior concentration rate. We discuss some implications of our results and outlines possible future directions Section <ref>.



§ NOTATIONS


Throughout the article, Y_1,…, Y_n are independent and identically
distributed with density f_0 ∈ℱ, the set of all densities on ℝ absolutely continuous with respect to the Lebesgue measure λ. The supremum and _1-norm are denoted by ·_∞ and ·_1, respectively. We let ·_p, ν denote the norm of
L_p(ν), the space of measurable functions with ν-integrable pth absolute power. For two density functions f, g ∈ℱ, let h denote the Hellinger distance defined as h^2(f, g) = √(f)-√(g)_2, λ=∫ (f^1/2 - g^1/2)^2 dλ, K(f,g) the Kullback-Leibler 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] →ℝ endowed with the supremum norm. For β >0, we let C^β[0, 1] denote the Hö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 write “≾” for inequality up to a constant multiple. Let ϕ(x) = (2π)^-1/2exp(-x^2/2) denote the standard normal density, and let ϕ_σ(x) = (1/σ) ϕ(x/σ). Let an asterisk denote a convolution e.g., (ϕ_σ * f)(y) = ∫ϕ_σ(y - x)f(x)dx, and let ϕ_σ^(i) * 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 <cit.>. A Borel measurable random element W with values in a separable Banach space (𝔹, ·) (e.g., C[0,1]) is called Gaussian if the random variable b^*W is normally distributed for any element b^* ∈𝔹^*, the dual space of 𝔹. The reproducing kernel
Hilbert space (RKHS) ℍ attached to a zero-mean Gaussian process W is
defined as the completion of the linear space of functions t ↦ EW(t)H
relative to the inner product

    ⟨ W(·)H_1;  W(· )H_2⟩_ℍ =  H_1H_2,

where H, H_1 and H_2 are finite linear combinations of the form ∑_ia_iW(s_i)
with a_i ∈ℝ and s_i in the index set of W.

Let W = (W_t: t ∈ℝ) be a Gaussian process with squared exponential covariance kernel. The spectral measure m_w of W is absolutely continuous with respect to the Lebesgue measure λ on ℝ with the Radon-Nikodym derivative given by

    dm_w/dλ(x) = 1/2π^1/2e^-x^2/4.


Define a scaled Gaussian process W^a=(W_at: t ∈ [0,1]), viewed as a map in C[0,1]. Let ℍ^a denote the RKHS of W^a, with the corresponding norm ·_ℍ^a. The unit ball in the RKHS is denoted ℍ^a_1. 

Throughout the paper, C, C_1, C_2, … denote global constants whose value may change one line to another.



§ ASSUMPTIONS ON THE TRUE DENSITY


It has been widely recognized that one needs certain smoothness assumptions and tail conditions on the true density f_0 to derive posterior convergence rates. We make the following assumptions in our case:

Assume logf_0∈ C^β[0, 1]. Let l_j(x)=d^j/dx^jlog f_0(x) be the derivatives for j=1, …, r with r= ⌊β⌋. For any β >0, assume that there exists a constant L > 0 such that

    |l_r(x) - l_r(y)| ≤ L|x-y|^β - r,

for all x ≠ y.



Assume f_0 is compactly supported on [a_0,b_0], for -∞ < a_0 <b_0 < ∞, and that there exists some interval [a,b] ⊂ [a_0, b_0] such that f_0 is nondecreasing on [a_0, a], bounded away from 0 on [a,b] and non-increasing on [b, b_0].

Assumption <ref> is useful in simplifying expressions for f_0 as convolutions with a given density, providing a key piece in our theoretical developments. Similar assumption on the local smoothness appeared in <cit.>, while in our case the global smoothness assumption is sufficient since f_0 is assumed to be compactly supported. 
 Assumption <ref> guarantees that for every δ > 0, there exists a constant C > 0 such that f_0 * ϕ_σ≥ Cf_0 for every σ < δ.  




§ THE NL-LVM MODEL


Consider the nonlinear latent variable model,

    y_i    = μ(η_i) + ϵ _i,   ϵ_i ∼(0, σ^2),   (i=1, …, n)
    μ   ∼Π_μ,   σ∼Π_σ,  η_i ∼(0,1),

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 the transfer function μ and scale σ is obtained on marginalizing out the latent variable as

    f(y; μ, σ) def= f_μ, σ(y)= ∫_0^1ϕ_σ(y-μ(x))dx.

Define a map g: C[0,1] × [0,∞) →ℱ with g(μ, σ) = f_μ, σ.
One can induce a prior Π on ℱ via the mapping g by placing independent priors
Π_μ and Π_σ on C[0,1] and [0, ∞) respectively, with Π = (Π_μ⊗Π_σ) ∘ g^-1. <cit.> 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 f_0 be a continuous density with cumulative distribution function F_0(t) = ∫_-∞^t f_0(x) dx. Assume f_0 to be non-zero almost everywhere within its support, so that F_0 : (f_0) → [0,1] is strictly monotone and hence has an inverse F_0^-1 : [0,1] →(f_0) satisfying F_0{F_0^-1(t)} = t for all t ∈(f_0).
Letting μ_0(x) = F_0^-1(x), one obtains

    f_μ_0, σ(y) = ∫_0^1ϕ_σ(y-F_0^-1(x))dx = ∫_-∞^∞ϕ_σ(y-t) f_0(t) dt,

where the second equality follows from the change of variable theorem. Thus, f_μ_0,σ(y) = ϕ_σ*f_0, i.e., f_μ_0, σ is the convolution of f_0 with a normal density having mean 0 and standard deviation σ. It is well known that the convolution ϕ_σ*f_0 can approximate f_0 arbitrary closely as the bandwidth σ→ 0 in the sense that for f_0 ∈ L_p(λ) with p ≥1, ϕ_σ*f_0 - f_0_p, λ→ 0 as σ→ 0. For Holder-smooth functions, the order of approximation can be characterized in terms of the smoothness. If f_0 ∈ C^β[0, 1] with β≤ 2, it follows from standard Taylor series expansion that || ϕ_σ*f_0 - f_0||_∞ = O(σ^β). For β > 2, a similar Taylor series expansion yields a sub-optimal error || ϕ_σ*f_0 - f_0||_∞ = O(σ^2). In this case, we can refine the approximation by convoluting with a sequence of functions f_j constructed by the procedure,

    f_j+1 = f_0 - _σ f_j,   _σf_j = ϕ_σ*f_j - f_j,    j ≥ 1.
 
For f_0 ∈ C^β[0,1] with β∈ (2j, 2j+2] we have ϕ_σ*f_j - f_0_∞ = O(σ^β) <cit.>. Although the f_js need not be non-negative in general, we show that they are non-negative on (f_0) when f_0 is compactly supported. It can be additionally shown that the normalizing constant is 1 + O(σ^β); let h_j denote the density obtained by normalizing f_j. We then approximate f_0 by ϕ_σ * h_β, where h_β = h_j for β∈ (2j, 2j+2].



Let λ̃ denote the Lebesgue measure on [0,1], or equivalently, the [0,1] distribution. For any measurable function μ : [0,1] →ℝ, let ν_μ denote the induced measure on (ℝ, ℬ), with ℬ denoting the Borel sigma-field on ℝ. 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,

    ∫_0^1ϕ_σ(y-μ(x))dx = ∫ϕ_σ(y-t) dν_μ(t),

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 <ref> assures us that this specification leads to large L_1 support on the space of densities.


Suppose the prior Π_μ on μ has full sup-norm support on C[0,1] so that (μ - μ^*_∞ < ϵ) > 0 for any ϵ > 0 and μ^* ∈ C[0,1], and the prior Π_σ on σ has full support on [0, ∞). If f_0 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 L_1 neighborhoods of any density f_0. We summarize the above discussion in the following theorem, with a proof provided in the appendix.





If Π_μ has full sup-norm support on C[0,1] and Π_σ has full support on [0, ∞), then the L_1 support of the induced prior Π on ℱ contains all densities f_0 which have a finite first moment and are non-zero almost everywhere on their support.


The conditions of Theorem <ref> are satisfied for a wide range of Gaussian process priors on μ (for example, a GP with a squared exponential or Matérn covariance kernel).


When f_0 has full support on ℝ, the quantile function μ_0 is unbounded near 0 and 1, so that μ_0_∞ = ∞. However, ∫_0^1μ_0(t) dt = ∫_ℝx f_0(x) dx, which implies that μ_0 can be identified as an element of L_1[0,1] if f_0 has finite first moment. Since C[0,1] is dense in L_1[0,1], the previous conclusion regarding L_1 support can be shown to hold in the non-compact case too.






§ THE MAIN THEOREM


We consider the case where f_0 satisfies Assumption <ref> and Assumption <ref>. For β-Hölder density f_0, we consider density h_β as defined after expression (<ref>).  Denote μ_0 the quantile function F_h_β^-1: [0,1] → [a_0, b_0], a continuous monotone function inheriting the smoothness of h_β. Note that h_β has the same smoothness of f_0 based on the construction of f_j, 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 Π_σ.

We assume μ follows a centered and rescaled Gaussian process denoted by (0, c^A), where A denotes the rescaled parameter, and assume A is a density satisfying for a>0,

    C_1a^pexp(-D_1a log^q a)≤ g(a) ≤ C_2a^pexp(-D_2a log^q a)

for positive constants C_1, C_2, D_1, D_2, nonnegative constant p and q, and every sufficiently large a>0.


We assume σ∼(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 <cit.> 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.
 
If f_0 satisfies Assumption <ref> and the priors Π_μ and Π_σ are as in
Assumptions <ref> and <ref> respectively, the best obtainable rate of posterior convergence relative to Hellinger metric h is

    ϵ_n = max(ϵ̃_n, ϵ̅_n),

where ϵ̃_n=n^-β/2β+1(log n)^t_1, ϵ̅_n= n^-β/2β+1(log n)^t_2, with nonnegative constants t_1=β(2∨ q)/(2β +1), t_2= t_1+1.


Unlike the treatment in discrete mixture models <cit.> where a compactly supported density is approximated with a discrete mixture of normals, the main trick here is to approximate the true density f_0 by the convolution ϕ_σ*f_0 and allow the prior on the transfer function to appropriately concentrate around the true quantile function μ_0 ∈ C[0,1].



§ 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 <ref>.

For μ_1, μ_2 ∈ C[0, 1] and σ_1, σ_2 > 0,

    h^2(f_μ_1,σ_1, f_μ_2, σ_2) ≤ 1- √(2σ_1σ_2/σ_1^2 + σ_2^2)exp{-μ_1 - μ_2_∞^2/4(σ_1^2 + σ_2^2)}.



Note that by Hölder's inequality,

    f_μ_1, σ_1(y)f_μ_2, σ_2(y) ≥{∫_0^1√(ϕ_σ_1(y - μ_1(x)))√(ϕ_σ_2(y - μ_2(x)))dx}^2.

Hence,

    h^2(f_μ_1,σ_1, f_μ_2, σ_2)    ≤   ∫[∫_0^1ϕ_σ_1(y-μ_1(x))dx  + ∫_0^1ϕ_σ_2(y-μ_2(x))dx 
       -   2∫_0^1√(ϕ_σ_1(y - μ_1(x)))√(ϕ_σ_2(y - μ_2(x)))dx]dy.

By changing the order of integration (applying Fubini's theorem since the function within the integral is jointly integrable) we get,

    h^2(f_μ_1,σ_1, f_μ_2, σ_2)    ≤   ∫_0^1h^2(f_μ_1(x),σ_1, f_μ_2(x), σ_2)dx 
       =   ∫_0^1[1- √(2σ_1σ_2/σ_1^2 + σ_2^2)exp{-(μ_1(x) - μ_2(x))^2/4(σ_1^2 + σ_2^2)}]dx
       ≤    1- √(2σ_1σ_2/σ_1^2 + σ_2^2)exp{-μ_1 - μ_2_∞^2/4(σ_1^2 + σ_2^2)}.





When σ_1 = σ_2  = σ, h^2(f_μ_1,σ, f_μ_2, σ) ≤ 1 - exp{μ_1 -μ_2_∞^2/ 8 σ^2}, which implies that h^2(f_μ_1,σ, f_μ_2, σ) ≾μ_1 -μ_2_∞^2/σ^2.



The standard inequality h^2(f_μ_1,σ_1, f_μ_2, σ_2) ≤f_μ_1,σ_1- f_μ_2, σ_2_1 relating the Hellinger distance to the total variation distance leads to the cruder bound

    h^2(f_μ_1,σ_1, f_μ_2, σ_2) ≤ C_1 μ_1 -μ_2_∞/(σ_1 ∧σ_2) + C_2|σ_2 - σ_1|/(σ_1 ∧σ_2),

which is linear in μ_1 -μ_2_∞. This bound is less sharp than what is obtained in Lemma <ref> and does not suffice for obtaining the optimal rate of convergence.


To control the Kullback-Leibler divergence between the true density f_0 and the model f_μ, σ, we derive an upper bound for logf_0/f_μ, σ_∞ in Lemma <ref>.

If f_0 satisfies Assumption <ref>,

    logf_0/f_μ, σ_∞≤ C + μ - μ_0_∞^2/σ^2

for some constant C > 0.


Note that

    f_μ, σ(y)    =   1/√(2π)σ∫_0^1exp{-(y-μ(x))^2/2σ^2}dx
       ≥   1/√(2π)σ∫_0^1exp{-(y-μ(x))^2/σ^2}dx exp{-μ-μ_0_∞^2/σ^2}
       ≥    C ϕ_σ/√(2) * f_0 (y) exp{-μ-μ_0_∞^2/σ^2}
       ≥    C f_0(y) exp{-μ-μ_0_∞^2/σ^2},

where the last inequality follows from Lemma 6 of <cit.> since f_0 is compactly supported by Assumption <ref>. This provides the desired inequality. 










For β∈ (2j, 2j+2], j ≥ 0 and f_j constructed by <ref>, we have the expression f_j = ∑_i=0^j(-1)^i j+1  i+1ϕ_σ^(i) f_0.


The proof can be found in Appendix <ref>. The expression of f_j as a linear combination of true density and the folded convolutions indicates that f_j is as smooth as f_0. 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 f_j to ensure the approximation error remains of order O(σ^β).    



For any β > 0, let f_0 satisfy Assumption <ref> and <ref>, integer j be such that for β∈ (2j, 2j+2], f_j constructed by (<ref>). For any constant L and all x ∈ [a_0, b_0], we have 

    ϕ_σ * f_β (x) = f_0(x) (1 + O(D(x)σ^β)),

where

    D(x) = ∑_i=1^r c_i |l_j(x)|^β/i + c_r+1,

for nonnegative constants c_i, i = 1, …, r, and c_r+1 a multiple of L.


Following the proof of Lemma 1 in <cit.>, for any x, y ∈ [a_0, b_0],

    logf_0(y)≤logf_0(x) + ∑_i=1^rl_j(x)/j!(y-x)^j + L|y-x|^β,


    logf_0(y)≥logf_0(x) + ∑_i=1^rl_j(x)/j!(y-x)^j - L|y-x|^β.

Define

    B^u_f_0,r(x,y) = ∑_i=1^rl_j(x)/j!(y-x)^j + L|y-x|^β,


    B^l_f_0,r(x,y) = ∑_i=1^rl_j(x)/j!(y-x)^j - L|y-x|^β.

Then we have

    e^B^u_f_0,r≤ 1 + B^u_f_0,r + 1/2!(B^u_f_0,r)^2 + … + M |B^u_f_0,r|^r+1,


    e^B^l_f_0,r≥ 1 + B^l_f_0,r + 1/2!(B^l_f_0,r)^2 + … - M |B^l_f_0,r|^r+1.

where 

    M = 1/(r+1)!exp{sup_x, y  ∈ [a_0, b_0], x ≠ y (|∑_j=1^rl_j(x)/j! (y-x)^j| + L|y-x|^β)}.

Note that f_0 is bounded on [a_0, b_0], we consider the convolution on the whole real line by extending f_0 analytically outside [a_0, b_0]. For β∈ (1, 2], r = 1 and x ∈ (a_0, b_0), 

    ϕ_σ* f_0(x)    ≤    f_0(x) ∫ e^B^u_f_0,r(x, y)ϕ_σ(y-x) dy
       ≤    f_0(x) ∫_ℝϕ_σ(y-x) [ 1 + L|y-x|^β + M ( l^2_1(x)(y-x)^2 + L^2|y-x|^2β) ] dy.

Since l_j(x)'s are all continuous on [a_0, b_0], there exist finite constants M_j such that |l_j| ≤ M_j and |y-x| ≤ |b_0 - a_0|. The integral in the last inequality can be bounded by 

    ∫_ℝϕ_σ(y-x) [ 1 + L|y-x|^β + M ((M_1|b_0 - a_0|)^2-β |l_1(x)(y-x)|^β + (L^2|b_0 - a_0|^β )|y-x|^β) ] dy

Therefore,

    ϕ_σ* f_0(x) ≤ f_0(x) ( 1 + (r_1|l_1(x)|^β +r_2) σ^β),

where r_1 = M (M_1|b_0 - a_0|)^2-βμ_β,   r_2 =  (L + ML^2)μ_β. 

In the other direction, 

    ϕ_σ* f_0(x) ≥  f_0(x) ∫ϕ_σ(y-x) [ 1 - L|y-x|^β - M ( l^2_1(x)(y-x)^2 + L^2|y-x|^2β) ] dy.

Thus we achieve (<ref>).

For any β > 2, suppose β∈ (2j, 2j +2], j > 1. First we calculate ϕ_σ* f_0, ϕ^(2)_σ* f_0, …, ϕ^(j)_σ* f_0(x), by Lemma <ref> to get ϕ_σ* f_β(x). The calculation of ϕ^(i)_σ* f_0(x) is the same as ϕ_σ* f_0(x) except taking the convolution with ϕ_√(i)σ. The terms σ^2, σ^4, …, σ^2j caused by the factors containing |y-x|^k, k < β in ϕ^(i)_σ f_0 can be canceled out by Lemma <ref>.
For terms containing |y-x|^k, k ≥β, we take out |y-x|^β and bound the rest by a certain power of |l_j(x)| or some constant. 




Let f_0 satisfy Assumption <ref> and <ref>. With A_σ = {x: f_0(x) ≥σ^H} , we have

    ∫_A^c_σf_0(x)dx = O(σ^2β),   ∫_A^c_σϕ_σ*f_j(x) dx = O(σ^2β),

for all non-negative integer j, sufficiently small σ and sufficiently Large H.


Under Assumption <ref> there exists (a, b) ⊂ [a_0, b_0] such that A^c_σ⊂ [a_0, a) ∪ (b, b_0] if we choose σ sufficiently small, so that f_0 (x) ≤σ^H for x ∈A^c_σ. Therefore, ∫_A^c_σf_0(x) ≤σ^H |b_0 - a_0| ≤ O(σ^2β) if we choose H ≥ 2β. Using Lemma <ref>, 

    ∫_A^c_σϕ_σ* f_j(x) dx =  ∫_A^c_σ f_0(x) (1 + O(D(x)σ^β)) ≤ O(σ^H),

with bounded D(x) and H ≥ 2β it is easy to bound the second integral by O(σ^2β).



Suppose f_0(x) satisfies Assumption <ref> and <ref>. For β >2 and j such that β∈ (2j,2j+2], we can construct the density h_β from (<ref>) and show that h_β satisfies Lemma <ref> and Lemma <ref>.


To get the density function we first show that f_j is nonnegative and compute the normalizing constant ∫ f_j(x) = 1 + O(σ^β). Following the proof of Lemma 2 in <cit.>, we treat log f_0 as a function in C^2 [0,1] and obtain the same form of ϕ_σ*f_0 as (<ref>). For small enough σ we can find  ρ_1 ∈ (0,1) very close to 0 such that

    ϕ_σ* f_0(x) = f_0(x) (1 + O(D^(2)(x)σ^2)) < f_0(x)(1 + ρ_1),

 where D^(2) contains |l_1(x)| and |l_2(x)| with certain power, so D^(2) is bounded. Then we have
 
    f_1(x) = 2f_0(x) - K_σ f_0(x) > 2f_0(x) - f_0(x) (1 + ρ_1) = f_0(x) (1 - ρ_1).

 Then we treat logf_0 as a function with β = 4, j = 1. Similarly, we can get
 
    ϕ_σ* f_1(x) = f_0(x) (1 + O(D^(4)(x)σ^4)),

 where D^(4) contains |l_1(x)|, …. |l_4(x)|. We can find 0 < ρ_2 < ρ_1 such that ϕ_σ* f_1(x) < f_0(x) (1 + ρ_2), then can get
 
    f_2 (x) = f_0(x) - (ϕ_σ* f_1(x) - f_1(x)) > f_0(x)(1 - ρ_1 - ρ_2) > f_0(x)(1 - 2ρ_1).

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

Now we calculate the normalizing constant for f_j. When β < 2, with Lemma <ref>,

    ∫ f_1(x) = ∫ f_0 - (ϕ_σ* f_0 - f_0) ≤∫ f_0 + |∫ (ϕ_σ*f_0 - f_0)| ≤ 1 + O(σ^β).

For β∈ (2,4], 

    ∫ f_2(x) = ∫ f_0 -  (ϕ_σ* f_1 - f_1) ≤∫ f_1 + |∫ (ϕ_σ*f_1 - f_0)| ≤ 1 + O(σ^β).

 Then by induction, we have ∫ f_j = 1 + O(σ^β), so that we have the density

    h_β= f_j/1 + O(σ^β),   β∈ (2j, 2j+2].

 Now to show h_β satisfying (<ref>), note that 
 
    ϕ_σ* h_β = ϕ_σ* f_j/∫ f_j = f_0(x)(1+ O(D(x)σ^β))/1+ O(σ^β).

 Since D(x) is bounded and for sufficiently small σ we consider term 
 
    (1+ O(D(x)σ^β))/1+ O(σ^β)≾(1+ O((D(x)+1)σ^β)+ O(D(x)σ^2β))/1+ O(σ^β) = 1+ O(D(x)σ^β),

which directly leads to the same form as (<ref>), and obviously Lemma <ref> is satisfied.



 Let f_0 satisfy Assumption <ref> and <ref>, and integer j be such that β∈ (2j, 2j+2]. Then we can show that the density h_β defined by (<ref>) satisfies, 
    ∫ f_0(x) logf_0(x)/ϕ_σ* h_β(x) = O(σ^2β),

for sufficiently small σ and all x ∈ [a_0, b_0].


Again consider the set A_σ = {x: f_0(x) ≥σ^H} with arbitrarily large H. We separate the Kullback-Leibler divergence into

    ∫_[a_0, b_0] f_0 logf_0/ϕ_σ* h_β≤∫_A_σ(f_0 - ϕ_σ*h_β)^2/ϕ_σ* h_β +  ∫_A^c_σ f_0 logf_0/ϕ_σ* h_β +  ∫_A^c_σ(ϕ_σ*h_β - f_0).

Under Assumption <ref> and by Remark 3 in <cit.>, for small enough σ there exists constant C such that for all x ∈ [a_0, b_0], ϕ_σ*f_0 ≥ Cf_0, especially on set A^c_σ f_0 satisfies ϕ_σ*f_0 ≥ f_0/3. Also in the proof of Lemma <ref> we can find ρ∈ (0,1) such that f_β > ρ f_0. Then we have on set A_σ with sufficiently small σ

    ϕ_σ*h_β = ϕ_σ*f_β/1+O(σ^β)≥ρϕ_σ*f_0/1 + O(σ^β)≥ K f_0,

for some positive and finite constant K. Applying Lemma <ref>, the first integral on the right side of <ref> can be bounded by 

    ∫_A_σ(f_0 - ϕ_σ*h_β)^2/ϕ_σ* h_β   ≤∫_A_σ[f_0(x) - f_0(x)(1 + O(D(x)σ^β))]^2/K_1 f_0(x)
       ≾∫_A_σ f_0(x) O(D^2(x)σ^2β) = O(σ^2β).

To bound the second integral of r.h.s again by Remark 3 in <cit.> we get ϕ_σ*h_β≥ρ/3(1+ O(σ^β))f_0, so easily we can find a constant C < 1 such that ϕ_σ*h_β≥ C f_0. With Lemma <ref> clearly the second and third term can be bounded by O(σ^2β).



Let ℍ_1^a denote the unit ball of RKHS of the Gaussian process with rescaled parameter a and 𝔹_1 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]ℍ_1^a, ·_∞) ≤ Kr( log1/ϵ)^2.




 Since we can write any element of ℍ_1^a as a function of Re(z) by  <cit.> which can be analytically extended to some interval containing Ω^a = {z∈ℂ: |Re(z)| ≤ R} with R = δ/(6max(a,1)), so for any h ∈Ω^a, |2aRe(z)| ≤ |δ/6max(a,1)· 2a| = δ/3. Consider any b with |b-a|≤1, we can show that any element of ℍ_1^b can be extended analytically to Ω^a noting that for z ∈Ω^a related to the maximum norm, 

    2bRe(z)≤2aRe(z) + 2(b-a)Re(z)≤δ/3+2Re(z)≤2/3δ.
  
Therefore, ℱ^a forms one ϵ-net over ℍ_1^b. We find one set Γ = {a_i, i = 1, …, k} with k = ⌊ r⌋ +1 and a_k = r, such that for any b∈ [0,r] there exists some a_i satisfying |b-a_i|≤ 1, so that ∪_i≤ kℱ^a_i forms an ϵ-net over ∪_a ≤ rℍ_1^a. Since the covering number of ∪_i≤ kℱ^a_i is bounded by summation of covering number of ℱ^a_i, we obtain

    log N(ϵ, ∪_a ∈ [0,r]ℍ_1^a, ·_∞) ≤log(∑_i=1^k #(ℱ^a_i)) ≤log(k ·#(ℱ^r)) ≤  Kr( log1/ϵ)^2.
 
To complete the proof, note that the piecewise polynomials are constructed on the partition of [0,1], ∪_i≤ mB_i, where B_i's are disjoint interval with length shorter than R= δ/(6max(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 B_i, both the approximation error and tail estimate are invariant to interval length R, therefore we have #(ℱ^a) ≤#(ℱ^b) if a ≤ b, for a,b ∈ [0,r].  


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/R^n, where |h|^2 ≤ C^2. 






§ PROOF OF THE MAIN THEOREM

Proof of Theorem <ref>:
Following <cit.>, we need to find sequences ϵ̅_n,ϵ̃_n → 0 with
nmin{ϵ̅_n^2,ϵ̃_n^2}→∞  such that there exist constants C_1, C_2, C_3, C_4> 0 and sets ℱ_n ⊂ℱ so that,

    log N(ϵ_n, ℱ_n, d) ≤ C_1nϵ̅_n^2 
       Π(ℱ_n^c) ≤ C_3exp{-nϵ̃_n^2(C_2+4)}
       Π( f_μ, σ:  ∫ f_0 logf_0/f_μ, σ≤ϵ̃_n^2, ∫ f_0 log(f_0/f_μ, σ)^2 ≤ϵ̃_n^2 ) ≥ C_4exp{-C_2nϵ̃_n^2}.

Then we can conclude that for ϵ_n = max{ϵ̅_n, ϵ̃_n} and sufficiently large M > 0, the posterior probability

    Π_n(f_μ, σ: d(f_μ, σ, f_0) > Mϵ_n | Y_1, …, Y_n) → 0   a.s.  P_f_0.


We consider the Gaussian process μ∼ W^A given A, with A satisfying Assumption <ref>.

We will first verify (<ref>) along the lines of <cit.>.
Note that

    h^2(f_0, f_μ, σ) ≾ h^2(f_0, f_μ_0, σ) + h^2(f_μ_0, σ, f_μ, σ).

Since f_μ_0, σ = ϕ_σ*h_β, by Lemma <ref>, one obtains under Assumptions <ref> and <ref>, 
    h^2(f_0, f_μ_0, σ) ≤∫ f_0 log(f_0/f_μ_0, σ) ≾ O(σ^2β).


From Lemma <ref> and the following remark, we obtain

    h^2(f_μ_0, σ, f_μ, σ) ≾μ- μ_0_∞^2/σ^2.

From Lemma 8 of <cit.>, one has

    ∫ f_0 log(f_0/f_μ, σ)^i≤ h^2(f_0, f_μ,σ)(1 + logf_0/f_μ, σ_∞)^i

for i=1,2.

From (<ref>)-(<ref>), for any b ≥ 1 and ϵ̃_n^2 = σ_n^2β,

    {σ∈ [σ_n, σ_n + σ_n^b], μ - μ_0_∞≾σ_n^β+1}⊂
    {∫ f_0 logf_0/f_μ, σ≾σ_n^2β, ∫ f_0 log(f_0/f_μ, σ)^2 ≾σ_n^2β}.


Since μ_0 ∈ C^β+1[0,1], from Section 5.1 of <cit.>,

    (μ - μ_0_∞≤ 2δ_n) ≥ C_4exp{-C_5(1/δ_n)^1/β+1log(1/δ_n)^2∨ q}(C_6/δ_n)^p+1/β+1,

for δ_n → 0 and constants C_4, C_5, C_6 > 0.  Letting δ_n = σ_n^3, we obtain

    (μ - μ_0_∞≤ 2δ_n) ≥exp{-C_7(1/σ_n)log^2∨ q(1/σ_n^β+1)},

for some constant C_7 > 0.
Since σ∼ IG(a_σ, b_σ),
we have

    (σ∈ [σ_n, 2σ_n ])    =   b_σ^a_σ/Γ(a_σ)∫_σ_n^2σ_nx^-(a_σ+1) e^-b_σ/xdx 
       ≥   b_σ^a_σ/Γ(a_σ)∫_σ_n^2σ_n e^-2b_σ/xdx 
       ≥   b_σ^a_σ/Γ(a_σ)σ_nexp{-b_σ/σ_n }
       ≥   exp{-C_8/σ_n },

for some constant C_8> 0. Hence

    {σ∈ [σ_n, 2σ_n ], μ - μ_0_∞≾σ_n^β+1}   ≥   exp{-C_7(1/σ_n)log^2∨ q(1/σ_n^β+1)}exp{-C_8/σ_n }
       ≥   exp{-2C_7(1/σ_n)log^2∨ q(1/σ_n^β+1)}.

Then (<ref>) will be satisfied with ϵ̃_n = n^-β/2β+1log^t_1(n), where t_1=β(2∨ q)/2β+1 and some C_9 > 0.

Next we construct a sequence of subsets ℱ_n such that <ref> and <ref> are satisfied with
ϵ̅_n = n^-β/2β+1log ^t_2n and ϵ̃_n for some global constant t_2 > 0.

Letting 𝔹_1 denote the unit ball of
C[0,1] and given positive sequences M_n, r_n, define

    B_n =  ∪_a < r_n(M_n ℍ^a_1) + δ̅_n𝔹_1

as in <cit.>, with δ̅_n =ϵ̅_nl_n/K_1, K_1 = 2(2/π)^1/2 and let

    ℱ_n = {f_μ, σ: μ∈ B_n, l_n < σ < h_n }.

First we need to calculate N(ϵ̅_n, ℱ_n, ·_1). Observe that for σ_2 >σ_1 > σ_2/2,

    f_μ_1, σ_1 - f_μ_2, σ_2_1 ≤(2/π)^1/2μ_1 - μ_2_∞/σ_1 + 3(σ_2 - σ_1)/σ_1.


Taking κ_n =min{ϵ̅_n/6, 1} and σ_m^n = l_n (1+ κ_n)^m, m ≥ 0, we obtain a partition of [l_n, h_n] as l_n=σ_0^n < σ_1^n < ⋯ < σ_m_n-1^n < h_n ≤σ_m_n^n with

    m_n= (logh_n/l_n) 1/log( 1+κ_n) +1.

One can show that 3(σ_m^n -σ_m-1^n)/σ_m-1^n = 3κ_n ≤ϵ̅_n/2. Let
{μ̃_k^n, k=1, …, N(δ̅_n, B_n, ·_∞)}  be
a δ̅_n-net of B_n. Now consider the set

    {(μ̃_k^n, σ_m^n): k=1, …, N(δ̅_n, B_n, ·_∞), 0≤ m≤ m_n }.

Then for any f= f_μ, σ∈ℱ_n, we can find (μ̃_k^n, σ_m^n) such that
μ - μ̃_k^n_∞ < δ̅_n. In addition, if one has σ∈ (σ_m-1^n, σ_m^n],
then

    f_μ, σ - f_μ_k^n, σ^n_m_1 ≤ϵ̅_n.

Hence the set in (<ref>) is an ϵ̅_n-net of ℱ_n and its covering number is given by

    m_nN(δ̅_n, B_n, ·_∞).

From the proof of Theorem 3.1 in <cit.>, for any M_n, r_n with r_n > a_0, we obtain

    log N(2δ̅_n, B_n, ·_∞) ≤ K_2 r_n ( log(M_n/δ̅_n))^2.

Again from the proof of Theorem 3.1 in <cit.>, for r_n > 1 and for M_n^2 > 16K_3r_n (log (r_n / δ̅_n))^2, we have

    (W^A ∉ B_n) ≤K_4r_n^p e^-K_5r_nlog^q r_n/K_5log^q r_n + exp{-M_n^2/8}

for constants K_3, K_4, K_5 > 0.

Next we calculate P(σ∉ [l_n, h_n]). Observe that

    (σ∉ [l_n, h_n ])    =    P(σ^-1 < h_n^-1) + P(σ^-1 > l_n^-1)
       ≤   ∑_k=α_σ^∞e^-b_σh_n^-1(b_σh_n^-1)^k/k! + b_σ^a_σ/Γ(a_σ)∫_l_n^-1^∞ e^-b_σx/2dx 
       ≤    e^-a_σlog(h_n) + b_σ^a_σ/Γ(a_σ)e^-b_σl_n^-1/2.


Thus with h_n = O(exp{n^1/2β+1(log n)^2t_1}), l_n = O(n^-1/2β+1(log n)^-2t_1), r_n =O(n^1/2β+1(log n)^2t_1), M_n = O(n^1/2β+1(log n)^t_1+1),
(<ref>) and (<ref>) implies

    Π(ℱ_n^c)= exp{-K_6nϵ̃^2_n}

for some constant K_6 > 0 guaranteeing that (<ref>) is satisfied with ϵ̃_n = n^-β/2β+1(log n)^t_1.

Also with ϵ̅_n = n^-β/2β+1(log n)^t_1+1,  it follows from
(<ref>) and (<ref>) that

    log N(ϵ̅_n, ℱ_n, ·_1) ≤ K_7 n^1/2β+1(log n)^2t_1+2

for some constant K_7 > 0.

Hence max{ϵ̅_n, ϵ̃_n} = n^-β/2β+1(log n)^t_1+1.



§ 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 <cit.> 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 <cit.>. However, for densities supported on ℝ, the results fail as the  corresponding 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. 







§ APPENDIX



 §.§ Proof of Theorem <ref>


Let f_0 be a density with quantile function μ_0 that satisfies the conditions of Theorem <ref>. Observe that μ_0_1 = ∫_t=0^1 μ_0(t) dt = ∫_-∞^∞z f_0(z) dz < ∞ since f_0 has a finite first moment, and thus μ_0 ∈_1[0, 1]. Fix ϵ > 0. We want to show that Π{ B_ϵ(f_0) } > 0, where B_ϵ(f_0) = {f  : f - f_0_1 < ϵ}.

Note that μ_0 ∉ C[0, 1], so that ( μ - μ_0_∞ < ϵ) can be zero for small enough ϵ. The main idea is to find a continuous function μ̃_̃0̃ close to μ_0 in L_1 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 ϕ_σ*f_0 - f_0_1 → 0 as σ→ 0, find σ_1 such that ϕ_σ*f_0 - f_0_1 < ϵ/2 for σ < σ_1. Pick any σ_0 < σ_1.
Since C[0, 1] is dense in _1[0, 1], for any δ > 0, we can find a continuous function μ̃_̃0̃ such that μ_0 - μ̃_̃0̃_1 < δ. Now, f_μ, σ - f_μ̃_̃0̃, σ_1 ≤ C μ - μ̃_̃0̃_1/σ for a global constant C. Thus, for δ = ϵσ_0/4,

    {f_μ, σ : σ_0 < σ < σ_1, μ - μ̃_̃0̃_∞ < δ}⊂{f_μ, σ : f_0 - f_μ, σ_1 < ϵ},

since f_0 - f_μ, σ_1 < f_0 - f_μ_0, σ_1 + f_μ_0, σ - f_μ̃_̃0̃, σ_1 + f_μ̃_̃0̃, σ - f_μ, σ_1 and f_μ_0, σ = ϕ_σ*f_0. Thus, Π{ B_ϵ(f_0) } > (μ - μ̃_̃0̃_∞ < δ)  (σ_0 < σ < σ_1) > 0, since Π_μ has full sup-norm support and Π_σ has full support on [0, ∞).





 §.§ Proof of Lemma <ref>


Consider f_j constructed by (<ref>). When i=1, f_1 = 2f_0 - ϕ_σ*f_0, so the form holds. 
By induction, suppose this form holds for j > 1, then 

    f_j+1   =    f_0 - (ϕ_σ*f_j - f_j) 
       =    f_0 +  ∑_i=0^j(-1)^i+1j+1  i+1ϕ_σ^(i+1)* f_0 + ∑_i=0^j(-1)^i j+1  i+1ϕ_σ^(i)* f_0 
       =    (j+2)f_0 +  ∑_i=1^j+1(-1)^ij+1  i+1ϕ_σ^(i)* f_0 +  ∑_i=1^j(-1)^ij+1  iϕ_σ^(i)* f_0
       =    (j+2)f_0  + ∑_i=1^j (-1)^i( j+1  i+1 + j+1  i) ϕ_σ^(i)*f_0 + (-1)^j+1ϕ_σ^(i +1)*f_0 
       =    (j+2)f_0 + ∑_i=1^j(-1)^ij+2  i+1ϕ_σ^(i)* f_0 + (-1)^j+1ϕ_σ^(i +1)*f_0
       =   ∑_i=0^j+1(-1)^ij+2  i+1ϕ_σ^(i)* f_0.

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




plain

