# Credible intervals and bootstrap confidence intervals in monotone regression

Piet Groeneboom and Geurt Jongbloed  Delft Institute of Applied Mathematics,
Delft University of Technology
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In the recent paper [5], a Bayesian approach for constructing confidence
intervals in monotone regression problems is proposed, based on credible
intervals. We view this method from a frequentist point of view, and show that
it corresponds to a percentile bootstrap method of which we give two versions.
It is shown that a (non-percentile) smoothed bootstrap method has better
behavior and does not need correction for over- or undercoverage. The proofs
use martingale methods.

62G05,

62N01,

62-04,

credible intervals,

confidence intervals,

Chernoff’s distribution,

smoothed bootstrap,

monotone regression,

###### keywords:

[class=AMS]

###### keywords:

##  1 Introduction

In many fields of application, inference on a monotone function is both
natural and needed. Many examples of applications can be found in the books on
shape constrained statistical inference [1], [16], [19] and [8]. Natural
estimators of monotone functions exist. In order to assess the accuracy of
these, there is a need for uncertainty quantification. In the frequentist
setting, this can be done by (pointwise) confidence intervals for the monotone
function of interest. From a Bayesian perspective, credible intervals are a
common method to quantify uncertainty.

We consider the following model for our observations
$D_{n}=\\{(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\\}$, also considered in [5].

$\displaystyle Y_{i}=f_{0}(X_{i})+\epsilon_{i},\qquad i=1,\dots,n.$

Here $f_{0}:[0,1]\mapsto\mathbb{R}$ is monotone nondecreasing, the
$\varepsilon_{i}$ are i.i.d. sub-Gaussian with expectation $0$ and variance
$\sigma_{0}^{2}$, independent of the $X_{i}$’s, and the $X_{i}$ are i.i.d.
with non-vanishing density $g$ on $[0,1]$. The classical least squares
estimator (LSE) $\hat{f}_{n}$ of $f_{0}$, under the condition that
$\hat{f}_{n}$ is nondecreasing minimizes

$\displaystyle\sum_{i=1}^{n}\left(Y_{i}-f(X_{i})\right)^{2}$ (1.1)

over all nondecreasing functions $f:[0,1]\mapsto\mathbb{R}$. This estimator,
which is also the maximum likelihood estimator if we assume the $\epsilon_{i}$
to be i.i.d. centered normally distributed, can be explicitly constructed
based on the data. Denoting by $0<x_{1}<x_{2}<\cdots<x_{n}<1$ the observed
ordered $X_{i}$’s and by $y_{i}$ the corresponding $Y$ values (so
$(x_{1},y_{1}),\ldots,(x_{n},y_{n})$ represent the original data pairs, but
with first coordinate sorted in increasing order), Lemma 2.1 in [8] shows that
$\hat{f}_{n}$ can be taken piecewise constant on the intervals
$(x_{i-1},x_{i}]$, $1\leq i\leq n$ and that $\hat{f}_{n}(x_{i})$ can be
obtained as the left derivative of the greatest convex minorant of the diagram
consisting of the points $\\{P_{i}\,:\,0\leq i\leq n\\}$ with

$P_{0}=(0,0)\mbox{ and for }1\leq i\leq
n,\,\,\,P_{i}=\left(\frac{i}{n},\frac{1}{n}\sum_{j=1}^{i}y_{j}\right),$

evaluated at the point $P_{i}$.

As is clear from this characterization, $\hat{f}_{n}$ will be a nondecreasing
step function with its jumps concentrated on a data-dependent subset of the
observed points $x_{i}$. Fixing a number of points in $[0,1]$, say
$0=\tau_{0}<\tau_{1}<\tau_{2}<\cdots<\tau_{m}=1$, one can also minimize (1.1)
over all nondecreasing functions, piecewise constant on the intervals
$I_{j}=(\tau_{j-1},\tau_{j}]$. Then, writing $n_{j}=|\\{i\,:\,i\in I_{j}\\}|$
and $\bar{y}_{j}=(\sum_{x_{i}\in I_{j}}y_{i})/n_{j}$, we have

$\displaystyle\sum_{i=1}^{n}\left(y_{i}-f(x_{i})\right)^{2}$
$\displaystyle=\sum_{j=1}^{m}\sum_{x_{i}\in
I_{j}}\left(y_{i}-\bar{y}_{j}+\bar{y}_{j}-f(\tau_{j})\right)^{2}=$
$\displaystyle=\sum_{j=1}^{m}\sum_{x_{i}\in
I_{j}}\left(y_{i}-\bar{y}_{j}\right)^{2}+\sum_{j=1}^{m}\left(\bar{y}_{j}-f(\tau_{j})\right)^{2}n_{j},$
(1.2)

where we use that for the current function class, $f(x_{i})=f(\tau_{j})$ for
$x_{i}\in I_{j}$. As the first term in (1.2) does not involve $f$, minimizing
it boils down to a weighted isotonic regression. The solution to this
minimization problem also allows for a graphical construction. The optimal
value of $f(\tau_{j})$ is the left derivative, taken at the point $P_{j}$, of
the greatest convex minorant of the diagram of points consisting of

$\displaystyle
P_{0}=(0,0),\,\,P_{j}=\left(\frac{1}{n}\sum_{k=1}^{j}n_{k},\frac{1}{n}\sum_{k=1}^{j}n_{k}\bar{y}_{k}\right),\,\,\,1\leq
j\leq m.$

From decomposition (1.2) it also follows that the piecewise constant function
$f$ defined by $f(x_{i})=\bar{y}_{j}$ for $x_{i}\in I_{j}$ is the least
squares estimator over the class of piecewise constant functions without
imposing the restriction of monotonicity.

In [5], a Bayesian method is proposed for constructing pointwise confidence
intervals for a monotone regression function, based on credible intervals. The
method is proved to give overcoverage for large sample sizes, but a correction
table is given in [5] to correct for the overcoverage. Purely based on the
algorithm that results in the credible intervals, the approach can be seen as
a particular percentile bootstrap method.

In Section 2 we describe the approach in [5] to construct confidence intervals
via credible intervals. In Section 3 we give the interpretation of the
credible intervals as percentile bootstrap intervals and in particular Theorem
3.1 for the bootstrap procedure, corresponding to the key Theorem 3.3 in [5]
for the construction of the credible intervals. In proving Theorem 3.1 we use
a martingale method.

In analogy with the Bayesian procedure, we construct the bootstrap intervals
by generating normal noise variables (following [5]), using the empirical
Bayes method for estimating the variance of these variables, defined in
Section 3. In subsection 3.2 we define a classical bootstrap procedure, where
we resample with replacement from the original data, and do not have to
estimate the variance. These two methods correspond, respectively, to the
“regression method” (holding the regressors $X_{i}$ in the regression model
fixed), and the “correlation model” (where we consider the $X_{i}$ as random)
in the terminology of [12]. The results of the three methods are highly
similar.

It has been proved by several authors that the straightforward bootstrap is
inconsistent in this situation (see, e.g., [13], [18] and [17] for results
related to this phenomenon).This straightforward bootstrap uses resampling
with replacement from the pairs $(X_{i},Y_{i})$ and computes the monotone
least squares estimator $\hat{f}_{n}^{*}$ based bootstrap samples and
approximates the distribution of
$n^{1/3}\left(\hat{f}_{n}(t_{0})-f_{0}(t_{0})\right)$ by that of the analogous
‘bootstrap quantity’
$n^{1/3}\left(\hat{f}_{n}^{*}(t_{0})-\hat{f}_{n}(t_{0})\right)$. The Bayesian
approach and the percentile bootstrap approach circumvent this difficulty by
using the convergence in distribution of the random variable (as a function of
$D_{n}=\\{(X_{1},Y_{n}),\dots,(X_{n},Y_{n})$})

$\displaystyle{\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
x\Bigm{|}D_{n}\right)$ (1.3)

to

$\displaystyle{\mathbb{P}}\left[\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{\rm
argmin}_{t\in\mathbb{R}}\left[W_{1}(t)+W_{2}(t)+t^{2}\right]\leq
x\Bigm{|}W_{1}\right]$ (1.4)

see Theorem 3.3 in [5] and Theorem 3.1, where $W_{1}$ and $W_{2}$ are two
independent standard two-sided Brownian motions, originating from zero. Here
$\hat{f}_{n}^{*}$ is either the projection-posterior Bayes estimate (to be
described in Section 2), in which case we would write

$\displaystyle\Pi\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
x\Bigm{|}D_{n}\right)$

instead of (1.3), or the percentile bootstrap estimate $\hat{f}_{n}^{*}$. The
limit (1.4) leads to credible intervals which asymptotically give
overcoverage, which can be corrected for as described in [5].

In Section 4 an altogether different method for constructing the confidence
intervals is given, where we use the smoothed (non-percentile) bootstrap. Here
we keep the regressors fixed again, and resample with replacement residuals
w.r.t. a smooth estimate of the regression function: the Smoothed Least
Squares Estimator (SLSE). In this way, using theory from [9], consistent
confidence intervals are constructed.

In fact, instead of
$n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}=n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\hat{f}_{n}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}}$
we can now consider

$\displaystyle
n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}},$

where $\hat{f}_{n}^{*}$ is based on sampling with replacement from the
residuals w.r.t. the SLSE $\tilde{f}_{nh}$ with bandwidth $h$ of order
$n^{-1/5}$. In contrast with Theorem 3.3 in [5] and Theorem 3.1 in the present
paper, we now have convergence of

$\displaystyle{\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
0\Bigm{|}D_{n}\right)$

to the uniform distribution on $[0,1]$ (using the symmetry of the limit
distribution of $n^{1/3}\\{\hat{f}_{n}(t_{0})-f_{0}(t_{0})\\}$), see Theorem
4.2.

For the non-percentile bootstrap, however, it is more natural to consider

$\displaystyle
n^{1/3}\bigl{\\{}\hat{f}_{n}(t_{0})-\\{\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\\}-f_{0}(t_{0})\bigr{\\}}$

(avoiding “looking up the wrong tables, backwards”, see the discussion on p.
938 of [11]), for which we also get convergence to the the uniform
distribution of

$\displaystyle{\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}(t_{0})-\\{\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\\}-f_{0}(t_{0})\bigr{\\}}\leq
0\Bigm{|}D_{n}\right),$

implying the consistency of the smoothed bootstrap method. This method of
constructing confidence intervals seems superior in comparison to the Bayesian
method and the percentile bootstrap intervals, as is suggested by our
simulations of the coverage of the different methods.

##  2 Credible intervals

In [5], a Bayesian approach to construct confidence intervals for a monotone
regression function is proposed. A prior distribution is defined on the class
of functions on $[0,1]$, supported on a sieve of piecewise constant functions.
More specifically, the interval $[0,1]$ is partitioned into $J$ intervals
$I_{j}=((j-1)/J,j/J]$, $1\leq j\leq J$. In the notation of the previous
section, $\tau_{j}=j/n$. A draw from the prior distribution is then
represented by

$\displaystyle
f_{\bm{\theta}}=\sum_{j=1}^{J}\theta_{j}1_{I_{j}},\qquad\bm{\theta}=(\theta_{1},\dots,\theta_{J}),$
(2.1)

where the $\theta_{j}$ are independent normal random variables with
expectation $\zeta_{j}$ and variance $\sigma_{0}^{2}\lambda_{j}^{2}$, where
$0<\lambda_{j}<\infty$ (including noise variance $\sigma_{0}^{2}$ as a factor
is only done for convenience in formulas to follow). Note that function (2.1)
will not automatically be monotone, a requirement that would seem natural in
this setting. The main reason not to impose this, is that with this prior
distribution, the posterior distribution can be conveniently analytically
computed. Indeed, as seen in the Appendix, the posterior distribution of
$\bm{\theta}$ has independent coordinates, where $\theta_{j}$ has distribution

$\displaystyle\theta_{j}\sim
N\left(\frac{n_{j}\bar{y}_{j}+\zeta_{j}/\lambda_{j}^{2}}{n_{j}+1/\lambda_{j}^{2}},\frac{\sigma_{0}^{2}}{n_{j}+1/\lambda_{j}^{2}}\right),\qquad
n_{j}=\sum_{i=1}^{n}1_{I_{j}}(x_{i}).$ (2.2)

Here, as before, $\bar{y}_{j}$ is the mean of the $y_{i}$ for the $x_{i}$
belonging to the $j$-th interval $I_{j}$. As mentioned in the previous
section, this corresponds to the MLE of $f$ over the (nonrestricted) class of
functions which are constant on the intervals $I_{j}$. A draw from the
posterior on the set of piecewise constant functions on $[0,1]$ proceeds via
(2.1), based on a draw from the posterior of ${\bm{\theta}}$. The resulting
function will in general not be monotone, so the support of the posterior
extends outside the set of monotone functions on $[0,1]$.

Following ideas of [14], [2] and [3], in [5] a draw $f_{\bm{\theta}}$ from the
‘raw posterior’ is subsequently projected on the set of nondecreasing
functions on $[0,1]$, piecewise constant on the intervals $I_{j}$, $1\leq
j\leq J$, via weighted isotonic regression. This projection $f_{\theta}^{*}$
is computed using Lemma 2.1 of [8]. This boils down to computing the left
derivative of the greatest convex minorant of the cusum diagram consisting of
the points $P_{j}$, for $0\leq j\leq J$ with $P_{0}=(0,0)$ and

$\displaystyle
P_{j}=\left(\frac{1}{n}\sum_{k=1}^{j}n_{k},\frac{1}{n}\sum_{k=1}^{j}n_{k}\theta_{k}\right),\,\,1\leq
j\leq J$

if all $n_{j}>0$, see (3.2) in [5] (note that Lemma 2.1 in [8] has the
condition that all weights are strictly positive). It is clear that computing
the isotonic regression can be restricted to those $j$ with $n_{j}>0$ and that
for those $j$ with $n_{j}=0$, $f_{\theta}^{*}$ can be given any value such
that monotonicity is not violated.

In this procedure, various choices need to be made. One is the number of
intervals $J$. In [5], the asymptotic bounds

$\displaystyle n^{1/3}\ll J\ll n^{2/3}$ (2.3)

are given and the closest integer to $n^{1/3}\log n$ is chosen in the
simulations. Here the symbol “$\ll$” means “is of lower order than”, as
$n\to\infty$. Also the noise variance $\sigma_{0}^{2}$ needs to be dealt with.
For this, [5] choose the natural empirical Bayes estimate (fixing $\zeta$ and
$\Lambda$), given by

$\displaystyle\hat{\sigma}_{n}^{2}=n^{-1}(\bm{y}-\bm{B}\bm{\zeta})^{T}(\bm{B}\bm{\Lambda}\bm{B}^{T}+\bm{I})^{-1}(\bm{y}-\bm{B}\bm{\zeta})\qquad\bm{\Lambda}=\text{diag}(\lambda_{1}^{2},\dots,\lambda_{J}^{2}),$
(2.4)

where ${\bm{B}}=(b_{ij})$ the $n\times J$ ‘design matrix’ with entries
$b_{ij}=1_{I_{j}}(x_{i})$, corresponding to the regression model
$\bm{y}=\bm{B}\bm{\theta}+\bm{\epsilon}$ following from the representation
$f_{\bm{\theta}}(x_{i})=({\bm{B}}{\bm{\theta}})_{i}$; see the Appendix. As
also shown in the Appendix, this estimate can be rewritten as

$\displaystyle\hat{\sigma}_{n}^{2}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J}\sum_{x_{i}\in
I_{j}}\left(y_{i}-\bar{y}_{j}\right)^{2}+\frac{1}{n}\sum_{j=1}^{J}\frac{n_{j}(\bar{y}_{j}-\zeta_{j})^{2}}{1+n_{j}\lambda_{j}^{2}}$
(2.5)
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\bar{f}(x_{i})\right)^{2}+\frac{1}{n}\sum_{j=1}^{J}\frac{n_{j}(\bar{y}_{j}-\zeta_{j})^{2}}{1+n_{j}\lambda_{j}^{2}}$

where $\bar{f}(x)=\sum_{j=1}^{J}\bar{y}_{j}1_{I_{j}}(x)$ is the aforementioned
maximum likelihood estimate of $f_{0}$ over all piecewise constant (not
necessarily monotone) functions on the intervals $I_{j}$. The first term in
this expression is the mean of the squared residuals of the observations with
respect to $\bar{f}$. This is a quite natural estimator of the variance. The
influence of the hyper parameters $\bm{\zeta}$ and $\bm{\Lambda}$ on the
estimate of $\sigma_{0}^{2}$ can be inferred from the second term in the
expression.

As shown in the Appendix, the empirical Bayes estimate for $\zeta$, not taking
into account monotonicity is given by

$\displaystyle\bar{\zeta}=(\bar{y}_{1},\bar{y}_{2},\ldots,\bar{y}_{J})^{T}.$
(2.6)

Substituting this in (2.5) makes the second term vanish. Using the empirical
Bayes estimate over the monotone vectors $\zeta$, being the isotonic
regression of $\bar{\zeta}$ with weights $n_{j}/(1+n_{j}\lambda_{j}^{2})$
increases the empirical Bayes estimate for $\sigma_{0}^{2}$.

With the choices $\bm{\zeta}=0$ and $\lambda_{j}\equiv\lambda$ made in [5],
the empirical Bayes estimate for $\sigma_{0}^{2}$ becomes

$\displaystyle\hat{\sigma}_{n}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\bar{f}(x_{i})\right)^{2}+\frac{1}{n}\sum_{j=1}^{J}\frac{n_{j}\bar{y}_{j}^{2}}{1+n_{j}\lambda^{2}}$
(2.7)

For relatively large values of $\lambda^{2}n_{j}$, the second term becomes
negligible to the first.

Because the density $g$ generating the $X_{i}$’s is nonvanishing on $[0,1]$,
the (random) number $N_{j}$ of points in intervals of length of the order
$1/J=1/J_{n}$ is (in the setting of [5]) of the order $n/J_{n}$. With the
restriction $n^{1/3}<<J_{n}<<n^{2/3}$, this means that $N_{j}$ will be of
bigger order than $n^{1/3}$; taking $J_{n}\approx n^{1/3}\log n$, $N_{j}$ will
be of order $n^{2/3}/\log n$. Therefore, for reasonable choice of $\lambda$,
$\lambda^{2}N_{j}>>1$ with high probability when $n$ is large.

Considering (2.2) with $\zeta_{j}\equiv 0$ and $\lambda_{i}\equiv\lambda$,
fixed, as chosen in [5], a draw from the raw posterior of $\theta_{j}$ can be
viewed as

$\theta_{j}=\bar{f}(j/J)+\tilde{\epsilon}_{j}$

where

$\bar{f}(j/J)=\frac{\bar{y}_{j}}{1+1/(n_{j}\lambda^{2})}\mbox{ and
}\tilde{\epsilon}_{j}\sim^{\rm
indep}N\left(0,\frac{\sigma_{0}^{2}/n_{j}}{1+1/(n_{j}\lambda^{2})}\right),$
(2.8)

where $\sigma_{0}^{2}$ is estimated by its Empirical Bayes estimate (2.7).
Again due to the restriction $n^{1/3}<<J_{n}<<n^{2/3}$, $\bar{f}$ in (2.8) is
a (generally non-monotone) local average estimator of $f_{0}$. The added noise
$(\tilde{\epsilon}_{j})$ is normal and reflects the variance of the original
$\bar{Y}_{j}$. This means that the draw from the projected posterior is
computed as left derivative of the cumulative sum diagram consisting of the
points $P_{j}$, $0\leq j\leq J$ with $P_{0}=(0,0)$ and

$\displaystyle
P_{j}=\left(\sum_{k=1}^{j}n_{k},\sum_{k=1}^{j}n_{k}\left(\frac{\bar{y}_{k}}{1+1/(n_{k}\lambda^{2})}+\tilde{\epsilon}_{k}\right)\right)$
(2.9)

In [5], the following example is considered:

$\displaystyle f_{0}(x)=x^{2}+x/5,\qquad x\in[0,1].$

Here the $X_{i}$ are independently uniformly distributed on $[0,1]$ and the
$\varepsilon_{i}$ have a normal $N(0,0.01)$ distribution. The choices for $J$,
$\zeta_{j}$ and $\lambda_{j}$ are $J=\lfloor n^{1/3}\log n\rfloor$,
$\zeta_{j}=0$ and $\lambda_{j}=10$, following the parametrization in the R
code, kindly sent to us by Moumita Chakraborty. A picture of a single draw
$f_{\bm{\theta}}$ from the raw posterior and its isotonic projection
$f_{\bm{\theta}}^{*}$, for a sample of size $n=500$ is shown in Figure 1.

(a)

(b)

Figure 1: (a) A single draw $f_{\bm{\theta}}$ (blue) from the raw posterior
and (b) the corresponding projection $f_{\bm{\theta}}^{*}$ (blue), for a
sample of size $n=1000$ and $f_{0}(x)=x^{2}+x/5$ (red, dashed).

Now, one can generate $1000$ posterior samples of
$\bm{\theta}=(\theta_{1},\dots,\theta_{J})$ from the posterior normal
distribution, specified in (2.2), and consider the $\tfrac{1}{2}\alpha$th and
$(1-\tfrac{1}{2}\alpha)$th percentiles of the isotonic projections
$\hat{f}_{\bm{\theta}}^{*}(t)$ at a fixed point $t$. Would this give us, at
least asymptotically, valid $95\%$ confidence intervals for $f_{0}(t)$?

The question is answered in [5] by Theorem 3.3 on p. 1017:
$\Pi\left(n^{1/3}\\{\hat{f}_{\bm{\theta}}^{*}(t)-f_{0}(t)\\}\leq
z|D_{n}\right)$ converges to a limit distribution, leading to wider intervals
than in the situation in which we have the Chernoff distribution as limit. The
fraction by which they become wider is given in [5].

##  3 Credible intervals as bootstrap percentile confidence intervals

###  3.1 The percentile bootstrap for the regression model

In the Bayes approach, we considered random parameters $\theta_{j}$, with
(posterior) distribution given in (2.2). In the simulations, accompanying the
paper [5], the prior parameter $\bm{\zeta}$ was taken $\bm{\zeta}=\bm{0}$ and
$\lambda_{j}\equiv\lambda>0$. Moreover, the empirical Bayes estimator
$\hat{\sigma}_{n}^{2}$ was taken as estimator for $\sigma_{0}^{2}$.

With these choices, we get:

$\displaystyle\theta_{j}\sim
N\left(\frac{\bar{y}_{j}}{1+1/(n_{j}\lambda^{2})},\frac{\hat{\sigma}_{n}^{2}}{n_{j}\bigl{\\{}1+1/(n_{j}\lambda^{2})\bigr{\\}}}\right),1\leq
j\leq J,\,\,{\rm independently}.$

This means asymptotically, in first order:

$\displaystyle\theta_{j}\sim
N\left(\bar{y}_{j},\frac{\hat{\sigma}_{n}^{2}}{n_{j}}\right),\,\,1\leq j\leq
J$ (3.1)

if we keep $\lambda$ bounded away from zero ($\lambda^{2}=100$ was taken in
the simulations with paper [5]). Next the confidence intervals were determined
by taking the percentiles of simulated values of the (weighted) monotonic
projections of the $\theta_{j}$’s with distribution given by (3.1).

Algorithmically, this can be viewed as a percentile bootstrap method, where a
bootstrap sample is generated by adding noise to an estimate of the regression
function. The regression estimate in this setting is the (weighted) least
squares estimate of $f_{0}$, piecewise constant on intervals $I_{j}$ and not
taking into account the monotonicity constraint (so: $\bar{y}_{j}$ on
$I_{j}$). The noise is sampled from a centered normal distribution with
estimated variance. Then, $\hat{f}^{*}_{n}$ is determined by computing the
(weighted) isotonic regression based on the bootstrap dataset. Adopting the
”bootstrap notation” rather than the ”Bayesian notation” $\theta_{j}$, define
where

$\displaystyle Y_{j}^{*}\sim N(\bar{Y}_{j},\hat{\sigma}_{n}^{2}/N_{j}),\qquad
j=1,\dots,J$ (3.2)

and note that given the original data,
$(Y_{1}^{*},\dots,Y_{J}^{*})=^{D}(\theta_{1},\ldots,\theta_{J})$, in view of
(3.1) and (3.2). Using this notation, $\hat{f}_{n}^{*}$ is found by taking the
left derivative of the convex minorant of the cusum diagram, running through
the points

$\displaystyle
P_{j}^{*}=\left(\sum_{i=1}^{j}N_{i},\sum_{i=1}^{j}N_{i}Y_{i}^{*}\right),\qquad
j=1,\dots,J.$ (3.3)

To study the asymptotic behavior of $\hat{f}^{*}_{n}$, we define the (local)
“bootstrap” process

$\displaystyle\widetilde{W}_{n}^{*}(t)=n^{-1/3}\Biggl{\\{}\sum_{j:j/J\in[0,t_{0}+n^{-1/3}t]}N_{j}\left\\{Y_{j}^{*}-\bar{Y}_{j}\right\\}-\sum_{j:j/J\in[0,t_{0}]}N_{j}\left\\{Y_{j}^{*}-\bar{Y}_{j}\right\\}\Biggr{\\}}$
(3.4)

and the (local) “sample” process $\widetilde{W}_{n}$ by:

$\displaystyle\widetilde{W}_{n}(t)$
$\displaystyle=n^{-1/3}\Biggl{\\{}\sum_{j:j/J\in[0,t_{0}+n^{-1/3}t]}1_{\\{N_{j}>0\\}}N_{j}\Bigl{\\{}\bar{Y}_{j}-\frac{n}{N_{j}}\int_{I_{j}}f_{0}(u)\,d{\mathbb{G}}_{n}(u)\Bigr{\\}}$
(3.5)
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\sum_{j:j/J\in[0,t_{0}]}1_{\\{N_{j}>0\\}}N_{j}\Bigl{\\{}\bar{Y}_{j}-\frac{n}{N_{j}}\int_{I_{j}}f_{0}(u)\,d{\mathbb{G}}_{n}(u)\Bigr{\\}}\Biggr{\\}}.$
(3.6)

With these definitions we have the following theorem, similar to Theorem 3.3
in [5].

######  Theorem 3.1.

Let $D_{n}=\\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\\}$ and let $\hat{f}_{n}^{*}$
be a draw generated according to the bootstrap procedure described above.
Then, for each fixed $x\in\mathbb{R}$, as $n\to\infty$,

$\displaystyle{\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
x\Bigm{|}D_{n}\right)$
$\displaystyle\stackrel{{\scriptstyle{\mathscr{D}}}}{{\longrightarrow}}{\mathbb{P}}\left[\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{\rm
argmin}_{t\in\mathbb{R}}\left[W_{1}(t)+W_{2}(t)+t^{2}\right]\leq
x\Bigm{|}W_{1}\right]$ (3.7)

where $W_{1}$ and $W_{2}$ are independent standard two-sided Brownian motions.

Note that, for $t>0$,

$\displaystyle
n^{-1/3}\Biggl{\\{}\sum_{j:j/J\in[0,t_{0}+n^{-1/3}t]}N_{j}\bigl{\\{}Y_{j}^{*}-f_{0}(t_{0})\bigr{\\}}-\sum_{j:j/J\in[0,t_{0}]}N_{j}\bigl{\\{}Y_{j}^{*}-f_{0}(t_{0})\bigr{\\}}\Biggr{\\}}$
$\displaystyle=\widetilde{W}_{n}^{*}(t)+\widetilde{W}_{n}(t)+n^{-1/3}\sum_{j:j/J\in[t_{0},t_{0}+n^{-1/3}t]}n\,\int_{I_{j}}\\{f_{0}(u)-f_{0}(t_{0})\\}\,d{\mathbb{G}}_{n}(u)$
$\displaystyle\sim\widetilde{W}_{n}^{*}(t)+\widetilde{W}_{n}(t)+\tfrac{1}{2}f_{0}^{\prime}(t_{0})g(t_{0})t^{2},$
(3.8)

with a similar expansion for $t\leq 0$.

So the percentile bootstrap estimates have the same behavior as the Bayes
estimates in [5]. Histograms of estimates of the posterior probabilities for
the Bayesian procedure in [5] and the corresponding conditional probabilities
of the percentile bootstrap in Lemma 3.1 for varying $D_{n}$ of size $n=2000$
and $t_{0}=0.5$ are shown in Figure 2. The estimates are the relative
frequencies in $1000$ posterior, resp. percentile bootstrap samples for each
of the original (1000) samples.

(a)

(b)

Figure 2: (a) Histogram of estimate of
$\Pi\\{n^{1/3}\bigl{\\{}f_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
0\bigm{|}D_{n}\\}$ for the projection-posterior samples in [5], (b) Histogram
of estimate of
${\mathbb{P}}\\{n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
0\bigm{|}D_{n}\\}$ for the percentile bootstrap. Both histograms are based on
$1000$ samples $D_{n}$ of size $n=2000$ and $t_{0}=0.5$.

######  Remark 3.1.

Let $\Delta_{n}=\Pi\\{n^{1/3}\\{f_{n}^{*}(t_{0})-\hat{f}_{n}(t_{0})\\}\leq
0|D_{n}\\}$. In Figure 1 on p. 1017 of [5] three pictures of $\Delta_{n}$ are
shown for three different sets of simulated data, where $\hat{f}_{n}$ is the
LSE. Is is not completely clear to us how $\Delta_{n}$ is sampled here. Since
we do not have an explict expression for $\Delta_{n}$, it seems that an
estimate of $\Delta_{n}$ has to be based on a sample of posterior draws
$f^{*}(t_{0})$. If we use such a procedure and consider the fluctuation of
$\Delta_{n}$ as a function of $D_{n}$, we get a histogram similar to the
histograms in Figure 1 of [5]. The estimates are relative frequencies in
$1000$ samples of size $2000$. See Figure 3.

Figure 3: Histogram of estimate of
$\Pi\\{n^{1/3}\bigl{\\{}f_{n}^{*}(t_{0})-\hat{f}_{n}(t_{0})\bigr{\\}}\leq
0\bigm{|}D_{n}\\}$, for varying $D_{n}$ of size $n=2000$ and $t_{0}=0.5$. The
estimates are based on relative frequencies in $1000$ draws.

Theorem 3.1 is the consequence of the following two lemmas.

######  Lemma 3.1.

Let $W$ be standard two-sided Brownian motion on $\mathbb{R}$, originating
from zero. Let $D(\mathbb{R})$ be the space of right continuous functions,
with left limits (cadlag functions) on $\mathbb{R}$, equipped with the metric
of uniform convergence on compact sets, and let $t_{0}\in(0,1)$. Let
$\widetilde{W}_{n}^{*}$ be defined by (3.4). Then, along almost all sequences
$(X_{1},Y_{1}),(X_{2},Y_{2}),\dots$, the process $\widetilde{W}_{n}^{*}$
defined by (3.4) converges in $D(\mathbb{R})$ in distribution conditionally to
the process $V$, defined by

$V(t)=\sigma_{0}\sqrt{g(t_{0})}\,W(t),\,t\in\mathbb{R}.$ (3.9)

Here $W$ is standard two-sided Brownian motion, originating from zero.

###### Proof.

We consider the case $t\geq 0$. It is clear that
$t\mapsto\widetilde{W}_{n}(t)$ is a martingale with respect to the family of
$\sigma$-algebras ${\cal F}^{*}_{n,t},\,t\geq 0$, defined by:

$\displaystyle{\cal
F}^{*}_{n,t}=\sigma\left\\{(j/J,\bar{y}_{j}):j/J\in(t_{0},t_{0}+n^{-1/3}t]\right\\}.\qquad
t\geq 0.$

The quadratic variation process is, for $t\geq 0$ given by:

$\displaystyle\left[\widetilde{W}_{n}^{*}\right](t)=n^{-2/3}\sum_{j:j/J\in(t_{0},t_{0}+n^{-1/3}t]}N_{j}^{2}\left\\{\theta_{j}^{*}-\bar{y}_{j}\right\\}^{2}\sim
n^{-2/3}\sum_{j:j/J\in(t_{0},t_{0}+n^{-1/3}t]}N_{j}\sigma_{0}^{2}\,.$

If, for example as in ([5]), $J\sim n^{1/3}\log n$, we get:

$\displaystyle N_{j}\sim n^{2/3}g(t_{0})/\log n,$

and

$\displaystyle
n^{-2/3}\sum_{j:j/J\in(t_{0},t_{0}+n^{-1/3}t]}N_{j}\sigma_{0}^{2}\sim(\log
n)^{-1}\sum_{j:j/J\in(t_{0},t_{0}+n^{-1/3}t]}g(t_{0})\sigma_{0}^{2}\sim\sigma_{0}^{2}g(t_{0})\,t.$

The case $t<0$ is treated similarly. The result now follows from Rebolledo’s
theorem, see Theorem 3.6, p. 68 of [8] and [15]. ∎

######  Lemma 3.2.

Let $W$, $t_{0}$ and $V$ be as defined in Lemma 3.1 and $\widetilde{W}_{n}$ by
(3.5). Then the process $\widetilde{W}_{n}$ converges in $D(\mathbb{R})$ in
distribution, conditionally on the sequence $X_{1},X_{2},\dots$, to the
process $V$.

###### Proof of Lemma 3.2.

This is proved in the same way as Lemma 3.1, using that (3.5) is a martingale.
∎

###### Proof of Theorem 3.1.

We use the “switch relation” (see, e.g., Section 3.8 of [8] and section 5.1 of
[10]; the terminology is due to Iain Johnstone to denote a construction
introduced in a course given by the first author in Stanford, 1990). The
bootstrap estimate $f_{n}^{*}$ is computed as left derivative of the greatest
convex minorant of cumulative sum diagram (3.3). Let the processes $G_{n}$ and
$V_{n}^{*}$ be defined by

$\displaystyle G_{n}(t)=\sum_{j/J\leq t}N_{j}/n,\qquad
V_{n}^{*}(t)=\sum_{j/J\leq t}N_{j}Y_{j}^{*}/n,\qquad t\in[0,1].$

Moreover, let $U_{n}^{*}$ be defined by

$\displaystyle
U_{n}^{*}(a)=\text{argmin}\\{t\in[0,1]:V_{n}^{*}(t)-aG_{n}(t)\\},$

for $a$ in the range of $f_{0}$. Then we have the “switch relation”:

$\displaystyle\hat{f}_{n}^{*}(t)\geq a\iff G_{n}(t)\geq
G_{n}(U_{n}^{*}(a))\iff t\geq U_{n}^{*}(a),$

(compare with (3.35), p. 69 of [8]). So we get if $a=f_{0}(t_{0})$,

$\displaystyle{\mathbb{P}}\left\\{n^{1/3}\\{\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\\}\geq
x|D_{n}\right\\}={\mathbb{P}}\left\\{\hat{f}_{n}^{*}(t_{0})\geq
a+n^{1/3}x|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{U_{n}^{*}(a+n^{-1/3}x)\leq
t_{0}|D_{n}\right\\}={\mathbb{P}}\left\\{n^{1/3}\bigl{\\{}U_{n}^{*}(a+n^{-1/3}x)-t_{0}\bigr{\\}}\leq
0|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\text{argmin}\left[t\in[0,1]:V_{n}^{*}(t)-(a+n^{-1/3}x)\,G_{n}(t)\right]\leq
0|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\text{argmin}\left[t\in[0,1]:V_{n}^{*}(t)-V_{n}^{*}(t_{0})-(a+n^{-1/3}x)\\{G_{n}(t)-G_{n}(t_{0})\\}\right]\leq
0|D_{n}\right\\},$

where the last equality holds since the values of the argmin function do not
change if we add constants to the function for which we determine the argmin.

By Lemmas 3.1 and 3.2 we get the local expansion:

$\displaystyle{\mathbb{P}}\left\\{n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\bigr{\\}}\geq
x\Bigm{|}D_{n}\right\\}$
$\displaystyle\sim{\mathbb{P}}\left[\text{argmin}_{t}\left[\tilde{W}_{n}^{*}(t)+\tilde{W}_{n}(t)+\tfrac{1}{2}f_{0}^{\prime}(t_{0})g(t_{0})t^{2}-xg(t_{0})\right]\leq
0\bigm{|}D_{n}\right],$

which (using Brownian scaling) converges in distribution to

$\displaystyle{\mathbb{P}}\left[\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{\rm
argmin}_{t\in\mathbb{R}}\left[W_{1}(t)+W_{2}(t)+t^{2}\right]\leq-x\Bigm{|}W_{1}\right]$
$\displaystyle={\mathbb{P}}\left[\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{\rm
argmin}_{t\in\mathbb{R}}\left[W_{1}(t)+W_{2}(t)+t^{2}\right]\geq
x\Bigm{|}W_{1}\right].$

The last line uses the type of symmetry used in the proof of Theorem 5.2 of
[10]. ∎

In the proof of Theorem 3.1 we use the tightness of
$n^{1/3}\\{U_{n}^{*}(a+n^{-1/3}x)-t_{0}\\}$, which can be proved along
entirely similar lines as the proof of Lemma 3.5 in [8].

###  3.2 Convergence of a classical percentile bootstrap

It is of interest to investigate what happens if we perform a classical
empirical bootstrap, where we resample with replacement from the pairs
$(X_{i},Y_{i})$. This situation, where we also treat the $X_{i}$ as random
from the start instead of keeping them fixed, is called the “correlation
model” in [12]. In this case we compute the local means

$\displaystyle\bar{Y}^{*}_{j}=(N_{j}^{*})^{-1}\sum_{X_{i}^{*}\in
I_{j}}Y_{i}^{*},\qquad I_{j}=((j-1)/J,j/J],\qquad
N_{j}^{*}=\\#\\{i:X_{i}^{*}\in I_{j}\\},$ (3.10)

where the $(X_{i}^{*},Y_{i}^{*})$ are (discretely) uniformly (re-)sampled with
replacement from the set $D_{n}=\\{(X_{1},Y_{i}),\dots,(X_{n},Y_{n})\\}$. If
$N_{j}^{*}=0$ we define $\bar{Y}_{j}^{*}=0$, these values play no role in the
isotonization step.

Note that we can write alternatively, if $N_{j}^{*}>0$:

$\displaystyle\bar{Y}_{j}^{*}=(N_{j}^{*})^{-1}\sum_{X_{i}\in
I_{j}}M_{in}^{*}Y_{i},\qquad N_{j}^{*}=\sum_{i:X_{i}\in I_{j}}M_{in}^{*},$
(3.11)

where

$\displaystyle\left(M_{1n}^{*},\dots,M_{nn}^{*}\right)\sim\text{Multinomial}\left(n;n^{-1},\dots,n^{-1}\right).$

This means that

$\displaystyle{\mathbb{E}}\left\\{N_{j}^{*}\bar{Y}_{j}^{*}\bigm{|}(X_{1},Y_{1}),\dots,(X_{1},Y_{n})\right\\}=N_{j}\bar{Y}_{j},\qquad
j=1,\dots,J.$

The points of the cusum diagram needed to compute the bootstrap realization of
the LSE are given by

$\displaystyle
P_{j}^{*}=\left(\sum_{i=1}^{j}N_{i}^{*},\sum_{i=1}^{j}N_{i}^{*}\bar{Y}_{i}^{*}\right),\qquad
j=1,\dots,J.$

In order to study the local asymptotics of the greatest convex minorant of
this diagram, we consider the process

$\displaystyle\widetilde{W}_{n}^{*}(t)$
$\displaystyle=n^{-1/3}\Biggl{\\{}\sum_{j:j/J\in[0,t_{0}+n^{-1/3}t]}\sum_{X_{i}\in
I_{j}}(M_{in}^{*}-1)(Y_{i}-a_{0})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\sum_{j:j/J\in[0,t_{0}]}\sum_{X_{i}\in
I_{j}}(M_{in}^{*}-1)(Y_{i}-a_{0})\Biggr{\\}},$ (3.12)

where $a_{0}=f_{0}(t_{0})$, and

$\displaystyle\widetilde{W}_{n}(t)$
$\displaystyle=n^{-1/3}\Biggl{\\{}\sum_{j:j/J\in[0,t_{0}+n^{-1/3}t]}\sum_{X_{i}\in
I_{j}}(Y_{i}-a_{0})-\sum_{j:j/J\in[0,t_{0}]}\sum_{X_{i}\in
I_{j}}(Y_{i}-a_{0})\Biggr{\\}}.$

Defining

$\displaystyle U_{n}(a)=\text{argmin}\Bigl{\\{}t\in[0,1]:n^{-1}\sum_{j:j/J\leq
t}\bigl{\\{}N_{j}^{*}\bar{Y}_{j}^{*}-a\,N_{j}^{*}\bigr{\\}}\Bigr{\\}},$

results analogous to Lemmas 3.1 and 3.2 hold. For example we get Lemma 3.3
(the analogue to Lemma 3.1) which is proved in the Appendix.

######  Lemma 3.3.

Let $W$, $t_{0}$ and $V$ be as defined in Lemma 3.1 and
$\widetilde{W}_{n}^{*}$ be defined by (3.2). Then, along almost all sequences
$(X_{1},Y_{1}),(X_{2},Y_{2}),\dots$, the process $\widetilde{W}_{n}^{*}$
converges in $D(\mathbb{R})$ in distribution conditionally to the process $V$

So we get the same behavior as in subsection 3.1, but the present approach has
the interesting feature that we do not have to estimate the variance of the
errors separately. We can just resample with replacement from the original
sample $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$ and compute the estimator
$\hat{f}_{n}^{*}$ in the bootstrap samples.

The simulations, based on the regression function $f(x)=x^{2}+x/5$ with normal
noise with expectation $0$ and variance $0.01$ show almost no difference
beween the three methods if $n=20,000$, see Figure 4. At smaller sample size,
like, e.g., $n=1000$, the overcoverage is still not reached, as can be seen in
Figure 5. So the phenomenon of overcoverage also here only seems to occur with
very large sample sizes.

(a)

(b)

(c)

Figure 4: Coverage percentages for credible intervals and percentile
confidence intervals, for $n=20,000$. (a) credible intervals, (b) percentile
confidence intervals of section 3.1, (c) percentile confidence intervals of
section 3.2. The red line is at level 0.96324, but the intervals are based on
the $0.025$ and $0.975$ quantiles of the credible, respectively percentile
bootstrap simulations. The level 0.96324 was determined from the values of the
function $A$ in [5].

(a)

(b)

(c)

Figure 5: Coverage percentages for credible intervals and percentile
confidence intervals, for $n=1000$. (a) credible intervals, (b) percentile
confidence intervals of section 3.1, (c) percentile confidence intervals of
section 3.2. The red line is at level 0.96324, but the intervals are based on
the $0.025$ and $0.975$ quantiles of the credible, respectively percentile
bootstrap simulations. The level 0.96324 was determined from the values of the
function $A$ in [5].

##  4 Cube root $n$ consistent smoothed bootstrap confidence intervals

In [18] it was shown for interval censoring models that cube root $n$
convergent bootstrap confidence intervals can be computed for the distribution
function at a fixed point with the right asymptotic coverage. Key in this, is
the convergence of the nonparametric maximum likelihood estimator to
Chernoff’s distribution. We show that a similar approach is possible in the
present context.

In the regression context this means that we use, as in [9], the smoothed
least squares estimator (the SLSE) $\tilde{f}_{nh}$, for a bandwidth
$h=cn^{-1/5}$. To define $\tilde{f}_{nh}$, let $K$ be a symmetric twice
continuously differentiable nonnegative kernel with support $[-1,1]$ such that
$\int K(u)\,du=1$. Let $h>0$ be a bandwidth and define the scaled kernel
$K_{h}$ by

$\displaystyle
K_{h}(u)=\frac{1}{h}K\left(\frac{u}{h}\right),\,\,\,u\in\mathbb{R}.$ (4.1)

The SLSE $\tilde{f}_{nh}$ is then for $t\in[h,1-h]$ defined by

$\displaystyle\tilde{f}_{nh}(t)=\int K_{h}(t-x)\,\hat{f}_{n}(x)\,dx.$ (4.2)

For $t\notin[h,1-h]$ we use the boundary correction, defined in [9] (see (2.6)
and (2.7) in [9]).

We now generate residuals $E_{i}$ with respect to (4.2), defined by

$\displaystyle E_{i}=Y_{i}-\tilde{f}_{nh}(X_{i}),\qquad i=1,\dots,n,$

and compute the centered residuals $\tilde{E}_{i}$,

$\displaystyle\tilde{E}_{i}=E_{i}-n^{-1}\sum_{j=1}^{n}E_{j},\qquad
i=1,\dots,n.$ (4.3)

From these residuals, we generate bootstrap samples

$\displaystyle(X_{i},Y_{i}^{*}),\qquad
Y_{i}^{*}=\tilde{f}_{nh}(X_{i})+\tilde{E}_{i}^{*},\qquad i=1,\dots,n,$ (4.4)

where the $\tilde{E}_{i}^{*}$ are drawn uniformly with replacement from the
residuals $\tilde{E}_{i}$ defined by (4.3). For the bootstrap samples (4.4),
we compute the monotone (non-smoothed) LSE $\hat{f}_{n}^{*}$ and consider the
differences

$\displaystyle\hat{f}_{n}^{*}(t)-\tilde{f}_{nh}(t),$ (4.5)

and the $95\%$ bootstrap confidence intervals, given by

$\displaystyle\left(\hat{f}_{n}(t)-Q_{0.975}^{*},\hat{f}_{n}(t)-Q_{0.025}^{*}\right),$
(4.6)

where $Q_{0.025}^{*}$ and $Q_{0.975}^{*}$ are the $2.5$th and $97.5$th
percentiles of $1000$ bootstrap samples of (4.5) and $\hat{f}_{n}$ is the LSE
in the original sample. Note that this is the more conventional bootstrap
approach, rather than the percentile method.

(a)

(b)

(c)

Figure 6: Confidence intervals, based on (4.6), for $n=100,500$ and $1000$.
The blue curve is the nonparametric monotone LSE, and the red dashed curve the
real regression function $f_{0}(x)=x^{2}+x/5$. (a) $n=100$, (b) $n=500$, (c)
$n=1000$.

The SLSE $\tilde{f}_{nh}$. which plays a central role in the construction of
these confidence intervals, has the limit behavior specified in the following
theorem, which is Theorem 1 in [9].

######  Theorem 4.1.

Let $f_{0}$ be a nondecreasing continuous function on $[0,1]$. Let
$X_{1},X_{2},\ldots$ be i.i.d random variables with continuous density $g$,
staying away from zero on $[0,1]$, where the derivative $g^{\prime}$ is
continuous and bounded on $(0,1)$. Furthermore, let
$\varepsilon_{1},\varepsilon_{2}\ldots$ be i.i.d. random variables distributed
according to a sub-Gaussian distribution with expectation zero and variance
$0<\sigma_{0}^{2}<\infty$, independent of the $X_{i}$’s. Then consider
$Y_{i}$, defined by

$\displaystyle Y_{i}=f_{0}(X_{i})+\varepsilon_{i},\,\,\,\,i=1,2,\ldots$

Suppose $t_{0}\in(0,1)$ such that $f_{0}$ has a strictly positive derivative
and a continuous second derivative $f_{0}^{\prime\prime}(t_{0})\neq 0$ at
$t_{0}$. Then, for the SLSE $\tilde{f}_{nh}$ defined by (4.2) based on the
pairs $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, and $h\sim cn^{-1/5}$ for $c>0$,

$\displaystyle
n^{2/5}\left\\{\tilde{f}_{nh}(t_{0})-f_{0}(t_{0})\right\\}\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(\beta,\sigma^{2}).$

Here

$\displaystyle\beta=\tfrac{1}{2}c^{2}f_{0}^{\prime\prime}(t_{0})\int
u^{2}K(u)\,du\,\,\,\mbox{ and
}\,\,\,\sigma^{2}=\frac{\sigma^{2}_{0}}{cg(t)}\int K(u)^{2}\,du.$ (4.7)

The asymptotically Mean Squared Error optimal constant $c$ is given by:

$\displaystyle c=\left\\{\frac{\sigma^{2}_{0}}{g(t_{0})}\int
K(u)^{2}\,du\Bigm{/}\left\\{f_{0}^{\prime\prime}(t_{0})\int
u^{2}K(u)\,du\right\\}^{2}\right\\}^{1/5}.$

We have the following lemma of which the proof is given in the Appendix.

######  Lemma 4.1.

Let $W$, $t_{0}$ and $V$ be as defined in Lemma 3.1 and
$\widetilde{W}_{n}^{*}$ be defined by:

$\displaystyle\widetilde{W}_{n}^{*}(t)=n^{-1/3}\Biggl{\\{}\sum_{i:X_{i}\in[0,t_{0}+n^{-1/3}t]}\tilde{E}_{i}^{*}-\sum_{i:X_{i}\in[0,t_{0}]}\tilde{E}_{i}^{*}\Biggr{\\}},\qquad
t\in\mathbb{R},$

where $\tilde{E}_{i}^{*}$ is drawn with replacement from the residuals
$\tilde{E}_{i}$, defined by (4.3). Then, along almost all sequences
$(X_{1},Y_{1}),(X_{2},Y_{2}),\dots$, the process $\widetilde{W}_{n}^{*}$
converges in $D(\mathbb{R})$ in distribution conditionally to the process $V$

(a)

(b)

(c)

Figure 7: Coverage percentages for confidence intervals, based on (4.6), for
$n=100,500$ and $1000$. (a) $n=100$, (b) $n=500$, (c) $n=1000$.

This leads to the following corollary.

######  Corollary 4.1.

Let the bootstrap LSE $f_{n}^{*}$ be constructed as in (4.5) and let $t_{0}$
and the SLSE $\tilde{f}_{nh}$ be defined as in Theorem 4.1. Then, along almost
all sequences $(X_{1},Y_{1}),\dots$, we have, under the condtions of Theorem
4.1,

$\displaystyle
n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\bigr{\\}}\stackrel{{\scriptstyle{\mathscr{D}}}}{{\longrightarrow}}Z,$

where

$\displaystyle
Z=\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{\rm
argmin}_{t\in\mathbb{R}}\left[W(t)+t^{2}\right],$

and $W$ is standard two-sided Brownian motion, originating from zero.

###### Proof.

We use the “switch relation” again (see the proof of Theorem 3.1). Let $G_{n}$
be the empirical distribution function of the $X_{i}$ and let $V_{n}^{*}$ be
defined by

$\displaystyle V_{n}^{*}(t)=\sum_{X_{i}\leq t}Y_{i}^{*},$

where $Y_{i}^{*}$ is defined by (4.4). Let $U_{n}^{*}$ be defined by

$\displaystyle
U_{n}^{*}(a)=\text{argmin}\\{t\in[0,1]:V_{n}^{*}(t)-aG_{n}(t)\\},$

for $a$ in the range of $f_{0}$. Then we have the “switch relation”:

$\displaystyle\hat{f}_{n}^{*}(t)\geq a\iff G_{n}(t)\geq
G_{n}(U_{n}^{*}(a))\iff t\geq U_{n}^{*}(a).$

Hence we get if $a_{n}=\tilde{f}_{nh}(t_{0})$ and
$D_{n}=\\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\\}$,

$\displaystyle{\mathbb{P}}\left\\{n^{1/3}\\{\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\\}\geq
x|D_{n}\right\\}={\mathbb{P}}\left\\{\hat{f}_{n}^{*}(t_{0})\geq
a_{n}+n^{-1/3}x|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{U_{n}^{*}(a_{n}+n^{-1/3}x)\leq
t_{0}|D_{n}\right\\}={\mathbb{P}}\left\\{n^{1/3}\bigl{\\{}U_{n}^{*}(a_{n}+n^{-1/3}x)-t_{0}\bigr{\\}}\leq
0|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\text{argmin}\left[t\in[0,1]:V_{n}^{*}(t)-(a_{n}+n^{-1/3}x)\,G_{n}(t)\right]\leq
0|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\text{argmin}\left[t\in[0,1]:V_{n}^{*}(t)-V_{n}^{*}(t_{0})-(a_{n}+n^{-1/3}x)\\{G_{n}(t)-G_{n}(t_{0})\\}\right]\leq
0|D_{n}\right\\},$

where the last equality holds since the values of the argmin function do not
change if we add constants to the function for which we determine the argmin.

We have:

$\displaystyle{\mathbb{P}}\left\\{\text{argmin}\left[t\in[0,1]:V_{n}^{*}(t)-V_{n}^{*}(t_{0})-(a_{n}+n^{-1/3}x)\\{G_{n}(t)-G_{n}(t_{0})\\}\right]\leq
0|D_{n}\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\text{argmin}_{t}\left[W_{n}^{*}(t)+n^{1/3}\int_{(t_{0},t_{0}+n^{-1/3}t]}\tilde{f}_{nh}(u)\,dG_{n}(u)\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\left.\left.-n^{1/3}(\tilde{f}_{nh}(t_{0})+n^{-1/3}x)\\{G_{n}(t_{0}+n^{-1/3})-G_{n}(t_{0})\\}\right]\leq
0|D_{n}\right\\}$
$\displaystyle\sim{\mathbb{P}}\left\\{\text{argmin}_{t}\left[W_{n}^{*}(t)+\tfrac{1}{2}f_{nh}^{\prime}(t_{0})g_{0}(t_{0})t^{2}-xg(t_{0})t\right]\leq
0|D_{n}\right\\}$
$\displaystyle\sim{\mathbb{P}}\left\\{\text{argmin}_{t}\left[W_{n}^{*}(t)+\tfrac{1}{2}f_{0}^{\prime}(t_{0})g_{0}(t_{0})t^{2}-xg(t_{0})t\right]\leq
0|D_{n}\right\\}$
$\displaystyle\longrightarrow{\mathbb{P}}\left\\{\text{argmin}_{t}\left[\sigma_{0}\sqrt{g(t_{0})}\,W(t)+\tfrac{1}{2}f_{0}^{\prime}(t_{0})g(t_{0})t^{2}-xg(t_{0})t\right]\leq
0\right\\},$

where we use Lemma 4.1 in the last step. By Brownian scaling, we can write

$\displaystyle{\mathbb{P}}\left\\{\text{argmin}_{t}\left[\sigma_{0}\sqrt{g(t_{0})}\,W(t)+\tfrac{1}{2}f_{0}^{\prime}(t_{0})g_{0}(t_{0})t^{2}-xg(t_{0})t\right]\leq
0\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{argmin}_{t}\left[W(t)+t^{2}\right]\leq-x\right\\}$
$\displaystyle={\mathbb{P}}\left\\{\left(\frac{4\sigma_{0}^{2}f_{0}^{\prime}(t_{0})}{g(t_{0})}\right)^{1/3}\text{argmin}_{t}\left[W(t)+t^{2}\right]\geq
x\right\\}.$

∎

######  Remark 4.1.

Corollary 4.1 shows that the limit distribution of
$n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\bigr{\\}}$,
along almost all sequences $D_{n}$, is the same as the limit distributionof
$n^{1/3}\bigl{\\{}\hat{f}_{n}(t_{0})-\tilde{f}_{0}(t_{0})\bigr{\\}}$. In the
latter case the limit distribution was first derived in [4].

Using this corollary it is clear that in using this method the (ordinary, not
percentile) smoothed bootstrap simulations recreate the actual asymptotic
distribution correctly, and that we do not have to use a correction for over-
or undercoverage. It is also clear from Figure 7 that its behavior is much
better than the behavior of the confidence intervals in the preceding section.
Even for sample size $n=100$ the confidence intervals are more or less “on
target”. In comparison, the credible intervals are still far off the target
for these sample sizes, and the overcoverage has still not set in for these
sample sizes, see Figure 8.

(a)

(b)

(c)

Figure 8: Coverage percentages for the confidence intervals, for (a) $n=100$,
(b) $n=500$, (c) $n=1000$.

We can also prove the result below, illustrating the fact that there is no
need for correction for over-or under coverage.

######  Theorem 4.2.

Let $D_{n}=\\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\\}$ and $h\asymp n^{-1/5}$,
and let $t_{0}$ and $\tilde{f}_{nh}$ be defined as in Corollary 4.1. Let the
conditions of Theorem 4.1 be satisfied. Then, for $z\in(0,1)$, as
$n\to\infty$,

$\displaystyle{\mathbb{P}}\left\\{{\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
0\Bigm{|}D_{n}\right)\leq z\right\\}\longrightarrow z.$ (4.8)

###### Proof.

Let

$\displaystyle
U_{n}=n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\bigr{\\}},\qquad
V_{n}=n^{1/3}\bigl{\\{}\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}},$

and let $\Phi$ be the distribution function of $Z$, defined in Corollary 4.1.
We have

$\displaystyle{\mathbb{P}}\left\\{{\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
0\Bigm{|}D_{n}\right)\leq z\right\\}$
$\displaystyle={\mathbb{P}}\left\\{{\mathbb{P}}\\{V_{n}\leq-U_{n}|D_{n}\\}\leq
z\right\\}$ $\displaystyle\sim{\mathbb{P}}\left\\{\Phi(-U_{n})\\}\leq
z\right\\}={\mathbb{P}}\left\\{-U_{n}\leq\Phi^{-1}(z)\right\\}\longrightarrow\Phi\circ\Phi^{-1}(z)=z,$

using the symmetry around zero of the distribution of $Z$ (the limit in
distribution of $U_{n}$). ∎

Figure 9: Histogram of $1000$ estimates of
${\mathbb{P}}\left(n^{1/3}\bigl{\\{}\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\bigr{\\}}\leq
0\bigm{|}D_{n}\right)$ for the smoothed bootstrap, for $1000$ samples $D_{n}$
of size $n=2000$ and $t_{0}=0.5$.

Note that $\hat{f}_{n}^{*}(t_{0})$ is now centered by $\tilde{f}_{nh}(t_{0})$
instead of $\hat{f}_{n}$, and that
$n^{1/3}\\{\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{n,h}(t_{0})\\}$ tends to the
right limit distribution, in contrast with
$n^{1/3}\\{\hat{f}_{n}^{*}(t_{0})-\hat{f}_{n}(t_{0})\\}$. The histogram of
estimates of
${\mathbb{P}}\bigl{(}n^{1/3}\\{\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})+\hat{f}_{n}(t_{0})-f_{0}(t_{0})\\}\leq
0\bigm{|}D_{n}\bigr{)}$, based on relative frequencies, is shown in Figure 9.

For the ordinary (non-percentile) bootstrap we get by an entirely similar
proof, in which we do not need the symmetry of the limit distribution:

######  Theorem 4.3.

Let $D_{n}=\\{((X_{1},Y_{1}),\dots,(X_{n},Y_{n})\\}$ and $h\asymp n^{-1/5}$.
Let the conditions of Theorem 4.1 be satisfied. Then, for $z\in(0,1)$, as
$n\to\infty$,

$\displaystyle{\mathbb{P}}\left\\{{\mathbb{P}}\left(\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\leq\hat{f}_{n}(t_{0})-f_{0}(t_{0})\Bigm{|}D_{n}\right)\leq
z\right\\}\longrightarrow z.$ (4.9)

The result gives an interesting consequence of what it means to say that the
”bootstrap works”. This phenomenon also occurs in the simple bootstrap setting
where one resamples with replacement from samples $U_{1},\ldots,U_{n}$ from a
normal $N(\mu,\sigma^{2})$ distribution with the aim to construct a confidence
set for the mean. Then also, for all $z\in(0,1)$,

$\displaystyle{\mathbb{P}}\left\\{{\mathbb{P}}\left(\bar{U}_{n}^{*}-\bar{U}_{n}\leq\bar{U}_{n}-\mu\Bigm{|}U_{1},\ldots
U_{n}\right)\leq z\right\\}\longrightarrow z.$ (4.10)

##  5 Concluding remarks

We showed that the construction of pointwise credible intervals for the
monotone regression function, as proposed in [5], has an interpretation as the
construction of percentile bootstrap intervals. The overcoverage, as explained
by Theorem 3.1, only sets in for very large sample sizes, like $n=20,000$; for
smaller sample sizes we have observed undercoverage.

Because the confidence intervals are based on piecewise constant estimates of
the regression function, on intervals of equal length, the effect of bias is
very pronounced, which does not hold for the confidence intervals, based on
the smoothed bootstrap in Section 4. The latter confidence intervals have the
further advantages of being on target also for smaller sample sizes, and not
needing correction for overcoverage or undercoverage, since the estimates are
consistent.

The consistency is also borne out by Theorem 4.3, showing convergence to the
uniform distribution of
${\mathbb{P}}(\hat{f}_{n}^{*}(t_{0})-\tilde{f}_{nh}(t_{0})\leq\hat{f}_{n}(t_{0})-f_{0}(t_{0})|D_{n})$,
$D_{n}=\\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\\}$, for the smoothed bootstrap
estimates $\hat{f}_{n}^{*}$, in contrast with the situation for the credible
intervals and the percentile bootstrap intervals, where we need a correction
for convergence of ${\mathbb{P}}(\hat{f}_{n}^{*}(t_{0})-f_{0}(t_{0})\leq
0\bigm{|}D_{n})$ to a distribution different from the uniform distribution.

As shown in [9], it is also possible to use the smoothed least squares
estimator (SLSE) directly as the basis for confidence intervals. In this case,
resampling is done from residuals w.r.t. an oversmoothed estimate of the
regression function to treat the bias in the right way. The bias is in this
case much more of a problem because the variance and squared bias of the SLSE
are of the same order if the bandwidth is of order $n^{-1/5}$. A picture of
confidence intervals of this type is given in Figure 10. For more details, see
[9].

All simulations in our paper can be recreated using the ${\tt R}$ scripts in
[7].

Figure 10: The SLSE (blue, solid) and $95\%$ confidence intervals, using the
theory in [9], for sample size $n=100$ and $f_{0}(x)=x^{2}+x/5$ (red dashed
curve).

## Appendix

Derivation of raw posterior distribution (2.2)
Considering $\bm{X}$ fixed, we have $\bm{Y}=\bm{B}\bm{\theta}+\bm{\epsilon}$,
where $\bm{B}=(b_{ij})$ with $b_{ij}=1_{I_{j}}(X_{i})$ giving
$\bm{Y}|\bm{\theta}\sim N_{n}(\bm{B}\bm{\theta},\sigma_{0}^{2}I_{n})$. Writing
$\bm{\Lambda}={\rm Diag}(\lambda_{j}^{2})$, this can be combined with the
prior $\bm{\theta}\sim N_{J}(\bm{\zeta},\sigma_{0}^{2}\Lambda)$ yielding

$\displaystyle f(\bm{\theta}|\bm{Y})\propto
f(\bm{Y}|\bm{\theta})f(\bm{\theta})\propto\exp\left(-\frac{1}{2\sigma_{0}^{2}}\left[\bm{\theta}^{T}(\bf\Lambda^{-1}+\bm{B}^{T}\bm{B})\bm{\theta}-2\bm{\theta}^{T}(\bm{B}^{T}\bm{Y}+\bm{\Lambda}^{-1}\bm{\zeta})\right]\right)$

By ‘completing the square’, this function can be seen to be proportional to
the normal density with covariance matrix

$\sigma_{0}^{2}\left(\bm{\Lambda}^{-1}+\bm{B}^{T}\bm{B}\right)^{-1}=\sigma_{0}^{2}{\rm
Diag}(1/\lambda_{j}^{2}+n_{j})^{-1}={\rm
Diag}\left(\frac{\sigma_{0}^{2}}{1/\lambda_{j}^{2}+n_{j}}\right)$

where we use that both $\bm{\Lambda}$ and $\bm{B}^{T}\bm{B}$ are diagonal
matrices with diagonal entries $\lambda_{j}^{2}$ and
$n_{j}=\sum_{i}b_{ij}=\\#\\{i\,:x_{i}\in I_{j}\\}$ respectively. The
expectation is given by

$\left(\bm{\Lambda}^{-1}+\bm{B}^{T}\bm{B}\right)^{-1}(\bm{B}^{T}\bm{Y}+\bm{\Lambda}^{-1}\bm{\zeta})={\rm
Diag}\left(\frac{1}{1/\lambda_{j}^{2}+n_{j}}\right)\left(\begin{array}[]{c}\zeta_{1}/\lambda_{1}^{2}+\sum_{\\{i\,:\,x_{i}\in
I_{1}\\}}y_{i}\\\ \vdots\\\ \zeta_{J}/\lambda_{J}^{2}+\sum_{\\{i\,:\,x_{i}\in
I_{J}\\}}y_{i}\end{array}\right),$

boiling down to the expression in (2.2).

Derivation of (2.5)
First note that $I+B\Lambda B^{T}$ is a block diagonal matrix, where block
$j$, has size $n_{j}\times n_{j}$, diagonal elements $1+\lambda_{j}^{2}$ and
off-diagonal elements $\lambda_{j}^{2}$ for $1\leq j\leq J$. The $j$-th block
can be written as

$A_{j}=I_{n_{j}\times
n_{j}}+\lambda_{j}^{2}{\mathbf{1}}_{n_{j}}{\mathbf{1}}_{n_{j}}^{T},$

where $I_{k\times k}$ is the $k\times k$ identity matrix and
${\mathbf{1}}_{k}$ is the column vector of length $k$ with all elements equal
to one.

This means that $(I+B\Lambda B^{T})^{-1}$ is also a block diagonal matrix,
with $j$-th block

$A_{j}^{-1}=I_{n_{j}\times
n_{j}}-\frac{\lambda_{j}^{2}}{1+n_{j}\lambda_{j}^{2}}{\mathbf{1}}_{n_{j}}{\mathbf{1}}_{n_{j}}^{T},$

by the Sherman Morrison formula. For convenience, write $z=y-B\zeta$ and
denote by subscript $[j]$ the part of a vector corresponding to the $j$-th
block (so $i$ for which $x_{i}\in I_{j}$; $z_{[j]}$ has length $n_{j}$). Then

$\left[(I+B\Lambda
B^{T})^{-1}z\right]_{[j]}=A_{j}^{-1}z_{[j]}=z_{[j]}-\frac{\lambda_{j}^{2}}{1+n_{j}\lambda_{j}^{2}}{\mathbf{1}}_{n_{j}}{\mathbf{1}}_{n_{j}}^{T}z_{[j]}=z_{[j]}-\frac{\lambda_{j}^{2}n_{j}\bar{z}_{[j]}}{1+n_{j}\lambda_{j}^{2}}{\mathbf{1}}_{n_{j}},$

where $\bar{z}_{[j]}$ denotes the average of the entries of $z_{[j]}$.
Therefore

$z_{[j]}^{T}\left[(I+B\Lambda
B^{T})^{-1}z\right]_{[j]}=z_{[j]}^{T}z_{[j]}-\frac{\lambda_{j}^{2}n_{j}^{2}\bar{z}_{[j]}^{2}}{1+n_{j}\lambda_{j}^{2}}.$

Hence,

$\displaystyle z^{T}(I+B\Lambda
B^{T})^{-1}z=\sum_{j=1}^{J}z_{[j]}^{T}\left[(I+B\Lambda
B^{T})^{-1}z\right]_{[j]}=z^{T}z-\sum_{j=1}^{J}\frac{n_{j}\bar{z}_{[j]}^{2}}{1+1/(n_{j}\lambda_{j}^{2})}.$
(5.1)

Now write $(1+1/(n_{j}\lambda_{j}^{2}))^{-1}=1-\delta_{j}$, so
$\delta_{j}=(1+n_{j}\lambda_{j}^{2})^{-1}$. Then (5.1) can be further
rewritten as

$\displaystyle\,\,\,z^{T}(I+B\Lambda
B^{T})^{-1}z=z^{T}z-\sum_{j=1}^{J}(1-\delta_{j})n_{j}\bar{z}_{[j]}^{2}=z^{T}z-\sum_{j=1}^{J}n_{j}\bar{z}_{[j]}^{2}+\sum_{j=1}^{J}\delta_{j}n_{j}\bar{z}_{[j]}^{2}$
(5.2) $\displaystyle\,\,\,=\sum_{j=1}^{J}\sum_{x_{i}\in
I_{j}}\left((y_{i}-\zeta_{j})^{2}-(\bar{y}_{j}-\zeta_{j})^{2}\right)+\sum_{j=1}^{J}\delta_{j}n_{j}(\bar{y}_{j}-\zeta_{j})^{2}$
$\displaystyle\,\,\,=\sum_{j=1}^{J}\sum_{x_{i}\in
I_{j}}\left(y_{i}^{2}-\bar{y}_{j}^{2}\right)+\sum_{j=1}^{J}\delta_{j}n_{j}(\bar{y}_{j}-\zeta_{j})^{2}=\sum_{j=1}^{J}\sum_{x_{i}\in
I_{j}}\left(y_{i}-\bar{y}_{j}\right)^{2}+\sum_{j=1}^{J}\delta_{j}n_{j}(\bar{y}_{j}-\zeta_{j})^{2}.$

Substituting $\delta_{j}=(1+n_{j}\lambda_{j}^{2})^{-1}$ and writing the first
term as one sum, yields (2.5).

Derivation of Empirical Bayes estimators (2.6) and (2.4) The distribution of
the observed $\bm{Y}$ can be expressed in terms of the parameters
$\sigma_{0}^{2}$, $\Lambda$ and $\zeta$,

${\bm{Y}}={\bm{B}}{\bm{\theta}}+\epsilon={\bm{B}}({\bm{\zeta}}+\sigma_{0}\tilde{\epsilon})+\sigma_{0}\epsilon={\bm{B}}\zeta+\sigma_{0}({\bm{B}}\tilde{\epsilon}+\epsilon),$

where $\epsilon\sim N_{n}(0,I_{n\times n})$ and $\tilde{\epsilon}\sim
N_{J}(0,\Lambda)$ are independent. Therefore,

${\bm{Y}}\sim N_{n}\left({\bm{B}}\zeta,\sigma_{0}^{2}\left(I_{n\times
n}+{\bm{B}}{\bm{\Lambda}}{\bm{B}}^{T}\right)\right).$

Maximizing the likelihood in $\zeta$, for fixed values of $\sigma_{0}^{2}$ and
$\Lambda$ entails minimizing

$(y-B\zeta)^{T}\left(I_{n\times
n}+{\bm{B}}{\bm{\Lambda}}{\bm{B}}^{T}\right)^{-1}(y-B\zeta).$

Recognizing (5.2) in this expression, it is clear that the empirical Bayes
estimate of $\zeta$ is either given by the vector
$(\bar{y}_{1},\ldots,\bar{y}_{J})^{T}$ if the likelihood is maximized over
$\mathbb{R}^{J}$ or its isotonic regression with weights
$n_{j}\delta_{j}=n_{j}/(1+n_{j}\lambda_{j}^{2})$ if monotonicity is taken into
account.

For any fixed value of $\zeta$, maximizing the the log likelihood of
$\sigma_{0}$ corresponds to minimizing

$\frac{n}{2}\log\sigma_{0}^{2}+(y-B\zeta)^{T}\left(I_{n\times
n}+{\bm{B}}{\bm{\Lambda}}{\bm{B}}^{T}\right)^{-1}(y-B\zeta)/(2\sigma_{0}^{2}),$

yielding (2.4).

###### Proof of Lemma 3.3.

We employ a construction, used in the proof of Lemma 2.2 in [6]. Let $A_{n}$
be the interval $[t_{0}-n^{-1/3}\log n,t_{0}+n^{-1/3}\log n]$ and let
$(U_{1}^{*},V_{1}^{*}),(U_{2}^{*},V_{2}^{*}),\dots$ be an i.i.d sequence of
points, (discretely) uniformly distributed on the set of points
$(X_{i},Y_{i})$ such that $X_{i}\in A_{n}$.

Let $M_{n}$ be the number of points $X_{i}\in A_{n}$. The number of bootstrap
draws such that the first component belongs to $A_{n}$ has distribution

$\displaystyle M_{n}^{*}\sim\text{Binom}(n,M_{n}/n),$ (5.3)

so, taking the random variable $M_{n}^{*}$ defined by (5.3), independent of
the sequence $(U_{1}^{*},V_{1}^{*})$, $(U_{2}^{*},V_{2}^{*}),\dots$, we can
represent the bootstrap variables such that the first component belongs to
$A_{n}$ by

$\displaystyle\sum_{i=1}^{M_{n}^{*}}\delta_{\\{(U_{i}^{*},V_{i}^{*})\\}},$

where $\delta_{x}$ denotes Dirac measure.

We can couple this process with a Poisson process

$\displaystyle\sum_{i=1}^{N_{n}}\delta_{\\{(U_{i}^{*},V_{i}^{*})\\}},$

where

$\displaystyle N_{n}\sim\text{Poisson}\,(M_{n}),$

independent of the $(U_{i}^{*},V_{i}^{*})$, using the construction with a
Uniform(0,1) random variable $U$ as in the construction in the proof of Lemma
2.2 in [6]. We find in this way:

$\displaystyle{\mathbb{P}}\left\\{\sum_{i=1}^{M_{n}^{*}}\delta_{\\{(U_{i}^{*},V_{i}^{*})\\}}\neq\sum_{i=1}^{N_{n}}\delta_{\\{(U_{i}^{*},V_{i}^{*})\\}}\right\\}\leq
2M_{n}/n,$

where $M_{n}/n$ tends to zero almost surely, using an inequality from [20].

This means that we can replace $M_{in}^{*}$ by $N_{in}$ in (3.2), where the
$N_{in}$ are independent Poisson$(1)$ random variables and can replace
$\widetilde{W}_{n}^{*}$ by its Poissonized version

$\displaystyle W_{n}^{(P)}(t)=$ $\displaystyle
n^{-1/3}\sum_{j:j/J\in[0,t_{0}+n^{-1/3}t]}\sum_{X_{i}\in
I_{j}}(N_{in}-1)(Y_{i}-a_{0}-n^{-1/3}x)$ (5.4)
$\displaystyle\qquad\qquad\qquad\qquad-n^{-1/3}\sum_{j:j/J\in[0,t_{0}]}\sum_{X_{i}\in
I_{j}}(N_{in}-1)(Y_{i}-a_{0}-n^{-1/3}x).$ (5.5)

For the latter process we have the martingale structure again, and the
quadratic variation process $[W_{n}^{(P)}](t),\,t\geq 0$ satisfies

$\displaystyle\left[W_{n}^{(P)}\right](t)=n^{-2/3}\sum_{j:j/J\in(t_{0},t_{0}+n^{-1/3}t]}\Bigl{\\{}\sum_{X_{i}\in
I_{j}}(N_{in}-1)(Y_{i}-a_{0}-n^{-1/3}x)\Bigr{\\}}^{2}$
$\displaystyle\longrightarrow\sigma_{0}^{2}g(t_{0})t.$

for almost all sequences $(X_{1},Y_{1}),\dots$. A similar reation holds for
$t<0$. So the result follows in the same way as in the proof of Lemma 3.1. ∎

###### Proof of Lemma 4.1.

We consider the case $t\geq 0$. It is clear that, conditionally on
$(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$, $t\mapsto\widetilde{W}_{n}^{*}(t)$ is a
martingale with respect to the family of $\sigma$-algebras ${\cal
F}^{*}_{n,t},\,t\geq 0$, defined by:

$\displaystyle{\cal
F}^{*}_{n,t}=\sigma\left\\{(X_{i},E_{i}^{*}):X_{i}\in[t_{0}+n^{-1/3}t]\right\\}.\qquad
t\geq 0.$

The quadratic variation process is given by:

$\displaystyle\left[\widetilde{W}_{n}^{*}\right](t)=n^{-2/3}\sum_{i:X_{i}\in[t_{0},t_{0}+n^{-1/3}t]}\left(E_{i}^{*}\right)^{2}$

We have:

$\displaystyle{\mathbb{E}}\left\\{\left[\widetilde{W}_{n}^{*}\right](t)\Bigm{|}(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\right\\}$
$\displaystyle=n^{-2/3}\sum_{i:X_{i}\in[t_{0},t_{0}+n^{-1/3}t]}{\mathbb{E}}\left\\{\left(E_{i}^{*}\right)^{2}\Bigm{|}(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\right\\}$
$\displaystyle=n^{-2/3}\sum_{i:X_{i}\in[t_{0},t_{0}+n^{-1/3}t]}n^{-1}\sum_{j=1}^{n}\tilde{E}_{j}^{2}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}\sigma_{0}^{2}g(t_{0})t.$

Note that

$\displaystyle\bigl{\\{}Y_{j}-\tilde{f}_{nh}(X_{j})\bigr{\\}}^{2}$
$\displaystyle=\bigl{\\{}Y_{j}-f_{0}(X_{j})\bigr{\\}}^{2}+\bigl{\\{}f_{0}(X_{j})-\tilde{f}_{nh}(X_{j})\bigr{\\}}^{2}+2\bigl{\\{}Y_{j}-f_{0}(X_{j})\bigr{\\}}\bigl{\\{}f_{0}(X_{j})-\tilde{f}_{nh}(X_{j})\bigr{\\}},$

and that therefore

$\displaystyle n^{-1}\sum_{j=1}^{n}\tilde{E}_{j}^{2}\sim
n^{-1}\sum_{j=1}^{n}\bigl{\\{}Y_{j}-f_{0}(X_{j})\bigr{\\}}^{2}$

almost surely, as $n\to\infty$, by the properties of $\tilde{f}_{nh}$.

We can treat the case $t<0$ in a similar way. This means that we can apply
Rebolledo’s theorem, see Theorem 3.6, p. 68 of [8] and [15]. The conclusion of
the lemma now follows. ∎

## References

  * Barlow et al. [1972] [author] Barlow, R. E.R. E., Bartholomew, D. J.D. J., Bremner, J. M.J. M. and Brunk, H. D.H. D. (1972). Statistical inference under order restrictions. The theory and application of isotonic regression. John Wiley & Sons, London-New York-Sydney Wiley Series in Probability and Mathematical Statistics. 0326887 (48 ##5229)
  * Bhaumik and Ghosal [2015] [author] Bhaumik, PrithwishP. and Ghosal, SubhashisS. (2015). Bayesian two-step estimation in differential equation models.
  * Bhaumik and Ghosal [2017] [author] Bhaumik, PrithwishP. and Ghosal, SubhashisS. (2017). Efficient Bayesian estimation and uncertainty quantification in ordinary differential equation models.
  * Brunk [1970] [author] Brunk, H. D.H. D. (1970). Estimation of isotonic regression. In Nonparametric Techniques in Statistical Inference (Proc. Sympos., Indiana Univ., Bloomington, Ind., 1969) 177–197. Cambridge Univ. Press, London. 0277070
  * Chakraborty and Ghosal [2021] [author] Chakraborty, MoumitaM. and Ghosal, SubhashisS. (2021). Coverage of credible intervals in nonparametric monotone regression. Ann. Statist. 49 1011–1028. 10.1214/20-aos1989 4255117
  * Groeneboom [1988] [author] Groeneboom, PietP. (1988). Limit theorems for convex hulls. Probab. Theory Related Fields 79 327–368. 10.1007/BF00342231 959514
  * Groeneboom [2021] [author] Groeneboom, PietP. (2021). Monotone Regression. https://github.com/pietg/monotone-regression.
  * Groeneboom and Jongbloed [2014] [author] Groeneboom, PietP. and Jongbloed, GeurtG. (2014). Nonparametric Estimation under Shape Constraints. Cambridge Univ. Press, Cambridge.
  * Groeneboom and Jongbloed [2023] [author] Groeneboom, PietP. and Jongbloed, GeurtG. (2023). Confidence intervals in monotone regression. Submitted.
  * Groeneboom and Wellner [1992] [author] Groeneboom, P.P. and Wellner, J. A.J. A. (1992). Information bounds and nonparametric maximum likelihood estimation. DMV Seminar 19\. Birkhäuser Verlag, Basel. 1180321 (94k:62056)
  * Hall [1988] [author] Hall, PeterP. (1988). Theoretical comparison of bootstrap confidence intervals. Ann. Statist. 16 927–985. With a discussion and a reply by the author. 10.1214/aos/1176350933 959185
  * Hall [1992] [author] Hall, PeterP. (1992). The bootstrap and Edgeworth expansion. Springer Series in Statistics. Springer.
  * Kosorok [2008] [author] Kosorok, M. R.M. R. (2008). Bootstrapping the Grenander estimator. In Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen. Inst. Math. Stat. Collect. 1 282–292. Inst. Math. Statist., Beachwood, OH.
  * Lin and Dunson [2014] [author] Lin, LizhenL. and Dunson, David BD. B. (2014). Bayesian monotone regression using Gaussian process projection. Biometrika 101 303–317.
  * Rebolledo [1980] [author] Rebolledo, R.R. (1980). Central limit theorems for local martingales. Z. Wahrsch. Verw. Gebiete 51 269–286. 10.1007/BF00587353 566321 (81g:60023)
  * Robertson, Wright and Dykstra [1988] [author] Robertson, T.T., Wright, F. T.F. T. and Dykstra, R. L.R. L. (1988). Order restricted statistical inference. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons Ltd., Chichester. 961262 (90b:62001)
  * Sen, Banerjee and Woodroofe [2010] [author] Sen, B.B., Banerjee, M.M. and Woodroofe, M. B.M. B. (2010). Inconsistency of bootstrap: the Grenander estimator. Ann. Statist. 38 1953–1977. 10.1214/09-AOS777 2676880 (2011f:62046)
  * Sen and Xu [2015] [author] Sen, BodhisattvaB. and Xu, GongjunG. (2015). Model based bootstrap methods for interval censored data. Comput. Statist. Data Anal. 81 121–129. 10.1016/j.csda.2014.07.007 3257405
  * Silvapulle and Sen [2005] [author] Silvapulle, Mervyn JM. J. and Sen, Pranab KumarP. K. (2005). Constrained statistical inference: Inequality, order and shape restrictions. John Wiley & Sons.
  * Vervaat [1969] [author] Vervaat, W.W. (1969). Upper bounds for the distance in total variation between the binomial or negative binomial and the Poisson distribution. Statistica Neerlandica 23 79–86. 10.1111/j.1467-9574.1969.tb00075.x 242235