Confidence intervals by bootstrap

This package demonstrates how to calculate bootstrap confidence intervals (BCa and extreme percentile) for univariate continuous data.

Author:  Sergej V. Aksenov

Department of Microbiology and Immunology
University of Michigan
Ann Arbor, MI  48109, USA
http://aksenov.freeshell.org/

Last updated 29 April 2002

Installation and initialization

The package file Bootstrap.m should be located inside the Bootstrap directory.  
The Bootstrap directory should be put in one of the canonical places for AddOns.   
Two recommended places are first inside the Mathematica distribution in:

ToFileName[ {$TopDirectory,  "AddOns",  "Applications"}]

a second in the preferences directory:

\(ToFileName[\ {$PreferencesDirectory, \ "AddOns", \ "Applications"}]\)

When one of these locations is used there is no need to set $Path and the loading instructions are identical for all versions of Mathematica.

Important note: 

this package requires package RootSearch by Ted Ersek to be installed in its default location: create directory Enhancements in one of the above canonical locations and put file RootSearch.m there.   
Please download RootSearch.m from MathSource at
  http://www.mathsource.com/cgi-bin/msitem?0211-734

The  Bootstrap` package can  be loaded with

<<Bootstrap`Bootstrap`

Introduction

Suppose that we have data F→x=((Subscript[x, 1]), …, (Subscript[x, n])) that are obtained by random sampling from an unknown distribution.  Suppose further that there is a parameter θ=t(F) that we are interested in estimating.  A point estimate of this parameter is obtained by using e.g. empirical distribution function based on data, \(θ\&^\)=t(\(F\&^\)).

A point estimate is usually supplemented by a confidence interval \(TraditionalForm\`I = \((\(θ\_lo\)\&^, \(θ\_up\)\&^)\)\) with approximate coverage 1-2α that is saying that the probability of finding the true value of the parameter inside this interval is exactly 1-2α.  This can be equivalently formulated as Prob\(TraditionalForm\`{θ < \(θ\_lo\)\&^} = α\) and Prob\(TraditionalForm\`{θ < \(θ\_lo\)\&^} = \(\(α\)\(.\)\)\)

Bootstrap is used to generate a distribution of \(TraditionalForm\`\(θ\^*\)\&^\) that approximates the distribution of the true parameter θ by taking B resamples of the data.

Bootstrap percentile intervals are constructed by taking αth and 1-αth percentiles of the bootstrap distribution of  \(TraditionalForm\`\(θ\^*\)\&^\): \(TraditionalForm\`I = \((θ\&^\^\(\(*\)\((α)\)\), θ\&^\^\(\(*\)\((1 - α)\)\))\)\), where \(TraditionalForm\`θ\&^\^\(\(*\)\((α)\)\)\)is the Bαth element of the ordered list of B replications of \(TraditionalForm\`\(θ\^*\)\&^\).  Likewise one determines \(TraditionalForm\`θ\&^\^\(\(*\)\((1 - α)\)\)\).

In this package, one can calculate two important modifications of the bootstrap percentile intervals: bias corrected and accelerated (BCa) and extreme percentiles intervals.

BCa interval is given by \(TraditionalForm\`I\_BCa = \((θ\&^\^\(\(*\)\((α\_1)\)\), θ\&^\^\(\(*\)\((α\_2)\)\))\)\) where the modified percentile points are calculated as

\(α\_1 = Φ \((z\&^\_0 + \(z\&^\_0 + z\^\((α)\)\)\/\(1 - \(a\&^\) \((z\&^\_0 + z\^\((α)\))\)\))\)\)
\(α\_2 = Φ \((z\&^\_0 + \(z\&^\_0 + z\^\((1 - α)\)\)\/\(1 - \(a\&^\) \((z\&^\_0 + z\^\((1 - α)\))\)\))\)\)

Here Φ is the standard normal cumulative distribution function and (z^((α))) is 100αth percentile point of a standard normal distribution.

Bias-correction constant \(TraditionalForm\`z\&^\_0\) is obtained directly from the proportion of bootstrap replications less or equal to the original estimate \(TraditionalForm\`θ\&^\), and the acceleration constant can be most conveniently estimated with jackknife values of the statistic \(TraditionalForm\`θ\&^\_\((i)\)\)(obtained with the original sample with the ith point deleted):

\(z\&^\_0 = \(Φ\^\(-1\)\) \((\(# \((\(θ\&^\^*\) < θ\&^)\)\)\/B)\)\)
\(a\&^ = \(∑\+\(i = 1\)\%n\((θ\&^\_\((·)\) - θ\&^\_\((i)\))\)\^3\)\/\(6 {∑\+\(i = 1\)\%n\((θ\&^\_\((·)\) - θ\&^\_\((i)\))\)\^2}\^\(3/2\)\)\)

where \(θ\&^\_\((·)\) = ∑\+\(i = 1\)\%n θ\&^\_\((i)\)/n\).

Extreme percentiles intervals \(TraditionalForm\`I\_EX = \((\(\(θ\&^\)\_\((1)\)\%*\), \(\(θ\&^\)\_\((B)\)\%*\))\)\)are obtained by first and last elements in the ordered list of bootstrap replications, where the size of bootstrap effort B is determined by solving extreme coverage equations for the optimal B:

\(\((1)\)\ \ \ \ Prob {θ ≤ \(\(θ\&^\)\_\((B)\)\%*\)} = 
      1 - \((B + 1)\)\^\(-1\) - \((1/6)\) \((b\^3/B)\)\ \(A\_1/
              n\^\(1/2\)\)/σ\^3\)
\(\((2)\)\ \ \ \ Prob {θ ≥ \(\(θ\&^\)\_\((1)\)\%*\)} = 
      1 - \((B + 1)\)\^\(-1\) + \((1/6)\) \((b\^3/B)\)\ \(A\_1/
              n\^\(1/2\)\)/σ\^3\)
\(\((3)\)\ \ \ \ Prob {\(\(θ\&^\)\_\((1)\)\%*\) ≤ θ ≤ \(\(θ\&^\)\_\((B)\)\%*\)} = 
      1 - 2 \((B + 1)\)\^\(-1\) - \((1/36)\) \((b\^6/B)\)\ \(A\_1\%2/
              n\)/σ\^6\)

where b is the maximum positive solution of

(4)    B ϕ (b - (b^(-1))) = b

Here, Eqs. (1), (2), and (3) are solved for one-sided upper, one-sided lower, and two-sided intervals, respectively.  To obtain an equi-tailed two-sided interval of coverage α, one has to solve two equations (1) and (2) equating each to (1+α)/2 and then follow the algorithm described in Lee's paper (the purpose of that special calculation is to reduce the variability of endpoints).

References

B. Efron, R.J. Tibshirani, "An introduction to the Bootstrap", Chapman and Hall, Inc., New York, London, 1993.

S.M.S. Lee, "Nonparametric confidence intervals based on extreme bootstrap percentiles", Statistica Sinica, 10: 475, 2000.

Examples

Example 1

In Example 1 we calculate two-sided BCa percentile interval for the sample mean.

In general, a statistic (function of a random variable) to be bootstrapped should be defined as a Mathematica function.  In the case of the mean, however, there exists a standard Mathematica function Mean that can be given as an argument to bootstrapping functions.

Generate a random sample from the standard normal distribution (first fix random seed to reproduce the sample later).

SeedRandom[228930]

data1 = RandomArray[NormalDistribution[0, 1], 30]

Here's the observed value of the statistic (sample mean):

Mean[data1]

Here are the jackknife replications of the statistic:

?JackknifeReplications

jack1 = JackknifeReplications[Mean, data1]

Here are the nonparametric bootstrap replications of the statistic:

?BootstrapReplications

SeedRandom[338901]

boot1 = BootstrapReplications[Mean, 2000, data1];

Here are the BCa confidence intervals:

?BootstrapCI

BootstrapCI[Mean[data1], boot1[[1]], jack1, 
      BootstrapCIMethod -> BCaPercentiles]

Here are the variances of endpoints of the above BCa interval:

?VarianceNonParametricBCaCIEndpoints

VarianceNonParametricBCaCIEndpoints[Mean[data1], boot1[[1]], jack1, 
      boot1[[2]]]

Standard deviations..

Sqrt[%]

Now let's calculate confidence interval for the mean based on the normal theory using standard Mathematica package Statistics`ConfidenceIntervals`

<< Statistics`ConfidenceIntervals`

Notice that the standard interval is a little wider than the bootstrap interval..
However, given rather large variances of the endpoints, the results are equivalent.

MeanCI[data1]

Example 2

In Example 2 we calculate minimum number of bootstrap resamples needed for the extreme percentile interval for the sample variance.

First we generate a random sample from the standard normal distribution..

SeedRandom[932648]

data2 = RandomArray[NormalDistribution[0, 1], 20];

We will need jackknife replications..

jack2 = JackknifeReplications[Variance, data2]

Now we calculate the minimum numbers of resamples for the one-sided upper 0.99 confidence limit for the sample variance.

?BootstrapResamples

BootstrapResamples[Variance[data2], jack2, 
      ConfidenceLevel -> 0.99, 
      ConfidenceInterval -> OneSidedLower]

Now we calculate minimum numbers for different type of intervals and different coverage levels..
These results can be compared with the Table 1 in S.M.S Lee's paper in Stat. Sinica, 10 (2000), 475.
The discrepancies are due to estimation of A_1 and sigma by jackknife in our function, whereas Lee used analytical expressions.

Please note the following about the resulting table:

1) columns correspond to one-sided upper, one-sided lower, two-sided asymmetrical, two-sided equi-tailed intervals;

2) rows correspond to confidence levels 0.8,0.85,0.9,0.925,0.95,0.975,0.99 that refer to two-sided coverage;  therefore, for comparison, to obtain one-sided and two-sided equitailed intervals of coverage, that would be equivalent to two-sided coverage α, one has to specify ConfidenceLevel -> (1+α)/2.

Map[
        Map[Apply[
                BootstrapResamples[Variance[data2], jack2, 
                    ConfidenceLevel -> #1, 
                    ConfidenceInterval -> #2] &, #] &, {{(1 + #)/2, 
                OneSidedUpper}, {(1 + #)/2, OneSidedLower}, {#, 
                TwoSidedAsymmetrical}, {(1 + #)/2, 
                TwoSidedEquiTailed}}] &, {0.8, 0.85, 0.9, 0.925, 0.95, 0.975, 
          0.99}] // TableForm

Example 3

In Example 3 we calculate extreme percentile interval for a goodness-of-fit statistic.

Generate sample data from standard normal distribution..

SeedRandom[390013]

data3 = RandomArray[NormalDistribution[0, 1], 50];

CDF of the normal distribution that we will fit to the sample..

ncdf[xlist_, mu_, sigma_] := 
      CDF[NormalDistribution[mu, sigma], xlist]

Attributes[ncdf] = {Listable}

Root mean square error is our goodness-of-fit statistic..

rmserror[m_, s_] := 
Sqrt[
        Apply[Plus, (ylist - ncdf[xlist, m, s])^2]]

This function takes in a list of random data, fits the distribution, and outputs the value of a goodness-of-fit statistic.
It is defined in a form suitable for input into functions of the Bootstrap package.

fun[data_] := (xlist = Transpose[QuantileForm[data]][[2]]; 
        ylist = Transpose[QuantileForm[data]][[1]] // N; 
        fitres = 
          FindMinimum[rmserror[m, s], {m, {0.0, 0.1}}, {s, {1.0, 1.1}}]; 
        fitres[[1]])

We will need observed value and jackknife replications..

obs = fun[data3]

jack3 = JackknifeReplications[fun, data3]

Calculate minimum number of resamples needed for construction of one-sided upper 0.95 extreme percentile confidence interval.
Notice that this number is much smaller than 1000-2000 usually recommended for BCa intervals..

BootstrapResamples[obs, jack3, 
      ConfidenceInterval -> OneSidedUpper]

Calculate bootstrap replications of the goodness-of-fit statistic..

SeedRandom[590246]

boot3 = BootstrapReplications[fun, 17, data3];

Calculate the upper confidence limit..
The observed rms value is well below the 0.95 upper confidence limit so the hypothesis that fitted distribution generated the sample is not rejected at the significance level 0.05.

BootstrapCI[obs, boot3[[1]], jack3, 
      ConfidenceInterval -> OneSidedUpper]

Bootstrap options

The following lists what options are available for functions in the package:  Options and their values are explained with examples below.

Options[BootstrapResamples]

Options[BootstrapReplications]

Options[VarianceNonParametricBCaCIEndpoints]

Options[BootstrapCI]

Since functions BootstrapResamples and BootstrapCI use RootSearch, one can pass RootSearch options to them..  Please refer to the RootSearch documentation for more information on values of these options.

Options[RootSearch]

BootstrapCIMethod

The following explains option BootstrapCIMethod.

?BootstrapCIMethod

The following explains possible values for option BootstrapCIMethod

?BCaPercentiles

?ExtremePercentiles

The default value is

BootstrapCIMethod /. Options[BootstrapCI]

Uses of BCaPercentiles value was demonstrated in the Example 1 above; uses of ExtremePercentiles value were demonstrated in the Examples 2 and 3 above.

BootstrapSamplingMethod

The following explains option BootstrapSamplingMethod.

?BootstrapSamplingMethod

The following explains possible values for option BootstrapSamplingMethod

?NonParametric

?Parametric

The default value is

BootstrapSamplingMethod /. Options[BootstrapReplications]

Assume now that sample from Example 1 is generated (and stored in data1) and jackknife replications are calculated (and stored in jack1)..

Calculate parametric bootstrap resamples

SeedRandom[119082]

boot1par = 
        BootstrapReplications[Mean, 2000, NormalDistribution[0, 1], 30, 
          BootstrapSamplingMethod -> Parametric];

Notice that the second element of the list boot1par is empty; it is nonempty only in the nonparametric mode.

boot1par[[2]]

Now let's calculate BCa interval again..

BootstrapCI[Mean[data1], boot1par[[1]], jack1, 
      BootstrapCIMethod -> BCaPercentiles]

To appreciate the difference, let's look at some parameters of intervals.
Parametric intervals are shorter and more asymmetric..

Length of the parametric interval:

(Part[#, 2] - Part[#, 1]) &[{-0.2788466530017277`, 
        0.44859708391329184`}]

Length of the nonparametric interval:

(Part[#, 2] - Part[#, 1]) &[{-0.3407003507280921`, 
        0.39938215012236294`}]

Shape of the parametric interval:

(Part[#, 2] - Mean[data1])/(Mean[data1] - 
              Part[#, 1]) &[{-0.2788466530017277`, 
        0.44859708391329184`}]

Shape of the nonparametric interval:

(Part[#, 2] - Mean[data1])/(Mean[data1] - 
              Part[#, 1]) &[{-0.3407003507280921`, 
        0.39938215012236294`}]

ConfidenceInterval

The following explains option ConfidenceInterval.

?ConfidenceInterval

The following explains possible values for option ConfidenceInterval.

?OneSidedUpper

?OneSidedLower

?TwoSidedAsymmetrical

?TwoSidedEquiTailed

The default value is

ConfidenceInterval /. Options[BootstrapCI]

Uses of these value were demonstrated in Example 2 above.

ConfidenceLevel

Option ConfidenceLevel is a standard Mathematica option..

?ConfidenceLevel

Default value is

ConfidenceLevel /. Options[BootstrapCI]

Notes

How to specify option ConfidenceLevel for equi-tailed intervals

Suppose that alpha is the desired coverage of an interval.

Then for upper and lower one-sided and for the two-sided asymmetrical intervals one should use option

ConfidenceLevel -> alpha

However, for the two-sided equitailed interval, one should use option

ConfidenceLevel -> (1 + alpha)/2

How equations are solved in BootstrapResamples

In function BootstrapResamples, Ersek's function RootSearch is instructed to search for all roots of the equation (4) in the Introduction in the interval (0,5), because it was found empirically that b=5 corresponds to some millions of B, a practically not achievable number of resamples.

Moreover, only the larger one from the two positive roots of equation (4) in the Introduction is taken, because the asymtotic theory of extreme intervals assumes b→∞.

How to use options with user functions

You can use options with the function you want bootstrapped/jackknived; just pass them to functions BootstrapReplications and JackknifeReplications and they will be passed on to all your function calls inside.

Use the syntax below.

JackknifeReplications[<funname>, <datalist>, OptionName1 -> OptionValue1,  ... ]

Note however, that you should have a list of default options assigned to your function <funname>, e.g. with the built-in function Options, because Bootstrap.m functions use FilterOptions package and will use only those options that are assigned to the function <funname>.

For example, since <funname> should take only one argument, data list, and the calculation might involve setting other parameters, you might be forced to write a wrapper for your custom function and use options with that wrapper to pass all necessary parameters for your calculation.  In this case, don't forget to assign all the options that you intend to use to the wrapper function.

For example, if your calculation involves standard Mathematica function FindMinimum, use the assignment below.

Options[<funname>] = Options[FindMinimum]

How to implement computation failures in user functions

User functions that calculate statistics of random data and are used in functions BootstrapReplications and JackknifeReplications can fail to compute.  BootstrapReplications and JackknifeReplications can handle these situations but user functions has to be written appropriately.

A user function should evaluate its intrinsic error conditions and return symbol $Failed in such cases.

Function BootstrapReplications detects the return values $Failed and ignores the bootstrap resample that caused failure.  A new resample is drawn until the requested number of replications is accumulated.

Function JackknifeReplications detects the return value $Failed and aborts the calculation; its return value is a list {$Failed, <ind>}, where <ind> is the index of the jackknife resample that caused failure.  Index of the jackknife resample is simply the index of the point from the original (as supplied to JackknifeReplications) data list that is omitted in that particular jackknife sample.