- Email: [email protected]
Fast approximately balanced bootstrap without construction
ARTICLE IN PRESS
Statistics & Probability Letters 76 (2006) 1861–1872 www.elsevier.com/locate/stapro
Fast approximately balanced bootstrap without construction$ C. Devon Lina, Wilson W. Lub, R.R. Sittera, a
Department of Statistics and Actuarial Science, Simon Fraser University, Burnaby, BC, Canada V5A 1S6 b Department of Mathematics and Statistics, Acadia University, Wolfville, NS, Canada B4P 2R6 Received 3 March 2006; accepted 24 April 2006 Available online 5 June 2006
Abstract We develop second-order balanced and approximately balanced resamples that retain the nice properties demonstrated by Graham et al. (1990. Balanced design of bootstrap simulations. J. Roy. Statist. Soc. Ser. B 52, 185–202) but require essentially no design construction. r 2006 Elsevier B.V. All rights reserved. Keywords: Balanced incomplete block design; Balanced repeated replications; Orthogonal array; Stratified sample
1. Introduction Bootstrap simulations (resamples) are often used to approximate aspects of the sampling distribution, such as variance, of estimators when the form of the underlying distribution is unknown (see Efron and Tibshirani, 1993). The number of resamples required so as to obtain an accurate approximation can be large. This can be an important issue for complex estimation situations. There have been various methods proposed for reducing the required number of resamples needed to obtain a fixed level of accuracy. These include Davison et al. (1986), Graham et al. (1990), Efron (1990) and Nigam and Rao (1996). Balanced bootstrapping uses constructions of balanced arrays, balanced incomplete block designs and orthogonal latin squares, to create second-order balanced bootstrap samples that reduce the Monte Carlo (simulation) error and thus require fewer resamples to attain the same level of accuracy. The main impedence on their general usage has been the complexity of construction of such arrays. In this article, we give three strategies for quick and automatic construction of (nearly) second-order balanced arrays that are simple to use and attain significant reductions in Monte Carlo error and thus significantly reduce the number of required resamples to achieve the same level of accuracy. We will use similar notation as and a similar setup to Graham et al. (1990) to aid the reader. In general, a random sample x ¼ ðx1 ; . . . ; xn Þ0 is drawn from a distribution F ðÞ and the estimate T ¼ tðxÞ of y is obtained. The interest lies in approximating the distributional properties of T or some related quantities. Of particular interest are V ðTÞ and PðTpaÞ for fixed a. This can be done by substituting some estimate F^ for F. When no $
Supported by a grant from the Natural Sciences and Engineering Research Council of Canada.
Corresponding author. Tel.: +1 604 291 3334; fax: +1 604 291 4368.
E-mail address: [email protected] (R.R. Sitter). 0167-7152/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2006.04.035
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1862
assumptions on F are made, this is typically the empirical distribution function F^ ðxÞ ¼
Pn
j¼1 I ½xj px =n,
where ^ I ½xj px ¼ 1 if xj px and ¼ 0 otherwise. Since exact calculation of properties of T under F is typically not possible, the bootstrap uses simulation to approximate them. This is done by generating B random resamples : xb ¼ ðxb1 ; . . . ; xbn Þ0 drawn from F^ and calculating T b ¼ tðxb Þ for b ¼ 1; . . . ; B. Since PrðTpajF Þ¼PrðTpajF^ Þ P : : and V ðTjF Þ¼V ðTjF^ Þ, we can approximate via simulation PrðTpajF^ Þ¼ b I ½T b pa =B and :P V ðTjF^ Þ¼ b ðT b TÞ2 =B. The issue in this article is the Monte Carlo error in using these approximations. We will restrict to estimators T ¼ tðF^ Þ of y ¼ tðF Þ that can be expanded as 1 T ¼ tðF^ Þ ¼ tðF Þ þ n
n X
Lðxj ; F Þ þ
j¼1
n X n 1 X Qðxj ; xk ; F Þ þ , 2n2 j¼1 k¼1
where L and Q are, respectively, the first and second functional derivatives of t (see Hinkley and Wei (1984)). Note that L is typically called the influence function. If x is obtained by random sampling from F^ , a similar expansion is available for T ¼ tðx Þ by replacing F^ and F by F^ and F^ , respectively, where F^ is the empirical distribution function of x . That is, T ¼ T þ
n n X n 1X 1 X f j L^ j þ 2 f f Q^ þ , n j¼1 2n j¼1 k¼1 j k jk
P where L^ j ¼ Lðxj ; F^ Þ, Q^ jk ¼ Qðxj ; xk ; F^ Þ, and f j is the frequency of xj in x . Here j L^ j ¼ 0. It has been observed that the differences between various bootstrap simulation methods are embedded in the properties of the frequencies f bj . 2. Balanced bootstrap Efron (1979) notes that f ¼ ðf 1 ; . . . ; f n Þ0 follows the n-category uniform multinomial P distribution. Davison et al. (1986) introduces the notion of a first-order balanced bootstrap by forcing Bb¼1 f bj ¼ B via use of the permutation bootstrap, whereby the B n matrix of f bj follow a hypergeometric distribution with row sums n and column sums B. Graham et al. (1990) observe that if T is symmetric in the x’s, the above expansions show that the mean and variance of T can be approximated without simulation error up to terms of order 1=n provided the f are second-order balanced. They go on to use balanced incomplete block designs to construct B ¼ kn sets f such that B f X bj ¼ 1, B b¼1
(2.1)
B f 2 X bj ¼ ð2n 1Þ=n ¼ 2 1=n B b¼1 B f f X bj bk ¼ ðn 1Þ=n ¼ 1 1=n B b¼1
: ð¼2Þ,
: ð¼1Þ,
(2.2)
(2.3)
where by approximately equal we mean to order n1 . These are in agreement with the multinomial ^ ^ expectations, Eð f j jF^ Þ ¼ 1, Eð f 2 j jF Þ and Eð f j f k jF Þ. To see this, note that the balanced bootstrap variance estimator of T is 2 #2 , B B n X X X : 2 4 ðf 1ÞL^ j =n vbb ðTÞ ¼ ðT b TÞ =B¼ B bj b¼1
¼ n2
b¼1 n X j¼1
B X
j¼1
n B X 1 1X L^ 2j L^ j L^ k ðf bj 1Þ2 þ n2 ðf bj 1Þðf bk 1Þ. B b¼1 B jak b¼1
ð2:4Þ
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1863
P P Under conditions (2.2) and (2.3), b ðf bj 1Þ2 =B ¼ 1 1=n and b ðf bj 1Þðf bk 1Þ=B ¼ 1=n, and thus P 2 : vbb ðTÞ¼ j L^ j =n2 , which is the leading term of Eðvboot ðTÞjF^ Þ under multinomial resampling. This implies that the Monte Carlo variance of vbb ðTÞ is reduced to lower order. P The simplest example is T ¼ x¯ ¼ nj¼1 xj =n. In this case L^ j ¼ xj x¯ and under conditions (2.2) and (2.3), P : vbb ðxÞ ¯ ¼ ðn 1Þs2 =n2 ¼s2 =n, where s2 ¼ j ðxj xÞ ¯ 2 =ðn 1Þ. Graham et al. (1990) also show via simulation that the balanced bootstrap greatly reduces the simulation variability for percentile estimates when 0:05ppp0:95. Unfortunately, balanced bootstrapping is not routinely used in practice. This is almost certainly due to the difficulty in constructing balanced sets of f bj . In the next section, we will introduce three very simple and automatic methods for obtaining such f bj that satisfy the balance conditions or approximately satisfy them so as to easily reduce the simulation error in the bootstrap without any knowledge of design construction techniques. 3. Fast approximately balanced bootstrap Before introducing the three methods for creating fast approximately balanced bootstrap (FABB) resamples, let us consider Hadamard matrices, which will form the basis for all of the methods. A Hadamard matrix H ¼ fdjk g is a B B matrix of 1’s and 1’s, where B is an integer multiple of 4. It can always be rewritten column of 1’s and a row of 1’s, termed standard form. The elements, djk satisfy Pto have aP P d ¼ d ¼ 0, j jk k jk j djk djl ¼ 0 for kal. Table 1 gives an 8 8 example where 1’s are represented by þ and . Hadamard matrices have been well-studied and there exist many algorithms for their construction. Hedayat et al. (1999) give eight algorithms that can be used to obtain B B Hadamard matrices for all B a multiple of 4 up to size 200 200. A simple Fortran program that uses these eight algorithms to yield a B B Hadamard matrix (in standard form) for every B, a multiple of 4 up to 200 can be obtained from the authors. 3.1. The fold-over method It will turn out to be convenient to extend the definitions for the balancing conditions (2.1)–(2.3) to allow for varying resample sizes. To do so, let nb denote the size of the bth resample, and alter the balance conditions to be B f X n bj ¼ 1, (3.1) n B b b¼1 B X f 2 bj ¼ 2, ðn=nb Þ2 B b¼1
(3.2)
Table 1 An 8 8 Hadamard matrix
1 2 3 4 5 6 7 8
y1
y2
y3
y4
y5
y6
y7
y8
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ þ þ þ þ
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1864
B X f bj f bk ¼ 1 1=ðn 1Þ ðn=nb Þ2 B b¼1
: ð¼1Þ.
(3.3)
These are related to the BOMA conditions of Sitter (1993) for obtaining balanced half-samples in stratified random sampling. It is not difficult to see that a similar development to that of (2.4) is available under conditions (3.1)–(3.3) to :P yield vbb ðTÞ¼ j L^ 2j =½nðn 1Þ. In the simplest example where T ¼ x, ¯ vbb ðxÞ ¯ ¼ s2 =n, instead of ðn 1Þs2 =n2 as with the usual balanced bootstrap. To develop a fast method for constructing a set of B resamples that satisfy (3.1)–(3.3), consider the matrix in Table 2. This matrix was obtained by taking the 8-run Hadamard matrix in Table 1, folding it over (adding eight rows by changing signs) and then removing the row of +’s and the row of ’s. Note that the 14 remaining sets represented by keeping the column numbers corresponding to +’s in each row form a resolvable BIBD. Now, if we let each row define a resample by keeping two copies of each yj which has a þ below it, then each row defines a resample of size n. For example, row 2 defines the resample fy1 ; y1 ; y6 ; y6 ; y7 ; y7 ; y8 ; y8 g. Thus, it is not difficult to show that this matrix defines 14 resamples such that the corresponding frequencies satisfy (3.1)–(3.3). That is, if we let dbj denote the bjth element of the matrix, define f bj ¼ dbj þ 1. It turns out that this strategy is quite general. In fact, one can generalize it to approximately satisfy conditions (3.1)–(3.3) for any n via the following theorem. Theorem 1. For n þ t ¼ n0 ¼ 4m, 0ptp3, B ¼ 2n0 2, the two-way array X ðB; nÞ, built by folding-over an n0 n0 Hadamard matrix in standard form, and then removing t columns and the row of þ1’s and the row of 1’s, satisfies (3.1)–(3.3) exactly for t ¼ 0 and for t ¼ 1, 2 or 3, (1) (2) (3)
PB
2 b¼1 ðn=nb Þf bj =B ¼ 1 þ Oð1=n Þ (exact for t ¼ 1 2 2 2 b¼1 ðn=nb Þ f bj =B ¼ 2 þ Oð1=n Þ, PB 2 2 b¼1 ðn=nb Þ f bj f bk =B ¼ 1 1=ðn 1Þ þ Oð1=n Þ.
PB
and 2);
Proof. See appendix. Thus, using the fast simple algorithms for construction of Hadamard matrices combined with Theorem 1 : yields a simple and easy way to construct balanced and approximately balanced bootstraps with B¼2n. Table 2 Construction of a 14 8 array
1 2 3 4 5 6 7 8 9 10 11 12 13 14
y1
y2
y3
y4
y5
y6
y7
y8
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1865
3.2. The Hadamard method It turns out that one can gain flexibility in the choice of B by sacrificing a bit in the approximations of Theorem 1. It is clear from conditions (3.2) and (3.3) and approximation (2.4) that the 1=ðn 1Þ term in condition (3.3) is not necessary to capture the leading term of the variance. It is merely kept so that the approximation will be exact for linear statistics such as x. ¯ If we are willing to sacrifice this property, it turns out there is no need to fold-over, stated as follows. Theorem 2. For n þ t ¼ n0 ¼ 4m, 1ptp4, B ¼ n0 1, the two-way array X ðB; nÞ, built by choosing an n0 n0 Hadamard matrix in standard form, and then removing t columns and the row of þ1’s, satisfies the following conditions for t ¼ 1, 2, 3 or 4: (1) (2) (3)
PB
b¼1 ðn=nb Þf bj =B ¼ 1 þ Oð1=nÞ; 2 2 b¼1 ðn=nb Þ f bj =B ¼ 2 þ Oð1=nÞ; PB 2 b¼1 ðn=nb Þ f bj f bk =B ¼ 1 þ Oð1=nÞ.
PB
Proof. See appendix. 3.3. The random-pairing method The fold-over method requires approximately 2n resamples. The Hadamard method requires approximately n resamples. This makes them particularly useful when n is small to moderate. Once n becomes larger, one would likely prefer to use a regular bootstrap with Bon. In this section, we develop a third method that requires approximately n=2 replicates. It is related to some random methods for balanced repeated replication (BRR) variance estimates used in surveys (see Wolter, 1985). For the moment, assume n is even. In our context, the idea would reduce to randomly ordering y1 ; . . . ; yn , denoted y1 ; . . . ; yn and then forming successive pairs ðy1 ; y2 Þ; ðy3 ; y4 Þ; . . . ; ðyn1 ; yn Þ. We now use something akin to the BRR method treating these pairs as strata with two units in each stratum. That is, we form a B ðn=2Þ matrix, H ¼ fdbk g, by removing columns : from a B B Hadamard matrix, where B is the smallest multiple of 4 strictly larger than n=2, so B¼n=2. The rows of this matrix define the f bj as follows: f bð2k1Þ ¼ 1 þ dbk and f bð2kÞ ¼ 1 dbk , for k ¼ 1; . . . ; n=2. If this method is used, nb ¼ nPand it is easily seen that conditions (3.1) and (3.2) are satisfied. However, condition (3.3) is not satisfied. Instead, Bb¼1 ðn=nb Þ2 f bj f bk =B ¼ 0 and 1 for j and k paired together or not, respectively. Thus, this method does not remove the leading term from the Monte Carlo variance of vbb ðTÞ. However, the following theorem implies that the leading term of the Monte Carlo variance is reduced substantially. Theorem 3. Using the random-pairing method, the leading term of V ðvbb Þ=V ðvb Þ, where B ¼ n=2 þ t, satisfies V ðvbb Þ nðn þ 2tÞ½n^g þ ðn2 3n þ 3Þ nðn 2Þ2 ðn þ 2tÞ p , ¼ V ðvb Þ ðn 1Þðn 3Þ½n^g þ ð2n 3Þðn 1Þ 2ðn 1Þ3 ðn 3Þ
(3.4)
where V P ðvbb Þ and V ðvb Þ are the Monte Carlo variances of vbb and vb , respectively, and g^ ¼ ðn 1Þ P ð i L^ 4i Þ=ð i L^ 2i Þ2 is an estimate of the kurtosis. Proof. See appendix. For all t ¼ 1; . . . ; 4, the right-hand side of (3.4) is decreasing in n to a limit of 0.5, summarized in Table 3. Of course, this is an upper bound and depending upon the sample value of g^ the random pairing method can reduce the Monte Carlo variance much more. For example, if g^ 4ðn 1Þ=2 the upper bound in (3.4) becomes nðn 2Þðn þ 2tÞ=½ðn 1Þ2 ð5n 6Þ and for moderate sized n is nearer 0:2 than 0.5. We explore the actual performance through some limited simulations in the next section. Extending to n odd can be accomplished fairly easily by randomly pairing n 3 units, and then using the remaining three units to randomly create two strata that have one unit in common. This adequately approximates the even n case.
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1866
Table 3 The bound of ratio of Monte Carlo variance of random pairing balanced bootstrap relative to the regular bootstrap n
20
30
40
50
t¼1 t¼2 t¼3 t¼4
0.61 0.67 0.72 0.78
0.57 0.61 0.64 0.68
0.55 0.58 0.61 0.63
0.54 0.56 0.58 0.60
4. Simulations 4.1. Percentiles for the sample mean and t-statistic In this section, we investigate the Monte Carlo error when estimating percentiles rather than moments of the sample average in a similar fashion to that used by Graham et al. (1990, Example 9). Consider the skew sample of n ¼ 22 x values 8:0; 9:6; 10:0; 10:4; 11:7; 13:0; 14:0; 15:0; 15:8; 16:6; 16:9; 17:2; 17:25; 17:3; 19:55; 21:8; 22:9; 24:0; 25:45; 26:9; 30:35; 33:8: To evaluate the bootstrap percentile estimates for x¯ and the Student t-statistic t ¼ n1=2 ðx¯ mÞ=s, we obtain x¯ 1 px¯ 2 p px¯ B for each of the three FABBs. Since the three FABBs use differing resample sizes, we will compare them to the usual bootstrap using the same resample size. Let x¯ p be the 100pth percentile of x¯ 1 px¯ 2 p px¯ B . The simulation error is then ep ¼ x¯ p x¯ p ,
(4.1)
where P ðx¯ px¯ p Þ ¼ p and P denotes probability under the resampling distribution. The values of x¯ p and tp are obtained via an independent simulation taking K ¼ 10; 000 resamples. For each of the three FABBs, we repeat the resampling 1000 times, independently, and similarly for the usual bootstrap. The resulting 1000 errors ep in (4.1) are summarized by their mean and variance for p ¼ i=ðB þ 1Þ for i ¼ 1; . . . ; B and plotted against the corresponding exact percentiles in Fig. 1 (for the mean). The corresponding figure for the t-statistic (not shown) is qualitatively the same. Clearly the fold-over balanced bootstrap yields significant reductions in simulation variability for percentile estimates and maintains near unbiasedness. These results parallel those in Graham et al. (1990), as do similar simulations for the correlation coefficient (not shown). The Hadamard method performs nearly as well. The random pairing method also gives significant gains. 4.2. Variance estimates for the sample mean and sample correlation To study the simulation variance of the FABB variance estimators relative to the usual bootstrap we conducted a limited simulation study employing the following model: 1=2
xi ¼ bzi þ zi i ,
(4.2)
where i Nð0; s2 Þ and independent of zi gammaðg; hÞ. Thus, the mean and variance of zi are given by mz ¼ gh and s2z ¼ gh2 , respectively. Further, the mean and variance of xi are mx ¼ bmz and s2x ¼ b2 s2z þ mz s2 , and corrðxi ; zi Þ ¼ r ¼ bsz =sx . We fixed b ¼ 1:0, mz ¼ 100 and s ¼ 30. For various settings of s2z to vary r, and various values of n, we generated N 1 ¼ 1000 independent random samples from model (4.2) with sample size n. For each such sample we repeatedly (N 2 ¼ 1000 times) obtained the fold-over balanced bootstrap variance estimator, vfobb (resample : size B ¼ 2n0 2¼2n), the Hadamard balanced bootstrap variance estimator, vhbb (resample size
ARTICLE IN PRESS standard deviation of error
C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
average error
0.6 0.4 0.2 0.0 −0.2 −0.6 0.0
0.2
(a)
0.4
0.6
0.8
0.2 0.0 −0.2 −0.6
(c)
0.4
0.6
0.8
0.2 0.0 −0.2 −0.6
(e)
0.6
exact percentle
0.0 0.2
0.8
1.0
0.6
0.8
1.0
0.8
1.0
0.8
1.0
1.0 0.8 0.6 0.4 0.2 0.0 0.2
0.4
0.6
exact percentle
1.0 0.8 0.6 0.4 0.2 0.0 0.0
(f)
0.4
exact percentle
0.0
standard deviation of error
average error
0.4
0.4
0.2
(d)
0.6
0.2
0.4
1.0
exact percentle
0.0
0.6
0.0
standard deviation of error
average error
0.4
0.2
0.8
(b)
0.6
0.0
1.0
1.0
exact percentle
1867
0.2
0.4
0.6
exact percentle
Fig. 1. Averages and (b) standard deviations of 1000 simulation errors in bootstrap percentile estimates for x: ¯ , regular bootstrap; , fold-over balanced bootstrap. (c) and (d), and (e) and (f), are similar plots for the Hadamard and random-pairing balanced bootstraps, respectively.
: : B ¼ n0 1¼n), the random-pairing balanced bootstrap variance estimator, vrpbb (resample size B ¼ n0 =2¼n=2) and the usual bootstrap variance estimator, vb (resample size B chosen to be the same as that of the balanced ^ bootstrap being compared to), for T 1 ¼ x, ¯ T 2 ¼ z¯ and T 3 ¼ r. If we let E 1 and V 1 denote the expectation and variance under the sampling from the model, and let E 2 and V 2 denote the expectation and variance under the resampling conditional on each sample, respectively, V ðvÞ ¼ V 1 E 2 ðvÞ þ E 1 V 2 ðvÞ for each variance estimator v, where the second term is the average Monte Carlo variance. To estimate the average Monte Carlo variance for the three approximately balanced bootstrap variance estimators and for the P usual bootstrap, for each of the N 1 samples we estimate the conditional Monte Carlo ~ ðvÞ ¼ N 2 ðvl v¯ Þ2 =N 2 , where v denotes one of the four bootstrap variance estimators and variance using V l¼1 P 2 v¯ ¼ N l¼1 vl =N 2 . This yields N 1 estimates of the conditional simulation variance, V 2 , of each of the three bootstrap variance estimators for T 1 , T 2 and T 3 . We then take the average of these N 1 values as an estimate of the average Monte Carlo variance, denoted EVðvÞ. Table 4 gives REVfobb ¼ EV½vfobb ðT i Þ=EV½vb ðT i Þ, REVhbb ¼ EV½vhbb ðT i Þ=EV½vb ðT i Þ and REVrpbb ¼ EV½vrpbb ðT i Þ=EV½vb ðT i Þ for i ¼ 1, 2 and 3, n ¼ 22, 24, 50, and 52, and r ¼ 0:2, 0.5, and 0.8, where the B in the usual bootstrap was chosen to match the balanced bootstrap it is being compared to in each ratio. Viewing the results for T 1 ¼ x¯ and T 2 ¼ z¯ in Table 4, we see that REVðvfobb Þ ¼ 0 for n ¼ 24 and n ¼ 52 verifying Case 1 of Theorem 1, and are extremely small for n ¼ 22 and 50 verifying Case 3 of Theorem 1 (Note, the n ¼ 50 case yielded non-zero results if viewed to six decimal places). The Hadamard method performed nearly as well. Viewing the results for the random-pairing method, we see that the Monte Carlo variance is almost cut in half in the worst cases, and is reduced far more in many cases. Note that Eð^gÞ was
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1868
Table 4 Ratio of average Monte Carlo variances of balanced bootstraps relative to the regular bootstrap for x, ¯ z¯ and r^ n
r ¼ 0:2
r ¼ 0:5
r ¼ 0:8
T1
T2
T3
T1
T2
T3
T1
T2
T3
REVfobb
22 24 50 52
0.0003 0.0000 0.0000 0.0000
0.0030 0.0000 0.0000 0.0000
0.1129 0.0909 0.0462 0.0420
0.0002 0.0000 0.0000 0.0000
0.0002 0.0000 0.0000 0.0000
0.1120 0.1251 0.1174 0.1179
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.1202 0.1268 0.1969 0.1949
REVhbb
22 24 50 52
0.0113 0.0262 0.0018 0.0048
0.0108 0.0256 0.0018 0.0047
0.3925 0.3767 0.2110 0.2025
0.0071 0.0179 0.0014 0.0038
0.0056 0.0152 0.0013 0.0035
0.3114 0.3592 0.3245 0.3308
0.0016 0.0065 0.0006 0.0013
0.0012 0.0065 0.0006 0.0012
0.2857 0.3039 0.4150 0.4140
REVrpbb
22 24 50 52
0.4927 0.5979 0.5332 0.5115
0.4700 0.5861 0.5188 0.5011
0.7981 0.8615 0.7416 0.7034
0.3083 0.4022 0.4118 0.4001
0.2452 0.3424 0.3726 0.3659
0.5776 0.6984 0.7111 0.7083
0.0699 0.1417 0.1661 0.1330
0.0538 0.1389 0.1700 0.1194
0.5001 0.5531 0.7300 0.7131
chosen quite large for r ¼ 0.5 and 0.8 for z¯ and r ¼ 0.8 for x¯ (overall ranging from 3 to 20). Roughly speaking, one would need to use a resample size about 2–20 times larger with vb to attain the same Monte Carlo variance as that of vrpbb in these situations. ^ we see that the gains are not restricted to linear estimators but are almost as Viewing results for T 3 ¼ r, dramatic for this nonlinear statistic. For r^ it does appear that obtaining the closer approximation by folding over has a larger relative impact. In this case, one would need to use a resample size about 8–25 times larger with vb to attain the same simulation variance as that of vfobb , about 3–10 times larger for vhbb and 1.2 or 2 times larger for vrpbb . We repeated the process using n ¼ 21, 23, 49 and 51 by removing one observation and using the method described at the end of Section 3 for the odd n case for the random-pairing method. The results were qualitatively the same (not shown). 5. Concluding remarks We develop three methods for obtaining approximately balanced bootstrap resamples that retain the advantages in performance of second-order balancing without a need for tedious and difficult design construction techniques. Acknowledgement We thank D. Bingham for helpful comments. Appendix Proof of Theorem 1. For H n0 ¼ fdjk gn0 n0 a Hadamard matrix in standard form: (H1) Each column (row) of H n0 , not including the column (row) of 1’s, consists of n0 =2 ¼ 2m, þ1’s and 2m; 1’s; (H2) For two distinct columns (rows) of H n0 , u ¼ ðu1 ; . . . ; un0 Þ0 and v ¼ ðv1 ; . . . ; vn0 Þ0 , each of the pairs ðui ; vi Þ ¼ ðþ1; þ1Þ, ðþ1; 1Þ, ð1; þ1Þ and ð1; 1Þ appears m times.
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1869
(H3) Since H n0 is a strength 2 orthogonal array, its fold-over, H~ n0 will be strength 3, that is, the rows of any three columns of H~ n0 ¼ ðH Tn0 ; H Tn0 ÞT will consist of an equal number of copies of the eight possible combinations of ð1; 1; 1Þ. From (H1)–(H3), B X
f bj
0
¼ 2n 2;
b¼1
n0 1 X
f bj f bk ¼ n0 4 and
b¼1
In Pnaddition, for each replicate b, nb ¼ x are j¼1 bj
8 0 > > > < 1 Db ¼ > 0; 2 > > : 1; 3
Pn
j¼1 f bj
B X
f bj f bk ¼ n0 .
(A.1)
b¼n0
¼nþ
Pn
j¼1 xbj
¼ n þ Db , where the possible values of Db ¼
for t ¼ 0; for t ¼ 1; for t ¼ 2;
(A.2)
for t ¼ 3:
Therefore, we have the approximation, n=nb ¼ nðn þ Db Þ1 ¼ 1 Db =n þ Oðn2 Þ, and similarly, ðn=nb Þ2 ¼ 1 2Db =n þ Oðn2 Þ. Due to this, together with (A.1), the left-hand sides of conditions (3.1), (3.2) and (3.3) can be expressed as B f X n bj b¼1
nb
B B 1X 1 X 1 ¼ f bj Db f bj þ O 2 B b¼1 nB b¼1 n B B 1 X 1 ¼1 Db f bj þ O 2 , nB b¼1 n
B 2 f 2 X n bj b¼1
nb
B B 2X 4 X 1 ¼ f bj Db f bj þ O 2 B b¼1 nB b¼1 n B B 4 X 1 ¼2 Db f bj þ O 2 , nB b¼1 n
ðA:3Þ
ðA:4Þ
and B 2 f f X n bj bk b¼1
nb
B
B B 1X 2 X 1 ¼ f f Db f bj f bk þ O 2 B b¼1 bj bk nB b¼1 n B 1 2 X 1 ¼1 Db f bj f bk þ O 2 , n 1 nB b¼1 n
ðA:5Þ
respectively. P P It only remains to prove that both Bb¼1 Db f bj and Bb¼1 Db f bj f bk are terms of Oð1Þ for each of t ¼ 0, 1, 2 and 3. P P For t ¼ 0, we have Bb¼1 Db f bj ¼ Bb¼1 Db f bj f bk ¼ 0, since Db ¼ 0 by (A.2), and conditions (3.1)–(3.3) are thus exactly satisfied. To consider t ¼ 1, 2 and 3, we first note the relationships between any replicate b and its fold-over b0 ; Db ¼ Db0 and xbj ¼ xb0 j , so that f bj f bk ¼ ð1 þ xbj Þð1 þ xbk Þ ¼ ð1 xb0 j Þð1 xb0 k Þ ¼ ð1 þ xb0 j Þð1 þ xb0 k Þ 2ðxb0 j þ xb0 k Þ ¼ f b0 j f b0 k 2ðxb0 j þ xb0 k Þ.
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1870
Thus, B X
Db f bj ¼
b¼1
n0 1 X
Db ð1 þ xbj Þ þ
n0 1 X
ðDb þ Db xbj Þ þ
b¼1
¼2
Db0 ð1 þ xb0 j Þ
b0 ¼n0
b¼1
¼
B X
n0 1 X
½ðDb Þ þ ðDb Þðxbj Þ
b¼1
n0 1 X
Db xbj ,
b¼1
and B X
Db f bj f bk ¼
b¼1
n0 1 X
Db f bj f bk þ
¼
Db f bj f bk þ
b¼1
¼2
Db0 f b0 j f b0 k
b0 ¼n0
b¼1 n0 1 X
B X
n0 1 X
n0 1 X
ðDb Þ½f bj f bk 2ðxbj þ xbk Þ
b¼1
Db xbj þ 2
b¼1
n0 1 X
Db xbk .
b¼1
P 0 1 Therefore, we need only prove that nb¼1 Db xbj ¼ Oð1Þ. Without loss of generality, we assume that in all cases, t ¼ 1, 2 and 3, the column consisting of þ1’s for its first n0 1Pelements and 1’s Pn0 1for the rest, is deleted. 0 1 Case 1: t ¼ 1. Since Db 1 for b ¼ 1; . . . ; n0 1, we have nb¼1 Db xbj ¼ b¼1 xbj ¼ 1. Case 2: t ¼ 2. D ¼ 2 if ‘‘þþ’’ is deleted from replicate b; or ¼ 0 if ‘‘þ’’ is deleted. Thus, P b Pn0 1 D x ¼ 2 x ¼ 2ð1Þ ¼ 2, where the second to last equality follows from (H2). b¼1 b bj b2fb:Db ¼2g bj Case 3: t ¼ 3. Db ¼ 3; 1; 1; 1, corresponding to the four patterns of deletions, ‘‘þ þ þ’’, ‘‘ þ þ’’, ‘‘þ þ’’ and ‘‘ þ’’. DenoteP A1 , A2 , A3 and A4 as collections of replicates associated with those four patterns, respectively. Let S ij ¼ b2Ai xbj for i ¼ 1; 2; 3; 4 and Pn0j1¼ 1; . . . ; n. From (H2) it is easy to see that S1j þ S 2j ¼ 1; S 1j þ S3j ¼ 1; S 2j þ S 4j ¼ 0. Therefore, b¼1 Db xbj ¼ 2ðS 1j þ S 2j Þ ðS 1j þ S 3j Þ þ ðS 2j þ S4j Þ ¼ 3. & PB PB Proof Theorem 2. For the Hadamard method, we have b¼1 xbj =B ¼ 1, b¼1 f bj ¼ 1 1=B, and PB of f f =B ¼ 1 3=B. For b ¼ 1; . . . ; B, and t ¼ 1; 2; 3; 4, note that jD jp4. Thus, b b¼1 bj bk 4 1 X 4 1 X B B 1 D f p f o ¼ O , nB b¼1 b bj n B b¼1 bj n n
(A.6)
4 1 X 4 1 X B B 1 D f f p f f o ¼ O . nB b¼1 b bj bk n B b¼1 bj bk n n
(A.7)
and
Theorem 2 follows by using (A.6) and (A.7) together with (A.3)–(A.5). & P P P P P P 4 2 2 Proof of Theorem 3. For simplicity, let C 1 ¼ ni¼1 L^ i , C 2 ¼ ni¼1 j4i L^ i L^ j , C 3 ¼ ni¼1 j4i ni0 Xi;i0 ¼1 Pn P Pn P P ^ ^ ^ ^ L^ i L^ j L^ i0 L^ j0 .Davison et al. (1986, Eq. (2.8)) gives the 0 0 Li Lj Li 0 Lj 0 and C 4 ¼ 0 0 0 0 0 j 4i
i¼1
j4i
i 4i;i aj
j 4i ;j aj
leading term of the variance of the usual bootstrap for the bootstrap variance estimator defined as P ¯ 2 b ðT b T Þ =ðB 1Þ. Though the results are qualitatively the same for this formulation, we choose to P instead use vb ¼ b ðT b TÞ2 =B as it simplifies some of the calculations. Using a parallel approach to
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1871
Davison et al., we get 8 !2 9 X n n < = X 3 2 : 1 L^ 4i þ 2 L^ i V ðvb ðTÞÞ¼ 4 ; n Bn : i¼1 i¼1 C 1 þ 2C 2 n^g þ ð2n 3Þðn 1Þ ¼ , nðn 1Þ Bn4 P P where g^ ¼ ðn 1ÞC 1 =ðC 1 þ 2C 2 Þ ¼ ½ i L^ 4i =ðn 1Þ=½ i L^ 2i =ðn 1Þ2 . The fact that the random pairing method satisfies conditions (3.1) and (3.2) and ( B 1 ajk ¼ 1; 1X ðf bj 1Þðf bk 1Þ ¼ 0 otherwise; B b¼1
ðA:8Þ
where ajk ¼ 1 if units j and k are paired together and 0 otherwise (note that Pðajk ¼ 1Þ ¼ 1=ðn 1Þ), together P P P 2 : with (2.4) implies that vbb ðTÞ¼n2 ð nj¼1 L^ j nj¼1 kaj L^ j L^ k ajk Þ. Thus, the leading term of the Monte Carlo variance of vbb ðTÞ under the random-pairing method reduces to ! n X X : 2 ^ ^ Lj Lk ajk V ðvbb Þ¼V 1=n 2
¼
j¼1 kaj n X X
44 n4 j¼1
2 2 L^ j L^ k V ðajk Þ þ 2
n X X n X X
3 L^ j L^ k L^ j 0 L^ k0 Cov ðajk ; aj 0 k0 Þ5.
ðA:9Þ
j¼1 k4j j 0 ¼1;j 0 Xj k0 4j 0
k4j
Since V ðajk Þ ¼ Pðajk ¼ 1Þ½1 Pðajk ¼ 1Þ ¼ ðn 1Þ1 ½1 ðn 1Þ1 ¼ ðn 2Þ=ðn 1Þ2 and Cov ðajk ; aj0 k0 Þ ¼ Pðajk ¼ 1; aj0 k0 ¼ 1Þ Pðajk ¼ 1ÞPðaj0 k0 ¼ 1Þ 8 4ðn=2Þ½ðn 2Þ=2ðn 4Þ! 1 2 > > ¼ > > 2 < n! ðn 1Þ ðn 1Þ2 ðn 3Þ ¼ > 1 > > > : ðn 1Þ2
if jakaj 0 ak0 ; otherwise;
(A.9) reduces to 16 8 : 1 4ðn 2Þ V ðvbb Þ¼ 4 C4 C2 þ ðC 3 C 4 Þ . n ðn 1Þ2 ðn 1Þ2 ðn 3Þ ðn 1Þ2 P P Next, consider relations among the C i ’s, specifically, ð ni¼1 j4i L^ i L^ j Þ2 ¼ C 2 þ 2C 3 and n X X i¼1 j4i
!2 L^ i L^ j
n X n n X 1 X 2 ¼ L^ i L^ j L^ i 4 i¼1 j¼1 i¼1
!2
n 1 X 2 ¼ L^ i 4 i¼1
where the second equality follows from the fact that On the other hand, 8C 4 ¼
n X X n X
X
i¼1 jai i0 ajai j 0 ai0 ajai
L^ i L^ j L^ i0 L^ j0 ¼ 3
n X i¼1
P
2 L^ i
!2
1 ¼ ðC 1 þ 2C 2 Þ, 4
^ ¼ 0, and thus C 1 ¼ 2C 2 þ 8C 3 .
j Lj
!2 6
n X
4 L^ i ¼ 6C 2 3C 1 .
i¼1
Thus, it is straightforward to show that the leading term of V ðvbb Þ reduces to : C 1 þ 2C 2 2n^g þ 2ðn2 3n þ 3Þ V ðvbb Þ¼ . n4 ðn 1Þ2 ðn 3Þ
(A.10)
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1872
We consider the leading term of V ðvbb Þ=V ðvb Þ, where V ðvb Þ is given in (A.8) replacing B with n=2 þ t ¼ 4m, where t ¼ 1; 2; 3; 4, that is V ðvbb Þ : n^g þ ðn2 3n þ 3Þ nðn þ 2tÞ ¼ . (A.11) V ðvb Þ ðn 1Þðn 3Þ n^g þ ð2n 3Þðn 1Þ The result in (3.4) follows immediately from two facts: that the right-hand side of (A.11) is monotonically decreasing in g^ ; and that g^ Xðn 1Þ=n. The first of these follows from the fact that the derivative with respect to g^ is negative, while the second follows from 1 X ^2 1 X ^4 n 1 1 X ^2 2 ðLi L^ 2 =nÞ2 ¼ Li Li , 0p n1 n1 n n1 P 2 2 where L^ ¼ L^ . &
i
i
References Davison, A.C., Hinkley, D.V., Schechtman, E., 1986. Efficient bootstrap simulation. Biometrika 73, 555–566. Efron, B., 1979. Bootstrap methods: another look at the jackknife. Ann. Statist. 7, 1–26. Efron, B., 1990. More efficient bootstrap computations. J. Amer. Statist. Assoc. 85, 79–89. Efron, B., Tibshirani, R.J., 1993. An Introduction to Bootstrap. Chapman & Hall, London, UK. Graham, R.L., Hinkley, D.V., John, P.W.M., Shi, S., 1990. Balanced design of bootstrap simulations. J. Roy. Statist. Soc. Ser. B 52, 185–202. Hedayat, A.S., Sloan, N.J.A., Stufken, J., 1999. Orthogonal Arrays: Theory and Applications. Springer, New York. Hinkley, D.V., Wei, B.C., 1984. Improvements of jackknife confidence limit method. Biometrika 71, 331–339. Nigam, A.K., Rao, J.N.K., 1996. On balanced bootstrap for stratified multistage samples. Statist. Sinica 6, 199–214. Sitter, R.R., 1993. Balanced repeated replications based on orthogonal multi-arrays. Biometrika 80, 211–221. Wolter, K.M., 1985. Introduction to Variance Estimation. Springer, New York.
Statistics & Probability Letters 76 (2006) 1861–1872 www.elsevier.com/locate/stapro
Fast approximately balanced bootstrap without construction$ C. Devon Lina, Wilson W. Lub, R.R. Sittera, a
Department of Statistics and Actuarial Science, Simon Fraser University, Burnaby, BC, Canada V5A 1S6 b Department of Mathematics and Statistics, Acadia University, Wolfville, NS, Canada B4P 2R6 Received 3 March 2006; accepted 24 April 2006 Available online 5 June 2006
Abstract We develop second-order balanced and approximately balanced resamples that retain the nice properties demonstrated by Graham et al. (1990. Balanced design of bootstrap simulations. J. Roy. Statist. Soc. Ser. B 52, 185–202) but require essentially no design construction. r 2006 Elsevier B.V. All rights reserved. Keywords: Balanced incomplete block design; Balanced repeated replications; Orthogonal array; Stratified sample
1. Introduction Bootstrap simulations (resamples) are often used to approximate aspects of the sampling distribution, such as variance, of estimators when the form of the underlying distribution is unknown (see Efron and Tibshirani, 1993). The number of resamples required so as to obtain an accurate approximation can be large. This can be an important issue for complex estimation situations. There have been various methods proposed for reducing the required number of resamples needed to obtain a fixed level of accuracy. These include Davison et al. (1986), Graham et al. (1990), Efron (1990) and Nigam and Rao (1996). Balanced bootstrapping uses constructions of balanced arrays, balanced incomplete block designs and orthogonal latin squares, to create second-order balanced bootstrap samples that reduce the Monte Carlo (simulation) error and thus require fewer resamples to attain the same level of accuracy. The main impedence on their general usage has been the complexity of construction of such arrays. In this article, we give three strategies for quick and automatic construction of (nearly) second-order balanced arrays that are simple to use and attain significant reductions in Monte Carlo error and thus significantly reduce the number of required resamples to achieve the same level of accuracy. We will use similar notation as and a similar setup to Graham et al. (1990) to aid the reader. In general, a random sample x ¼ ðx1 ; . . . ; xn Þ0 is drawn from a distribution F ðÞ and the estimate T ¼ tðxÞ of y is obtained. The interest lies in approximating the distributional properties of T or some related quantities. Of particular interest are V ðTÞ and PðTpaÞ for fixed a. This can be done by substituting some estimate F^ for F. When no $
Supported by a grant from the Natural Sciences and Engineering Research Council of Canada.
Corresponding author. Tel.: +1 604 291 3334; fax: +1 604 291 4368.
E-mail address: [email protected] (R.R. Sitter). 0167-7152/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2006.04.035
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1862
assumptions on F are made, this is typically the empirical distribution function F^ ðxÞ ¼
Pn
j¼1 I ½xj px =n,
where ^ I ½xj px ¼ 1 if xj px and ¼ 0 otherwise. Since exact calculation of properties of T under F is typically not possible, the bootstrap uses simulation to approximate them. This is done by generating B random resamples : xb ¼ ðxb1 ; . . . ; xbn Þ0 drawn from F^ and calculating T b ¼ tðxb Þ for b ¼ 1; . . . ; B. Since PrðTpajF Þ¼PrðTpajF^ Þ P : : and V ðTjF Þ¼V ðTjF^ Þ, we can approximate via simulation PrðTpajF^ Þ¼ b I ½T b pa =B and :P V ðTjF^ Þ¼ b ðT b TÞ2 =B. The issue in this article is the Monte Carlo error in using these approximations. We will restrict to estimators T ¼ tðF^ Þ of y ¼ tðF Þ that can be expanded as 1 T ¼ tðF^ Þ ¼ tðF Þ þ n
n X
Lðxj ; F Þ þ
j¼1
n X n 1 X Qðxj ; xk ; F Þ þ , 2n2 j¼1 k¼1
where L and Q are, respectively, the first and second functional derivatives of t (see Hinkley and Wei (1984)). Note that L is typically called the influence function. If x is obtained by random sampling from F^ , a similar expansion is available for T ¼ tðx Þ by replacing F^ and F by F^ and F^ , respectively, where F^ is the empirical distribution function of x . That is, T ¼ T þ
n n X n 1X 1 X f j L^ j þ 2 f f Q^ þ , n j¼1 2n j¼1 k¼1 j k jk
P where L^ j ¼ Lðxj ; F^ Þ, Q^ jk ¼ Qðxj ; xk ; F^ Þ, and f j is the frequency of xj in x . Here j L^ j ¼ 0. It has been observed that the differences between various bootstrap simulation methods are embedded in the properties of the frequencies f bj . 2. Balanced bootstrap Efron (1979) notes that f ¼ ðf 1 ; . . . ; f n Þ0 follows the n-category uniform multinomial P distribution. Davison et al. (1986) introduces the notion of a first-order balanced bootstrap by forcing Bb¼1 f bj ¼ B via use of the permutation bootstrap, whereby the B n matrix of f bj follow a hypergeometric distribution with row sums n and column sums B. Graham et al. (1990) observe that if T is symmetric in the x’s, the above expansions show that the mean and variance of T can be approximated without simulation error up to terms of order 1=n provided the f are second-order balanced. They go on to use balanced incomplete block designs to construct B ¼ kn sets f such that B f X bj ¼ 1, B b¼1
(2.1)
B f 2 X bj ¼ ð2n 1Þ=n ¼ 2 1=n B b¼1 B f f X bj bk ¼ ðn 1Þ=n ¼ 1 1=n B b¼1
: ð¼2Þ,
: ð¼1Þ,
(2.2)
(2.3)
where by approximately equal we mean to order n1 . These are in agreement with the multinomial ^ ^ expectations, Eð f j jF^ Þ ¼ 1, Eð f 2 j jF Þ and Eð f j f k jF Þ. To see this, note that the balanced bootstrap variance estimator of T is 2 #2 , B B n X X X : 2 4 ðf 1ÞL^ j =n vbb ðTÞ ¼ ðT b TÞ =B¼ B bj b¼1
¼ n2
b¼1 n X j¼1
B X
j¼1
n B X 1 1X L^ 2j L^ j L^ k ðf bj 1Þ2 þ n2 ðf bj 1Þðf bk 1Þ. B b¼1 B jak b¼1
ð2:4Þ
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1863
P P Under conditions (2.2) and (2.3), b ðf bj 1Þ2 =B ¼ 1 1=n and b ðf bj 1Þðf bk 1Þ=B ¼ 1=n, and thus P 2 : vbb ðTÞ¼ j L^ j =n2 , which is the leading term of Eðvboot ðTÞjF^ Þ under multinomial resampling. This implies that the Monte Carlo variance of vbb ðTÞ is reduced to lower order. P The simplest example is T ¼ x¯ ¼ nj¼1 xj =n. In this case L^ j ¼ xj x¯ and under conditions (2.2) and (2.3), P : vbb ðxÞ ¯ ¼ ðn 1Þs2 =n2 ¼s2 =n, where s2 ¼ j ðxj xÞ ¯ 2 =ðn 1Þ. Graham et al. (1990) also show via simulation that the balanced bootstrap greatly reduces the simulation variability for percentile estimates when 0:05ppp0:95. Unfortunately, balanced bootstrapping is not routinely used in practice. This is almost certainly due to the difficulty in constructing balanced sets of f bj . In the next section, we will introduce three very simple and automatic methods for obtaining such f bj that satisfy the balance conditions or approximately satisfy them so as to easily reduce the simulation error in the bootstrap without any knowledge of design construction techniques. 3. Fast approximately balanced bootstrap Before introducing the three methods for creating fast approximately balanced bootstrap (FABB) resamples, let us consider Hadamard matrices, which will form the basis for all of the methods. A Hadamard matrix H ¼ fdjk g is a B B matrix of 1’s and 1’s, where B is an integer multiple of 4. It can always be rewritten column of 1’s and a row of 1’s, termed standard form. The elements, djk satisfy Pto have aP P d ¼ d ¼ 0, j jk k jk j djk djl ¼ 0 for kal. Table 1 gives an 8 8 example where 1’s are represented by þ and . Hadamard matrices have been well-studied and there exist many algorithms for their construction. Hedayat et al. (1999) give eight algorithms that can be used to obtain B B Hadamard matrices for all B a multiple of 4 up to size 200 200. A simple Fortran program that uses these eight algorithms to yield a B B Hadamard matrix (in standard form) for every B, a multiple of 4 up to 200 can be obtained from the authors. 3.1. The fold-over method It will turn out to be convenient to extend the definitions for the balancing conditions (2.1)–(2.3) to allow for varying resample sizes. To do so, let nb denote the size of the bth resample, and alter the balance conditions to be B f X n bj ¼ 1, (3.1) n B b b¼1 B X f 2 bj ¼ 2, ðn=nb Þ2 B b¼1
(3.2)
Table 1 An 8 8 Hadamard matrix
1 2 3 4 5 6 7 8
y1
y2
y3
y4
y5
y6
y7
y8
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ
þ þ þ þ þ þ þ þ
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1864
B X f bj f bk ¼ 1 1=ðn 1Þ ðn=nb Þ2 B b¼1
: ð¼1Þ.
(3.3)
These are related to the BOMA conditions of Sitter (1993) for obtaining balanced half-samples in stratified random sampling. It is not difficult to see that a similar development to that of (2.4) is available under conditions (3.1)–(3.3) to :P yield vbb ðTÞ¼ j L^ 2j =½nðn 1Þ. In the simplest example where T ¼ x, ¯ vbb ðxÞ ¯ ¼ s2 =n, instead of ðn 1Þs2 =n2 as with the usual balanced bootstrap. To develop a fast method for constructing a set of B resamples that satisfy (3.1)–(3.3), consider the matrix in Table 2. This matrix was obtained by taking the 8-run Hadamard matrix in Table 1, folding it over (adding eight rows by changing signs) and then removing the row of +’s and the row of ’s. Note that the 14 remaining sets represented by keeping the column numbers corresponding to +’s in each row form a resolvable BIBD. Now, if we let each row define a resample by keeping two copies of each yj which has a þ below it, then each row defines a resample of size n. For example, row 2 defines the resample fy1 ; y1 ; y6 ; y6 ; y7 ; y7 ; y8 ; y8 g. Thus, it is not difficult to show that this matrix defines 14 resamples such that the corresponding frequencies satisfy (3.1)–(3.3). That is, if we let dbj denote the bjth element of the matrix, define f bj ¼ dbj þ 1. It turns out that this strategy is quite general. In fact, one can generalize it to approximately satisfy conditions (3.1)–(3.3) for any n via the following theorem. Theorem 1. For n þ t ¼ n0 ¼ 4m, 0ptp3, B ¼ 2n0 2, the two-way array X ðB; nÞ, built by folding-over an n0 n0 Hadamard matrix in standard form, and then removing t columns and the row of þ1’s and the row of 1’s, satisfies (3.1)–(3.3) exactly for t ¼ 0 and for t ¼ 1, 2 or 3, (1) (2) (3)
PB
2 b¼1 ðn=nb Þf bj =B ¼ 1 þ Oð1=n Þ (exact for t ¼ 1 2 2 2 b¼1 ðn=nb Þ f bj =B ¼ 2 þ Oð1=n Þ, PB 2 2 b¼1 ðn=nb Þ f bj f bk =B ¼ 1 1=ðn 1Þ þ Oð1=n Þ.
PB
and 2);
Proof. See appendix. Thus, using the fast simple algorithms for construction of Hadamard matrices combined with Theorem 1 : yields a simple and easy way to construct balanced and approximately balanced bootstraps with B¼2n. Table 2 Construction of a 14 8 array
1 2 3 4 5 6 7 8 9 10 11 12 13 14
y1
y2
y3
y4
y5
y6
y7
y8
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
þ þ þ þ þ þ þ
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1865
3.2. The Hadamard method It turns out that one can gain flexibility in the choice of B by sacrificing a bit in the approximations of Theorem 1. It is clear from conditions (3.2) and (3.3) and approximation (2.4) that the 1=ðn 1Þ term in condition (3.3) is not necessary to capture the leading term of the variance. It is merely kept so that the approximation will be exact for linear statistics such as x. ¯ If we are willing to sacrifice this property, it turns out there is no need to fold-over, stated as follows. Theorem 2. For n þ t ¼ n0 ¼ 4m, 1ptp4, B ¼ n0 1, the two-way array X ðB; nÞ, built by choosing an n0 n0 Hadamard matrix in standard form, and then removing t columns and the row of þ1’s, satisfies the following conditions for t ¼ 1, 2, 3 or 4: (1) (2) (3)
PB
b¼1 ðn=nb Þf bj =B ¼ 1 þ Oð1=nÞ; 2 2 b¼1 ðn=nb Þ f bj =B ¼ 2 þ Oð1=nÞ; PB 2 b¼1 ðn=nb Þ f bj f bk =B ¼ 1 þ Oð1=nÞ.
PB
Proof. See appendix. 3.3. The random-pairing method The fold-over method requires approximately 2n resamples. The Hadamard method requires approximately n resamples. This makes them particularly useful when n is small to moderate. Once n becomes larger, one would likely prefer to use a regular bootstrap with Bon. In this section, we develop a third method that requires approximately n=2 replicates. It is related to some random methods for balanced repeated replication (BRR) variance estimates used in surveys (see Wolter, 1985). For the moment, assume n is even. In our context, the idea would reduce to randomly ordering y1 ; . . . ; yn , denoted y1 ; . . . ; yn and then forming successive pairs ðy1 ; y2 Þ; ðy3 ; y4 Þ; . . . ; ðyn1 ; yn Þ. We now use something akin to the BRR method treating these pairs as strata with two units in each stratum. That is, we form a B ðn=2Þ matrix, H ¼ fdbk g, by removing columns : from a B B Hadamard matrix, where B is the smallest multiple of 4 strictly larger than n=2, so B¼n=2. The rows of this matrix define the f bj as follows: f bð2k1Þ ¼ 1 þ dbk and f bð2kÞ ¼ 1 dbk , for k ¼ 1; . . . ; n=2. If this method is used, nb ¼ nPand it is easily seen that conditions (3.1) and (3.2) are satisfied. However, condition (3.3) is not satisfied. Instead, Bb¼1 ðn=nb Þ2 f bj f bk =B ¼ 0 and 1 for j and k paired together or not, respectively. Thus, this method does not remove the leading term from the Monte Carlo variance of vbb ðTÞ. However, the following theorem implies that the leading term of the Monte Carlo variance is reduced substantially. Theorem 3. Using the random-pairing method, the leading term of V ðvbb Þ=V ðvb Þ, where B ¼ n=2 þ t, satisfies V ðvbb Þ nðn þ 2tÞ½n^g þ ðn2 3n þ 3Þ nðn 2Þ2 ðn þ 2tÞ p , ¼ V ðvb Þ ðn 1Þðn 3Þ½n^g þ ð2n 3Þðn 1Þ 2ðn 1Þ3 ðn 3Þ
(3.4)
where V P ðvbb Þ and V ðvb Þ are the Monte Carlo variances of vbb and vb , respectively, and g^ ¼ ðn 1Þ P ð i L^ 4i Þ=ð i L^ 2i Þ2 is an estimate of the kurtosis. Proof. See appendix. For all t ¼ 1; . . . ; 4, the right-hand side of (3.4) is decreasing in n to a limit of 0.5, summarized in Table 3. Of course, this is an upper bound and depending upon the sample value of g^ the random pairing method can reduce the Monte Carlo variance much more. For example, if g^ 4ðn 1Þ=2 the upper bound in (3.4) becomes nðn 2Þðn þ 2tÞ=½ðn 1Þ2 ð5n 6Þ and for moderate sized n is nearer 0:2 than 0.5. We explore the actual performance through some limited simulations in the next section. Extending to n odd can be accomplished fairly easily by randomly pairing n 3 units, and then using the remaining three units to randomly create two strata that have one unit in common. This adequately approximates the even n case.
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1866
Table 3 The bound of ratio of Monte Carlo variance of random pairing balanced bootstrap relative to the regular bootstrap n
20
30
40
50
t¼1 t¼2 t¼3 t¼4
0.61 0.67 0.72 0.78
0.57 0.61 0.64 0.68
0.55 0.58 0.61 0.63
0.54 0.56 0.58 0.60
4. Simulations 4.1. Percentiles for the sample mean and t-statistic In this section, we investigate the Monte Carlo error when estimating percentiles rather than moments of the sample average in a similar fashion to that used by Graham et al. (1990, Example 9). Consider the skew sample of n ¼ 22 x values 8:0; 9:6; 10:0; 10:4; 11:7; 13:0; 14:0; 15:0; 15:8; 16:6; 16:9; 17:2; 17:25; 17:3; 19:55; 21:8; 22:9; 24:0; 25:45; 26:9; 30:35; 33:8: To evaluate the bootstrap percentile estimates for x¯ and the Student t-statistic t ¼ n1=2 ðx¯ mÞ=s, we obtain x¯ 1 px¯ 2 p px¯ B for each of the three FABBs. Since the three FABBs use differing resample sizes, we will compare them to the usual bootstrap using the same resample size. Let x¯ p be the 100pth percentile of x¯ 1 px¯ 2 p px¯ B . The simulation error is then ep ¼ x¯ p x¯ p ,
(4.1)
where P ðx¯ px¯ p Þ ¼ p and P denotes probability under the resampling distribution. The values of x¯ p and tp are obtained via an independent simulation taking K ¼ 10; 000 resamples. For each of the three FABBs, we repeat the resampling 1000 times, independently, and similarly for the usual bootstrap. The resulting 1000 errors ep in (4.1) are summarized by their mean and variance for p ¼ i=ðB þ 1Þ for i ¼ 1; . . . ; B and plotted against the corresponding exact percentiles in Fig. 1 (for the mean). The corresponding figure for the t-statistic (not shown) is qualitatively the same. Clearly the fold-over balanced bootstrap yields significant reductions in simulation variability for percentile estimates and maintains near unbiasedness. These results parallel those in Graham et al. (1990), as do similar simulations for the correlation coefficient (not shown). The Hadamard method performs nearly as well. The random pairing method also gives significant gains. 4.2. Variance estimates for the sample mean and sample correlation To study the simulation variance of the FABB variance estimators relative to the usual bootstrap we conducted a limited simulation study employing the following model: 1=2
xi ¼ bzi þ zi i ,
(4.2)
where i Nð0; s2 Þ and independent of zi gammaðg; hÞ. Thus, the mean and variance of zi are given by mz ¼ gh and s2z ¼ gh2 , respectively. Further, the mean and variance of xi are mx ¼ bmz and s2x ¼ b2 s2z þ mz s2 , and corrðxi ; zi Þ ¼ r ¼ bsz =sx . We fixed b ¼ 1:0, mz ¼ 100 and s ¼ 30. For various settings of s2z to vary r, and various values of n, we generated N 1 ¼ 1000 independent random samples from model (4.2) with sample size n. For each such sample we repeatedly (N 2 ¼ 1000 times) obtained the fold-over balanced bootstrap variance estimator, vfobb (resample : size B ¼ 2n0 2¼2n), the Hadamard balanced bootstrap variance estimator, vhbb (resample size
ARTICLE IN PRESS standard deviation of error
C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
average error
0.6 0.4 0.2 0.0 −0.2 −0.6 0.0
0.2
(a)
0.4
0.6
0.8
0.2 0.0 −0.2 −0.6
(c)
0.4
0.6
0.8
0.2 0.0 −0.2 −0.6
(e)
0.6
exact percentle
0.0 0.2
0.8
1.0
0.6
0.8
1.0
0.8
1.0
0.8
1.0
1.0 0.8 0.6 0.4 0.2 0.0 0.2
0.4
0.6
exact percentle
1.0 0.8 0.6 0.4 0.2 0.0 0.0
(f)
0.4
exact percentle
0.0
standard deviation of error
average error
0.4
0.4
0.2
(d)
0.6
0.2
0.4
1.0
exact percentle
0.0
0.6
0.0
standard deviation of error
average error
0.4
0.2
0.8
(b)
0.6
0.0
1.0
1.0
exact percentle
1867
0.2
0.4
0.6
exact percentle
Fig. 1. Averages and (b) standard deviations of 1000 simulation errors in bootstrap percentile estimates for x: ¯ , regular bootstrap; , fold-over balanced bootstrap. (c) and (d), and (e) and (f), are similar plots for the Hadamard and random-pairing balanced bootstraps, respectively.
: : B ¼ n0 1¼n), the random-pairing balanced bootstrap variance estimator, vrpbb (resample size B ¼ n0 =2¼n=2) and the usual bootstrap variance estimator, vb (resample size B chosen to be the same as that of the balanced ^ bootstrap being compared to), for T 1 ¼ x, ¯ T 2 ¼ z¯ and T 3 ¼ r. If we let E 1 and V 1 denote the expectation and variance under the sampling from the model, and let E 2 and V 2 denote the expectation and variance under the resampling conditional on each sample, respectively, V ðvÞ ¼ V 1 E 2 ðvÞ þ E 1 V 2 ðvÞ for each variance estimator v, where the second term is the average Monte Carlo variance. To estimate the average Monte Carlo variance for the three approximately balanced bootstrap variance estimators and for the P usual bootstrap, for each of the N 1 samples we estimate the conditional Monte Carlo ~ ðvÞ ¼ N 2 ðvl v¯ Þ2 =N 2 , where v denotes one of the four bootstrap variance estimators and variance using V l¼1 P 2 v¯ ¼ N l¼1 vl =N 2 . This yields N 1 estimates of the conditional simulation variance, V 2 , of each of the three bootstrap variance estimators for T 1 , T 2 and T 3 . We then take the average of these N 1 values as an estimate of the average Monte Carlo variance, denoted EVðvÞ. Table 4 gives REVfobb ¼ EV½vfobb ðT i Þ=EV½vb ðT i Þ, REVhbb ¼ EV½vhbb ðT i Þ=EV½vb ðT i Þ and REVrpbb ¼ EV½vrpbb ðT i Þ=EV½vb ðT i Þ for i ¼ 1, 2 and 3, n ¼ 22, 24, 50, and 52, and r ¼ 0:2, 0.5, and 0.8, where the B in the usual bootstrap was chosen to match the balanced bootstrap it is being compared to in each ratio. Viewing the results for T 1 ¼ x¯ and T 2 ¼ z¯ in Table 4, we see that REVðvfobb Þ ¼ 0 for n ¼ 24 and n ¼ 52 verifying Case 1 of Theorem 1, and are extremely small for n ¼ 22 and 50 verifying Case 3 of Theorem 1 (Note, the n ¼ 50 case yielded non-zero results if viewed to six decimal places). The Hadamard method performed nearly as well. Viewing the results for the random-pairing method, we see that the Monte Carlo variance is almost cut in half in the worst cases, and is reduced far more in many cases. Note that Eð^gÞ was
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1868
Table 4 Ratio of average Monte Carlo variances of balanced bootstraps relative to the regular bootstrap for x, ¯ z¯ and r^ n
r ¼ 0:2
r ¼ 0:5
r ¼ 0:8
T1
T2
T3
T1
T2
T3
T1
T2
T3
REVfobb
22 24 50 52
0.0003 0.0000 0.0000 0.0000
0.0030 0.0000 0.0000 0.0000
0.1129 0.0909 0.0462 0.0420
0.0002 0.0000 0.0000 0.0000
0.0002 0.0000 0.0000 0.0000
0.1120 0.1251 0.1174 0.1179
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.1202 0.1268 0.1969 0.1949
REVhbb
22 24 50 52
0.0113 0.0262 0.0018 0.0048
0.0108 0.0256 0.0018 0.0047
0.3925 0.3767 0.2110 0.2025
0.0071 0.0179 0.0014 0.0038
0.0056 0.0152 0.0013 0.0035
0.3114 0.3592 0.3245 0.3308
0.0016 0.0065 0.0006 0.0013
0.0012 0.0065 0.0006 0.0012
0.2857 0.3039 0.4150 0.4140
REVrpbb
22 24 50 52
0.4927 0.5979 0.5332 0.5115
0.4700 0.5861 0.5188 0.5011
0.7981 0.8615 0.7416 0.7034
0.3083 0.4022 0.4118 0.4001
0.2452 0.3424 0.3726 0.3659
0.5776 0.6984 0.7111 0.7083
0.0699 0.1417 0.1661 0.1330
0.0538 0.1389 0.1700 0.1194
0.5001 0.5531 0.7300 0.7131
chosen quite large for r ¼ 0.5 and 0.8 for z¯ and r ¼ 0.8 for x¯ (overall ranging from 3 to 20). Roughly speaking, one would need to use a resample size about 2–20 times larger with vb to attain the same Monte Carlo variance as that of vrpbb in these situations. ^ we see that the gains are not restricted to linear estimators but are almost as Viewing results for T 3 ¼ r, dramatic for this nonlinear statistic. For r^ it does appear that obtaining the closer approximation by folding over has a larger relative impact. In this case, one would need to use a resample size about 8–25 times larger with vb to attain the same simulation variance as that of vfobb , about 3–10 times larger for vhbb and 1.2 or 2 times larger for vrpbb . We repeated the process using n ¼ 21, 23, 49 and 51 by removing one observation and using the method described at the end of Section 3 for the odd n case for the random-pairing method. The results were qualitatively the same (not shown). 5. Concluding remarks We develop three methods for obtaining approximately balanced bootstrap resamples that retain the advantages in performance of second-order balancing without a need for tedious and difficult design construction techniques. Acknowledgement We thank D. Bingham for helpful comments. Appendix Proof of Theorem 1. For H n0 ¼ fdjk gn0 n0 a Hadamard matrix in standard form: (H1) Each column (row) of H n0 , not including the column (row) of 1’s, consists of n0 =2 ¼ 2m, þ1’s and 2m; 1’s; (H2) For two distinct columns (rows) of H n0 , u ¼ ðu1 ; . . . ; un0 Þ0 and v ¼ ðv1 ; . . . ; vn0 Þ0 , each of the pairs ðui ; vi Þ ¼ ðþ1; þ1Þ, ðþ1; 1Þ, ð1; þ1Þ and ð1; 1Þ appears m times.
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1869
(H3) Since H n0 is a strength 2 orthogonal array, its fold-over, H~ n0 will be strength 3, that is, the rows of any three columns of H~ n0 ¼ ðH Tn0 ; H Tn0 ÞT will consist of an equal number of copies of the eight possible combinations of ð1; 1; 1Þ. From (H1)–(H3), B X
f bj
0
¼ 2n 2;
b¼1
n0 1 X
f bj f bk ¼ n0 4 and
b¼1
In Pnaddition, for each replicate b, nb ¼ x are j¼1 bj
8 0 > > > < 1 Db ¼ > 0; 2 > > : 1; 3
Pn
j¼1 f bj
B X
f bj f bk ¼ n0 .
(A.1)
b¼n0
¼nþ
Pn
j¼1 xbj
¼ n þ Db , where the possible values of Db ¼
for t ¼ 0; for t ¼ 1; for t ¼ 2;
(A.2)
for t ¼ 3:
Therefore, we have the approximation, n=nb ¼ nðn þ Db Þ1 ¼ 1 Db =n þ Oðn2 Þ, and similarly, ðn=nb Þ2 ¼ 1 2Db =n þ Oðn2 Þ. Due to this, together with (A.1), the left-hand sides of conditions (3.1), (3.2) and (3.3) can be expressed as B f X n bj b¼1
nb
B B 1X 1 X 1 ¼ f bj Db f bj þ O 2 B b¼1 nB b¼1 n B B 1 X 1 ¼1 Db f bj þ O 2 , nB b¼1 n
B 2 f 2 X n bj b¼1
nb
B B 2X 4 X 1 ¼ f bj Db f bj þ O 2 B b¼1 nB b¼1 n B B 4 X 1 ¼2 Db f bj þ O 2 , nB b¼1 n
ðA:3Þ
ðA:4Þ
and B 2 f f X n bj bk b¼1
nb
B
B B 1X 2 X 1 ¼ f f Db f bj f bk þ O 2 B b¼1 bj bk nB b¼1 n B 1 2 X 1 ¼1 Db f bj f bk þ O 2 , n 1 nB b¼1 n
ðA:5Þ
respectively. P P It only remains to prove that both Bb¼1 Db f bj and Bb¼1 Db f bj f bk are terms of Oð1Þ for each of t ¼ 0, 1, 2 and 3. P P For t ¼ 0, we have Bb¼1 Db f bj ¼ Bb¼1 Db f bj f bk ¼ 0, since Db ¼ 0 by (A.2), and conditions (3.1)–(3.3) are thus exactly satisfied. To consider t ¼ 1, 2 and 3, we first note the relationships between any replicate b and its fold-over b0 ; Db ¼ Db0 and xbj ¼ xb0 j , so that f bj f bk ¼ ð1 þ xbj Þð1 þ xbk Þ ¼ ð1 xb0 j Þð1 xb0 k Þ ¼ ð1 þ xb0 j Þð1 þ xb0 k Þ 2ðxb0 j þ xb0 k Þ ¼ f b0 j f b0 k 2ðxb0 j þ xb0 k Þ.
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1870
Thus, B X
Db f bj ¼
b¼1
n0 1 X
Db ð1 þ xbj Þ þ
n0 1 X
ðDb þ Db xbj Þ þ
b¼1
¼2
Db0 ð1 þ xb0 j Þ
b0 ¼n0
b¼1
¼
B X
n0 1 X
½ðDb Þ þ ðDb Þðxbj Þ
b¼1
n0 1 X
Db xbj ,
b¼1
and B X
Db f bj f bk ¼
b¼1
n0 1 X
Db f bj f bk þ
¼
Db f bj f bk þ
b¼1
¼2
Db0 f b0 j f b0 k
b0 ¼n0
b¼1 n0 1 X
B X
n0 1 X
n0 1 X
ðDb Þ½f bj f bk 2ðxbj þ xbk Þ
b¼1
Db xbj þ 2
b¼1
n0 1 X
Db xbk .
b¼1
P 0 1 Therefore, we need only prove that nb¼1 Db xbj ¼ Oð1Þ. Without loss of generality, we assume that in all cases, t ¼ 1, 2 and 3, the column consisting of þ1’s for its first n0 1Pelements and 1’s Pn0 1for the rest, is deleted. 0 1 Case 1: t ¼ 1. Since Db 1 for b ¼ 1; . . . ; n0 1, we have nb¼1 Db xbj ¼ b¼1 xbj ¼ 1. Case 2: t ¼ 2. D ¼ 2 if ‘‘þþ’’ is deleted from replicate b; or ¼ 0 if ‘‘þ’’ is deleted. Thus, P b Pn0 1 D x ¼ 2 x ¼ 2ð1Þ ¼ 2, where the second to last equality follows from (H2). b¼1 b bj b2fb:Db ¼2g bj Case 3: t ¼ 3. Db ¼ 3; 1; 1; 1, corresponding to the four patterns of deletions, ‘‘þ þ þ’’, ‘‘ þ þ’’, ‘‘þ þ’’ and ‘‘ þ’’. DenoteP A1 , A2 , A3 and A4 as collections of replicates associated with those four patterns, respectively. Let S ij ¼ b2Ai xbj for i ¼ 1; 2; 3; 4 and Pn0j1¼ 1; . . . ; n. From (H2) it is easy to see that S1j þ S 2j ¼ 1; S 1j þ S3j ¼ 1; S 2j þ S 4j ¼ 0. Therefore, b¼1 Db xbj ¼ 2ðS 1j þ S 2j Þ ðS 1j þ S 3j Þ þ ðS 2j þ S4j Þ ¼ 3. & PB PB Proof Theorem 2. For the Hadamard method, we have b¼1 xbj =B ¼ 1, b¼1 f bj ¼ 1 1=B, and PB of f f =B ¼ 1 3=B. For b ¼ 1; . . . ; B, and t ¼ 1; 2; 3; 4, note that jD jp4. Thus, b b¼1 bj bk 4 1 X 4 1 X B B 1 D f p f o ¼ O , nB b¼1 b bj n B b¼1 bj n n
(A.6)
4 1 X 4 1 X B B 1 D f f p f f o ¼ O . nB b¼1 b bj bk n B b¼1 bj bk n n
(A.7)
and
Theorem 2 follows by using (A.6) and (A.7) together with (A.3)–(A.5). & P P P P P P 4 2 2 Proof of Theorem 3. For simplicity, let C 1 ¼ ni¼1 L^ i , C 2 ¼ ni¼1 j4i L^ i L^ j , C 3 ¼ ni¼1 j4i ni0 Xi;i0 ¼1 Pn P Pn P P ^ ^ ^ ^ L^ i L^ j L^ i0 L^ j0 .Davison et al. (1986, Eq. (2.8)) gives the 0 0 Li Lj Li 0 Lj 0 and C 4 ¼ 0 0 0 0 0 j 4i
i¼1
j4i
i 4i;i aj
j 4i ;j aj
leading term of the variance of the usual bootstrap for the bootstrap variance estimator defined as P ¯ 2 b ðT b T Þ =ðB 1Þ. Though the results are qualitatively the same for this formulation, we choose to P instead use vb ¼ b ðT b TÞ2 =B as it simplifies some of the calculations. Using a parallel approach to
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1871
Davison et al., we get 8 !2 9 X n n < = X 3 2 : 1 L^ 4i þ 2 L^ i V ðvb ðTÞÞ¼ 4 ; n Bn : i¼1 i¼1 C 1 þ 2C 2 n^g þ ð2n 3Þðn 1Þ ¼ , nðn 1Þ Bn4 P P where g^ ¼ ðn 1ÞC 1 =ðC 1 þ 2C 2 Þ ¼ ½ i L^ 4i =ðn 1Þ=½ i L^ 2i =ðn 1Þ2 . The fact that the random pairing method satisfies conditions (3.1) and (3.2) and ( B 1 ajk ¼ 1; 1X ðf bj 1Þðf bk 1Þ ¼ 0 otherwise; B b¼1
ðA:8Þ
where ajk ¼ 1 if units j and k are paired together and 0 otherwise (note that Pðajk ¼ 1Þ ¼ 1=ðn 1Þ), together P P P 2 : with (2.4) implies that vbb ðTÞ¼n2 ð nj¼1 L^ j nj¼1 kaj L^ j L^ k ajk Þ. Thus, the leading term of the Monte Carlo variance of vbb ðTÞ under the random-pairing method reduces to ! n X X : 2 ^ ^ Lj Lk ajk V ðvbb Þ¼V 1=n 2
¼
j¼1 kaj n X X
44 n4 j¼1
2 2 L^ j L^ k V ðajk Þ þ 2
n X X n X X
3 L^ j L^ k L^ j 0 L^ k0 Cov ðajk ; aj 0 k0 Þ5.
ðA:9Þ
j¼1 k4j j 0 ¼1;j 0 Xj k0 4j 0
k4j
Since V ðajk Þ ¼ Pðajk ¼ 1Þ½1 Pðajk ¼ 1Þ ¼ ðn 1Þ1 ½1 ðn 1Þ1 ¼ ðn 2Þ=ðn 1Þ2 and Cov ðajk ; aj0 k0 Þ ¼ Pðajk ¼ 1; aj0 k0 ¼ 1Þ Pðajk ¼ 1ÞPðaj0 k0 ¼ 1Þ 8 4ðn=2Þ½ðn 2Þ=2ðn 4Þ! 1 2 > > ¼ > > 2 < n! ðn 1Þ ðn 1Þ2 ðn 3Þ ¼ > 1 > > > : ðn 1Þ2
if jakaj 0 ak0 ; otherwise;
(A.9) reduces to 16 8 : 1 4ðn 2Þ V ðvbb Þ¼ 4 C4 C2 þ ðC 3 C 4 Þ . n ðn 1Þ2 ðn 1Þ2 ðn 3Þ ðn 1Þ2 P P Next, consider relations among the C i ’s, specifically, ð ni¼1 j4i L^ i L^ j Þ2 ¼ C 2 þ 2C 3 and n X X i¼1 j4i
!2 L^ i L^ j
n X n n X 1 X 2 ¼ L^ i L^ j L^ i 4 i¼1 j¼1 i¼1
!2
n 1 X 2 ¼ L^ i 4 i¼1
where the second equality follows from the fact that On the other hand, 8C 4 ¼
n X X n X
X
i¼1 jai i0 ajai j 0 ai0 ajai
L^ i L^ j L^ i0 L^ j0 ¼ 3
n X i¼1
P
2 L^ i
!2
1 ¼ ðC 1 þ 2C 2 Þ, 4
^ ¼ 0, and thus C 1 ¼ 2C 2 þ 8C 3 .
j Lj
!2 6
n X
4 L^ i ¼ 6C 2 3C 1 .
i¼1
Thus, it is straightforward to show that the leading term of V ðvbb Þ reduces to : C 1 þ 2C 2 2n^g þ 2ðn2 3n þ 3Þ V ðvbb Þ¼ . n4 ðn 1Þ2 ðn 3Þ
(A.10)
ARTICLE IN PRESS C.D. Lin et al. / Statistics & Probability Letters 76 (2006) 1861–1872
1872
We consider the leading term of V ðvbb Þ=V ðvb Þ, where V ðvb Þ is given in (A.8) replacing B with n=2 þ t ¼ 4m, where t ¼ 1; 2; 3; 4, that is V ðvbb Þ : n^g þ ðn2 3n þ 3Þ nðn þ 2tÞ ¼ . (A.11) V ðvb Þ ðn 1Þðn 3Þ n^g þ ð2n 3Þðn 1Þ The result in (3.4) follows immediately from two facts: that the right-hand side of (A.11) is monotonically decreasing in g^ ; and that g^ Xðn 1Þ=n. The first of these follows from the fact that the derivative with respect to g^ is negative, while the second follows from 1 X ^2 1 X ^4 n 1 1 X ^2 2 ðLi L^ 2 =nÞ2 ¼ Li Li , 0p n1 n1 n n1 P 2 2 where L^ ¼ L^ . &
i
i
References Davison, A.C., Hinkley, D.V., Schechtman, E., 1986. Efficient bootstrap simulation. Biometrika 73, 555–566. Efron, B., 1979. Bootstrap methods: another look at the jackknife. Ann. Statist. 7, 1–26. Efron, B., 1990. More efficient bootstrap computations. J. Amer. Statist. Assoc. 85, 79–89. Efron, B., Tibshirani, R.J., 1993. An Introduction to Bootstrap. Chapman & Hall, London, UK. Graham, R.L., Hinkley, D.V., John, P.W.M., Shi, S., 1990. Balanced design of bootstrap simulations. J. Roy. Statist. Soc. Ser. B 52, 185–202. Hedayat, A.S., Sloan, N.J.A., Stufken, J., 1999. Orthogonal Arrays: Theory and Applications. Springer, New York. Hinkley, D.V., Wei, B.C., 1984. Improvements of jackknife confidence limit method. Biometrika 71, 331–339. Nigam, A.K., Rao, J.N.K., 1996. On balanced bootstrap for stratified multistage samples. Statist. Sinica 6, 199–214. Sitter, R.R., 1993. Balanced repeated replications based on orthogonal multi-arrays. Biometrika 80, 211–221. Wolter, K.M., 1985. Introduction to Variance Estimation. Springer, New York.