Bootstrap

Vera

1 Introduction

The estimation of parameters in probability distributions is a basic problem in statistics that one tends to encounter already during the very first course on the subject. Along with the estimate itself we are usually interested in its accuracy, which is can be described in terms of the bias and the variance of the estimate, as well as confidence intervals around it. Sometimes, such measures of accuracy can be derived analytically. Often, they can not. Bootstrap is a technique that can be used to try and estimate them numerically from a single data set.

2 The general idea

Let X_{1}, ..., X_{n} be a i.i.d. sample from distribution F, that is P(X_{i} < x) = F (x), and let X_{(1)}, ..., X_{(n)} be the corresponding ordered sample. We are interested in some parameter θ which is associated with the distribution (mean, median, variance etc). There is also an estimator θ̂ = t({X_{1}, ..., X_{n}}), with t denoting some function, that we can use to estimate θ from data. It is the deviation of θ̂ from θ that is of primarily interest, or, more precisely, the distribution of θ̂.

Ideally, to get an approximation of the estimator distribution we would like to repeat the data-generating experiment, say, N times, calculating θ̂ for each of the N data sets. That is, we would draw N samples of size n from the true distribution F. In practice, this is impossible.

The bootstrap method is based on the following simple idea: Even if we do not know F we can approximate it from data and use this approximation, F̂, instead of F itself.

This idea leads to several flavours of bootstrap that deviate in how, exactly, the approximation F̂ is obtained. Two broad areas are the parametric and non-parametric bootstrap.

The non-parametric estimate is the so called empirical distribution (you will see the corresponding pdf if you simply do a histogram of the data), that can be formally defined as follows:

With # denoting the number of members of a set, the empirical distribution function F̃ is given by

F̃(x) = \frac{#\{i : X_{i} ≤ x\}}{n} = \{ \begin{matrix}0 & if & x ∈ (−∞, X_{(1)}), \\ i/n & if & x ∈ [X_{(i)}, X_{(i+1)}) & for i ∈ \{1, …, n − 1\}, \\ 1 & if & x ∈ [X_{(n)}, ∞). \end{matrix}

That is, it is a discrete distribution that puts mass 1/n on each data point in your sample. It can be shown that F̃ converges to F as n → ∞.

The parametric estimate assumes that the data comes from a certain distribution family (Normal, Gamma etc). That is, we say that we know the general functional form of the pdf, but not the exact parameters. Those parameters can then be estimated from the data (typically with Maximum Likelihood) and plugged in the pdf to get F̂. This estimation method leads to more accurate inference if we guessed the distribution family correctly but, on the other hand, F̂ will not approach F if the family assumption is wrong.

2.1 Algorithms

The two algorithms below implement the above. To concretize, let us say that X_{i} ∼ N(0, 1), is the median and it is estimated by θ̂ = X_{(n/2)}.

Non-parametric bootstrap

Assuming a data set x = (x_{1}, ..., x_{n}) is available.

1. Fix the number of bootstrap re-samples N. Often N ∈ [1000, 2000].

2. Sample a new data set x^{∗} set of size n from x with replacement (this is equivalent to sampling from the empirical cdf F̂).

3. Estimate θ from x^{∗}. Call the estimate it θ̂_{i}^{∗}. Store.

4. Repeat 2 and 3 N times.

5. Consider the empirical distribution of (θ̂_{1}^{∗}, ..., θ̂_{N}^{∗}) as an approximate of the true cdf of θ̂.

Parametric bootstrap

Assuming a data set x = (x_{1}, ..., x_{n}) is available.

1. Assume that the data comes from a known distribution family F_{ψ} described by a set of parameters ψ (for a Normal distribution ψ = (μ, σ) with μ being the expected value and σ the standard deviation).

2. Estimate ψ with, for example, Maximum likelihood, obtaining the estimate ψ̂.

3. Fix the number of bootstrap re-samples N. Often N ∈ [1000, 2000].

4. Sample a new data set x^{∗} set of size n from Fψ̂.

5. Estimate θ from x^{∗}. Call the estimate it θ̂_{i}^{∗}. Store.

6. Repeat 4 and 5 N times.

7. Consider the empirical distribution of (θ̂_{1}^{∗}, ..., θ̂_{N}^{∗}) as an approximate of the true cdf of θ̂.

In the Figure 1 the distribution of θ̂ (the median) approximated with the non-parametric and parametric bootstrap is plotted. Note that the parametric distribution is smoother than the non-parametric one, since the samples were drawn from a continuous distribution. As stated earlier, using the parametric version could be more effective, but also potentially more difficult to set up. In particular, we have to sample from Fψ̂ which is not easily done if F is not one of the common distributions (such as Normal) for which a pre-programmed sample function is usually available.

Figure 1: The data distribution. The non-parametric sample. The parametric sample.

3 Bias and variance estimation

The theoretical bias and variance of an estimator θ̂ are defined as

Bias(θ̂) = E[θ̂ − θ] = E[θ̂] − θ
Var(θ̂) = E [(θ̂ − E[θ̂])^{2}] ≈ E[(θ̂ − θ)^{2}]

with the approximation in the variance definition being valid if the bias is small. In words, the bias is a measure of a systematic error (θ̂ tends to be either smaller of larger than θ) while the variance is a measure of random error.

In order to obtain the bootstrap estimates of bias and variance we plug in the original estimate t̂heta (which is a constant given data) in place of θ and θ̂^{∗} (the distribution of which we get from bootstrap) in place of θ̂. This leads us to the following approximations:

Bias(θ̂) ≈ \frac{1}{N} ∑_{i} θ̂_{i}^{∗} − θ̂ = θ̂_{.}^{∗} − θ̂

Var(θ̂) ≈ \frac{1}{N − 1} ∑_{i}(θ̂_{i}^{∗} − θ̂_{.}^{∗})^{2}

That is, the variance is, as usual, estimated by the sample variance (only for the bootstrap sample of θ̂) and bias is estimated by how much the original θ̂ deviates from the average of the bootstrap sample.

4 Confidence intervals

There are several methods for CI construction with Bootstrap, the most popular being ”normal”, ”basic” and ”percentile”. Let θ̂_{(i)}^{∗} be the ordered bootstrap estimates, with i = 1, ..., N indicating the different samples. Then the two-sided CIs at the significance level α are defined as

Basic  [2θ̂ − θ̂_{(N+1)(1−α/2)}^{∗}, 2θ̂ − θ̂_{(N+1)α/2}^{∗}]
Normal  [θ̂ − z_{α/2}ŝe, θ̂ + z_{α/2}ŝe]
Percentile  [θ̂_{(N+1)α/2}^{∗}, θ̂_{(N+1)(1−α/2)}^{∗}]

with z_{α} denoting an α quantile from a Normal distribution and ŝe the estimated standard deviation of θ̂ calculated from the bootstrap sample.

Basic CI

To obtain this confidence interval we start with θ̂ − θ and, like in the variance and bias calculations, approximate with θ̂^{∗} − θ̂. Let us call the p-quantile of the distribution of θ̂^{∗} − θ̂ for q_{p}. We then have that an α-level CI for θ̂ − θ can be approximated with

q_{α/2} ≤ θ̂ − θ ≤ q_{1−α/2} ⇔ θ̂ − q_{1−α/2} ≤ θ ≤ θ̂ − q_{α/2}

We will get the formula for the Basic CI by noting that q_{p} = θ̂_{p}^{∗} − θ̂ and substituting it in the expression above.

This interval construction relies on an assumption about the distribution of θ̂−θ, namely that it is independent of θ. This assumption, called pivotality, does not necessarily hold in most cases. However, the interval gives acceptable results even if θ̂ − θ is close to pivotal, that is its dependence on θ is weak.

Normal CI

This confidence interval probably looks familiar since it is almost an exact replica of the commonly used confidence interval for a population mean. Indeed, similarly to that familiar CI, we can get it by assuming that Z = (θ̂−θ)/ŝe ∼ N(0, 1). Observe that this also implicitly assumes pivotality.

Percentile CI

May seem natural and obvious but actually requires a quite convoluted argument to motivate. Without going into details, you get this CI by assuming that there exist a transformation h such that the distribution of h(θ̂) − h(θ) is pivotal, symmetric and centered around 0. You then construct a Basic CI for h(θ) rather than θ itself, get the α quantiles and transform those back to the original scale by applying h^{−1}. Observe that although we do not heed to know h explicitly, the existence of such a transformation is a must, which is not always the case.

5 The restrictions of bootstrap

Yes, those do exist. Bootstrap tends to give results easily (too much so, in fact), but it is possible that those results are completely wrong. More than that, they can be completely wrong without being obvious about it. Let take a look at a few classics.

5.1 Infinite variance

In the ”general idea” of bootstrapping we plugged in an estimate of the density F, F̂, and then sampled from it. This works only if F̂ actually is a good estimate of F, that is it captures the essential features of F despite being based on a finite sample. This may not work well if F is very heavy tailed, e. g. has infinite variance. The intuition is that in this case extremely large or small values can occur, and when they do they have a great effect on the bootstrap estimate of the distribution of θ̂, making it unstable. As a consequence, the measures of accuracy of θ̂, such as CI, will be unreliable.

The classic example of this is the mean estimator θ̂ = X̄ with the data generated from a Cauchy distribution. For it, both the first and the second moments (e.g. mean and variance) are infinite, leading to nonsensical confidence intervals even for large sample sizes. A less obvious example is a non-central Student t distribution with 2 degrees of freedom. This distribution has pdf

f_{X}(x) = \frac{1}{2 (1 + (x − μ)^{2})^{3/2}} for x ∈ ℝ.

where μ is the location parameter, defined as E[X]. So, the first moment is finite and its estimator X̄ is consistent. The second moment, however, is infinite, and the right tail of the distribution grows heavier with increasing μ. This leads to the approximate 95% CI coverage probabilities displayed in Figure 2 when the ”basic” and ”percentile” methods are used on a non-parametric bootstrap sample (the ”Normal” method is, of course, inapplicable due to the infinite variance). The actual coverage is quite far from the the supposed 95% even when the sample size is as large as 500.

Figure 2: Left: The approximate coverage probabilities for the ”basic” (triangles) and ”percentile” (circles) CI when μ = 10. Right: histogram of 1000 non-parametric bootstrap samples of θ̂^{∗}, the data from Uni(0, θ) distribution. θ = 1.

This potential problem can be circumvented with the sub-sampling methods, such as Jack-Knife. This method, however, is not without its own issues, especially clear if θ̂ is not affected much by sub-sampling, for example if it is a median.

5.2 Parameter on the boundary

The classical example here is the Uni(0, θ) distribution. The maximum likelihood estimate θ̂ is simply max(x_{1}, ..., x_{n}), which will always be biased (θ̂ < θ). In this case the non-parametric bootstrap leads to a very discrete distribution, with more than half of the bootstrap estimates θ̂^{∗} equal to the original θ̂ (Figure 2). Clearly, if the quantiles used in CIs are taken from this distribution the results will be far from accurate. However, the parametric bootstrap will give a much smoother distribution and more reliable results.

5.3 Lack of pivotality

In all the CI descriptions above the word ”pivotality” shows up. So we can guess that it is a bad thing not to have (see Lab, task 3). To circumvent this, something called ”studentized bootstrap” can be used.

The idea behind the method is simple and can be seen as an extrapolation of the Normal bootstrap CI. There, we assumed that the distribution of the standardized test statistic Z is known, namely that it is N(0, 1). In studentized bootstrap we instead say that the distribution of Z is not known, but can be approximated from data. Through, again, bootstrap.

So, we proceed in the similar way as before and approximate Z with Z^{∗} = (θ̂^{∗} − θ̂)/ŝe^{∗}, where ŝe^{∗} is the estimate of the standard deviation of θ̂^{∗}. This quantity, Z^{∗}, is calculated for each bootstrap sample, so that N values are obtained. Those values can then be viewed as a sample from the distribution of Z. They are ordered to get the ordered sample (Z_{(1)}^{∗}, ..., Z_{(N)}^{∗}), and we use Z_{(N+1)α/2}^{∗} and Z_{(N+1)(1−α/2)}^{∗} as quantiles to put in the formula for the Normal CI instead of z_{α/2}. That is, the studentized CI is defined as

θ̂ − Z_{(N+1)α/2}^{∗}ŝe ≤ θ ≤ θ̂ − Z_{(N+1)(1−α/2)}^{∗}ŝe

This CI, unlike the three previous ones, is generally considered to be very reliable. However, there is a catch, and this catch is the estimate of the standard deviation of θ̂^{∗}, ŝe^{∗}. This estimate is not easily obtained. You can get it either parametrically (but then you need a model which you probably don’t have) or through re-sampling. This means that you have to do an extra bootstrap for each of the original bootstrap samples. That is, you will have a loop within a loop, and it can become very heavy computationally.

6 software

Will be done in R. Functions you may find useful:

sample, replicate, rgamma, boot, boot.ci, dgamma, nlm

7 Laboration: Bootstrap re-sampling.

The aim of this laboration is to study the performance of bootstrap estimates, as well as corresponding variance, bias and CI, using different bootstrap techniques.

7.1 Tasks

1. Non-parametric bootstrap. We start with the easier task of setting up a nonparametric bootstrap. First, generate some data to work with. Let us say that the data comes from a definitely non-normal, highly skewed distribution. Such as, for example, Gamma.

• Generate a data set of size n = 100 from a Gamma(k, 1/γ) distribution, where k = 1 is the shape parameter and γ the scale parameter. Do a histogram of the data. What are the theoretical mean and variance of a Gamma distribution? Observe that Gamma can be defined both through γ and 1/γ. Look out for which definition R uses.

• Let us say that our parameter of interest θ is the mean. Write a function that approximates the distribution of θ with non-parametric bootstrap. Present the distribution in a histogram. What is the estimated bias of θ̂? Variance?

• As with many other things nowadays, there already exist a package that does bootstrap for you. In R this package is boot. Use it to get the distribution of the mean and the estimates of the bias and variance.

• Construct CI for θ (normal, basic, percentile). Check your results with the preprogrammed function boot.ci.

• Bootstrap is known to give less than reliable results if the sample size n is small. Study this by doing the following. Fix an n = 10. Draw 1000 samples of this size. For each sample, calculate bootstrap CI, recording whether they cover the true value of θ (the theoretical mean, that is). Repeat this for n = 20, 30, ..., 100. Do the CI have the correct coverage probabilities? Present results in a graph. You can use either your own code or boot.ci.

2. Parametric bootstrap. (2 p) Repeat the previous task using parametric bootstrap instead of non-parametric. In order to do this you have to estimate the k and γ parameters in the Gamma distribution through, for example, Maximum Likelihood. A nice, closed form ML solution for the parameters does not exist, so you will have to write a function that finds the maximum of the Likelihood function L numerically.

Tip 1: most optimization packages look for minimum and not maximum, so work with -log(L) instead of L itself.

Tip 2: Look up nlm function.

3. Studentized CI. (2 p) Our parameter of interest θ is no longer the mean, it is the variance.

• Repeat the last question in task 1. What can be said about the performance of the standard confidence intervals?

• Construct Studentized confidence intervals for the variance. Check the coverage as in the above. Any improvement?

You can construct the studentized CI both from scratch with your own code and using boot and boot.ci. If you choose to do the latter, make sure that the ”statistic” function that you have to specify in boot returns two values: θ̂^{∗} and Var(θ̂^{∗}).

In order to pass the lab task 1 has to be completed.