Basic Bootstrap in Mathematica

Abstract

A Mathematica  notebook is presented which provided an overviews of Chapter 2 of Davison and Hinkley 's book  "Bootstrap Methods and their Application".  Mathematica  programs are developed to solve many of  the problem in this chapter, such as the bootstrap methods for standard errors, confidence intervals, etc.

Introduction

Data Assumption

We discuss techniques which are applicable to a single, homogeneous sample of data, denoted by (Subscript[Y, 1,   ...  , ]) (Subscript[Y, n]).  The sample values are thought of  as the outcomes of independent and identically distributed random variables (Subscript[Y, 1,   ...  , ]) (Subscript[Y, n]) whose probability density function (PDF) and cumulative distribution function (CDF) is denoted by f and F, respectively.  The sample is to be used to make inference about a population characteristic, denoted by θ, using a statistic T whose value in the sample is t.

Parametric and Nonparametric Method

When there is a particular mathematical model, with adjustable constants or parameters φ that fully determine  f , such a model is called parametric and statistical methods based on this model are parametric methods; when no such mathematical model is used, the statistical analysis is nonparametric, only using that (Subscript[Y, j]) are independent and identically distributed.

Empirical Distribution

The corresponding estimate of F is the empirical distribution function (EDF) \(TraditionalForm\`F\&^\), which puts equal probabilities (n^(-1))at each sample value (Subscript[Y, j]) and is defined as the sample proportion

\(TraditionalForm\`\(F\&^\)(y) = \(# {y\_j ≤ y}\)\/n\).

More formally

\(TraditionalForm
    \`\(\(F\&^\)(y) = 
      \(1\/\(n\ \)\) \(\(\ ∑\_\(j = 1\)\)\&n\) \(H(y - y\_j)\), \)\)

Or

\(TraditionalForm
      \`\(F\&^\)(y) = 
        \(1\/n\) \(\(∑\_\(j = 1\)\)\&n I(y\_j ≤ y)\)\),

where H(u) is the unit step function which jumps from 0 to 1 at u=0; I {A} is the indicator of the event A = {(Subscript[Y, j]) | (Subscript[Y, j]) ≤ y} .

Statistical Functions

Many simply statistics can be thought of in terms of properties of the EDF.  Often this can be expressed simply as t=t{\(\(F\&^\)(y)\)}, where t(·) is a statistical function.  Such a statistical function is of central importance in the nonparametric case because it also defines the parameter of interest θ through the "algorithm" θ = t(F) .

Example 1

The sample average , \(TraditionalForm\`\(Y\&-\)\) , estimates the population mean

u = ∫ y d F(y).

We substitute for \(TextForm\`\(\(F\&^\) \((y)\)\ \)\)to obtain

\(t(F\&^) = ∫y\ d\)\(\(F\&^\)(y)\)=\(∫y\ 
                \(d(\(1\/n\) \(\(∑\_\(j = 1\)\)\&n H(y - y\_j)\))\) = 
            \(\(1\/n\) 
                \(\(∑\_\(j = 1\)\)\&n
                    ∫y\ d\ \(H(y - y\_j)\)\) = 
              \(\(1\/n\) \(\(∑\_\(j = 1\)\)\&n y\_j\) = 
                \(Y\&-\)\)\)\).

The following exercises show that \(\(F\&^\)(y)\) converges to F as n⟶∞.  Then if t(·) is continuous in an appropriate sense, the definition T = t(·) implies that T converges to θ as n⟶∞, which is the property of consistency.

Exercise 1 ( 2.10  Problem 1)

Let \(TraditionalForm\`F\&^\) denote the EDF. Show that E{(OverHat[F])(y)} = F(y) and that var{(OverHat[F])(y)}=F(y) {1  – F(y)}/n.  Hence deduce that provided 0  <  F(y)  <  1, (OverHat[F])(y) has a limiting normal distribution for large n, and that Pr (| (OverHat[F])(y)–F(y)| ≤ ϵ )→1 as n→∞ for any positive ϵ.

To show that, let

\(TraditionalForm
    \`E(\(F\&^\)(y)) = 
      \(E\ {\(1\/n\) \(\(∑\_\(j = 1\)\)\&n I(y\_j ≤ y)\)} = 
        \(\(1\/n\) \(\(∑\_\(j = 1\)\)\&n E {I(y\_j ≤ y)}\) = 
          \(1\/n\) 
            \(\(∑\_\(j = 1\)\)\&n P(y\_j ≤ y)\)\)\)\)

= P(Y ≤ y) = F(y),

since (Subscript[Y, j]) is identically independent distribution.  We also know that

\(TraditionalForm
      \`var {\(F\&^\)(y)} = 
        \(\(1\/n\^2\) 
            \(\(∑\_\(j = 1\)\)\&n var {I(y\_j ≤ y)}\) = 
          \(\(1\/n\^2\) 
              \(\(∑\_\(j = 1\)\)\&\(n\ \)\(F(y)\) {1 - F(y)}\) = 
            \(1\/n\) \(F(y)\) {1 - F(y)}\)\)\),

since I((Subscript[y, j])  ≤  y)  =  1 or 0 when (Subscript[y, j])  ≤  y or (Subscript[y, j]) >  y, which is binormal for fixed y, and \(TraditionalForm
      \`Pr\ \((I = 1)\)\  = \ \(F(y)\)\ {1  –\ \(F(y)\)}\), then

\(TraditionalForm
      \`var {I(y\_j ≤ y)} = \(F(y)\) {1 - F(y)}\).

By Central Limit Theorem, \(TraditionalForm\`\(F\&^\)(y)\) has a limiting normal distribution for large n, and by Chebyshev's Inequality

\(TraditionalForm
    \`Pr\ \(( | \(F\&^\)(y) - F(y) | \( ≤ ϵ\))\) -> 1
      \)

as n→∞ for any positive ϵ, using var{ (OverHat[F])(y) } = (1/n) F(y) {1  – F(y)}.

Parametric Simulation

Much of statistical theory is devoted to calculating approximate distributions for particular statistics T, on which to base inferences about their estimates θ.  But theoretical properties of T might be hard to determine with sufficient accuracy. One practical alternative of repeated simulation of data sets from a fitted parametric model, and empirical calculation of relevant properties of T.

Suppose that (Subscript[F, φ])(y) and (Subscript[f, φ]) (y) are the CDF and PDF of the data respectively.  When φ is estimated by \(TraditionalForm\`φ\&^\), its substitution in the model gives the fitted model, with CDF \(TraditionalForm\`\(\ \(\(F\&^\)(y) = \(F\_φ\&^\)(y), \)\ 
      \)\)which can be used to calculate properties of T.  We shall use \(TraditionalForm\`\(Y\^*\)\) to denote the random variable distributed according to the fitted model\(TraditionalForm\`\(\ \(F\&^, \)\)\) and the superscript * will be used with mean and var and so forth when these moments are calculated according to the fitted distribution.

Moment Estimates

We write a dataset as \(TraditionalForm\`\(\(Y\_1\^*\),  ... , \(Y\_n\^*\)\ \)\)where the \(TraditionalForm\`\(\(Y\_j\^*\)\ \)\)are independently sampled from the fitted distribution \(TraditionalForm\`\(F\&^\)\).  When the statistic of interest is calculated from a simulated dataset, we denote it by \(TraditionalForm\`\(T\^*\)\).  From R repetitions of the data simulation we obtain \(TraditionalForm\`\(T\_1\^*\),  ... \ , \ \(T\_n\^*\)\).  Properties of T– θ are then estimated from \(TraditionalForm\`\(T\_1\^*\),  ... , \ \(T\_n\^*\)\).  For example, the estimator of the bias b(F)  =  E(T | F) – θ 
       of T is

\(TraditionalForm
    \`\(B = \(b(F\&^) = \(E(T | F\&^) - t = \(E\^*\)(\(T\^*\)) - t\)\), 
    \)\)

and this in turn is estimated by

\(TraditionalForm
    \`B\_R = \(\(R\^\(-1\)\) \(\(∑\_\(j = 1\)\)\&R\( T\_r\^*\)\) - t = 
        \(\(T\&-\)\^*\) - \(t . \)\)\)

The corresponding estimator of the variance of T is

\(TraditionalForm
    \`V\_R = \(1\/\(R - 1\)\) 
        \(\(\ ∑\_\(j = 1\)\)\&n
          \(\((\ \(T\_r\^*\) - \(\(T\&-\)\^*\))\)\^2 . \)\)\)

Distribution and Quantile Estimates

We are approximating the distribution of T – θ  by that of \(TraditionalForm\`\(\(T\^*\) - t\ \)\)then cumulative probabilities are estimated simply by the empirical distribution function of the simulated values \(TraditionalForm\`\(\(t\^*\)\  - \ \(t . \)\ \ \)\)More formally, if \(TraditionalForm
      \`G(u)\  = \ Pr(T–\ θ\  ≤ \ u\ )\), then the simulation estimate of G(u) is

\(\(\(G\_R\)\&^\)(u) = 
              \(# {\(t\_r\^*\) - t\  ≤ \ u}\)\/R\)=\(\(1\/R\) 
            \(\(∑\_\(r = 1\)\)\&n 
                I {\(t\_r\^*\) - t ≤ u}\)\),

where I {A} is the indicator of the event A.

We are interested in quantiles of the distribution of T – θ , and these are approximated using ordered values of \(TraditionalForm\`\(t\^*\) - t\) .   So we estimate the p quantile of T – θ  by the (R + 1) pth ordered value of (t^(* )) – t, that is \(TraditionalForm\`\(t\^*\)\_\(\((R + 1)\) p\) - \(t . \)\) We assume that R is chosen so that (R + 1) p is an integer.

Nonparametric Simulation

Assuming that (Subscript[Y, 1,   ...  , ]) (Subscript[Y, n ])are independent and identically distributed according to an unknown distribution function F.  We use the EDF \(TraditionalForm\`F\&^\) to estimate the unknown CDF.

Example 2 (Average)

In the case of the average, exact moments under sampling form the EDF are found.  For example,

\(TraditionalForm
    \`\(E\^*\)(\(Y\&-\)\^\(*\ \)) = 
      \(\(E\^*\)(\(Y\^*\)) = 
        \(\(\(\ ∑\_\(j = 1\)\)\&n\) \(1\/n\) Y\_j = \(Y\&-\)\)\)\)

since Example 1. Similarly

\(TraditionalForm
    \`\(var\^*\)(\(\(Y\&-\)\^*\)) = 
      \(\(1\/n\) \(\(var\^*\)(\(Y\^*\))\) = 
        \(\(1\/n\) \(E\^*\) {\(Y\^*\) - \(E\^*\)(\(Y\^*\))}\^2 = 
          \(1\/n\) \(\(\ ∑\_\(j = 1\)\)\&n\) \(1\/n\) 
            \((Y\_j - \(Y\&-\))\)\^2\)\)\)

Exercise 2 ( 2.10 Problem 3)

Under nonparametric simulation from a random sample (Subscript[y, 1,   ...  , ]) (Subscript[y, n]) in which  \(TraditionalForm
      \`T = \(n\^\(-1\)\) \(∑\ \((Y\_j\  - \ \(Y\&-\)\ )\)\^2\)\) takes value t.  To show that

\(TraditionalForm\`\(\(E\^*\)(\(T\^*\)) = \(\((n - 1)\) t\)\/n, 
      \ \ \ \ \ \ \ 
      \(var\^*\)(\(T\^*\)) = 
        \(\((n - 1)\)\^2\/n\) {m\_4\/n + \(\((3 - n)\) t\^2\)\/\(n(n - 1)\)}
            \ \ \)\),

where (Subscript[m, 4])=(n^(-1))(∑_j)(SuperStar[((Subscript[y, j]) - y)]),

let

\(TraditionalForm
    \`\(E\^*\)(\(T\^*\)) = 
      \(\(n\^\(-1\)\) 
          \(\(E\^*\)(
            \(∑\_\(j = 1\)\)\&n\( y\_j\^*\)\^2 - n\ \(\(y\^*\)\&-\)\^2)
            \) = \(\(n\^\(-1\)\) 
              \(\(∑\_\(j = 1\)\)\&n\(\( E\^*\)(\(y\_j\^*\))\)\^2\) - 
            \(\(E\^*\)(\(\(y\&-\)\^*\))\)\^2 = 
          \(\(E\^*\)(\(y\_j\^*\))\)\^2 - 
            \(\(\(E\^*\)(\(\(y\&-\)\^*\))\)\^2 . \)\)\)\)

Since

\(TraditionalForm
    \`\(\(E\^*\)(\(\(y\&-\)\^*\))\)\^2 = 
      \(var(\(\(y\&-\)\^*\)) + \((\(E\^*\) \(\(y\&-\)\^*\))\)\^2 = 
        \(1\/n\^2\) \(\(∑\_\(j = 1\)\)\&n\((y\_j - \(y\&-\))\)\^2\) + 
          \(y\&-\)\^2\)\)

by Example 2, and

\(TraditionalForm
      \`\(E\^*\)(y\_j\^2) = 
        \(∫\(\((\(y\_j\^*\))\)\^2\) d\ F\&^\ \((\(y\_j\^*\))\) = 
          \(∫\(y\^2\) d\ \(\(F\&^\)(y)\) = 
            \(\(1\/n\) \(\(\ ∑\_\(j = 1\)\)\&n\) 
                \(∫\(y\^2\) d\ \(H(y - y\_j)\)\) = 
              \(1\/n\) \(\(∑\_\(j = 1\)\)\&n y\_j\^2\)\)\)\)\),

then

\(TraditionalForm
    \`\(E\^*\)(\(T\^*\)) = 
      \(\(1\/n\) \(\(∑\_\(j = 1\)\)\&n y\_j\^2\) - 
          \(1\/n\^2\) \(\(∑\_\(j = 1\)\)\&n\((y\_j - \(y\&-\))\)\^2\) - 
          \(y\&-\)\^2 = 
        \(\(n - 1\)\/n\^2\) 
          \((\(∑\_\(j = 1\ \)\)\&n y\_j\^2 - n\ \(y\&-\)\^2)\)\)\)

\(TraditionalForm
    \`\( = \(\(n - 1\)\/n\) {
          \(1\/n\) \(\(∑\_\(j = 1\)\)\&n\((y\_j - y\&\(-\ \))\)\^2\)}
          \) = \(\((n - 1)\) t\)\/\(n . \)\)

Let

\(TraditionalForm
      \`\(var\^*\)(\(T\^*\)) = 
        \(\(var\^*\) {
              \(n\^\(-1\)\) 
                \(\(∑\_\(j = 1\)\)\&n
                    \((\(y\_\(j\ \)\^*\) - \(y\&-\))\)\^2\)} = 
          \(\((n - 1)\)\^2\/n\^2\) 
            \(var\^*\) {
              \(1\/\(n - 1\)\) 
                \(\(∑\_\(j = 1\)\)\&n
                    \((\(y\_j\^*\) - \(\(y\&-\)\^*\))\)\^2\)}\)\).

We know that

\(TraditionalForm
      \`\(var\^*\) {
            \(1\/\(n - 1\)\) 
              \(\(∑\_\(j = 1\)\)\&n\((\(y\_j\^*\) - \(\(y\&-\)\^*\))\)\^2
                  \)} = \(k\_4\^*\)\/n + \(2 \( k\_2\^*\)\)\/\(n - 1\)\),

where

\(TraditionalForm
    \`\(k\_2\^*\) = 
      \(\(var\^*\)(\(y\^*\)) = 
        \(\(1\/n\) \(\(∑\_\(j = 1\)\)\&n\((y\_j - \(y\&-\))\)\^2\) = t
          \)\)\)

\(TraditionalForm
    \`\(k\_4\^*\) = 
      \(E\^*\) y\^\(*4\) - 4 \((\(E\^*\) \(y\^*\)\^3)\) . \(E\^*\) \(y\^*\) - 
        3 \((\(E\^*\) y\^\(*2\))\)\^2 + 
        12 \((\(E\^*\) y\^\(*2\))\) . \((\(E\^*\) \(y\^*\))\)\^2 - 
        6 \((\(E\^*\) \(y\^*\))\)\^4\)

After some algebra caculation

\(TraditionalForm
    \`\(var\^*\) {
          \(1\/\(n - 1\)\) 
            \(\(∑\_\(j = 1\)\)\&n\((\(y\_j\^*\) - \(\(y\&-\)\^*\))\)\^2
                \)} = \(\(E\^*\)(\(y\^*\) - \(\(y\&-\)\&\(\  \)\^*\))\)\^4 - 
        3 {\(E\^*\) \(y\^*\)\^2 - \((\(E\^*\) \(y\^*\))\)\^2}\^2\)

\(TraditionalForm
    \`\( = \(E\^\(*\ \)\)(\(y\^*\) - \(\(y\&-\)\^*\))\)\^4 - 
      3 \((\(var\^*\) \(y\^*\))\)\^2\)

Hence

\(TraditionalForm\`\(\(var\^*\)(\(T\^*\)) = \)\)\(TraditionalForm
      \`\(\((n - 1)\)\^2\/n\^2\) {
          \(k\_4\^*\)\/n + \(2 \((\(k\_2\^*\))\)\^2\)\/\(n - 1\)}\)\(TraditionalForm
      \`\( = \(\((n - 1)\)\^2\/n\^2\) {
          \(1\/n\) \(\(E\^*\)(\(y\^*\) - \(\(y\&-\)\^*\))\)\^4 - 
            \(3\/n\) \((\(var\^*\) \(y\^*\))\)\^2 + 
            \(2 \((\(var\^*\) \(y\^*\))\)\^2\)\/\(n - 1\)}\)\)

\(TraditionalForm
    \`\( = \(\((n - 1)\)\^2\/n\^2\) {
          \(1\/n\) 
              \(\(∑\_\(j = 1\)\)\&n\((\(y\^*\) - \(\(y\&-\)\^*\))\)\^4
                  \) - \(3  t\^2\)\/n + \(2  t\^2\)\/\(n - 1\)}\) = 
      \(\((n - 1)\)\^2\/n\^2\) {
          \(1\/n\) 
              \(\(∑\_\(j = 1\)\)\&n\((\(y\^*\) - \(\(y\&-\)\^*\))\)\^4
                  \) + \(\((3 - n)\) t\^2\)\/\(n(n - 1)\)}\)

= ((((n - 1))^2)/(n^2)) 
        {((Subscript[m, 4])/n) + (((3 - n) (t^2))/(n(n - 1)))} .

Comparison with parametric methods

Exercise 3 ( 2.11 Practical 1)

City population data  reports n=49 data pairs, each corresponding to a city in the United States of America, the pair being the 1920 and 1930 populations of the city, which we denote by u and x.  The data is plotted in Figure 3.1.  Interest here is in the ratio of means, because this would enable us to estimate the total population of the USA in 1930 from 1920 figure.  That is, if the cities form a random sample with (U, X ) denoting the pair of population values for a randomly selected city, then the ratio of expectations θ  =  E(X)/E(U) is the parameter of interest.

Figure 3.2

It is natural to estimate the ratio by T= \(TraditionalForm\`\(\(X\&-\)/\(\(U\&-\) . \)\ \ \)\)Frist let us  consider parametric  simulation to estimate the distribution of the ratio when a bivariate lognormal distribution is fitted to the city poulation data.

The left panel of Figure 3.2 shows a histogram of \(TraditionalForm\`\(\(t\^*\) = \(\(X\&-\)\^*\)/\(\(U\&-\)\^*\)\ 
      \)\)and right panel shows a histogram of studentized estimates \(TraditionalForm\`\(\(z\^*\)\ , \)\) under parametric simulation, R=999.  Note the skewness of both \(TraditionalForm\`\(t\^*\)\)and \(\(z\^*\)\ \).  The studentized values \(\(z\^*\) = \((\(t\^*\) - t)\)/v\_L\^\(*\ 1/2\)\), where \(\(v\_L\^*\)\ \)is the delta method variance estimate based on simulated data by result of 2.10 Problems 9.  That is,

\(TraditionalForm
      \`\(v\_L\^*\) = 
        \(n\^\(-2\)\) \(\(\ ∑\_\(j = 1\)\)\&n\) 
          \((\(x\_j\^*\) - \(t\^*\) \(u\_j\^*\))\)\^2/\(\(u\&-\)\^*\)\^2\).

Second, comparing parametric simulation with nonparametric resampling  for sample of size n=10, R=999 simulation .  Figure 3.3 plots the histograms of  ratio \(TraditionalForm\`\(t\^*\)\ for\ the\ resamples\ (\) left panel) and  studentized estimates \(TraditionalForm\`\(\(\(z\^*\)(\ right\ panel\ ) . \)\ \)\)  From Figure 3.2 and 3.3, we see that histograms of \(TraditionalForm\`\(t\^*\)\)and \(TraditionalForm\`\(z\^*\)\)are similar in both cases.

Figure 3.3

Simple Confidence Interval

The major application for distribution and quantiles of an estimator T is in the calculation of confidence limits.  Here we describe briefly basic methods.

Using a normal approximation for the distribution of T

It may be possible to show that T is approximately normal with mean θ + β and variance υ; here β is the bias of T .  If β and υ are both unknown, we can use bootstrap estimates of bias and variance, and if the α quantile of the standard normal distribution is (Subscript[z, α ]) = (Φ^(-1))(α), 
           then an approximate (1 - 2  α) confidence interval for θ has limits

\(TraditionalForm
      \`t - β\&^ - \(υ\&^\^\(1/2\)\) z\_\(1 - α\ \), 
      \ \ \ \ \ \ \ \ 
      t - β\&^ - \(υ\&^\^\(1/2\)\) z\_\(α\ \)\),

as follows from

\(TraditionalForm
      \`Pr(β\&^ + \(υ\&^\^\(1/2\)\) z\_α ≤ 
            T - θ ≤ 
            β\&^ + \(υ\&^\^\(1/2\)\) z\_\(1 - α\))
          \( = \& . \)1 - 2  α\) .

Using the bootstrap estimates of quantiles for T - θ

An equitailed (1 - 2  α) confidence interval will have limits

\(TraditionalForm
    \`t - \((\(t\^*\)\_\(\((R + 1)\) \((1 - α)\)\) - t)\), \ \ \ 
    t - \((\(t\^*\)\_\(\((R + 1)\) α\) - t)\), \)

which is as the basic bootstrap confidence limits.  This is based on the probability implication

Pr(a < T - θ < b) = 1 - 2  α.

Using the studentized bootstrap statistic

Studentized statistic  Z = ((T - θ)/(V^(1/2))),

where V  is an estimate of var(T  |  F) .

The studentized (1 - 2  α) confidence interval limits for θ have the form

\(TraditionalForm
    \`\(t - \(v\^\(1/2\)\) \(z\_\(\(1 - α\ , \)\ \ \ \ \ \ \ \)\) t - 
      \(v\^\(1/2\)\) z\_α\ , \)\)

where (Subscript[z, p]) denotes the p quantile of Z. It is more accurate to estimate the quantiles of Z from replicates of the studentized bootstrap statistic, \(TraditionalForm\`\(Z\^*\) = \((\(T\^*\) - t)\)/\(V\^*\)\^\(1/2\)\), where \(TraditionalForm\`\(T\^*\)\) and \(TraditionalForm\`\(V\^*\)\) are based on a simulated random sample, \(TraditionalForm\`\(Y\_1\^*\),  ... \ , \(\(Y\_n\^*\) . \)\)  If the model is parametric, the \(TraditionalForm\`\(Y\_j\^*\)\) are generated from the fitted parametric distribution, and if the model is nonparametric, they are generated from the EDF \(TraditionalForm\`\(F\&^ . \)\)  This studentized bootstrap method is most likely to be of use in nonparametric problems.  We use the (R + 1) αth order statistic of the simulated values \(TraditionalForm\`\(z\_1\^*\),  ... \ , \(z\_R\^*\)\ for\ Z, \ namely
      \) \(TraditionalForm\`\(z\^*\)\_\(\((R + 1)\) α\)\), to estimate (Subscript[z, α]) .  Then studentized bootstrap confidence interval for θ has limits

\(TraditionalForm
      \`t - \(υ\^\(1/2\)\) 
          \(\(z\^*\)\_\(\(\((R + 1)\) \((1 - α)\), \)\ \ \ \ \)\) t - 
        \(υ\^\(1/2\)\) \(z\^*\)\_\(\((R + 1)\) α\)\),

Under the nonparametric model for city population data as  Exercises 3, we can calculate the 95 % basic bootstrap confidence interval for the ratio is ( 0.924, 1.787); and the 95 % studentized bootstrap confidence interval for the ratio is (1.241, 2.098).

Exercise 4 ( 2.11 Practical 2 )

Co.transfer data contains data on the carbon monoxide transfer factor for seven smokers with chickenpox, measured on admission to hospital and after a stay of one week.  The aim is to estimate the average change in the factor.  That is, the expectation of average difference after a stay of one week is the parametric interest t.

co = ReadList["o:\yzhang\boot\d.txt", Table[Real, { 2}]];

entry = Transpose[co][[1]];

week = Transpose[co][[2]];

<<Statistics`HypothesisTests`

MeanTest[ week -  entry, 0,  FullReport -> True]

<< Statistics`ConfidenceIntervals`

MeanCI[week - entry]

The interval shown above is the 95% interval of  true difference calculated based on a normal distribution for difference.

\(<< q:\459\eiss.m\)

NPlot[week - entry]

Graph shown above is normal quantile plot of difference of week and entry, which does not indicates that the differences are close to normal, since the size of sample is small. So confidence interval based on normal distribution is questionable. Next we try to calculate the basic bootstrap (1-2α) confidence interval, using function " BasicCon".

BootStats[data_,  statistic_, n_] := 
      Table[statistic[
          data[[Table[
              Random[Integer, {1, Length[data]}], {Length[data]}]
              ]]], {n}];

mean[data_] := Apply[Plus,  data]/Length[data];

BasicCon[data_,  statistic_, n_, α_] := 
      Module[{t}, t = Sort[BootStats[data, statistic, n]]; 
        t = 2  statistic[data] - 
            t[[{Floor[(n + 1)*(1 - α)], 
                Floor[(n + 1) α]}]]]

BasicCon[week - entry, mean, 999,  0.025]

Let us calculate the studentized bootstrap (1-2α) confidence interval using function " StudentCon".

var[data_] := 
      Apply[Plus, (data - mean[data])^2]/(Length[data] 
          (Length[data] - 1))

new[data_] := {mean[data], var[data]}

temp = BootStats[week - entry, new, 999];

Stats = Transpose[temp];

Student = 
      Sort[((Stats[[1]] - 
            mean[week - entry])/Sqrt)[
            Stats[[2]]]];

StudentCon[data_, statistic_, n_, α_] := 
      statistic[data] - 
        Sqrt[var[data]]*
          Student[[{Floor[(n + 1)*(1 - α)], 
              Floor[(n + 1) α]}]];

In function " StudentCon", the v = var(T | F) is estimated by sample variance of T (the mean of difference of week and entry). That is,\(TraditionalForm\`\(\ 
      v\&^ = \(1\/\(n(n - 1)\)\) \(\(∑\_\(j = 1\)\)\&n y\_j\)\)\).

StudentCon[week - entry, mean, 999, 0.025]

we check the distribution of bootstap estimate (t^(* ))and studentized bootstrap statistics \(TraditionalForm\`\(z\^*\)\) based on nonparametric resampling R=999 and see why we use nonparametric quantile estimate to calculate the confidence interval.

NPlot[Stats[[1]]]

The graph shown above is the normal quantile plot of \(TraditionalForm\`\(t\^*\)\) values based on resampling R=999.  It indicates that \(TraditionalForm\`\(\(t\^*\)\ \)\)has a little bit  fat tail.

NPlot[Student]

The graph shown above is the normal quantile plot of \(TraditionalForm\`\(z\^*\)\) values based on resample R=999.  It indicates that there may be some simulation outliers in left tail and this may effect the goodness  of results from the studentized bootstrap.

Finally, let us compare the variance of the bootstrap estimate \(TraditionalForm\`\(t\^*\)\) with the estimated variance.

var[week - entry]

\(\(ListPlot[temp]
\)\)

The graph shown above plots the variance of the bootstrap estimate \(TraditionalForm\`\(t\^*\)\) versus \(TraditionalForm\`\(t\^*\)\).  The sample variance of t is 33.7823.  The graph shows us how varibility of the variance of the bootstrap estimate \(TraditionalForm\`\(\(t\^*\)\ \)\)comparing with sample variance of t. This implies that we should do transformation to stablize the variance and calculate the much more reasonable confidence interval.

Exercise 5 ( 2.11 Problem 4 )

Let t be the median of a random sample of size n = 2  m + 1 with ordered values (Subscript[y, (1)])  <  ...  < (Subscript[y, (n)]);  t = (Subscript[y, (m + 1) . ])  Then \(TraditionalForm\`\(\(T\^*\) > \)\) (Subscript[y, (l)]) if and only if fewer than m + 1of the \(TraditionalForm\`\(Y\_j\^*\)\) are less than or equal to (Subscript[y, (l) . ])  Otherwise, there is contradiction.

Furthermore, let j =  # of\(TraditionalForm\`\(\ \(Y\^*\)\)\), which are less than or equal to (Subscript[y, (l)]), then if \(TraditionalForm\`\(T\^*\) > y\_\((l)\)\) if and only if j < m + 1,  
      i . e .   j ∈ (0, 1,  ...  , n).  So

\(\(Pr\^*\)(\(T\^*\) > y\_\((l)\))\)=(Overscript[( (∑_(j = 0))), m])(({{n}, {j}}))((((l/n)))^j)(((1 - (l/n)))^(n - j)),

since \(TraditionalForm\`\(Pr(\(Y\^*\) ≤ y\_\((l)\)) = l\/n, \ \ 
      Pr(\(Y\^*\) > y\_\((l)\)) = 1 - l\/n, \ and\ \)\)random variable j is binormal.

Suppose n=11, then m=5

Prob=\(TraditionalForm
      \`\(Pr\^*\)(\(T\^*\) ≤ y\_\((3)\)) = 
        \(\(Pr\^*\)(\(T\^*\) ≥ y\_\((9)\)) = 
          1 - \(Pr\^*\)(\(T\^*\) > y\_\((3)\))\)\).

Using the resampling distribution shown above

Prob = 
      N[1 - Sum[
            ((11!)/((11 - j)! j!)) ((((3/11)))^j) 
              (((1 - (3/11)))^((11 - j))), {j, 0, 5}]]

so the 5% and 95% quantiles of \(TraditionalForm\`\(T\^*\)\) are (Subscript[y, (3)]) and (Subscript[y, (9)]).  The basic bootstrap 90% confidence interval for the population median θ is \(TraditionalForm
      \`\(\((2  y\_\((6)\) - y\_\((9)\), \ 2  y\_\((6)\) - y\_\((3)\))\) . 
        \)\)

Let examine the coverage of the confidence interval of population median for samples from normal and Cauchy distribution.

<< "Statistics`ContinuousDistributions`"

gdist1  =  NormalDistribution[0, 1];

gdist2  =  CauchyDistribution[0, 1];

int[f_] := 
      Module[{r}, r = Sort[RandomArray[f,  11]]; 
        r = {2  r[[6]] - r[[9]], 2  r[[6]] - r[[3]]}]

int[gdist1]

int[gdist2]

Reducing Error

Two types of errors in resampling methods are statistical error and simulation error.  The first of these is due to the difference between F and \(TraditionalForm\`F\&^\).  The simulation error is wholly due to use of empirical estimates of properties under sampling from \(TraditionalForm\`\(F\&^ . \)\)

Statistical Error

The basic bootstrap ideal is to approximate a quantity c(F) by the estimate \(TraditionalForm\`\(c(F\&^\)\)), where \(TraditionalForm\`F\&^\) is either a parametric or a nonparametric estimate of F based on data (Subscript[y, 1,   ...  , ]) (Subscript[y, n . ])  The statistical error is then the difference between \(TraditionalForm\`\(\(\(c(F\&^)\)\ and\ \(c(F)\), \)\ \)\)and as far as possible we wish to minimize this or remove it entirely.  This is sometimes possible by careful choice of c(·).

For confidence interval caculation, we assume that T is decided upon, the statistical problem is to choose possible quantity \(TraditionalForm\`Q = q(F\&^, F)\) so that the resulting Q is as nearly pivotal as possible.

Parametric problem

In parametirc problem \(TraditionalForm\`F\&^ = \(F\_ψ\&^\ and\ F = F\_ψ\)\) have the same form, differing only in parameter values. We define a pivot as a function Q = q(T, θ), that is for all q,

Pr {q(T,  θ) ≤ q  | φ},

is independent of ψ.

Exercise 6 ( 2.10 Problem 2 )

Suppose that (Subscript[Y, 1]),   ...  ,  (Subscript[Y, n]) are independent exponential with mean μ; their average is

\(TraditionalForm\`\(Y\&-\) = \(n\^\(-1\)\) Y\_j\) .

(a) Show that \(TraditionalForm\`\(Y\&-\)\) has the gamma density , so its mean and variance are μ and (μ^2)/n .

Since (Subscript[y, 1]),   ...  , (Subscript[y, n]) are independent exponential with mean μ, density function f(y) = (1/μ) (e^(-y/μ)),  0 ≤ y ≤ ∞, and moment generation function (Subscript[m, y])(t) = 1/(1 - μ t),  t < 1/n,   moment generation function of (Overscript[((∑_(j = 1))), n]) (Subscript[y, j]) is

\(TraditionalForm
      \`\(M\_\(y\_j\)\)(t) = 
        \(\((\(m\_y\)(t))\)\^n = 1\/\((1 - μ\ t)\)\^n\)\) .

Then the generation function of (1/n) (Subscript[y, j]) is (1/(((1 - (n^(-1)) μ t))^n)).  The average (1/n) ∑(Subscript[y, j]) has a distribution of Gamma (n,   (μ/n)).  Then

E((1/n) (Subscript[y, j])) = μ;   
      var((1/n) (Subscript[y, j])) = ((μ^2)/n).

(b) Show that \(TraditionalForm\`log\ \(Y\&-\)\) is approximately normal with mean log μ and variance (n^(-1)) .

\(TraditionalForm\`\(Y\&-\)\) approximates to normal with mean μ and variance ((μ^2)/n) .    By Taylor expansion,

Log \(\(Y\&-\) ≈ log\ μ\)+\(\(1\/n\) \((\(Y\&-\) - μ)\)\),

where \(TraditionalForm\`log\ μ + \(1\/n\) \((\(Y\&-\) - μ)\)\) is normal with mean log μ and variance (1/n).  Hence log \(TraditionalForm\`\(\(Y\&-\)\ \)\)is approximately normal with mean log μ and variance (n^(-1)) .

(c) Compare the normal approximations of \(TraditionalForm\`\(\(Y\&-\)\ \)\)and log \(TraditionalForm\`\(Y\&-\)\) in calculating 95% confidence interval for μ.

Under normal approximation for \(TraditionalForm\`\(Y\&-\)\), 95% confidence interval for μ is

\(TraditionalForm
      \`\(Y\&-\)/\((1 + 2\/\@n)\) < μ < \(Y\&-\)/\((1 - 2\/\@n)\)\).

Under normal approximation for log \(TraditionalForm\`\(\(\(Y\&-\), \)\ \)\)95% confidence interval for μ is

\(TraditionalForm
      \`log\ \(Y\&-\) - 2\/\@n < log\ μ < log\ \(Y\&-\) + 2\/\@n\),

i.e.

\(TraditionalForm
      \`exp[log\ \(Y\&-\) - 2\/\@n] < μ < 
        exp[log\ \(Y\&-\) + 2\/\@n]\).

Under Gamma (n,   (μ/n)) for \(TraditionalForm\`\(Y\&-\)\), we know \(TraditionalForm\`\(Y\&-\)\)/ μ has Gamma (n, (n^(-1))), then we can calculat the exact confidence limits.

Air.conditioning data

aircon = {3, 5, 7, 18, 43, 85, 91, 98, 100, 130, 230, 487};

Air.conditioning data gives n=12 times between failures of air-conditioning equipment, for which we wish to estimate the underlying mean or its reciprocal, the failure rate.  Frist we try to find the model of data.

EDF[data_,  x_] := 
      (Length[Select[data,  #1 ≤ x &]])/Length[data]

Function "ECDF[data_, x_]" calculates the empirical  distribution function (OverHat[F])(x).

Plot[EDF[aircon, x], {x, Min[aircon], Max[aircon]}, 
      AxesLabel -> {"Failure time", "EDF"}];

Let us compare the EDF of air.conditioning and CDF of exponential distribution with mean μ set equal to sample average.

<<Statistics`ContinuousDistributions`

gdist = ExponentialDistribution[1/par];

par = N[Mean[aircon]]

CDF[gdist, x]

Plot[{EDF[aircon, x], CDF[gdist, x]}, {x, Min[aircon], Max[aircon]}, 
      AxesLabel -> {"Failure time", "CDF/EDF"}]

Comparison of the two curves suggests that the exponential model fits reasonable well.  So under the normal approximation for \(TraditionalForm\`\(Y\&-\)\) ,  95% confidence interval for μ is

con1 = {par/(1 + (2/(Sqrt[12]))),  par/(1 - (2/(Sqrt[12])))}

Under normal approximation for log \(TraditionalForm\`\(Y\&-\)\), 95% confidence interval for μ is

con2 = {Exp[Log [par] - (2/((Sqrt[12]) ))],  Exp[Log [par] + (2/(Sqrt[12]))]}

Under  exponential model for Y,(Overscript[Y, -]) has   Gamma (n,   (μ/n)) ,  (Overscript[so Y, -])/μ has Gamma (n, (n^(-1))) .

q1 = Quantile[GammaDistribution[12, (12^(-1))], 0.025]

q2 = Quantile[GammaDistribution[12, (12^(-1))], 0.975]

95% confidence interval for μ is

con3 = {par/q2, par/q1}

The 95% confidence interval of μ with Gamma model for \(TraditionalForm\`\(Y\&-\)\) is closer to the 95% confidence interval of μ with normal approximation for Log transformation of \(TraditionalForm\`\(Y\&-\)\) than with normal approximation for \(TraditionalForm\`\(Y\&-\)\) itself.

Exercise 7 ( 2.10 Problem 6 )

Suppose the (Subscript[Y, 1 ,  , , ,  ,  ]) (Subscript[Y, n]) are independently distributed with a two-parameter density (Subscript[f, θ, λ])(y).  The following simulation experiment would perform to check whether or not \(TraditionalForm
      \`Q = q(\(Y\_\(1\ , \ , , , \ , \ \)\) Y\_n; \ θ)\) is a pivot.

Step1: fix parameters (Subscript[θ, j]), λ and simulate (Subscript[Y, 1]),  ...  ,  (Subscript[Y, n]) from (Subscript[f, (Subscript[θ, j]), λ])(y) .

Step2: calculate the value of statistic Q  by (Subscript[Y, 1]),   ...   ,  (Subscript[Y, n]) and denote as (Subscript[Q, i, j]) .

Step3: repeat Step1 and Step2 N times and compute N values of statistics (Subscript[Q, 1  j]) ,   ...  ,  (Subscript[((Subscript[Q,  N])), j]), for (Subscript[θ, j . ])

Step4: compute the EDF of sample (Subscript[Q, 1  j]) ,   ...  ,  (Subscript[((Subscript[Q, N])), j]).

Step5: change of value of θ and repeat Step 1up to  Step 4.

Step6: compare the results of Step 4 and Step 5.

For example, if f is the gamma density, its density function is as following:

<< Statistics`ContinuousDistributions`

PDF[GammaDistribution[κ, μ/κ],   y]

Let \(TraditionalForm\`μ\&^\) be the MLE of μ, let

\(TraditionalForm
      \`\(l\_p\)(, μ) = 
        max \(\(∑\_\(j = 1\ \)\)\&n log\ 
              \(\(f\_\(μ, κ\)\)(y\_j)\)\)\),

be the profile log likelihood for μ and let \(TraditionalForm
      \`Q = 2 \({\(l\_p\)(μ\&^) - \(l\_p\)(μ)} . \)\)  In theory Q should be approximately a ((Subscript[χ, 1])^2) variable for large n.  Use simulation to examine whether or not Q is approximately pivotal for n=10 when κ is in the range (0.5, 2).

In this case the average of simulated sample is the MLE of μ.

f = GammaDistribution;

\(why[f_, m_, k_, a_, q_] := 
      Module[{i, g}, 
        g[i_] := 
          Module[{t}, 
            t = 
	Table[
                RandomArray[f[\(4 + b\)\/10, \(10  a\)\/\(4 + b\)], m], {b, 
                  k}]; 	t = 
              \((\(-Max[
                      Table[Apply[Plus, 
                          Log[PDF[f[\(4 + b\)\/10, \(10  a\)\/\(4 + b\)], 
                              t[[b
                                ]]\ ]]], {b, k}]]\) + 
                  Max[Table[
                      Apply[Plus, 
                        Log[PDF[
                            f[\(4 + b\)\/10, 
                              \(10  Mean[
                                    t[[b
                                      ]]]\)\/\(4 + b\)], 
                            t[[b
                              ]]\ ]]], {b, k}]])\); 
	
            2  t
		]; Table[g[i], {i, q}]]\)

Function "why[f, m, k, a, q]" calculates the simulated statistic Q with distribution f and simulated sample size m for fixed parameter value μ = a .   Variable k is defined as the number of values taken in range  (0.5, 2) for another parameter κ. Variable q is defined as the number of repeated calculation for Q .   That is, using function "why" , we can simulate sample (Subscript[Q, 1]),  ...  , (Subscript[Q, q])\(TraditionalForm\`\(\^\(\  \)\) of\) Q with fixed parameter value μ.

EDF[data_,  x_] := 
      (Length[Select[data,  #1 ≤ x &]])/Length[data]

Function "EDF[data, x]" calculates the EDF for data at value x.  Next using function "why" to calculate three simulated samples of Q with (Subscript[μ, 1]) = 1,  (Subscript[μ, 2]) = 2,  (Subscript[μ, 3]) = 3 respectively. Then we plot EDFs of these three samples.

why1 = why[f, 10, 16, 1, 20];

why2 = why[f, 10, 16, 2, 20];

why3 = why[f, 10, 16, 3, 20];

dist2 = ChiSquareDistribution[1];

Plot[{EDF[why1, x], EDF[why2, x], EDF[why3, x]}, {x, 0, 
        Max[{why1, why2, why3}]}, AxesLabel -> {"Q", "EDF"}]

The graph shown above plots EDFs of simulated samples calculated above. The three EDFs are not close.   So approximately Q is not pivotal for n = 10 when κ is in the range (0.5, 2).

gadata1 = why[f, 100, 16, 1, 20];

gadata2 = why[f, 500, 16, 1, 20];

com1 = 
      Plot[{EDF[gadata1, x], CDF[dist2, x]}, {x, 0, Max[gadata1]}, 
        Frame -> True, FrameLabel -> {"Q", "EDF/CDF"}, 
        PlotLabel -> "n=100", DisplayFunction -> Identity];

com2 = 
      Plot[{EDF[gadata2, x], CDF[dist2, x]}, {x, 0, Max[gadata2]}, 
        Frame -> True, FrameLabel -> {"Q", "EDF/CDF"}, 
        PlotLabel -> "n=500", DisplayFunction -> Identity];

Show[GraphicsArray[{com1, com2}], Frame -> True, 
      DisplayFunction -> $DisplayFunction];

The left panel of the graph shown above compares EDF of Q when n =100, μ = 1 and κ in range (0.5, 2) with CDF of ((Subscript[χ, 1])^2).  The right panel compares EDF of Q when n =500, μ =1 and κ in the same range with CDF of ((Subscript[χ, 1])^2).  They indicates that Q should be approximately a (Subscript[((χ^2)), 1]) variable for large n.

nonparametric problem

Exercise 8 ( 2.11 Practical 3 )

Data "cd4" contains the CD4 counts in hundreds for 20 HIV-positive patients at baseline and after one year of treatment with an experimental anti-viral drug.  We attempt to set a confidence interval for the correlation between the baseline and later counts, using the nonparametric bootstrap.

Figure 8.1

The left panel of Figure 8.1 is a scatter plot of sample standare deviation of correlation using delta method versus sample correlation for nonparametric resampling.  This plot suggests that the true model has standard deviation proportional to its correlation coefficient. The variance is  depedent of correlation.  This implies that we might transform the correlation coefficient.  The right panel of Figure 8.1 is a scatter plot of sample standard deviation of correlation versus sample correlation coefficient after fisher transformation. There is no obvious patten in the right plot.  The variance of large vales of correlation coefficient looks chang much, which may are caused by simulation outliers. We calculate studentized bootstrap 95% confidence interval on transformed scale, which is (0.559, 1.272), comparing with interval (0.208, 1.018) obtained without the transformation.

Simulation error

Simulation error arises when Monte Carlo simulations are performed and properties of statistics are approximated by their empirical properties in these simulations.  The Monte Carlo variability can only be removed entirely by an infinite simulation, which seems both impossible and unnecessary in practice.  The practical question is, how large simulation does need to be to achieve reasonable accuracy.

Exercise 9

This exercise tries to tell how many simulations are required for quantile estimation.  To get some idea, we make four replicate plots with 39, 99, 399, and 999 simulations.

<< Statistics`ContinuousDistributions`

n = {39, 99, 399, 999};

\(\(error[i_] := 
      With[{t = 
            RandomArray[NormalDistribution[0, 1], 
              n[[i]]]}, 
		
        quantiles = 
          \(Sort[t]\)[[{Ceiling[0.025*\((Length[t] + 1)\)], 
              Ceiling[0.975*\((Length[t] + 1)\)]}]]; 

        				nq = 
          Table[Quantile[NormalDistribution[0, 1], \((m - 0.5)\)/Length[t]], {
              m, 1, Length[t]}]; 
L1 = Transpose[{nq, Sort[t]}]; 
	
        ListPlot[L1, 
          Epilog -> {\ 
              Line[{{Min[t], 
                    quantiles[[1]]}, {
                    Max[t], 
                    quantiles[[1]]}}], 
              Line[{{Min[t], 
                    quantiles[[2]]}, {
                    Max[t], 
                    quantiles[[2]]}}]}]]
        
	\)\)

re[i_] := {{error[i], error[i]}, {error[i], error[i]}}

re1 = re[1]

Show[GraphicsArray[re1], 
      PlotLabel -> "Four replicate plots with 39 simulations", 
      Frame -> True]

re2 = re[2]

Show[GraphicsArray[re2], 
      PlotLabel -> "Four replicate plots with 99 simulations", 
      Frame -> True]

re3 = re[3]

Show[GraphicsArray[re3], 
      PlotLabel -> "Four replicate plots with 399 simulations", 
      Frame -> True]

re4 = re[4]

Show[GraphicsArray[re4], 
      PlotLabel -> "Four replicate plots with 999 simulations", 
      Frame -> True]

Four plots are quantile normal plots.  The two lines represent the 0.025 and 0.975 quantile estimates by  simulated data.  Since the data is simulated from standard normal distribution, two lines should be approximately through point (0, 1.96) and point (0, -1.96) respectively if simulations are enough to estimate quantile.  These graphs indicate that 39, 99 simulations are not enough to do a good job, 399 simulations do well and 999 simulations do much better.

Nonparametric Approximations for Variance and Bias

Delta Methods

In parametric analysis it is often possible to represent estimators T in terms of fundamental statistics (Subscript[U, 1,   ...  , ]) (Subscript[U, m, ]) such as sample moments, for which exact or approximate distributional calculations are relatively easy.  Then we can take advantage of the delta method to obtain distributional approximations for T itself.  For example, consider a scalar estimator T which is a function of the scalar statistic U based on a sample of size n, say T=g(U) .   Suppose that it is known that U is a consistency estimator of scalar parameter ζ and

\(TraditionalForm
    \`U\( ∼ \& . \)\(N(ζ, 
          \(n\^\(-1\)\) \(\(σ\^2\)(ζ)\)) . \)\)

If g(·) is a smooth function, the delta method result shows that provided that \(TraditionalForm\`\(\(g\& . \)(ζ) ≠ 0, \)\)

\(TraditionalForm
    \`\(T\( ∼ \& . \)N(θ, 
        \(n\^\(-1\)\) \({\(g\& . \)(ζ)}\^2\) 
          \(\(σ\^2\)(ζ)\)), \)\)

where θ = g(ζ);  the dot indicates differentation with respect to ζ. The principal feature is the delta method variance approximation:

\(TraditionalForm
    \`var {g(U)}\( = \& . \){
        \(\(\(g\& . \)(ζ)\)\^2\) \(var(U) . \)\)

To calculate the bias of T , we can use

\(TraditionalForm
    \`\(E(T)\( = \& . \)g(ζ) + 
        \(\(g\& . \)(ζ)\) \(E(U - ζ)\) + 
        \(1\/\(2  n\)\) \(\(g\& .. \)(ζ)\) 
          \(\(σ\^2\)(ζ)\); \)\)

or, if U is unbiased for ζ , then

\(TraditionalForm
      \`E(T)\( = \& . \)θ + 
          \(1\/\(2  n\)\) \(\(g\& .. \)(ζ)\) 
            \(\(σ\^2\)(ζ)\)\).

The following Exercise 9 extends above results to the case of vector U and vector T.

Exercise 10 ( 2.10 Problem 9 )

Suppose that the random variables U = {(Subscript[U, 1,   ...  ,  ]) (Subscript[U, m])) have means (Subscript[ζ, 1,   ...  ,  ]) (Subscript[ζ, m]) and covariances \(cov(\(\(U\_k\)\_, \) U\_l)\)=(n^(-1))(Subscript[ω, k l])(ζ), and that (Subscript[T, 1])=(Subscript[g, 1])(U),... ,(Subscript[T, q]) = (Subscript[g, q])(U) .

Suppose \(TraditionalForm
      \`T\_i = \(\(g\_i\)(U)\)\ \((\ i = 1, \  ... \ , q)\)\) is a differentiable function.  By Multivariate Taylor' s Theorem, the second-order Taylor series expansion of g about ζ={(Subscript[ζ, 1,   ...  , ]) (Subscript[ζ, m])} is

(Subscript[T, i])=(Subscript[g, i])(U)(OverDot[=])g(ζ) + 
              (Overscript[((∑_(k = 1))), m])
                  (( δ (Subscript[g, i k])(ζ))/(δ 
                        (Subscript[ζ, k]))) ((Subscript[U, i]) - (Subscript[ζ, i])) + 
              (1/2) (Overscript[( (∑_(k l = 1))), m]) 
                (((δ^2) 
                      (Subscript[g, i])(ζ))/(δ(Subscript[ζ, k]) 
                      δ(Subscript[ζ, l]))) ((Subscript[U, k]) - (Subscript[ζ, k])) 
                ((Subscript[U, l]) - (Subscript[ζ, l])).

Take expectations on both sides of above formula to get

\(TraditionalForm
    \`E(T\_i)\( = \& . \)g(ζ) + 
        1\/\(2  n\)\ 
          \(\(\ ∑\_\(k\ l = 1\)\)\&m\(\( ω\_\(k\ l\)\)(ζ)\) 
            \(δ\^2\) 
            \(\(\(g\_i\)(ζ)\)\/\(\(δζ\_k\) 
                  δζ\_l\) . \)\)\)

Using the first-order Taylor series expansion for g about ζ={(Subscript[ζ, 1,   ...  , ]) (Subscript[ζ, m])), that is,

(Subscript[T, i])=(Subscript[g, i])(U)(OverDot[=])g(ζ) + 
              (Overscript[((∑_(k = 1))), m])
                  (( δ (Subscript[g, i k])(ζ))/(δ 
                        (Subscript[ζ, k]))) ((Subscript[U, i]) - (Subscript[ζ, i])),

we can approximate the covariances between (Subscript[T, i])  by

\(TraditionalForm
      \`cov(T\_\(i, \)\ T\_j)\( = \& . \)\(n\^\(-1\)\) 
          \(\(∑\_\(k\ l = 1\)\)\&m\( ω\_\(k\ l\)\)(ζ)
              \)\)\(TraditionalForm
      \`\(\(δ\ \(\(g\_k\)(ζ)\)\)\/\(δ\ ζ\_k\)\) 
        δ\ \(\(g\_l\)(ζ)\)\/\(δ\ ζ\_l\)\).

In practice, we can use these results to estimate the approximate variance and bais.  For instance, suppose \(TraditionalForm\`\(x\&-\)\) and \(TraditionalForm\`\(y\&-\)\) are random variables with means (Subscript[μ, x]) and (Subscript[μ, y]) respectively.  The parametric function to be estimated is 
      g((Subscript[μ, x,  ]) (Subscript[μ, y])) = (Subscript[μ, x])/(Subscript[μ, y]) .   Then

\(TraditionalForm
    \`\(δg(\(μ\_\(\(x, \)\ \)\) μ\_y)\)\/δμ\_x = 
      \(1\/μ\_y\ \ \ \ \ and\ \ \ \ 
          \(\ \ \ δg(\(μ\_\(\(x, \)\ \)\) μ\_y)
              \)\/δμ\_y = \(-μ\_x\)\/\((μ\_y)\)\^2\)\)

By the approximated covariances shown above, we get

\(TraditionalForm
    \`var(t) = 
      \(var(\(x\&-\)\/\(y\&-\))
          \( = \& . \)\(1\/\((μ\_x)\)\^2\) \(var(\(x\&-\))\) + 
          \(\((μ\_x)\)\^2\/\((μ\_y)\)\^4\) \(var(\(y\&-\))\) - 
          \(\(2  μ\_x\)\/\((μ\_y)\)\^3\) 
            \(cov(\(x\&-\), \ \(y\&-\))\)\)\)

\(TraditionalForm
      \`\(\( = \& . \)\(n\^\(-2\)\) 
          \(\(∑\+\(j = 1\)\)\&n
              \((x\_j - \(t\_\(\  \)\) y\_j)\)\^2\/\(y\&-\)\^2\)\)\).

which is based on independent pairs ((Subscript[x, 1,  ]) (Subscript[y, 1])),  ...  , 
      ((Subscript[x, n, ]) (Subscript[y, n])) .

Influence function and nonparametric delta method

There is a simple variance approximation for many statistics T with the representation \(TraditionalForm\`t(\(F\&^\))\).  An extension of the Taylor series expansion for statistical functions  t(G) about true distribution function of population is

\(TraditionalForm
    \`\(t(G)\( = \& . \)t(F) + ∫\(\(L\_t\)(y; F)\) d\ \(G(y)\), 
    \)\)

where (Subscript[L, t]) ,  the first derivative of t(·) at F, is defined by

\(TraditionalForm
      \`\(L\_t\)(y; F) = 
        \(lim\+\(ϵ -> 0
                \)\(t {\((1 - ϵ)\) F + \ ϵ\ H\_y} - 
                t(F)\)\/ϵ = 
          \(\(δ\ 
                t {\((1 - ϵ)\) F + \ ϵ\ H\_y}
                  \)\/δϵ\( | \_\(ϵ = 0\)\)\)\)\),

with (Subscript[H, y])(u) = H(u - y) the Heaviside or unit step function jumping from 0 to 1 at u = y.

Set G = F to obtain

∫(Subscript[L, t])(y; F) d F(y) = 0.

The function (Subscript[L, t])(y) = (Subscript[L, t])(y; F) is called the influence function of T and its empirical approximation \(TraditionalForm\`l(y) = \(L\_t\)(y; \(F\&^\))\) is called the empirical influence function.  The particular values (Subscript[l, j]) = l((Subscript[y, j])) are called the empirical influence values.

Let G = \(TraditionalForm\`\(F\&^\)\), the first-oder approximation of t \(TraditionalForm\`\((\(F\&^\)\)\) ) is

\(TraditionalForm
    \`t(\(F\&^\))
        \( = \& . \)\(t(F) + 
          ∫\(\(L\_t\)(y; F)\) d\ \(\(F\&^\)(y)\) = 
        t(F) + \(1\/n\) 
            \(\(∑\+\(j = 1\)\)\&n\(\( L\_t\)(y\_j; F) . \)\)\)\)

By the central limit theorem , we get

\(TraditionalForm
      \`T - \(θ\ \ \)\( ∼ \&·\)\ \ 
        N(\ 0, \ \(v\_L\)(F))\),

because ∫(Subscript[L, t])(y; F) d F(y) = 0, 
       where

\(TraditionalForm
      \`\(v\_L\)(F) = 
        \(\(n\^\(-1\)\) var {\(L\_t\)(Y)} = 
          \(n\^\(-1\)\) \(∫\(\(L\_y\^2\)(y)\) d\ \(F(y)\)\)\)\).

Then the nonparametric delta method variance estimate is

\(TraditionalForm
    \`v\_L = \(\(v\_l\)(F\&^) = 
        \(n\^\(-2\)\) \(\(\ ∑\_\(j = 1\)\)\&n\( l\_j\^2 . \)\)\)\)

If t(F) = a {(Subscript[t, 1])(F),  ...  , (Subscript[t, m])(F)}, then

\(TraditionalForm
    \`\(L\_t\)(y) = 
      \(\(\ ∑\_\(i = 1\)\)\&m\(\( δ\ a\)\/\(δ\ t\_i\)\) 
        \(\(L\_\(t\_i\)\)(y) . \)\)\)

Example (Average)

Let\(TraditionalForm\`\(\ t = \(y\&-\)\)\), corresponding to the statistical function t(F) = ∫y d F(y).  To calculate the influence function of T, we write

\(TraditionalForm
      \`t {\((1 - ϵ)\) F + ϵ\ H\_y} = 
        \(∫y\ d\ {\((1 - ϵ)\) F + ϵ\ H\_y} = 
          \(\((1 - ϵ)\) \(∫y\ d\ F\) + 
              ϵ \(∫y\ d\ H\_y\) = 
            \((1 - ϵ)\) μ + ϵ\ y\)\)\),

and differentiate at ϵ = 0 to obtain

(Subscript[L, t])(y) = y - μ.

Then the empirical influence function is \(TraditionalForm\`\(\(l(y) = y - \(y\&-\), \)\ \)\)and the empirical influence values is \(TraditionalForm\`\(l\_j = y\_j - \(\(y\&-\) . \)\ \)\) Thus the delta method variance approximation is

(Subscript[v, l]) = (n - 1) (s^2)/(n^2),

where (s^2) is the unbaised sample variance of the (Subscript[y, j]) .

Exercises 11 ( 2.10 Problem 10 (a) )

To show that the influence function for a linear statistic t(F) = ∫a(x) d F(x), we calculate

\(TraditionalForm
      \`t {\((1 - ϵ)\) F + ϵ\ H\_y} = 
        \(∫\(a(x)\) d
              \ {\((1 - ϵ)\) F + ϵ\ H\_y} = 
          \((1 - ϵ)\) \(t(F)\) + ϵ\ \(a(y)\)\)\),

so influence function of t(F) is

\(TraditionalForm
      \`\(L\_t\)(y) = 
        \(\(δ {\((1 - ϵ)\) \(t(F)\) + ϵ\ \(a(y)\)}
                \)\/δϵ = 
          \(a(y) - t(F) = a(y) - E {a(y)}\)\)\).

For a sample moments about zero \(TraditionalForm
      \`\(\(\(μ\_r\)(F) = ∫\(x\^r\) d\ \(F(x)\), \)\ \)\)the influence functions are

(Subscript[L, (Subscript[μ, r])])(y) = (y^r) - (Subscript[μ, r])(F).

For the variance (Subscript[μ, 2])(F) - {(Subscript[μ, 1])(F)}^2,  
      the influence function is

\(TraditionalForm
      \`\(L\_v\)(y) = 
        \(\(L\_\(μ\_2\)\)(y) - 
            2 \(\( μ\_1\)(F)\) \(\(L\_\(μ\_1\)\)(y)\) = 
          y\^2 - \(μ\_2\)(F) - 
            2 \(\( μ\_1\)(F)\) \((y - \(μ\_1\)(F))\)\)\).

For the product moment correlation (Subscript[μ, r s]) = E((U^r) (X^s)), the influence function is

\(TraditionalForm
      \`\(L\_\(μ\_\(r\ s\)\)\)(u, y) = 
        \(\(u\ \)\^r\) y\^s - μ\_\(r\ s\)\).

Example ( correlation )

The sample correlation is the sample version of the product moment correlation, which for the pair Y = (U, X) can be defined in terms of (Subscript[μ, r s]) = E((U^r) (X^s))  by

\(TraditionalForm
      \`ρ = 
        \(ρ(F) = 
          \(\(μ\_\(11 - \)\) \(μ\_10\) 
              μ\_01\)\/{
                \((μ\_20 - \((μ\_10)\)\^2)\) 
                  \((μ\_02 - \((μ\_01)\)\^2)\)}\^\(1/2\)\)\).

The results of Exercises 12 and the partial derivatives of ρ (·) with respect to the μs lead to the influence function for the correlation coefficient,

\(TraditionalForm
    \`\(\(L\_ρ\)(u, x) = 
      \(u\_s\) x\_s - \(1\/2\) \(ρ(u\_s\^2 + x\_s\^2)\), \)\)

where (Subscript[u, s]) = 
                (u - (Subscript[u, 10]))/((((Subscript[u, 20]) - ((((Subscript[u, 10])))^2)))^(1/2)),  and 
              (Subscript[x, s])=(x - (Subscript[μ, 01]))/((Subscript[μ, 02]) - ((((Subscript[μ, 01])))^2)) are standardized variances.

Jackknife estimates

Another approach to approximating the influence function, but only at the sample values (Subscript[y, 1]),  ...  , (Subscript[y, n]) themselves,  is the 
      jackknife. Here (Subscript[l, j]) is approximated by

(Subscript[l, jack,  j]) = (n - 1) (t - (Subscript[t, -j])),

where (Subscript[t, -j]) is the estimate calculated with (Subscript[y, j]) omitted from the data.

The jackknife approximations to the bias and variance of T are

\(TraditionalForm
      \`b\_jack = 
        \(\(-1\)\/n\) \(\(∑\+\(j = 1\)\)\&n l\_\(jack, \ j\)\), 
      \ \ \ \ \ \ \ \ \ \ 
      v\_jack = 
        \(1\/\(n(n - 1)\)\) 
          \((∑\(l\^2\)\_\(jack, j\) - n\ \(b\^2\)\_jack\)\)).

Empirical influence values via regression

The approximation \(TraditionalForm
      \`t(G)\( = \& . \)t(F) + ∫\(\(L\_t\)(y; F)\) d\ \(G(y)\)\) can be applied to the bootstrap estimate \(TraditionalForm\`\(\(T\^*\) . \)\)  If the EDF of the bootstrap sample is denoted by \(TraditionalForm\`\(F\&^\^*\), \ then\)

\(TraditionalForm
    \`\(t(\(F\&^\^*\))\( = \& . \)t(F\&^) + 
        \(1\/n\) \(\(\ ∑\_\(j = 1\)\)\&n\( L\_t\)(\(y\_j\^*\); F\&^)\), 
    \)\)

or in simpler notation

\(TraditionalForm
    \`\(\(t\^*\) = 
      \(t + \(1\/n\) \(\(∑\_\(j = 1\)\)\&n\( f\_j\^*\) l\_j\) = 
        \(t\_L\^*\)\), \)\)

say, where \(TraditionalForm\`\(\(f\_j\^*\)\ \)\)is the number of times that \(TraditionalForm\`\(\(y\_j\^*\)\ \)\)equals (Subscript[y, j]),  for j = 1,  ...  , n .

Suppose that we have simulated R samples from \(TraditionalForm\`F\&^\).  Define \(TraditionalForm\`\(\(f\_\(r\ j\)\)\&*\)\) to be the frequency with which the data value (Subscript[y, j]) occurs in the rth bootstrap sample.  Then we have

\(TraditionalForm
    \`\(t\_r\^*\)\( = \& . \)t + 
        \(1\/n\) \(\(\ ∑\_\(j = 1\)\)\&n\(\( f\_\(r\ j\)\)\&*\)\ l\_j\), 
    \ \ \ \ \ \ \ \ \ r = 1,  ... \ , R\)

Exercises 12 ( 2.11 Practical 5)

We compare variance approximations for the correlation in data cd4 from Exercises 10.

Figure 12.1

Left panel of the graph shown above plots Jackknife estimated influence values of correlation against the empirical influence values.  Right panel of this graph plots regression estimated influence values of correlation against the empirical influence values. They look like 45° line. This implies that there is no big difference between these estimated influence values.  Next let see the estimated variances by nonparametric delta method , Jackknife, regression, and bootstrap simulation. The nonparametric delta method variance equals to 0.006.  The Jackknife approach variance equals to 0.007.  Regression variance is 0.008.  The bootstrap variance is 0.008.