A Monte Carlo comparison of three consistent bootstrap procedures

Vol. 00, No. 00, Month 200x, 1–11 A Monte Carlo comparison of three consistent bootstrap procedures R. Pino Mej´ıas∗†, M.D. Jim´enez Gamero†, A. Enguix Gonz´alez† †Dpto. de Estad´ıstica e Investigaci´on Operativa, Universidad de Sevilla, Facultad de Matem´aticas, C/ Tarfia s.n., 41.012 Sevilla, Spain (Received 00 Month 200x; In final form 00 Month 200x) Since bootstrap samples are simple random samples with replacement from the original sample, the information content of some bootstrap samples can be very low. To avoid this fact, some authors have proposed several variants of the classical bootstrap. In this paper we consider two of them: the sequential or Poisson bootstrap and the reduced bootstrap. Both of them, like ordinary bootstrap, can yield second order accurate distribution estimators, that is, the three bootstrap procedures are asymptotically equivalent. The question that naturally arises is which of them should be used in a practical situation, in other words, which of them should be used for finite sample sizes. To try to answer this question, we have carried out a simulation study. Although no method was found to exhibit best performance in all the considered situations, some recommendations are given. Keywords: Bootstrap, Poisson bootstrap, reduced bootstrap, distribution estimation, finite sample performance. 1 Introduction The bootstrap, introduced by Efron [1], is a powerful tool for nonparametric estimation of sampling distributions and standard errors. It may be described as follows. Let X= (X1, X2, . . . , Xn) be a random sample from an unknown distribution F, and let Tn=Tn(X;F) be a statistic of interest. Let Fnbe the empirical distribution function of X1, X2, . . . , Xnand let X∗= (X∗ 1, X∗ 2, . . . , X∗ n) be a random sample drawn from Fn.X∗is called a bootstrap sample. The bootstrap method estimates the distribution of Tnthrough the conditional distribution of T∗ n=Tn(X∗;Fn), given X1, X2, . . . , Xn. Bootstrap samples are simple random samples of size nselected with replacement from the original sample. Thus, all bootstrap samples are not equally informative, when the information content of X∗is measured through the number of distinct original observations in it, νn. As Rao et al. [2] assert, this ∗Corresponding author. E-mail: [email protected] 2R. Pino Mej´ıas et al. variability is neither necessary nor desirable. To avoid it, these authors have introduced a sequential sampling method that keeps constant the information content of bootstrap samples, the sequential or Poisson bootstrap (originally, they named it sequential bootstrap; the name Poisson bootstrap was given latter by Babu et al. [3]). It consists of sampling sequentially (with replacement from the original sample) until mdistinct original observations appear. To keep the information constant, this method introduces another source of variability: the size of the sequential bootstrap samples. To reduce the variability of νn, keeping the size of bootstrap samples fixed, Jim´enez-Gamero et al. [4] have introduced a modified bootstrap that consists of only considering those bootstrap samples satisfying k1≤νn≤k2, for some 1≤k1≤k2≤n. They call it the reduced bootstrap, since it only uses a portion of the set of all possible bootstrap samples. The method is an extension of the one introduced by Mu˜noz-Garc´ıa et al. [5], that takes k2=n. Note that ordinary bootstrap is a particular case of reduced bootstrap with k1= 1 and k2=n. Another method to reduce the variability of νn, keeping the size of bootstrap samples fixed, can be found in Stromberg [6]. In order to reduce the breakdown point of the bootstrap variance estimator, which is 1/n regardless of the breakdown point of the estimator due to the fact that P(νn=n)>0, Stromberg [6] proposed to use a ”limited replacement bootstrap”, in which all bootstrap samples selecting the same element more than m∗times are excluded, that is, it only considers those bootstrap samples with Ni≤m∗, 1≤i≤n, where Ni=card{j:X∗ j=Xi}, 1 ≤i≤n. As before, note that ordinary bootstrap is a particular case of limited replacement bootstrap with m∗=n. Singh [7] has shown that the bootstrap estimator of the distribution of the sample mean is consistent and that it can be more accurate than the approximation given by the central limit theorem when higher-order population moments exist. Babu and Singh [8] have extended these results to functions of multivariate means and their studentized versions. Babu et al. [3] have shown that they are also true for the sequential bootstrap when m=E(νn). Jim´enez-Gamero et al. [9] have also proven them for the reduced bootstrap, for certain choices of k1and k2. Therefore, the results in Babu et al. [3] and Jim´enez-Gamero et al. [9] reveal that the usual bootstrap high variability of νnis not necessary. To our knowledge there is no similar result for the limited replacement bootstrap of Stromberg [6] so, from now on we will not consider this method. In summary, the above cited results tell us that we have three methods to consistently estimate the distribution function of statistics in an important class (the functions of multivariate means and their studentized versions), that is, all of them behave similarly for large samples. The question that naturally Comparison of bootstrap procedures 3 arises is which of them should be used in a practical situation, that is, for finite sample sizes. In order to try to answer this question, we have carried out a simulation study. Although it is not exhaustive (it is impossible to simulate all cases), we can conclude some general recommendations from the obtained results. The paper is organized as follows. In Section 2 we briefly review the methods to be compared. In Section 3 we describe the simulations and display the obtained results. Section 4 concludes giving some practical recommendations in the light of the obtained results. 2 The methods In this Section we first study the variability reduction properties of the considered methods, and then we review some properties of the corresponding distribution estimators. Here we only review the second order correctness properties because they just make the considered distribution estimators to have an edge over the approximation by the limiting normal distribution. 2.1 Variability reduction properties As in the Introduction, let νndenote the number of distinct observations in a bootstrap sample. Then under ordinary bootstrap, E(νn) = n·1−µ1−1 n¶n¸≃0.632n, var(νn) = nµ1−1 n¶n +n(n−1) µ1−2 n¶n −n2µ1−1 n¶2n ≃0.233n, To reduce the variability of νn, the sequential bootstrap draws observations from the original sample with replacement until there are m≃E(νn) distinct observations in the bootstrap sample. This way, the number of distinct observations is no longer a random variable. Nevertheless, the method yields bootstrap samples with random size N. The expectation and the variance of Nare (see Rao et al. [2]) E(N) = nµ1 n+1 n−1+... +1 n−m+ 1¶=n+O(1), 4R. Pino Mej´ıas et al. var(N) = m X j=1 n(j−1) (n−j+ 1)2=n(e−1) + O(1). Next result gives the expectation and the variance of νnfor the reduced bootstrap (see Appendix C for a proof). Proposition 2.1 Let X∗be a bootstrap sample, νn=νn(X∗)and let k1, k2 be two positive integers with 1≤k1≤k2≤n. Then, E(νn|k1≤νn≤k2) = k2−(k2−k1)µ1−1 k2¶n−k1 and var (νn|k1≤νn≤k2) = (k2−k1)(k2−k1−1) µ1−2 k2¶n−k1 + +(k2−k1)µ1−1 k2¶n−k1 −(k2−k1)2µ1−1 k2¶2(n−k1) . (1) From (1) we have that if k1=k1(n) and k2=k2(n) are such that (n−k1)/k2 converges to a real limit l6= 0, as n→ ∞, then var (νn|k1≤νn≤k2) = O(k2−k1). Hence, as var(νn)≃0.233n, if k2−k1=o(n), the variability of νnin the reduced bootstrap will be less than in the usual bootstrap. If (n−k1)/k2→ ∞, as n→ ∞, this implies that k2=o(n). As it is shown in Jim´enez-Gamero et al. [9], in this case the reduced bootstrap is not consistent and thus this choice of k1and k2is not advisable. A similar situation occurs when (n−k1)/k2→0, as n→ ∞. Hence, for adequate choices of k1and k2, the reduced bootstrap variability of νnis less than that of ordinary bootstrap. 2.2 Second order accuracy of the distribution estimators Let X1, X2, ..., Xnbe i.i.d. random variables having common distribution function F, mean µand variance σ2. Let ¯ Xn=1 n n X j=1 Xj, s2 n=1 n n X j=1 (Xj−¯ Xn)2and ¯ X∗ n=1 n n X j=1 X∗ j. Throughout this paper P∗will denote the bootstrap conditional probability law, given X1, X2, . . . , Xn. Comparison of bootstrap procedures 5 If E|X1|3<∞and Fsatisfies Cram´er’s condition, Theorem 1.D in Singh [7] shows that sup xn1/2¯¯¯P∗nn1/2(¯ X∗ n−¯ Xn)/sn≤xo− −Pnn1/2(¯ Xn−µ)/σ ≤xo¯¯¯=o(1). (2) This implies that the ordinary bootstrap approximation to P{n1/2(¯ Xn− µ)/σ ≤x}is better than the normal approximation, which has an error of order O(n−1/2). Theorem 3.1 in Babu et al. [3] gives the analogue of the above result for the Poisson bootstrap, with m=E(νn). Corollary 1 in Jim´enez- Gamero et al. [9] states conditions on k1and k2(1 ≤k1≤k2≤n) for (2) to hold for the reduced bootstrap distribution estimator. These conditions are the following: Condition C.1.a Φ(w2)−Φ(w1)≥α,∀n≥n0, for some n0∈Nand some fixed constant α > 0 and w1+w2=o(1) or w1→ −∞ and w2→+∞, Condition C.2.aw1,w2→0, where w1= (k1−1−np)n−1/2σ0−1,w2= (k2−np)n−1/2σ0−1,σ02=pq−q2,p= 1−e−1,q= 1 −pand Φ denotes the standard normal cumulative distribution function. Since P(νn≤k) = Φ(w) + o(1) with w= (k−np)/√nσ0(Johnson and Kotz [10, p. 318]), if k1and k2are such that some of the above conditions hold, then the proportion of bootstrap samples not to be considered because they have a small number of different elements (νn< k1) is (almost) the same as the proportion of bootstrap samples not to be considered because they have a big number of different elements (νn> k2). The above results can be extended to the multivariate case and to statistics which can be expressed as smooth functions of multivariate means. Let X1, X2, ..., Xnbe ni.i.d. d-dimensional random vectors with common distribution function F, mean µand variance matrix Σ, where dis a fixed positive integer. Let ¯ Xnbe the sample mean and let Σnbe the sample variance matrix. Let Hbe a real valued Borel measurable function on Rd. If His differentiable at y, we denote by h(y) the vector of first order partial derivatives of Hat y. Let τ2=h(µ)0Σh(µ) and τ2 n=h(¯ Xn)0Σnh(¯ Xn), where the prime denotes transpose. If EkX1k3<∞,Fsatisfies Cram´er’s condition, His three times continuously differentiable in a neighborhood of µ,h(µ)6= 0 and τ > 0, Corollary 2 in Babu and Singh [8] shows that sup xn1/2¯¯¯P∗hn1/2©H(¯ X∗ n)−H(¯ Xn)ª/τn≤xi− −Phn1/2©H(¯ Xn)−H(µ)ª/τ ≤xi¯¯¯=o(1). (3) 6R. Pino Mej´ıas et al. Theorem 4.1 in Babu et al. [3] and Theorem 4(b) in Jim´enez-Gamero et al. [9] give the analogues of the above result for the Poisson bootstrap (with m=E(νn)) and the reduced bootstrap (with k1and k2satisfying conditions C.1.a and C.2.a), respectively. A class of statistics frequently used as pivots or approximate pivots is the studentized version of the above considered statistics. The result in (3) is also true for this class of statistics. Let τ∗2 n=h(¯ X∗ n)0Σ∗ nh(¯ X∗ n), where Σ∗ n is the sample variance matrix of the bootstrap sample. For a d-dimensional random vector X0= (X(1), X(2), ..., X(d)), let Wbe the d(d+ 3)/2-vector obtained from Xby adjoining it the products X(i)X(j), 1 ≤i≤j≤d, that is, W0= (X(1), X(2), ..., X(d), X2 (1), X(1)X(2), ..., X(1)X(d), X2 (2), X(2)X(3), ..., X2 (d)). If EkW1k3<∞,Fsatisfies Cram´er’s condition, His three times continuously differentiable in a neighborhood of µ,h(µ)6= 0 and τ > 0, then Theorem 4 in Babu and Singh [8] shows that sup xn−1/2¯¯¯P∗hn1/2©H(¯ X∗ n)−H(¯ Xn)ª/τ∗ n≤xi− −Phn1/2©H(¯ Xn)−H(µ)ª/τn≤xi¯¯¯=o(1). (4) Theorem 4.3 in Babu et al. [3] and Theorem 5(b) in Jim´enez-Gamero et al. [9] give the analogues of the result in (4) for the Poisson bootstrap (with m=E(νn)) and the reduced bootstrap (with k1and k2satisfying conditions C.1.a and C.2.a), respectively. As we said in the Introduction, Rao et al. [2] asserted that the variability of νnis not desirable. An example of this fact emerges when one considers statistics which are a studentized version of a smooth function of sample means, as the ones involved in (4), since bootstrap samples with all their components equal yield a degenerate value of the considered statistics. This undesirable situation may be avoided by using the sequential bootstrap or the reduced bootstrap with some k1>1. 3 Simulations In Section 2 we have seen that the ordinary bootstrap, the sequential bootstrap with m=E(νn) and the reduced bootstrap for adequate choices of k1and k2 yield distribution estimators having all of them the same asymptotic accuracy. To study and compare the finite sample performance of the corresponding estimators we have carried out a simulation experiment. The considered methods are displayed in Table 1, where [x] denotes the greatest integer less or equal than x. The Appendices A and B describe the employed algorithms to generate reduced bootstrap samples and Poisson bootstrap samples, respectively. Comparison of bootstrap procedures 7 Table 1. Methods. method description 1 ordinary bootstrap 2 reduced bootstrap with k1= [np −2√npq] + 1, k2= [np + 2√npq] 3 reduced bootstrap with k1= [np −√npq] + 1, k2= [np +√npq] 4 reduced bootstrap with k1=k2= [np]+1 5 sequential bootstrap with m= [np]+1 To study the corresponding distribution estimators of the statistic t1(X) = √n¯ Xn−µ sn , we have generated M= 1000 samples of size n= 10 from a standard normal population, N(0,1). For each method in Table 1 and from each sample Xm, 1≤m≤M, we have generated B= 1000 bootstrap samples, X∗m,1, ..., X∗m,B. We have considered the following intervals (−∞,−3),[−3,−2.9),[−2.9,−2.8),[−2.8,−2.7), ..., [2.9,3),[3,∞) (5) and, for each method, we have calculated f(I, m), the fraction of t1(X∗m,b), 1≤b≤B, falling on interval I. As a global measure of the bias of each distribution estimator we have considered BS =X I∈I bias2(I), where Iis the set of the intervals in (5), bias(I) = 1 MPM m=1 f(I, m)−f(I) and f(I) gives the fraction of T(Xm), 1 ≤m≤M, falling on the interval I; and as a global measure of the mean squared error of each distribution estimator we have considered MS =X I∈I mse(I), where mse(I) = 1 MPM m=1 {f(I, m)−f(I)}2. We have repeated the above experiment for n= 20,50,100,200,500,1000 and also for the statistics t(X) = √nH(¯ Xn)−H(µ) |h(¯ Xn)|sn , with H(x) = exp(x) (statistic t2) and H(x) = x2(statistic t3). For statistic t3, 8R. Pino Mej´ıas et al. Table 2. Results for statistic t1with data from a N(0,1) population. sample size method 10 20 50 100 200 500 1000 1BS ×1030.61 0.98 0.63 1.26 1.11 0.66 1.05 MS ×1032.07 2.00 1.61 2.24 2.09 1.64 2.01 2BS ×1030.54 0.96 0.64 1.28 1.10 0.66 1.06 MS ×1032.07 2.00 1.63 2.23 2.07 1.64 2.04 3BS ×1030.53 0.96 0.64 1.28 1.10 0.65 1.03 MS ×1032.07 2.00 1.61 2.25 2.07 1.63 2.01 4BS ×1038.44 3.69 2.58 3.34 2.85 2.30 2.59 MS ×10310.96 4.68 3.56 4.30 3.80 3.26 3.55 5BS ×1030.64 0.98 0.63 1.26 1.12 0.66 1.06 MS ×1031.94 2.00 1.60 2.23 2.10 1.65 2.03 Table 3. Results for statistic t2with data from a N(0,1) population. sample size method 10 20 50 100 200 500 1000 1BS ×1031.00 0.92 0.99 0.92 0.99 1.26 0.53 MS ×1032.98 2.02 1.98 1.90 1.98 2.23 1.51 2BS ×1031.17 0.85 0.99 0.92 0.97 1.24 0.54 MS ×1033.29 2.01 1.99 1.90 1.94 2.21 1.51 3BS ×1031.19 0.85 0.99 0.91 0.98 1.27 0.53 MS ×1033.34 2.00 2.00 1.88 1.95 2.23 1.51 4BS ×1039.62 6.59 3.67 3.92 3.36 2.21 2.21 MS ×10312.53 7.61 4.65 4.87 4.31 3.17 3.18 5BS ×1030.66 0.99 0.98 0.94 0.98 1.24 0.53 MS ×1032.50 2.10 1.98 1.91 1.96 2.21 1.49 Table 4. Results for statistic t3with data from a N(5,1) population. sample size method 10 20 50 100 200 500 1000 1BS ×1031.02 0.63 1.45 0.69 0.86 0.84 0.92 MS ×1032.52 1.66 2.42 1.65 1.83 1.82 1.89 2BS ×1031.19 0.65 1.42 0.70 0.84 0.82 0.90 MS ×1032.76 1.68 2.40 1.67 1.83 1.79 1.88 3BS ×1031.17 0.65 1.45 0.71 0.86 0.84 0.89 MS ×1032.73 1.69 2.41 1.68 1.84 1.82 1.86 4BS ×1036.39 3.40 3.70 2.30 2.79 2.18 2.15 MS ×1038.93 4.38 4.67 3.27 3.75 3.14 3.11 5BS ×1030.96 0.62 1.45 0.69 0.83 0.82 0.89 MS ×1032.27 1.65 2.42 1.66 1.79 1.80 1.87 to avoid h(µ) = 0, the data comes from a N(5, 1) population. The obtained results are displayed in Tables 2, 3 and 4. Tables 5, 6 and 7 present the results of a similar experiment with data from a mixture of two normal populations, 0.2N(0,1) + 0.8N(5,22). All computations in this paper have been performed using programs written in the R language (R Development Core Team [11]). Comparison of bootstrap procedures 9 Table 5. Results for statistic t1with data from a 0.2N(0,1) + 0.8N(5,22) population. sample size method 10 20 50 100 200 500 1000 1BS ×1031.05 0.83 0.99 1.18 0.65 1.06 1.24 MS ×1032.51 1.83 1.96 2.16 1.62 2.02 2.21 2BS ×1030.91 0.86 0.92 1.17 0.66 1.02 1.22 MS ×1032.46 1.87 1.90 2.14 1.64 1.98 2.19 3BS ×1030.93 0.86 0.93 1.16 0.67 1.03 1.24 MS ×1032.47 1.87 1.92 2.12 1.63 2.01 2.21 4BS ×1038.95 3.22 3.80 3.47 2.05 3.12 3.10 MS ×10311.71 4.21 4.78 4.43 3.02 4.09 4.05 5BS ×1031.06 0.82 0.98 1.20 0.64 1.02 1.25 MS ×1032.34 1.83 1.94 2.18 1.61 2.00 2.22 Table 6. Results for statistic t2with data from a 0.2N(0,1) + 0.8N(5,22) population. sample size method 10 20 50 100 200 500 1000 1BS ×1030.67 0.92 0.72 0.99 1.00 1.11 0.74 MS ×1032.75 2.17 1.73 1.97 1.98 2.07 1.71 2BS ×1030.75 1.02 0.80 1.08 1.00 1.09 0.72 MS ×1033.04 2.31 1.82 2.06 1.95 2.05 1.69 3BS ×1030.78 0.98 0.82 1.07 1.00 1.10 0.73 MS ×1033.00 2.28 1.83 2.04 1.97 2.06 1.71 4BS ×10316.18 6.77 3.54 3.15 3.44 3.38 2.24 MS ×10319.51 7.89 4.54 4.11 4.39 4.36 3.21 5BS ×1030.95 0.90 0.72 0.98 1.02 1.09 0.73 MS ×1032.88 2.14 1.74 1.97 1.98 2.06 1.71 Table 7. Results for statistic t3with data from a 0.2N(0,1) + 0.8N(5,22) population. sample size method 10 20 50 100 200 500 1000 1BS ×1031.03 0.60 0.96 0.79 1.08 1.37 0.80 MS ×1032.87 1.71 1.96 1.76 2.03 2.35 1.75 2BS ×1031.07 0.61 0.98 0.79 1.07 1.34 0.81 MS ×1033.03 1.74 1.97 1.76 2.04 2.32 1.77 3BS ×1031.07 0.62 0.97 0.78 1.07 1.33 0.81 MS ×1033.03 1.73 1.96 1.74 2.05 2.30 1.77 4BS ×10311.66 4.19 3.51 3.34 3.05 4.01 3.07 MS ×10314.65 5.19 4.48 4.30 4.01 4.97 4.03 5BS ×1031.44 0.55 0.97 0.81 1.05 1.36 0.81 MS ×1033.11 1.65 1.97 1.79 2.01 2.34 1.78 4 Conclusions Looking at Tables 2–7 we see that method 4 has the biggest bias for all the considered cases. For the rest of the methods in Table 1, we observe that there is no method behaving much better than the others for all statistics, all sample sizes and all distributions. As expected, the differences in BS and MS decrease as the sample size increases. Although the performance of methods 1, 2, 3, 5