Supremum Norm Posterior Contraction and Credible Sets for Nonparametric Multivariate Regression

DOI:10.1214/15-AOS1398

©Institute of Mathematical Statistics, 2016

SUPREMUM NORM POSTERIOR CONTRACTION AND CREDIBLE SETS FOR NONPARAMETRIC MULTIVARIATE REGRESSION

B

W

W

Y

S

G

Université Paris Dauphine and North Carolina State University

In the setting of nonparametric multivariate regression with unknown er- ror variance σ², we study asymptotic properties of a Bayesian method for estimating a regression function f and its mixed partial derivatives. We use a random series of tensor product of B-splines with normal basis coefficients as a prior for f , and σ is either estimated using the empirical Bayes approach or is endowed with a suitable prior in a hierarchical Bayes approach. We es- tablish pointwise, L₂and L_∞-posterior contraction rates for f and its mixed partial derivatives, and show that they coincide with the minimax rates. Our results cover even the anisotropic situation, where the true regression func- tion may have different smoothness in different directions. Using the conver- gence bounds, we show that pointwise, L₂and L_∞-credible sets for f and its mixed partial derivatives have guaranteed frequentist coverage with opti- mal size. New results on tensor products of B-splines are also obtained in the course.

1. Introduction. Consider the nonparametric regression model Y

= f (X

) + ε

, i = 1, . . . , n,

(1.1)

where Y

is a response variable, X

is a d-dimensional covariate, and ε

, . . . , ε

are independent and identically distributed (i.i.d.) as N(0, σ

) with unknown 0 < σ < ∞. The covariates are deterministic or are sampled from some fixed dis- tribution independent of ε

. In both cases, each X

takes values in some rectangular region in R

, which is assumed to be [0, 1]

without loss of generality. We follow the Bayesian approach by representing f by a finite linear combination of tensor products of B-splines and endowing the coefficients with a multivariate normal prior. We consider both the empirical and the hierarchical Bayes approach for the variance σ

. For the latter approach, a conjugate inverse-gamma prior is particu- larly convenient.

We study frequentist behavior of the posterior distributions and the resulting credible sets for f and its mixed partial derivatives, in terms of pointwise, L

and L

(supremum) distances. We assume that the true regression function f

belongs to an anisotropic Hölder space (see Definition 2.1 below), and the errors under the true distribution are sub-Gaussian.

Received November 2014; revised September 2015.

MSC2010 subject classifications.Primary 62G08; secondary 62G05, 62G15, 62G20.

Key words and phrases. Tensor product B-splines, sup-norm posterior contraction, nonparametric multivariate regression, mixed partial derivatives, anisotropic smoothness.

1069

Posterior contraction rates for regression functions in the L

-norm are well stud- ied, but results for the stronger L

-norm are limited. Giné and Nickl [14] studied contraction rates in L

-metric, 1 ≤ r ≤ ∞, and obtained optimal rate using conju- gacy for the Gaussian white noise model and a rate for density estimation based on a random wavelet series and Dirichlet process mixture using a testing approach.

In the same context, Castillo [2] introduced techniques based on semiparametric Bernstein–von Misses (BvM) theorems to obtain optimal L

-contraction rates.

Hoffman et al. [17] derived adaptive optimal L

-contraction rate for the white noise model and also for density estimation. Scricciolo [25] applied the techniques of [14] to obtain L

-rates using Gaussian kernel mixtures prior for analytic true densities.

De Jonge and van Zanten [9] used finite random series based on tensor products of B-splines to construct a prior for nonparametric regression and derived adaptive L

-contraction rate for the regression function in the isotropic case. A BvM theo- rem for the posterior of σ is treated in [10]. Shen and Ghosal [28, 29] used tensor products of B-splines, respectively, for Bayesian multivariate density estimation and high dimensional density regression in the anisotropic case.

Nonparametric confidence bands for an unknown function were considered by [1, 30] and more recently by [5, 6, 13]. A Bayesian approaches the problem by constructing a credible set with a prescribed posterior probability. It is then natural to ask if the credible set has adequate frequentist coverage for large sample sizes.

For parametric problems, a BvM theorem concludes that Bayesian and frequentist measures of uncertainly are nearly the same in large samples. However, for the infinite dimensional normal mean model (equivalently the Gaussian white noise model), [7, 12] observed that for many true parameters in 

, credible regions can have inadequate coverage. Leahu [20] showed that if the prior variances are chosen very big so that the support of the prior extends beyond 

, then coverage can be obtained. Knapik et al. [18, 19] showed that for sequences with specific smoothness, by deliberately undersmoothing the prior, coverage of credible sets may be guaranteed. Sniekers and van der Vaart [31] obtained similar results for nonparametric regression using a scaled Brownian motion prior.

Castillo and Nickl [3] showed that for the Gaussian white noise model a BvM theorem can hold in weaker topologies for some natural priors, and the resulting credible sets appropriately modified will have asymptotically the correct coverage and optimal size. A similar result for the stronger L

-norm using this weak notion of BvM theorem is considered in [4]. Adaptive L

-credible regions with adequate frequentist coverage are constructed using the empirical Bayes approach in [34]

for the Gaussian white noise model and in [27] for the nonparametric regression model using smoothing splines. In the setting of the Gaussian white noise model, Ray [23] constructed adaptive L

-credible sets using a weak BvM theorem, and also adaptive L

-credible band using a spike and slab prior.

In this paper, we consider multivariate nonparametric regression with unknown

variance parameter and study posterior contraction rates and coverage of credible

sets in the pointwise, L

- and L

-senses, for the regression function f as well as its mixed partial derivatives. Study of posterior contraction rate in L

-norm is important for its natural interpretation and implications for other problems such as the convergence of the mode of a function. An L

-credible band is easier to visualize than a L

-credible set. We assume that the smoothness of the function is given but allow anisotropy, so the smoothness level may vary with the direc- tion. Anisotropic function has applications in estimating time-dependent spectral density of a locally stationary time series (see [22]), and variable selection (see [16]).

A prior on the regression function is constructed using a finite random series of tensor products of B-splines with normally distributed coefficients. Posterior conjugacy leads to explicit expression for the posterior distribution which is con- venient for computation as well as theoretical analysis. Although wavelets are also widely used to construct random series priors, B-splines have the added advan- tage in that mixed partial derivatives of f are expressible in terms of lower degree B-splines. This allows posterior analysis for mixed partial derivatives of f , a topic that is largely unaddressed in the literature, except implicitly as inverse problems in the Gaussian white noise model.

The paper is organized as follows. The next section introduces notation and assumptions. Section 3 describes the prior and the resulting posterior distribution.

Section 4 contains main results on pointwise and L

-contraction rates of f and its mixed partial derivatives. Section 5 presents results on coverage of the correspond- ing credible sets. Section 6 contains a simulation study of the proposed method.

Proofs are in Section 7. New results on tensor products of B-splines are presented in the Appendix.

2. Assumptions and preliminaries. We describe notation and assumptions used in this paper. Given two numerical sequences a

and b

, a

= O(b

) or a

 b

means a

/b

is bounded, while a

= o(b

) or a

 b

means a

/b

→ 0. Also, a

 b

means a

= O(b

) and b

= O(a

). For stochastic sequence Z

, Z

= O

(a

) means P( |Z

| ≤ Ca

) → 1 for some constant C > 0. Let N = {1, 2, . . .}

and N

= N ∪ {0}.

Define x

= (

|x

|

)

, 1 ≤ p < ∞, x

= max

|x

|, and write

x for x

, the Euclidean norm. We write x ≤ y if x

≤ y

, k = 1, . . . , d. For an m × m matrix A = ((a

)), let λ

(A) and λ

(A) be the smallest and largest eigenvalues, and the (r, s) matrix norm of A as A

= sup{Ax

: x

≤ 1}.

In particular, A

= |λ

(A) | and A

= max

j=1

|a

|. These norms are related by |a

| ≤ A

≤ A

for 1 ≤ i, j ≤ m. With another matrix B of the same size, A ≤ B means B − A is nonnegative definite. We denote by I

the m × m identity matrix and by 1

the d × 1 vector of ones.

For f : U → R on some bounded set U ⊆ R

with interior points, let f 

be the L

-norm, and f 

= sup

|f (x)|. For r = (r

, . . . , r

)

∈ N

, let D

be the partial derivative operator ∂

/∂x

· · · ∂x

, where |r| =

r

. If r = 0, we interpret D

f ≡ f . We say Z ∼ N

(ξ , ) if Z has a J -dimensional normal distribution with mean vector ξ and covariance matrix . For a random function {Z(t), t ∈ U}, write Z ∼ GP(ξ, ) if Z is a Gaussian process with EZ(t) = ξ(t) and Cov(Z(s), Z(t)) = (s, t).

D

2.1. The anisotropic Hölder space H

( [0, 1]

) of order α = (α

, . . . , α

)

consists of functions f : [0, 1]

→ R such that f 

< ∞, where

 · 

is the anisotropic Hölder norm max



D

f

+

d k=1



D

D

f

: r ∈ N

,

d k=1

r

/α

< 1



(2.1)

and e

∈ R

has 1 in the kth position and zero elsewhere.

Let α

be the harmonic mean of (α

, . . . , α

)

, that is, α

= d

α

. For x = (x

, . . . , x

)

, we define b

(x) = (B

(x

) · · · B

(x

), 1 ≤ j

≤ J

, k = 1, . . . , d) to be a collection of J =

J

tensor-product of B-splines, where B

(x

) is the kth component B-spline of fixed order q

≥ α

, with knot sequence 0 = t

< t

< · · · < t

< t

= 1, and let J

= q

+ N

and J = (J

, . . . , J

)

. In the prior construction, the knots depend on n and N

in- creases to infinity with n subject to

J

≤ n. At each k = 1, . . . , d, define δ

= t

− t

to be the one-step knot increment, and let 

= max

δ

be the mesh size. We assume that the knot sequence for each direction is quasi- uniform (Definition 6.4 of [24]), that is, 

/ min

δ

≤ C, for some C > 0.

This assumption is satisfied for the uniform and nested uniform partitions as spe- cial cases (Examples 6.6 and 6.7 of [24]) and we can choose a subset of knots from any given knot sequence to form a quasi-uniform sequence with C = 3 (Lemma 6.17 of [24]).

If the design points X

= (X

, . . . , X

)

for i = 1, . . . , n, are fixed, assume that there exists a cumulative distribution function G(x), with positive and contin- uous density on [0, 1]

such that

sup

x∈[0,1]^d



G

(x) − G(x)

= o

k=1

N

, (2.2)

where G

(x) = n

1

k=1[0,Xik]

(x) is the empirical distribution of {X

, i = 1, . . . , n }, with 1

( ·) the indicator function on U.

R

2.1. For example, let n = m

for some m ∈ N, the discrete uniform

design X

∈ {(j − 1)/(m − 1) : j = 1, . . . , m}

with i = 1, . . . , n, satisfies ( 2.2)

with G being the uniform distribution on [0, 1]

and N

 n

for k =

1, . . . , d.

For random design points, assume X

∼ G with a continuous density on [0, 1]

, then (2.2) holds with probability tending to one if N

 n

for k = 1, . . . , d, and α

> d/2 by Donsker’s theorem. In this paper, we shall prove re- sults on posterior contraction rates and credible sets based on fixed design points;

the random case can be treated by conditioning on X

, i = 1, . . . , n.

Let B = (b

(X

), . . . , b

(X

))

. Each entry of B

B is indexed by d -dimensional multi-indices, that is, for u = (u

, . . . , u

)

and v = (v

, . . . , v

)

with 1 ≤ u

, v

≤ J

, k = 1, . . . , d, the (u, v)th entry is (B

B)

=

_n i=1_d

k=1

B

(X

)B

(X

). The following generalization of matrix band- ing property will be useful.

D

2.2. Let A = ((a

)) be a matrix with rows and columns in- dexed by d-dimensional multi-indices 1

≤ u, v ≤ J, respectively, where arrange- ment of the elements are arbitrary. We say that A is h = (h

, . . . , h

)

banded if a

= 0 whenever |u

− v

| > h

for some k = 1, . . . , d.

Given X

= (X

, . . . , X

)

for i = 1, . . . , n, such that X

∈ [t

, t

], only q

adjacent basis functions (B

(X

), . . . , B

(X

))

will be nonzero for k = 1, . . . , d. Hence, if |u

− v

| > q

for some m = 1, . . . , d, then B

(X

)B

(X

) = 0, and we conclude (B

B)

= 0. It then follows that B

B is q = (q

, . . . , q

)

-banded.

Since approximation results for anisotropic functions by linear combinations of tensor-products of B-splines assume integer smoothness (see Chapter 12, Sec- tion 3 of [24]), we assume that α ∈ N

. For the isotropic case, the norm in (2.1) can be generalized (see Section 2.7.1 of [35]) and the approximation rate is ob- tained for all smoothness levels (Theorem 22 of Chapter XII in [8]). This allows generalization of posterior contraction results for arbitrary smoothness levels. We now describe the assumption on f

used in this paper.

A

1. Under the true distribution P

, Y

= f

(X

) + ε

, such that ε

are i.i.d. sub-Gaussian with mean 0 and variance σ

for i = 1, . . . , n. Also, f

∈ H

( [0, 1]

) with order α = (α

, . . . , α

)

∈ N

. If the design points are de- terministic, we assume that (2.2) holds. If the design points are random, we assume that α

> d/2.

Let E

( ·) and Var

( ·) be the expectation and variance operators taken with re- spect to P

. We write Y = (Y

, . . . , Y

)

, X = (X

, . . . , X

)

, F

= (f

(X

), . . . , f

(X

))

and ε = (ε

, . . . , ε

)

.

3. Prior and posterior conjugacy. We induce a prior on f by represent-

ing it as a tensor-product B-splines series, that is, f (x) = b

(x)

θ , where

θ = {θ

: 1 ≤ j

≤ J

, k = 1, . . . , d} are the basis coefficients. Then its r = (r

, . . . , r

)

mixed partial derivative is

D

f (x) =

J1

 j₁=1

· · ·

Jd

 j_d=1

θ

k=1

∂

∂x

B

(x

).

Define an operator D

k

acting on θ

such that D

θ

= θ

, and for r

≥ 1,

D

θ

(3.1)

= D

θ

− D

θ

(t

− t

)/(q

− r

) . Furthermore, let D

θ

= D

· · · D

θ

be the application of D

to θ

for all direction k = 1, . . . , d. Using equations (15) and (16) of Chapter X from [8], D

f (x) can be written as

J₁−r1

j₁=1

· · ·

J_d−rd

j_d=1

D

θ

k=1

B

(x

) = b

(x)

W

θ , (3.2)

where W

is a

(J

− r

) ×

J

matrix, with entries given by (A.1)–

(A.4) in Lemma A.2. These entries are coefficients associated with applying the weighted finite differencing operator of (3.1) iteratively on θ in all directions.

We represent the model in (1.1) by Y|(X, θ, σ) ∼ N

(Bθ , σ

I

). In this pa- per, we treat J = (J

, . . . , J

)

as deterministic and allow it to depend on n, d and α. On the basis coefficients, we assign θ |σ ∼ N

(η, σ

), where η

is uniformly bounded. The entries of  do not depend on n, and are indexed using d -dimensional multi-indices described above. We further assume that 

is a m = (m

, . . . , m

)

banded matrix with fixed m. Note that  depends on n only through its dimension J × J . Furthermore, as n → ∞, we assume that there exists constants 0 < c

≤ c

< ∞ such that

c

I

≤  ≤ c

I

. (3.3)

It follows that D

f |(Y, σ) ∼ GP(A

Y + c

η, σ

), where A

, c

and the covari- ance kernel are defined for x, y ∈ [0, 1]

by

A

(x) = b

(x)

W

B

B + 

B

, (3.4)

c

(x) = b

(x)

W

B

B + 



, (3.5)

(x, y) = b

(x)

W

B

B + 

W

b

(y).

(3.6)

Since the posterior mean of D

f is an affine transformation of Y, Assumption 1

implies that A

Y + c

η is a sub-Gaussian process under P

. If r = 0, defining

W

= I

, we obtain the conditional posterior distribution of f .

To deal with σ , observe that Y |σ ∼ N

[Bη, σ

(BB

+ I

) ]. Maximizing the corresponding log-likelihood with respect to σ leads to



σ

= n

(Y − Bη)

BB

+ I

(Y − Bη).

(3.7)

Empirical Bayes then entails substituting the maximum likelihood estimator

σ

for σ in the conditional posterior of D

f , that is,



D

f |Y, σ

|

= 

D

f |Y

∼ GP

A

Y + c

η,

σ

. (3.8)

In a hierarchical Bayes approach, we further endow σ with a continuous and pos- itive prior density. A conjugate inverse-gamma (IG) prior σ

∼ IG(β

/2, β

/2), with hyperparameters β

> 4 and β

> 0 is particularly convenient for both com- putation and theoretical analysis since by direct calculations, the posterior of σ

is

σ

|Y ∼ IG

(β

+ n)/2,

β

+ n

σ

/2

. (3.9)

Under the quasi-uniformity of the knots and (2.2), Lemma A.9 concludes that there exist constants 0 < C

≤ C

< ∞ such that

C

n

k=1

J

I

≤ B

B ≤ C

n

k=1

J

I

. (3.10)

In particular, B

B 

 n

J

. Combining the above with (3.3),

C

n

k=1

J

+ c

≤ λ

B

B + 

≤ λ

B

B + 

(3.11)

≤

C

n

J

+ c

.

4. Posterior contraction rates. To establish posterior contraction rates for f and its mixed partial derivatives with unknown σ , a key step is showing that the empirical Bayes estimator for σ in the empirical Bayes approach or the posterior distribution of σ in the hierarchical Bayes approach, are consistent, uniformly for the true regression function f

satisfying f



≤ R for any given R > 0.

P

4.1. Let J

 n

, k = 1, . . . , d. Then for any R >

0, the following assertions holds uniformly for the true regression f

satisfying

f



≤ R:

(a) the empirical Bayes estimator

σ

converges to the true σ

at the rate

max(n

, n

);

(b) if the inverse gamma prior IG(β

/2, β

/2) is used on σ

, then the posterior for σ contracts at σ

at the same rate;

(c) if the true distribution of the regression errors ε

, . . . , ε

is Gaussian, then for any prior on σ with positive and continuous density, the posterior distribution of σ is consistent.

For the rest of the paper, we shall treat f and its mixed partial derivatives in a unified framework by viewing f as D

f . Then the results on posterior con- traction and credible sets (Section 5) for f can be recovered by setting r = 0.

Since an explicit expression for the conditional posterior of D

f given σ is available due to the normal-normal conjugacy, we derive contraction rates by directly bounding posterior probabilities of deviations from the truth uniformly for σ in a shrinking neighborhood of σ

, which suffices in view of the con- sistency of the empirical Bayes estimator or that of the posterior distribution of σ . A decomposition of the posterior mean square error into posterior vari- ance, variance and squared bias of the posterior mean is used for pointwise contraction, and uniformized using maximal inequalities to establish contrac- tion with respect to the supremum distance. Contraction rates below are uni- form in f



≤ R. We write 

= n

and 

= (log n/n)

. Observe that for 

and 

to approach 0 as n → ∞, we will need

(r

/α

) < 1. For the hierarchical Bayes approach, we do not restrict to the inverse gamma prior for σ

but throughout assume that its posterior is consistent uniformly for f



≤ R for any R > 0.

T

4.2 (Pointwise contraction). If J

 n

for k = 1, . . . , d, then for any x ∈ [0, 1]

and M

→ ∞,

Empirical Bayes: E



D

f (x) − D

f

(x)

> M



|Y

→ 0.

Hierarchical Bayes: E



D

f (x) − D

f

(x)

> M



|Y

→ 0.

R

4.3. The above rate of contraction holds for the L

-distance as well under the same set of assumptions for both empirical and hierarchical Bayes ap- proaches. This follows since the posterior expectation of the squared L

-norm can be bounded by the integral of the corresponding uniform estimates of the pointwise case obtained in the proof of Theorem 4.2.

T

4.4 (L

-contraction). If J

 (n/ log n)

for k = 1, . . . , d, then for any M

→ ∞,

Empirical Bayes: E



D

f − D

f

∞

> M



|Y

→ 0.

Hierarchical Bayes: E



D

f − D

f

> M



|Y

→ 0.

Note that an extra logarithmic factor appears in the L

-rate in agreement with the corresponding minimax rate for the problem (see [32, 33]). A similar result for the white noise model using a prior based on wavelet basis expansion for the signal function is given by Theorem 1 of [14] for known variance. It is interesting to note that given any notion of posterior contraction and smoothness index, the same optimal J

, k = 1, . . . , d, applies to f and its mixed partial derivatives, so the Bayes procedure automatically adapts to the order of the derivative to be estimated.

5. Credible sets for f and its mixed partial derivatives. We begin by con- structing pointwise credible set for D

f (x) at x ∈ [0, 1]

, where r ∈ N

satisfies

_d

k=1

(r

/α

) < 1. Let γ

∈ [0, 1] be a sequence such that γ

→ 0 as n → ∞. De- fine z

to be the (1 − δ)-quantile of a standard normal. Since (D

f (x)|Y, σ) ∼ N(A

(x)Y + c

(x)η, σ

(x, x)), we can construct a (1 − γ

)-pointwise credible interval for D

f (x) from the relation



g :

g(x) − A

(x)Y − c

(x)η

≤ z

σ

(x, x) |Y, σ

= 1 − γ

. However, as σ is unknown, we use empirical Bayes by substituting σ by

σ

derived in (3.7), leading to the following empirical credible set:

C

(x) =

g :

g(x) − A

(x)Y − c

(x)η

≤ z

σ



(x, x)

.

For the hierarchical Bayes approach, the resulting credible region is given by C

(x) = {g : |g(x) − A

(x)Y − c

(x)η | ≤ R

(x) }, where R

(x) is the (1 − γ

)-quantile of the marginal posterior distribution of |D

f (x) − A

(x)Y − c

(x)η | after integrating out σ with respect to its posterior distribution. If the con- jugate inverse-gamma prior is used on σ

, then the cut-off may be expressed ex- plicitly in terms of quantiles of a generalized t-distribution. In general, the cut-off value R

(x) may be found by posterior sampling: generate σ from its marginal posterior distribution and D

f |(Y, σ) ∼ GP(A

Y + c

η, σ

).

T

5.1 (Pointwise credible intervals). If J

 n

, k = 1, . . . , d, then for γ

→ 0, the coverage of C

(x) tends to 1 and its radius is O

(

√ log (1/γ

)) at x ∈ [0, 1]

uniformly on f



≤ R.

If the posterior distribution of σ is consistent, then the same conclusion holds for the hierarchical Bayes credible set C

(x).

R

5.2. We can also define a (1 − γ

)-credible set in the L

-norm for

D

f given Y and σ as the set of all functions which differ from A

(x)Y + c

(x)η

in the L

-norm by σ h

, where h

is the 1 − γ

quantile of the L

-norm

of GP(0,

). Then the empirical Bayes credible set is obtained by substituting σ

by

σ

. The hierarchical Bayes credible set is obtained by replacing σ h

by

the 1 −γ

quantile of D

f −A

Y −c

η 

. Both credible regions have asymptotic

coverage 1 under the assumptions in Theorem 5.1.

The (1 −γ

)-empirical Bayes L

-credible set for D

f , can be expressed as {g :

g − A

Y − c

η 

≤ ρ

σ

h

}, where h

is the (1 − γ

)-quantile of the L

-norm of GP(0,

). It turns out that in order to obtain adequate frequentist coverage, this natural credible ball needs to be slightly inflated by a factor ρ

, leading to the inflated empirical Bayes credible region

C

=

g : g − A

Y − c

η 

≤ ρ

σ

h



. (5.1)

On the other hand, unlike in the pointwise or the L

-credible regions, we need not make γ

→ 0, but can allow any fixed γ < 1/2. In the hierarchical Bayes approach, we consider the analogous credible ball C

= {g : g − A

Y − c

η 

≤ ρ

R

}, where R

stands for the (1 − γ )-quantile of the marginal pos- terior distribution of D

f − A

Y − c

η 

integrating out σ with respect to its posterior distribution.

T

5.3 (L

-credible region). If J

 (n/ log n)

for k = 1, . . . , d, then for any ρ

→ ∞ and γ < 1/2, the coverage of C

tends to 1 and its radius is O

(

ρ

) uniformly in f



≤ R. Moreover, if the true distribution of the regression errors is Gaussian, then we can let ρ

= ρ for some sufficiently large constant ρ > 0.

If the posterior distribution of σ is consistent, then the same conclusion holds for the hierarchical Bayes L

-credible ball C

.

R

5.4. To control the size of C

and ensure guaranteed frequen- tist coverage, we can take ρ

to be a factor slowly tending to infinity, or a suffi- ciently large constant for the Gaussian situation. A similar correction factor was also used by [34] in the context of adaptive L

-credible region.

6. Simulation. We compare finite sample performance of pointwise cred- ible intervals and L

-credible bands for f in one dimension (i.e., d = 1, r = 0) with confidence intervals and L

-confidence bands proposed by The- orem 4.1 of [36]. Following [18], we consider the true function f

(x) =

√ 2

i

sin i cos {(i − 1/2)πx}, x ∈ [0, 1], which has smoothness α = 1.

We observed the signal f

with i.i.d. N(0, 0.1) errors at covariate values at X

= (i − 1)/(n − 1) for i = 1, . . . , n. We use cubic B-splines (i.e., q = 4) with uniform knot sequence, where we added 4 duplicate knots at 0 and 1. For the prior param- eters, we set η = 0 and  = I

. We construct (1 − γ

)-empirical credible intervals for γ

= 5/n with

σ

computed using (3.7). The corresponding confidence regions are constructed using Theorem 4.1 of [36] based on the least squares estimator f (x)

= b

(x)

θ for

θ = (B

B)

B

Y, and

σ

= (Y − B

θ )

(Y − B

θ )/(n − J ).

In the Bayesian context when the smoothing parameter J is to be determined from

the data, it is natural to use its posterior mode. However, for a fair comparison, we

used leave-one-out cross validation to determine J for both methods and also ob-

served that the posterior mode essentially chose the same values. We conduct our

FIG. 1. Pointwise coverage probabilities for credible and confidence intervals. The y-axis is the coverage probabilities and the x-axis is the covariate x.

experiment across sample sizes n = 100, 300, 500, 700, 1000, 2000. For pointwise credible and confidence intervals, we report the empirical coverage based on 1000 Monte Carlo runs for each n. All simulations were carried out in R using the bs function from the splines package.

The coverage probabilities of pointwise credible and confidence intervals are

shown in Figure 1. One distinguishing feature is the downward spike at around the

bump of f

at x = 0.3 in Figure 2, and the plots narrow down to this point as n

increases. Moreover, the pointwise coverage is 0 at this point for both Bayesian

and frequentist methods in all sample sizes considered. This phenomenon occurs

perhaps due to the fact that the true function at x = 0.3 has a sharp bend but the

function is much smoother elsewhere, so based on a limited sample the cross-

validation method oversmooths by choosing a smaller J than ideal. Both methods

yield almost the same pointwise coverage for large sample sizes, and are equivalent

in quantifying uncertainty of estimating f

. To cover the function at all points, we

consider the simultaneous (modified) credible band at the level 1 −γ = 0.95, given

by (A

(x)Y + c

(x)η) ± ρ

σ

h

.

FIG. 2. Dots: posterior mean (top) andf(bottom), Solid: true function, Dashes: 95% L_∞-credible (top) and confidence (bottom) bands.

The second assertion of Theorem 5.3 allows us to use a fixed ρ because our true errors are normally distributed which we choose as ρ = 0.5. To construct (1 − γ )- asymptotic confidence band, we use Theorem 4.2 of [36].

Table 1 shows the coverage of 95% simultaneous credible and confidence bands.

At n = 100, the apparent higher coverage of the confidence bands is due to the pos- itive bias of

σ

for small n. From n = 300 onward, the coverage of both credible and confidence bands steadily increase with n. The corresponding graphical rep-

TABLE1

95% simultaneous credible and confidence bands

n 100 300 500 700 1000 2000

Credible band coverage 0.852 0.896 0.954 0.945 0.964 0.972 Confidence band coverage 0.972 0.948 0.963 0.978 0.985 0.986

Credible band radius 0.235 0.155 0.148 0.127 0.121 0.098

Confidence band mean radius 0.27 0.165 0.147 0.132 0.129 0.101 Confidence band max radius 0.64 0.436 0.409 0.374 0.372 0.3

resentations of these bands are shown in Figure 2. The top panel corresponds to the proposed Bayesian method, where the dotted line stands for the posterior mean and dashed lines for the 95% credible band. The bottom panel corresponds to the frequentist method of [36], where the dotted line standing for the least squares esti- mator f

and the dashed lines for the 95% L

-confidence bands. In both panels, the solid line is the true function f

. Observe that the credible bands have fixed length, while the confidence bands have varying lengths. This is because the procedure of [36] is based on the supremum of the scaled absolute differences. Therefore, for the latter we present both average and maximum radius. The frequentist method has larger width at the endpoints due to the fact that there are fewer observations, and this results in larger maximum radius.

7. Proofs. We shall repeatedly use the following fact about approximation power of tensor product B-splines given by (12.37) of [24].

For any R > 0, if f



≤ R, there exists a θ

∈ R

such that for constant C > 0 depending only on α, q and d, we have



b

( ·)

θ

− f

≤ C

d

k=1

J

∂

∂x

f

∞



d

k=1

J

. (7.1)

Since b

( ·)

θ



≤ f



+ C

J

f



 R + d, sup

f0α,∞≤R

θ



 sup

f0α,∞≤R



b

(·)

θ

= O(1), (7.2)

by (12.25) of [24]. An extension of the approximation result for derivatives is given by the following lemma.

L

7.1. There exists C > 0 depending only on α, q and d such that for f

∈ H

( [0, 1]

),



b

( ·)

W

θ

− D

f

≤ C

 k=1

J

D

D

f

.

P

. Let I

=

[t

, t

]. Define a bounded linear opera- tor Qf (x) =

· · ·

(λ

f )

B

(x

) on H

([0, 1]

), where λ

=

λ

and λ

is the dual basis of B

( ·), that is, λ

is a linear functional such that λ

B

( ·) = 1

( ·) for k = 1, . . . , d (see Section 4.6 of [24]). Using Theorem 13.20 of [24], there exists a tensor-product Taylor’s polyno- mial p

(x) such that



D

(f

− p

)|

∞

≤ C

d

k=1

J

D

D

f

|

∞

,

where f |

is the restriction of f onto I

and C > 0 depends only on α, q and d. By equations (12.30) and (12.31) of Theorem 12.6 in [24], (D

f

− QD

f

) |



is bounded above by



D

(f

− p

)|

+

Q

D

f

− D

p

|

≤ C

D

(f

− p

) |

∞

≤ C

d k=1

J

D

D

f

|

∞

.

Since QD

f

= D

Qf

, identifying (θ

)

from (7.1) with λ

f

and applying equations (15) and (16) of Chapter X in [8], we see that QD

f

= b

( ·)

W

θ

. Now sum both sides over 1 ≤ j

≤ J

, k = 1, . . . , d. 

P

P

4.1. Define U = (BB

+ I

)

and J =

J

. By equation (33) of page 355 in [26],



E



σ

− σ

=

n

σ

tr(U) − σ

+ n

(F

− Bη)

U(F

− Bη)

 n

tr(I

− U) + (F

− Bθ

)

U(F

− Bθ

) (7.3)

+ (Bθ

− Bη)

U(Bθ

− Bη)

,

where we used (x + y)

D(x + y) ≤ 2x

Dx + 2y

Dy for any D ≥ 0. Let P

= B(B

B)

B

. Let A be an m × m matrix, C an m × r matrix, T an r × r matrix, and W an r × m matrix, with A and T invertible. Then by the binomial inverse theorem (Theorem 18.2.8 of [15])

(A + CTW)

= A

− A

C

T

+ WA

C

WA

. (7.4)

Therefore, two applications of (7.4) to U yield

BB

+ I

₋₁

= I

− B

B

B + 

B

= I

− P

+ V, (7.5)

where V = B(B

B)

[ + (B

B)

]

(B

B)

B

≥ 0. Hence, the first term in (7.3) is

n

tr(P

− V) ≤ n

tr(P

) = J/n.

(7.6)

Note U ≤ I

since BB

≥ 0, and the second term in (7.3) is bounded by n

U

F

− Bθ



≤ F

− Bθ





d

k=1

J

, (7.7)

in view of (7.1). By (7.5) and (I − P

)B = 0, the last term in ( 7.3) is n

(θ

− η)

[ + (B

B)

]

(θ

− η), which is bounded above by

n

c

+ C

J /n

J θ

− η

 J/n,

(7.8)