\date{Revised 25 February 2003}
\input{../courseDefs}

\begin{center}
Eleventh Set of Notes
\end{center}
\section{Resampling Methods}
For references, see
Beran, R., 1995. Stein confidence sets and the bootstrap,
{\em Stat. Sinica, 5}, 109--127;
Beran, R., 1990. Calibrating predictions regions,
{\em J. Amer. Stat. Assoc., 85}, 715--723;
Beran, R., 1990. Refining bootstrap simultaneous confidence sets,
{\em J. Amer. Stat. Assoc., 85}, 417--426;
Beran, R., 1987. Prepivoting to reduce level error of confidence sets,
{\em Biometrika, 74}, 457--468;
Feller, W., 1971. {\em An introduction to probability theory and its applications,
V. II}\/, 2nd edition, John Wiley and Sons, Inc., New York;
Lehmann, E.L., 1975. {\em Nonparametrics}\/, Holden-Day, Oakland;
Efron, B., 1982. {\em The Jackknife, the bootstrap, and other resampling plans,} SIAM,
Philadelphia;
Romano, J.P., 1988. A bootstrap revival of some nonparametric distance tests,
{\em J. Amer. Stat. Assoc., 83}, 698--708;
Romano, J.P., 1989. Bootstrap and randomization tests of some nonparametric hypotheses,
{\em Ann. Stat., 17}, 141--159.

\subsection{Permutation Tests}
Suppose we wish to test the efficacy of a purported cold remedy.
We select $n = n_c + n_t$ individuals for an experiment.
From these $n$, we randomly select $n_t$ to receive the remedy, and $n_c$ to receive
a placebo.
The selection is done in such a way that every subset of $n_t$ of the $n$ subjects
is equally likely to be the treatment group.
There are $\Choose{n}{n_t}$ such subsets.
We infect everyone with the cold virus.
Once they show symptoms, we give the remedy or
placebo.
Assume that everyone eventually recovers from a cold: we record the time interval
from the onset of symptoms to the complete remission of symptoms.
Let $\{ X_j \}_{j=1}^{n_c}$ be the times for the control individuals, and let
$\{Y_j\}_{j=1}^{n_t}$ be the times for the treated individuals.
Let $\bar{X} = \frac{1}{n_c} \sum_j X_j$ and $\bar{Y} = \frac{1}{n_t} \sum_j Y_j$.

Let the null hypothesis be that the remedy does not work. If that be true,
it is as if every individual's response were determined before the random assignment
to treatment or control. There would be (at most)
\beq
\Choose{n}{n_t} = \Choose{n}{n_c}
\eeq
possible pairs of values $(\bar{X}, \bar{Y})$, and so at most that many values
of the difference $\bar{X}-\bar{Y}$.
We would want to reject the null hypothesis if $\bar{X}$ is sufficiently larger than
$\bar{Y}$ (if the time-to-health with the treatment is shorter than the time-to-health with
the placebo).

Because, under the null hypothesis, the assignment to
treatment is random but individuals' responses are not, we can
obtain a critical value for testing
$\{H_0: \mbox{ treatment has no effect} \}$ at
level $\alpha$ by finding a value $c$ such that $\bar{X} - \bar{Y} > c$
for at most
a fraction $\alpha$ of the $\Choose{n}{n_t}$ partitions of the subjects into
treatment and control.
In principle, we could enumerate all of the $\Choose{n}{n_t}$ partitions,
calculate
$\bar{X}$ and $\bar{Y}$ for each of them, and find the critical
value exactly.
(Incidentally, asymptotically in $n$, $n_t$ and $n_c = n-n_t$ this test is equivalent to the
2-sample $t$-test.)
In practice, finding the exact
critical value is unwieldy even for fairly small values of $n$
(but computers get faster all the time).
For example, if we have 100 subjects of whom
50 are to be assigned to treatment, there are about 1.01e29 possible
assignments.
Each assignment requires on the order of 100 floating point
operations to compute $\bar{X} - \bar{Y}$.
Using a computer that can calculate at the rate 1GFlop/s ($10^9$
floating point
operations per second), it would take about $3.2e14$ years, to
compute all the values of $\bar{X} - \bar{Y}$.
The universe is about $1.5e10$ years old.
(For 50 subjects with 25 assigned to treatment, the figure would
be about 0.4 years.)

However, we might approximate a critical value for the test by
randomly partitioning
the $n$ subjects into groups of size $n_t$ and $n_c$,
computing $\bar{X} - \bar{Y}$
for them, and using the $1-\alpha$ quantile of the empirical
distribution of those
values among the random partitions.
Such an approach {\em resamples} the data: we repeatedly
draw samples of size
$n_t$ without replacement, compute the average response
$\bar{Y}$ of those $n_t$
subjects, the average response $\bar{X}$ of the $n_c = n - n_t$
not in that sample, and record $\bar{X} - \bar{Y}$.


\subsection{The Bootstrap}
The setting for the next few lectures is that we observe an iid sample
of size $n$,
$\{X_j \}_{j=1}^n$ iid $F$.
Each observation is real-valued.
We wish to estimate some parameter of the distribution of
$F$ that can be written as a functional of $F$, $T(F)$.
Examples include the mean, $T(F) = \int x dF(x)$, other moments,
{\em etc.}

The (unpenalized) nonparametric maximum likelihood estimator of $F$ from the data
$\{X_j \}$ is
just the empirical distribution $\hat{F}_n$,
which assigns mass $1/n$ to each observation:
\beq
\arg \max_{\mbox{probability distributions }G}
\Prob_G \{ X_j = x_j, \; j=1, \ldots, n \}
= \hat{F}_n.
\eeq
(Note, however, that the MLE of $F$ is not generally consistent in
problems with an
infinite number of parameters, such as estimating a density or a
distribution function.)

Using the general rule that the maximum likelihood estimator of a
function of a parameter is that function of the maximum likelihood
estimator of the parameter, we might be led to consider $T(\hat{F}_n)$
as an estimator of $T(F)$.

That is exactly what the sample mean does, as an estimator of the mean:
\beq
        T(\hat{F}_n) = \int x d\hat{F}_n(x) = \sum_{j=1}^n\frac{1}{n}X_j =
        \frac{1}{n} \sum_j X_j.
\eeq

Similarly, the maximum likelihood estimator of
\beq
    \Var(X) = T(F) = \int \left ( x - \int x dF \right )^2dF
\eeq
is
\beq
        T(\hat{F}_n) = \int \left ( x - \int x d\hat{F} \right )^2 d\hat{F}_n =
        \frac{1}{n}\sum_j \left (X_j - \frac{1}{n} \sum_k X_k \right )^2.
\eeq

In these cases, we get analytically tractable expressions for $T(\hat{F}_n)$.
What is often more interesting is estimating a property of the
distribution of $T(\hat{F}_n)$, for example the variance of the estimator
$T(\hat{F}_n)$.
The bootstrap approach is to approximate the distribution of
$T(\hat{F}_n)$ by the distribution of   $T(\hat{F}_n^*)$.
What this notation is intended to mean is that we approximate
the distribution of an estimator applied to the empirical
distribution $\hat{F}_n$ of a random sample of size
$n$ from a distribution $F$ by the distribution of that estimator
applied to a random sample $\hat{F}_n^*$ of size $n$ from a particular
realization $\hat{F}_n$ of
the empirical distribution of a sample of size $n$ from $F$.
In the case $T$ is the mean, so
$T(\hat{F}_n)$ is the sample mean, we could obtain
analytically the variance of the distribution of $T( \hat{F}_n^* )$:
Let $\{ X_j^* \}_{j=1}^n$ be an iid sample of size $n$ from $\hat{F}_n$.
Then
\beq
        \Var_{\hat{F}_n} \frac{1}{n}\sum_{j=1}^n X_j^*
        = \frac{1}{n^2} \sum_{j=1}^n (X_j - \bar{X})^2,
\eeq
where $\{X_j \}$ are the ``original'' data and $\bar{X}$ is their mean.
Generally, we would not get an analytically tractable espression
for the variance of an estimator under resampling from the empirical
distribution, but we could still approximate the distribution of
$T(\hat{F}_n)$ by
generating
size-$n$ data sets that are iid $\hat{F}_n$ (drawing samples of
size $n$ with replacement
from $\{x_j \}_{j=1}^n$), and applying $T$ to each.

The idea of the bootstrap is to approximate the distribution (under $F$)
of an estimator
$T(\hat{F}_n)$ by the distribution of the estimator under $\hat{F}_n$,
and to approximate {\em that} distribution by using a computer to
take pseudo-random samples of size
$n$ from  $\hat{F}_n$.

This basic idea is quite flexible, and can be applied to a wide variety of
testing and estimation problems, including finding confidence intervals
and sets.  It is related to some other ``resampling'' schemes in which
one re-weights the data to form other distributions.

\noindent
{\bf Exercise.}
Recall the normal linear regression model of the previous section.
Consider resampling as follows:
Let
\beq
        \hat{\epsilon} = A \hat{\theta}_{LS} - X,
\eeq
and let $\hat{F}_n$ be the empirical cdf of the components of $\hat{\epsilon}$.
Consider ``pseudo-data'' distributed as
\beq
        X^* = A \hat{\theta}_{LS} + \epsilon^*,
\eeq
where the components of $\epsilon^*$ are iid $\hat{F}_n$.
Consider applying least squares to the pseudo-data, and calculating the
variance of the components of $\theta$.
(If we were using Monte-Carlo, we would generate many independent sets of
pseudo-data, apply least-squares to each, and get an empirical joint distribution for
$\hat{\theta}_{\mbox{LS}}$ under $ \epsilon^*   \sim \hat{F}_n$.)
Give a formula for the resulting estimate of the variance of $\theta_j$.


Before doing more theory with the bootstrap, let's examine the jackknife.

\subsection{The Jackknife}
The idea behind the jackknife, which is originally due to Tukey and
Quenouille, is to
form from the data $\{ X_j \}_{j=1}^n$, $n$ sets of $n-1$ data,
leaving each datum out
in turn.
The ``distribution'' of $T$ applied to these $n$ sets is used to approximate
the distribution of $T{\hat{F}_n}$.
Let $\hat{F}_{(i)}$ denote  the empirical distribution of the data
set with the $i$th value deleted;
$T_{(i)} = T( \hat{F}_{(i)})$ is the corresponding
estimate of $T(F)$.
An estimate of the expected value of $T(\hat{F}_n)$ is
\beq
    \hat{T}_{(\cdot)} = \frac{1}{n} \sum_{i=1}^n T( \hat{F}_{(i)}) .
\eeq
Consider the bias of $T(\hat{F}_n)$:
\beq
    E_F T(\hat{F}_n) - T(F).
\eeq
Quenoille's jackknife estimate of the bias is
\beq
    \widehat{\mbox{BIAS}} = (n-1) (\hat{T}_{(\cdot)} - T(\hat{F}_n) ).
\eeq
It can be shown that if the bias of $T$ has a homogeneous
polynomial expansion
in $n^{-1}$ whose coefficients do not depend on $n$,
then the bias of the bias-corrected estimate
\beq
    \tilde{T} = nT(\hat{F}_n) - (n-1) T_{(\cdot)}
\eeq
is $O(n^{-2})$ instead of $O(n^{-1})$.

Applying the jackknife estimate of bias to correct
the plug-in estimate of
variance reproduces the formula for the sample variance (with
$1/(n-1)$) from the formula with $1/n$:
Define
\beq
    \bar{X} = \frac{1}{n} \sum_{j=1}^n X_j,
\eeq
\beq
    \bar{X}_{(i)} = \frac{1}{n-1} \sum_{j \ne i} X_j,
\eeq
\beq
    T(\hat{F}_n) = \hat{\sigma}^2 = \frac{1}{n} \sum_{j=1}^n (X_j - \bar{X})^2,
\eeq
\beq
    T(\hat{F}_{(i)}) = \frac{1}{n-1} \sum_{j \ne i} ( X_j - \bar{X}_{(i)})^2,
\eeq
\beq
    T(\hat{F}_{(\cdot)} = \frac{1}{n} \sum_{i=1}^n T(\hat{F}_{(i)}).
\eeq
Now
\beq
    \bar{X}_{(i)} = \frac{n\bar{X} - X_i}{n-1} = \bar{X} + \frac{1}{n-1} (\bar{X} - X_i),
\eeq
so
\begin{eqnarray}
     ( X_j - \bar{X}_{(i)})^2 &=&
        \left ( X_j - \bar{X} + \frac{1}{n-1} (\bar{X} - X_i) \right )^2 \nonumber \\
      &=& (X_j - \bar{X})^2 + \frac{2}{n-1} (X_j - \bar{X})(X_i - \bar{X}) +
      \frac{1}{(n-1)^2}(X_i - \bar{X})^2.
\end{eqnarray}
Note also that
\beq
    \sum_{j \ne i} (X_j - \bar{X}_{(i)})^2 =
    \sum_{j=1}^n (X_j - \bar{X}_{(i)})^2 - (X_i - \bar{X}_{(i)})^2.
\eeq
Thus
\begin{eqnarray}
    \sum_{i=1}^n \sum_{j \ne i} (X_j - \bar{X}_{(i)})^2 &=&
    \frac{1}{n-1} \sum_{i=1}^n \left [ \sum_{j=1}^n \left [
    (X_j - \bar{X})^2 + \frac{2}{n-1}(X_j - \bar{X})(X_i - \bar{X}) \right . \right . + \nonumber \\
    && + \left . \left . \frac{1}{(n-1)^2}(X_i - \bar{X})^2
    \right ]  - (X_i - \bar{X})^2 - \right . \nonumber \\
    &&  - \left . \left .
    \frac{2}{n-1}(X_i - \bar{X})^2 - \frac{1}{(n-1)^2}(X_i - \bar{X})^2 \right . \right ].
\end{eqnarray}
The last three terms all are multiples of $(X_i - \bar{X})^2$; the sum of the coefficients
is
\beq
    1 + 2/(n-1) + 1/(n-1)^2 = n^2/(n-1)^2.
\eeq
The middle term of the inner sum is a constant times $(X_j - \bar{X})$, which sums to zero over $j$.
Simplifying the previous displayed equation yields
\begin{eqnarray}
    \sum_{i=1}^n \sum_{j \ne i} (X_j - \bar{X}_{(i)})^2
    &=&
        \frac{1}{n-1} \sum_{i=1}^n
        \left ( n \hat{\sigma}^2 + \frac{n}{(n-1)^2}(X_i - \bar{X})^2 -
        \frac{n^2}{(n-1)^2} (X_i - \bar{X})^2
        \right ) \nonumber \\
    &=&
        \frac{1}{n-1} \sum_{i=1}^n (n \hat{\sigma}^2 - \frac{n}{n-1} (X_i - \bar{X})^2 ) \nonumber \\
    &=&
        \frac{1}{n-1} \left [ n^2 \hat{\sigma}^2 - \frac{n^2}{n-1} \hat{\sigma}^2 \right ] \nonumber \\
    &=&
        \frac{n(n-2)}{(n-1)^2}\hat{\sigma}^2.
\end{eqnarray}
The jackknife bias estimate is thus
\beq
    \widehat{\mbox{BIAS}} = (n-1)\left ( T(\hat{F}_{(\cdot)}) - T(\hat{F}_n) \right )
    = \hat{\sigma}^2 \frac{n(n-2) - (n-1)^2}{n-1} = \frac{-\hat{\sigma}^2}{n-1}.
\eeq
The bias-corrected MLE variance estimate is therefore
\beq
    \hat{\sigma}^2 \left ( 1 - \frac{1}{n-1} \right ) =
    \hat{\sigma}^2 \frac{n}{n-1} = \frac{1}{n-1} \sum_{j=1}^n (X_j - \bar{X})^2 = S^2,
\eeq
the usual sample variance.

The jackknife also can be used to estimate other properties of an
estimator, such as its variance.
The jackknife estimate of the variance of $T(\hat{F}_n)$ is
\beq
\hat{\Var}(T) = \frac{n-1}{n} \sum_{j=1}^n ( T_{(j)} - T_{(\cdot)} )^2.
\eeq

It is convenient to think of distributions on data sets to compare
the jackknife and the bootstrap.
We shall follow the notation in Efron (1982).
We condition on $(X_i = x_i)$ and treat the data as fixed in what
follows.
Let $\calS_n$ be the $n$-dimensional simplex
\beq
    \calS_n \equiv \{ \Prob^* = (P_i^*)_{i=1}^n \in \bfR^n
    : P_i^* \ge 0 \mbox{ and } \sum_{i=1}^n P_i^* = 1 \}.
\eeq
A {\em resampling vector}
$\bootProb = (P_k^*)_{k=1}^n $ is any element of $\calS_n$;
{\em i.e.}, an $n$-dimensional discrete probability vector.
To each $\bootProb = (P_k^*) \in \calS_n$ there corresponds
a re-weighted  empirical measure $\hat{F}^*(\bootProb)$ which
puts mass $P_k^* $ on $x_k$, and a value of the estimator
$T^* = T(\hat{F}(\bootProb)) = T(\bootProb)$.
The resampling vector $\Prob^0 = (1/n)_{j=1}^n$ corresponds to the
empirical distribution $\hat{F}_n$ (each datum $x_j$ has the same
mass).
The resampling vector
\beq
    \jackProb{i} = \frac{1}{n-1}(1, 1, \ldots, 0, 1, \ldots, 1),
\eeq
which has the zero in the $i$th place, is one of the $n$
resampling vectors the jackknife visits; denote the
corresponding value of the estimator $T$ by $T_{(i)}$.
The bootstrap visits all resampling vectors whose components are
multiples of $1/n$.

The bootstrap estimate of variance tends to be better than the
jackknife estimate of variance for nonlinear estimators because of
the distance between the empirical measure and the resampled measures:
\beq
    \| \bootProb - \Prob^0 \| = O_P(n^{-1/2}),
\eeq
while
\beq
    \| \jackProb{k} - \Prob^0 \| = O(n^{-1}).
\eeq
To see the former, recall that the difference between the
empirical distribution and the true distribution is $O_P(n^{-1/2})$:
For any two probability distributions $\Prob_1$, $\Prob_2$, on
$\bfR$, define the
Kolmogorov-Smirnov distance
\beq
    d_{KS}(\Prob_1, \Prob_2) \equiv
    \| \Prob_1 - \Prob_2 \|_{KS} \equiv \sup_{x \in \bfR} |
    \Prob_1\{(-\infty, x]\} - \Prob_2\{(-\infty, x]\} |.
\eeq
There exist universal constants $\chi_\alpha$
so that for every continuous (w.r.t. Lebesgue measure) distribution
$F$,
\beq
    \Prob_F
    \left \{ \| F - \hat{F}_n \|_{KS} \ge \chi_n(\alpha) \right \}
    = \alpha.
\eeq
This is the Dvoretzky-Kiefer-Wolfowitz inequality.
Massart ({\em Ann. Prob., 18}, 1269--1283, 1990;
stable URL
http://links.jstor.org/sici?sici=0091-1798\%28199007\%2918\%3A3\%3C1269\%3ATTCITD\%3E2.0.CO\%3B2-Q )
showed that the constant
\beq
    \chi_n(\alpha) \le \sqrt{\frac{\ln \frac{2}{\alpha}}{2n}}
\eeq
is {\em tight}.
Thinking of the bootstrap distribution (the empirical distribution $\hat{F}_n$) as the true
cdf and the resamples from it as the data gives the result that the distance between
the cdf of the bootstrap resample and the empirical cdf of the original data
is $O_P(n^{-1/2})$.

To see that the cdfs of the jackknife samples are $O(n^{-1})$ from the
empirical cdf $\hat{F}_n$, note that for univariate real-valued data, the difference
between $\hat{F}_n$ and the cdf of the jackknife data set that
leaves out the $j$th ranked observation $X_{(j)}$ is largest either at $X_{(j-1)}$ or
at $X_{(j)}$.
For $j = 1$ or $j = n$, the jackknife
samples that omit the smallest or largest observation, the $L_1$
distance between the jackknife measure and the empirical distribution
is exactly $1/n$.
Consider the jackknife cdf $\hat{F}_{n,(j)}$, the cdf of the sample without $X_{(j)}$,
$1 < j < n$.
\beq
    \hat{F}_{n,(j)}(X_{(j)}) = (j-1)/(n-1),
\eeq
while $\hat{F}_n((X_{(j)}) = j/n$; the difference is
\beq
    \frac{j}{n} - \frac{j-1}{n-1} = \frac{j(n-1) - n(j-1)}{n(n-1)} =
    \frac{n-j}{n(n-1)} = \frac{1}{n-1} - \frac{j}{n(n-1)}.
\eeq
On the other hand,
\beq
    \hat{F}_{n,(j-1)}(X_{(j-1)}) = (j-1)/(n-1),
\eeq
while $\hat{F}_n((X_{(j-1)})= (j-1)/n$; the difference is
\beq
    \frac{j-1}{n-1} - \frac{j-1}{n} = \frac{n(j-1) - (n-1)(j-1)}{n(n-1)} =
    \frac{j - 1}{n(n-1)}.
\eeq
Thus
\beq
    \| \hat{F}_{n,(j)} - \hat{F}_{n} \| = \frac{1}{n(n-1)} \max\{n-j, j-1\}.
\eeq
But $n/2 \le \max\{n-j, j-1\} \le n-1$, so
\beq
    \| \hat{F}_{n,(j)} - \hat{F}_n\| = O(n^{-1}).
\eeq

The neighborhood that the bootstrap samples is larger, and is
probabilistically of the right size to correspond to the uncertainty
of the empirical distribution function as an estimator of the
underlying distribution function $F$ (recall the
Kiefer-Dvoretzky-Wolfowitz inequality---a K-S ball of radius $O(n^{-1/2})$
has fixed coverage probability).
For linear functionals, this does not matter, but for strongly nonlinear
functionals, the bootstrap estimate of the variability tends to be more accurate than
the jackknife estimate of the variability.

Let us have a quick look at the distribution of the K-S distance between
a continuous distribution and the empirical distribution of a
sample $\{X_j\}_{j=1}^n$ iid $F$.
The discussion follows {\em Feller} (1971, pp.~36ff).
First we show that for continuous distributions $F$, the distribution of
$\| \hat{F}_n - F \|_{KS}$ does not depend on $F$.
To see this, note that $F(X_j) \sim U[0, 1]$:
Let $x_t \equiv \inf \{x \in \bfR : F(x_t) = t \}$.
Continuity of $F$ ensures that $x_t$ exists for all $t \in [0, 1]$.
Now the event $\{ X_j \le x_t \}$ is equivalent to the
event $\{F(X_j) \le F(x_t)\}$ up to a set of $F$-measure zero.
Thus
\beq
    t = \Prob_F \{ X_j \le x_t \} = \Prob_F \{ F(X_j) \le F(x_t) \} =
    \Prob_F \{ F(X_j) \le t \}, \,\, t \in [0, 1].
\eeq
Thus $\{ F(X_j) \}_{j=1}^n$ are iid $U[0, 1]$.
Let
\beq
    \hat{G}_n(t) = \#\{F(X_j) \le t\}/n =
    \#\{X_j \le x_t\}/n = \hat{F}_n (x_t)
\eeq
be the empirical cdf of the $F(X_j)$s.
Note that
\beq
    \sup_{x \in \bfR} | \hat{F}_n(x) - F(x) | =
    \sup_{t \in [0, 1]} | \hat{F}_n(x_t) - F(x_t) | =
    \sup_{t \in [0, 1]} | \hat{G}_n (t) - t |.
\eeq
The probability distribution of $\hat{G}_n$ is that of the cdf of $n$ iid
$U[0, 1]$ random variables (it does not depend on $F$), so the distribution
of the K-S distance between the empirical cdf and the true cdf is the
same for every continuous distribution.
It turns out that for distributions with atoms, the K-S distance between
the empirical and the true distribution functions is stochastically
smaller than it is for continuous distributions.

\subsection{Bootstrap and Randomization Tests}
This section is about Romano's papers. The set-up is as follows:
We observe $\{X_j \}_{j=1}^n$ iid $F$, where $F$ is a distribution on an
abstract sample space $\calS$.
The distribution $F \in \Omega$, a collection of distributions on $\calS$.
The null hypothesis is that $F \in \Omega_0 \subset \Omega$.  We assume that
$\Omega_0$ can be characterized as a set of distributions that are invariant under
a transformation on $\Omega$: let $\tau: \Omega \rightarrow \Omega_0$; we assume
that $\tau(G) = G$ for all $G \in \Omega_0$.

Let $\calV$ be a collection of subsets of a set $\calS$.
For a finite set $D \subset S$, let $\Delta^\calV(D)$ be the number of distinct sets
$\{ V \cap D: V \in \calV \}$.
For positive integers $n$, let
\beq
m^\calV(n) = \max_{D \subset \calS: \#\{D\} = n } \Delta^\calV(D).
\eeq
Let
\beq
    c(\calV) \equiv \inf \{ n : m^\calV(n) < 2^n \}.
\eeq
If $c(\calV) < \infty$, $\calV$ is a Vapnik-Cervonenkis (V-C) class.
That is, $\calV$ is a V-C class if
if the maximum number of distinct intersections of sets in $\calV$
with sets containing $n$ points grows sub-exponentially with $n$.
Intersections, finite unions, and Cartesian products of V-C classes are V-C classes.
In $\bfR^n$, the set of all ellipsoids, the set of all half-spaces, the set of all
lower-left quadrants, and the set of all convex sets with at most $p$ extreme
points are all V-C classes.

An alternative, equivalent definition of a V-C class is based on the following
\begin{Definition}
        Suppose $\calV$ is a collection of subsets of a set $\calX$, and that $D$ is a finite
        subset of $\calX$.
        We say $D$ is {\em shattered} by $\calV$ if every subset $d \subset D$
        can be written $d = V \cap D$ for some $V \in \calV$.
\end{Definition}

Because there are $2^n$ subsets of a set with $n$ elements, this is equivalent to saying
that there are $2^n$ different subsets of the form $D \cap V$ as $V$ ranges over $\calV$.

A collection $\calV$ is a V-C class if for some finite integer $n$, there exists a set
$D \subset \calX$ with $n$ elements that is not shattered by $\calV$.

\noindent
{\bf Exercise.}
Show that intersections and finite unions of V-C classes are V-C classes.
Give a (counter) example showing that a countable union of V-C classes need not
be a V-C class.

Let $\calV$ be a VC
class of subsets of $\calS$.
Define the pseudo-metric
\begin{eqnarray}
        \delta: \Omega \times \Omega &\rightarrow& \bfR^+
        \nonumber \\
        (G, H) & \rightarrow & \sup_{V \in \calV} | G(V) - H(V) | .
\end{eqnarray}
This is a generalization of the Kolmogorov-Smirnov distance for distributions on the
line.
In that case, the sets in $\calV$ are the half-lines $(-\infty, y]$ (which comprise a
V-C class).
Consider a set $D = \{x_j\}_{j=1}^n$ of points on the real line.
How many sets are there  of the form $V \cap D$, for $V$ a half-line
$(-\infty, y]$?
Just $n+1$.
Suppose the points are in increasing order.
Then the possibilities for $V \cap D$ are $\{\}$,
$\{ x_1 \}$, $\{ x_1, x_2 \}$, $\ldots$, $\{x_j \}_{j=1}^n$.
Thus $m^{\calV}(n) = n+1$, and $c(\calV) \equiv \inf \{ n : m^\calV(n) < 2^n \} = 2$
(for $n=0$, we have $0+1 = 2^0$, and for $n=1$, we have $1+1 = 2^1$, but for $n=2$,
we have $2+1 < 2^2$).


Assume that $\calV$ has been selected such that
$\delta(G, \tau G) = 0$ iff
$G \in \Omega_0$.
Romano proposes using the test statistic
\beq
    T_n = n^{1/2} \delta(\hat{P}_n, \tau \hat{P}_n ),
\eeq
where $\hat{P}_n$ is the empirical measure of $\{X_j \}_{j=1}^n$.
One rejects the hypothesis when  $\tau\hat{P}_n$ is far from $\hat{P}_n$.

One way to obtain a critical value from the test is by the bootstrap; one
resamples from $\hat{P}_n$, and looks at the distribution of the distance
between the empirical distribution of the bootstrap samples and $\tau$ applied to them.

Suppose that there is a (known) group $\calG_n$ of transformations of the sample space
$\calS_n$
such that under the null hypothesis, $F$ is invariant under $\calG_n$.
Then we can also construct a {\em randomization test} of the hypothesis $H_0$.
For simplicity, suppose that $\calG_n$ is finite, with $M_n$ elements
$\{ g_{nj} \}_{j=1}^{M_n}$.
Under the null hypothesis, conditional on $X = x$, the values $\{ g_{nj}x \}_{j=1}^{M_n}$
are equally likely.
The {\em orbit} of an element $x$ of a group $\calG$ is the set of all elements of
$\calG$ that can be obtained by (repeatedly) applying elements of $\calG$ to $x$.
For example, consider points in the plane and the group of rotations about the
origin.
Then the orbit of a point $x$ is the circle with radius $\|x\|$.
Compute the  test statistic for each $ g_{nj}x$ in the orbit of $x$.
Reject the null hypothesis if the statistic for $x$ exceeds the $\alpha$ quantile of
the test statistic for the set of values obtained from the orbit;
do not reject if it is less; reject with a given probability if the statistic equals
the $\alpha$ quantile in such a way as to get a level $\alpha$ test.
This is a randomization test.
Because the level of the randomization test is $\alpha$, conditional on the data,
integrating
over the distribution of the data shows that it is $\alpha$ unconditionally.

\subsubsection{Examples of hypotheses and functions $\tau$}.
Examples Romano gives include testing for independence of the components of each $X_j$,
testing for exchangeability of the components of each $X_j$, testing for spherical
symmetry of the distribution of $X_j$, testing for homogeneity among the $X_j$, and
testing for a change point.

In the example of testing for independence, the mapping $\tau$ takes the marginal
distributions of the joint distribution, then constructs a joint distribution that
is the product of the marginals.  For distributions with independent components,
this is the identity; otherwise, it maps a distribution into one with the same
marginals, but whose components are independent.
For testing for spherical symmetry, $\tau$ maps a distribution into one with the same  mass
at every distance from the origin, but that is uniform on spherical shells.
For testing for exchangability, Romano proposes looking at the largest difference between
$F$ and a permutation of the coordinates of $F$, over all permutations of the coordinates.
See his excellent and lucid paper for more details.

Romano shows that these tests are consistent against all alternatives, and that
the critical values given by the bootstrap and by randomization are asymptotically
equal with probability one.  Because the randomization tests are exact level $\alpha$
tests, they might be preferred.  Romano also discusses how to implement the tests
computationally.

\subsection{Bootstrap Confidence Sets}
Let $\calU$ be an index set (not necessarily countable).
Recall that a collection $\{ \calI_u \}_{u \in \calU}$ of confidence intervals for
parameters $\{\theta_u \}_{u \in \calU}$ has simultaneous $1-\alpha$
coverage probability if
\beq
\Prob_{\theta} \left \{ \cap_{u \in \calU} \{\calI_u \ni \theta_u \} \right \}
\ge 1-\alpha.
\eeq
If $\Prob \{ \calI_u \ni \theta_u\}$ does not depend on $u$, the confidence
intervals are said to be {\em balanced.}

Many of the procedures for forming joint confidence sets we have seen depend on
{\em pivots}, which are functions of the data and the parameter(s) whose
distribution is known (even though the parameter and the parent distribution are not).
For example, the Scheff\'{e} method relies on the fact that
the sum of squared differences between the data and the corresponding parameters,
divided by the variance estimate, has an $F$ distribution, regardless of the
parameter values.
Similarly, Tukey's maximum modulus method relies on the fact that the distribution
of the maximum of the studentized
absolute differences between the data and the corresponding
parameters   does not depend on the parameters.
Both of those examples are parametric, but the idea is more general:
the procedure we looked at for finding bounds on the density function subject
to shape restrictions just relied on the fact that there are uniform
bounds on the probability that the K-S distance between the empirical
distribution and the true distribution exceeds some threshold.

Even in cases where there is no known exact pivot, one can sometimes show that
some function of the data and parameters is asymptotically a pivot.
Working out the distributions of the functions involved is not typically
straightforward, and a general method of constructing (possibly
simultaneous) confidence sets  would be nice.

Efron gives several methods of basing confidence sets on the bootstrap.
Those methods are substantially improved (in theory, and in my experience)
by Beran's pre-pivoting approach, which leads to iterating the bootstrap.

Let $X_n$ denote a sample of size $n$ from $F$.
Let $R_n(\theta) = R_n(X_n, \theta)$ have cdf
$H_n$, and let $H_n^{-1}(\alpha)$
be the largest $\alpha$ quantile of the distribution of $R_n$.
Then
\beq
        \{ \gamma \in \Theta : R_n(\gamma) \le H_n^{-1}(1-\alpha) \}
\eeq
is a $1-\alpha$ confidence set for $\theta$.

\subsubsection{The Percentile Method}
The idea of the percentile method is to use the empirical bootstrap percentiles of
some quantity to approximate the true percentiles.
Consider constructing a confidence interval for a single real parameter $\theta = T(F)$.
We will estimate $\theta$ by $\hat{\theta} = T(\hat{F}_n)$.
We would like to know the distribution function $H_n = H_n(\cdot, F)$ of
$D_n (\theta) = T(\hat{F}_n) - \theta$.
Suppose we did.
Let $H_n^{-1}(\cdot) = H_n^{-1}(\cdot, F)$ be the inverse cdf of $D_n$.
Then
\beq
    \Prob_F \{ H_n^{-1}(\alpha/2) \le T(\hat{F}_n) - \theta  \le  H_n^{-1}(1- \alpha/2)
    \} = 1-\alpha,
\eeq
so
\beq
        \Prob_F \{ \theta \le T(\hat{F}_n) -  H_n^{-1}(\alpha/2)  \mbox{ and }
        \theta \ge T(\hat{F}_n) - H_n^{-1}(1-\alpha/2) \} = 1-\alpha,
\eeq
or, equivalently,
\beq
        \Prob_F \{ [T(\hat{F}_n) - H_n^{-1}(1-\alpha/2), T(\hat{F}_n) -  H_n^{-1}(\alpha/2)]
        \ni \theta \} = 1-\alpha,
\eeq
so the interval
$[T(\hat{F}_n) - H_n^{-1}(1-\alpha/2), T(\hat{F}_n) -  H_n^{-1}(\alpha/2)]$
would be a $1-\alpha$ confidence interval for $\theta$.

The idea behind the percentile method is to approximate $H_n(\cdot, F)$ by
$\hat{H}_n = H_n(\cdot, \hat{F}_n)$, the distribution of $D_n$
under resampling from $\hat{F}_n$ rather than $F$.
An alternative approach would be to take
$D_n (\theta) = |T(\hat{F}_n) - \theta|$; then
\beq
        \Prob_F \{ | T(\hat{F}_n) - \theta | \le  H_n^{-1}(1- \alpha)
        \} = 1-\alpha,
\eeq
so
\beq
        \Prob_F \{ [ T(\hat{F}_n -  H_n^{-1}(1-\alpha) , T(\hat{F}_n +  H_n^{-1}(1-\alpha)]
        \ni \theta  \} = 1-\alpha.
\eeq
In either case, the ``raw'' bootstrap approach is to approximate $H_n$ by resampling
under $\hat{F}_n$.

Beran proves a variety of results under the following condition:

\noindent
{\bf Condition 1.} (Beran, 1987)
For any sequence $\{F_n\}$ that converges to
$F$ in a metric $d$ on
cdfs, $H_n(\cdot, F_n)$ converges weakly to a continuous cdf
$H = H(\cdot, F)$ that depends only on $F$, and not the sequence $\{F_n\}$.

Suppose Condition 1 holds.
Then because $\hat{F}_n$ is consistent for $F$, the estimate $\hat{H}_n$ converges
in probability to $H$ in sup norm; moreover, the distribution of
$\hat{H}_n(R_n(\theta))$ converges to $U[0,1]$.

Instead of $D_n$, consider $R_n(\theta) = | T(\hat{F}_n) - \theta |$ or some
other (approximate) pivot.
Let $\hat{H}_n(\cdot, \hat{F}_n)$ be the bootstrap estimate of the cdf of $R_n$;
The set
\begin{eqnarray}
        B_n &=& \{ \gamma \in \Theta : \hat{H}_n (R_n(\gamma)) \le 1-\alpha \}
            \nonumber \\
        &=& \{ \gamma \in \Theta : R_n(\gamma) \le \hat{H}_n^{-1}(1-\alpha) \}
\label{eq:BnDef}
\end{eqnarray}
is (asymptotically) a $1-\alpha$ confidence set for $\theta$.

The level of this set for finite samples tends to be inaccurate.
It can be improved in the following way, due to Beran.

The original root, $R_n(\theta)$, whose limiting distribution depends on $F$,
was transformed into a new root $R_{n,1}(\theta) = \hat{H}_n(R_n(\theta) )$,
whose limiting distribution is $U[0,1]$.  The distribution of $R_{n,1}$
depends less strongly on $F$ than does that of $R_n$; Beran calls
mapping $R_n$ into $R_{n,1}$  {\em prepivoting}. The confidence set
\ref{eq:BnDef} acts as if the distribution of $R_{n,1}$ really is uniform,
which is not generally true.  One could instead treat $R_{n,1}$ itself as a root,
and pivot to reduce the dependence on $F$.

Let $H_{n,1} = H_{n,1}(\cdot, F)$ be the cdf of the new root $R_{n,1}(\theta)$,
estimate $H_{n,1}$ by $\hat{H}_{n,1} = H_{n,1}(\cdot, \hat{F}_n)$, and define
\begin{eqnarray}
        B_{n,1} &=& \{ \gamma \in \Theta : \hat{H}_{n,1}(R_{n,1}(\gamma)) \le 1-\alpha \}
             \nonumber \\
        &=& \{ \gamma \in \Theta : \hat{H}_{n,1}(\hat{H}_n(R_n(\gamma))) \le 1-\alpha \}
        \nonumber \\
        &=&   \{ \gamma \in \Theta :
                R_n(\gamma) \le \hat{H}_n^{-1}(\hat{H}_{n,1}^{-1}(1-\alpha)))
        \}.
\label{eq:Bn1Def}
\end{eqnarray}
Beran shows that this confidence set tends to have smaller error in its level than
does $B_n$.
The transformation can be iterated further, typically
resulting in additional reductions in the
level error.

\subsection{Approximating $B_{n,1}$ by Monte Carlo}
I'll follow Beran's (1987) notation (mostly).

Let $x_n$ denote the ``real'' sample of size $n$.
Let $x_n^*$ be a bootstrap sample of size $n$ drawn from the empirical cdf $\hat{F}_n$.
The components of $x_n^*$ are conditionally iid  given $x_n$.
Let $\hat{F}_n^*$ denote the ``empirical'' cdf of the bootstrap sample $x_n^*$.
Let $x_n^{**}$ denote a sample of size $n$ drawn from $\hat{F}_n^*$; the components of
$x_n^{**}$ are conditionally iid given $x_n$ and $x_n^*$.
Let $\hat{\theta}_n = T(\hat{F}_n)$, and $\hat{\theta}_n^* = T(\hat{F}_n^*)$.
Then
\beq
H_n(s, F) = \Prob_F \{ R_n (x_n, \theta) < s \}  ,
\eeq
and
\beq
H_{n,1}(s, F) = \Prob_F \left \{ \Prob_{\hat{F}_n} \{ R_n ( x_n^*, \hat{\theta}_n )
<  R_n(x_n, \theta) \} < s \right \}.
\eeq
The bootstrap estimates of these cdfs are
\beq
\hat{H}_n(s) = H_n(s, \hat{F}_n ) = \Prob_{\hat{F}_n} \{ R_n ( x_n^*, \hat{\theta}_n ) < s
\},
\eeq
and
\beq
\hat{H}_{n,1}(s) = H_{n,1}(s, \hat{F}_n) = \Prob_{\hat{F}_n}
\left \{ \Prob_{\hat{F}_n^*} \{ R_n(x_n^{**}, \hat{\theta}_n^* )
< R_n(x_n^*, \hat{\theta}_n) \} \right \}.
\eeq

The Monte Carlo approach is as follows:
\begin{enumerate}
\item
Draw $\{ y_k^* \}_{k=1}^M$ bootstrap samples of size $n$ from $\hat{F}_n$.
The ecdf of $\{ R_n(y_k^*, \hat{\theta}_n) \}_{k=1}^M $ is an approximation to $\hat{H}_n$.
\item
For $k = 1, \cdots, M$, let $\{ y_{k\ell}^{**} \}_{\ell=1}^N$ be
$N$ size $n$ bootstrap samples from the ecdf of $y_k^*$.
Let $\hat{\theta}_{n,k}^* = T(\hat{F}_{n,k}^*)$.
Let $Z_k$ be the fraction of the values
\beq
\{ R_n(y_{k,\ell}^{**}, \hat{\theta}_{n,k}^*  ) \}_{\ell=1}^N
\eeq
that are less than $R_n(y_k^*, \hat{\theta}_n)$.
The ecdf of $\{ Z_k \}$ is an approximation to $\hat{H}_{n,1}$ that improves
(in probability) as $M$ and $N$ grow.

\end{enumerate}

Note that this approach is extremely general.  Beran gives examples for
confidence sets for directions, {\em etc}.  The pivot can in principle be
a function of any number of parameters, which can yield simultaneous confidence
sets for parameters of any dimension.

\subsection{Other approaches to improving coverage probability }
There are other ways of iterating the bootstrap to improve the level accuracy of
bootstrap confidence sets.
Efron suggests trying to attain a different coverage probability so that
the coverage attained in the second generation samples is the nominal coverage probability.
That is, if one wants a 95\% confidence set, one tries different percentiles so that in
resampling from the sample, the attained coverage probability is 95\%.  Typically, the
percentile one uses in the second generation will be higher than 95\%.
Here is a sketch of the Monte-Carlo approach:
\begin{itemize}
\item Set a value of $\alpha^*$ (initially taking $\alpha^* = \alpha$ is reasonable)
\item From the sample, draw $M$ size-$n$ samples that are each iid
$\hat{F}_n$. Denote the ecdfs of the samples by $\{ \hat{F}_{n,j}^*\}$.
\item For each $j = 1, \ldots, M$, apply the percentile method to make a (nominal) level
$1-\alpha^*$ confidence
interval for $T(\hat{F}_n)$.
This gives $M$ confidence intervals; a fraction $1-\alpha'$ will cover
$T(\hat{F}_n)$. Typically,
$1- \alpha' \ne 1-\alpha$.
\item If $1-\alpha' < 1 - \alpha$, decrease $\alpha^*$ and return to the previous step.
If $1-\alpha' > 1 - \alpha$, increase $\alpha^*$ and return to the previous step.
If $1-\alpha' \approx 1-\alpha$ to the desired level of precision, go to the next
step.
\item
Report as a $1-\alpha$ confidence interval for $T(F)$ the (first generation)
bootstrap quantile confidence interval
that has nominal $1 - \alpha^*$ coverage probability.
\end{itemize}

An alternative approach to increasing coverage probability by iterating
the bootstrap is to use the same root, but to use a quantile
(among 2nd generation bootstrap samples) of
its $1-\alpha$ quantile rather than  the quantile observed in the first generation.
The heuristic justification is that we would ideally like to know the $1-\alpha$ quantile
of the pivot under sampling from the true distribution $F$.
We don't.  The percentile method estimates the $1-\alpha$ quantile of the
pivot under $F$ by the
$1-\alpha$ quantile of the pivot under $\Fhatn$, but this is subject to sampling variability.
To try to be conservative, we could use the bootstrap a second time  find an (approximate)
upper
$1-\alpha^*$ confidence interval for the $1-\alpha$ quantile of the pivot.

Here is a sketch of the Monte-Carlo approach:
\begin{itemize}
\item Pick a value $\alpha^* \in (0, 1/2)$ ({\em e.g.}, $\alpha^* = \alpha$).  This is a tuning
parameter.
\item From the sample, draw $M$ size-$n$ samples that are each iid
$\hat{F}_n$. Denote the ecdfs of the samples by $\{ \hat{F}_{n,j}^*\}$.
\item For each $j = 1, \ldots, M$, draw $N$ size-$n$ samples, each iid
$\hat{F}_{n,j}$. Find the $1-\alpha$ quantile of the pivot.
This gives $M$ values of the $1-\alpha$ quantile.
Let $c$ be the $1-\alpha^*$ quantile of the $1-\alpha$ quantiles.
\item
Report as a $1-\alpha$ confidence interval for  $T(F)$ the interval one gets
by taking $c$ to be the estimate of the $1-\alpha$ quantile of the pivot.
\end{itemize}
In a variety of simulations, this tends to be more conservative than Beran's method, and more
often attains at least the nominal coverage probability.

\noindent
{\bf Exercise.}
Consider forming a two-sided 95\% confidence interval for the mean $\theta$ of
a distribution $F$ based on the sample mean,
using $| \bar{X} - \theta |$ as a pivot.
Implement the three ``double-bootstrap'' approaches to finding a confidence interval
(Beran's pre-pivoting, Efron's  calibrated target percentile, and the
percentile-of-percentile).
Generate 100 synthetic samples of size 100 from the following distributions: normal,
lognormal, Cauchy,
mixtures of normals with the same mean but quite different variances (try different
mixture coefficients), and mixtures of normals with different means and different variances
(the means should differ enough that the result is bimodal).
Apply the three double bootstrap methods to each, resampling 1000 times from each of
1000 first-generation bootstrap samples.
Which method on the average has the lowest level error? Which method tends to be most
conservative?  Try to provide some intuition about the circumstances under
which each method fails, and the circumstances under which each method would be expected
to perform well.
How do you interpret coverage for the Cauchy?
{\bf Warning:} You might need to be clever in how you implement this to make it
a feasible calculation in S or Matlab.  If you try to store all the intermediate results,
the memory requirement is huge. On the other hand, if you use too many loops, the
execution time will be long.

\subsection{Bootstrap confidence sets based on Stein (shrinkage) estimates}
Beran (1995) discusses finding a confidence region for the mean vector $\theta \in \bfR^q$,
$q \ge 3$,
from data $X \sim N(\theta, I)$.
This is an example illustrating that {\em what} one bootstraps is important, and that
naive plug-in bootstrapping doesn't always work.

The sets are spheres centered at the shrinkage estimate
\beq
\hat{\theta}_S = \left ( 1 - \frac{q-2}{\|X\|^2} \right ) X,
\eeq
with random diameter $\hat{d}$.
That is, the confidence sets $C$ are of the form
\beq
C(\hat{\theta}_S, \hat{d}) =
\left \{ \gamma \in \bfR^q : \| \hat{\theta}_S - \gamma \| \le \hat{d}
\right \}.
\eeq
The problem is how to find $\hat{d} = \hat{d}(X; \alpha)$ such that
\beq
\Prob_\gamma \{ C(\hat{\theta}_S, \hat{d})  \ni \gamma \} \ge 1-\alpha
\eeq
whatever be $\gamma \in \bfR^q$.

This problem is parametric: $F$ is known up to the $q$-dimensional mean vector
$\theta$.
We can thus use a ``parametric bootstrap'' to generate data that are approximately from
$F$, instead of drawing directly
from $\hat{F}_n$: if we have an estimate $\hat{\theta}$ of $\theta$,
we can generate artificial data
distributed as $N( \hat{\theta}, I)$.
If $\hat{\theta}$ is a good estimator, the artificial data will  be distributed nearly
as $F$.  The issue is in what sense $\hat{\theta}$ needs to be good.

Beran shows (somewhat surprisingly) that resampling  from $N(\hat{\theta}_S,I)$
or from $N(X,I)$
do not tend to work well in calibrating $\hat{d}$.
The crucial  thing in using the bootstrap to calibrate the radius of the
confidence sphere seems to be to estimate $\| \theta \|$ well.

\begin{Definition}
The {\em geometrical risk} of a confidence set $C$ for the parameter $\theta \in \bfR^q$
is
\beq
G_q(C, \theta) \equiv q^{-1/2} E_\theta \sup_{\gamma \in C} \| \gamma - \theta \|.
\eeq
That is, the geometrical risk is the expected distance to the parameter from the
most distant point in the confidence set.
\end{Definition}

For confidence spheres
\beq
C = C(\hat{\theta}, \hat{d}) = \{ \gamma \in \bfR^q : \| \gamma - \hat{\theta} \| \le
\hat{r} \},
\eeq
the geometrical risk can be decomposed further: the distance from $\theta$
to the most distant
point in the confidence set is the distance from $\theta$ to the center of the sphere,
plus the radius of the sphere, so
\begin{eqnarray}
G_q(C(\hat{\theta}, \hat{r}), \theta) &=&
q^{-1/2} E_\theta \left (  \| \hat{\theta} - \theta \| + \hat{d} \right )
\nonumber \\
&=&
q^{-1/2} E_\theta \| \hat{\theta} - \theta \| + q^{-1/2} E_\theta \hat{d} .
\end{eqnarray}


\begin{Lemma}  \label{lemma:BeranLemma41}
(Beran, 1995, Lemma 4.1).
Define
\beq
W_q(X, \gamma) \equiv (q^{-1/2} ( \|X - \gamma \|^2 - q ), q^{-1/2} \gamma'(X - \gamma).
\eeq
Suppose $\{ \gamma_q \in \bfR^q \}$ is any sequence such that
\beq \label{eq:gammaqCond}
\frac{\| \gamma_q \|^2}{q}  \rightarrow a < \infty \mbox{ as } q \rightarrow \infty .
\eeq
Then
\beq
W_q(X, \gamma_q) \rightWarrow (\sqrt{2} Z_1, \sqrt{a} Z_2 )
\eeq
under $\Prob_{\gamma_q}$, where $Z_1$ and $Z_2$ are iid standard normal random variables.
(The symbol $\rightWarrow$ denotes weak convergence of distributions.)
\end{Lemma}

\noindent
{\bf Proof.}
Under $ \Prob_{\gamma_q}$, the distribution of $X - \gamma$ is rotationally invariant,
so the distribution of the components of $W_q$ depend on $\gamma$ only through
$\| \gamma \|$. Wlog, we may take each component of $\gamma_q$ to be
$q^{-1/2}\| \gamma_q\|$.
The distribution of the first component of $W_q$ is then that of the sum of squares
of $q$ iid standard normals (a chi-square rv with $q$ df),
minus the expected value of that sum, times $q^{-1/2}$.
The standard deviation of a chi-square random variable with $q$ df is $\sqrt{2q}$,
so the first component of $W_q$ is $\sqrt{2}$ times a standardized variable whose
distribution is asymptotically (in $q$) normal.
The second component of $W_q$ is a linear combination of iid standard normals; by
symmetry (as argued above), its distribution is that of
\begin{eqnarray}
q^{-1/2} \sum_{j=1}^q q^{-1/2}\| \gamma_q\| Z_j &=&
\| \gamma_q \|\sum_{j=1}^q Z_j
\nonumber \\
\rightarrow a^{1/2} Z_2.
\end{eqnarray}

Recall that the squared-error risk (normalized by $q^{-1/2}$)
of the James-Stein estimator is
$1 - q^{-1} E_\theta \{ (q-2)^2/\|X\|^2 \} < 1$.
The difference between the loss of $\hat{\theta}_S$ and an unbiased estimate of
its risk is
\beq
D_q(X, \theta) = q^{-1/2} \{ \| \hat{\theta}_S - \theta \|^2 -
[q - (q-2)^2/\|X\|^2] \}.
\eeq
By rotational invariance, the distribution of this quantity depends on $\theta$ only
through $\| \theta\|$; Beran writes the distribution as $H_q(\| \theta \|^2/q)$.
Beran shows that if  $\{ \gamma_q \in \bfR^q \}$ satisfies \ref{eq:gammaqCond},
then
\beq
H_q(\|\gamma_q\|^2/q) \rightWarrow N(0, \sigma^2(a)),
\eeq
where
\beq
\sigma^2(t) \equiv 2 - 4t/(1+t)^2 \ge 1.
\eeq
Define
\beq
\hat{\theta}_{\mbox{CL}} = [ 1 - (q-2)/\|X\|^2]_+^{1/2} X.
\eeq


\begin{Theorem}
(Beran, 1995, Theorem 3.1)
Suppose   $\{ \gamma_q \in \bfR^q \}$ satisfies    \ref{eq:gammaqCond}.
Then
\beq
H_q( \|\hat{\theta}_{\mbox{CL}}\|^2/q)   \rightWarrow N(0, \sigma^2(a))   ,
\eeq
\beq
H_q(\|X\|^2/q) \rightWarrow N(0, \sigma^2(1+a)),
\eeq
and
\beq
H_q(\| \hat{\theta}_S\|^2/q )\rightarrow N(0, \sigma^2(a^2/(1+a))),
\eeq
all in $P_{\gamma_q}$ probability.
\end{Theorem}

It follows that to estimate $H_q$ by the bootstrap consistently,
one should use
\beq
\hat{H}_B = H_q( \|\hat{\theta}_{\mbox{CL}}\|^2/q       )
\eeq
rather than estimating using either the norm of $X$ or the norm of the
James-Stein estimate $\hat{\theta}_{S}$ of $\theta$.

\noindent
{\bf Proof.}
Lemma \ref{lemma:BeranLemma41} implies that under the conditions of the theorem,
$\| \hat{\theta}_{\mbox{CL}}\|^2/q \rightarrow a$,
$\|X\|^2/q \rightarrow 1+a$, and $\| \hat{\theta}_S \|^2 /q \rightarrow a^2/(1+a)$.






\end{document}


