全文链接：http://tecdat.cn/?p=21625
自抽样统计是什么以及为什么使用它？
本文将自抽样方法与传统方法进行比较，并了解它为何有用。并在R语言软件中通过对汽车速度和制动距离数据（查看文末了解数据获取方式）进行线性回归预测来实践它。
统计学是从数据中学习的科学。统计知识有助于收集数据的正确方法，使用正确的方法分析数据，并有效地呈现从数据中得出的结果。这些方法对于做出决策和预测至关重要，无论是预测消费者对产品的需求；使用文本挖掘过滤垃圾邮件；或在自动驾驶汽车中做出实时决策。
大多数时候，在进行研究时，从总体中收集数据是不切实际的。这可能是因为预算和/或时间限制以及其他因素。取而代之的是，获取总体的一个子集，并从该子集中收集洞察力以了解有关该人口的更多信息。
这意味着可以从适当抽取的样本中高效且相对低成本地获得适当准确的信息。
然而，许多因素会影响样本反映总体的程度。因此，结论的有效性和可靠性如何？
因此，让我们谈谈自抽样统计。
Bootstrapping 是一种统计程序，它对单个数据集重新采样以创建许多模拟样本。_Bootstrap_又称自展法、自举法、_自助法_、靴带法 。
这个过程允许计算标准误差、置信区间和假设检验”（_Forst）。_自抽样方法是传统假设检验方法的一种非常有用的替代方法，因为它相当简单，并且可以减少传统方法中遇到的一些缺陷，这将在后面讨论。统计推断通常依赖于样本分布和感兴趣特征的标准误差。传统方法（或大样本方法）抽取一个大小为_n的样本_从总体中提取，该样本用于计算总体估计值，然后进行推断。现在，实际上，只观察到了一个样本。但是，存在抽样分布的概念，如果要对总体进行重新抽样，它是所有可能估计值的理论集合。该理论指出，在样本量大等特定条件下，抽样分布近似正态分布，分布的标准差等于标准误。但是如果样本量不够大怎么办？那么，就不一定能假设理论上的抽样分布是正态的。这使得很难确定估计的标准误差，也更难从数据中得出合理的结论。
点击标题查阅往期内容
数据分享|R语言Bootstrap、百分位Bootstrap法抽样参数估计置信区间分析通勤时间和学生锻炼数据
01
02
03
04
与传统方法一样，大小为_n_的样本是从自抽样方法中的总体中抽取的。让我们将此样本称为_S_。然后，不是使用理论来确定所有可能的估计值，而是通过对观测值进行重采样并从_S_中替换_m_次来创建采样分布，每个重采样集具有_n 个_观测值。现在，如果抽样得当，_S_应该代表总体。因此，通过放回重采样_S __m_次，就好像_m_样本是从原始总体中抽取的，得出的估计值将代表传统方法下的理论分布。
现在我们了解了自抽样方法，必须注意得出的结果与传统方法的结果基本相同。此外，自抽样方法将始终有效，因为它不假设数据的任何潜在分布。这与理论上假设数据是正态分布的传统方法形成对比。
传统的程序需要一个满足特定假设的检验统计量才能获得有效的结果，这在很大程度上取决于实验设计。传统方法也使用理论来说明抽样分布应该是什么样子，但如果不满足理论的假设，结果就会崩溃。
另一方面，自抽样方法获取原始样本数据，然后对其重新采样以创建许多 [模拟] 样本。这种方法不依赖于理论，因为可以简单地观察采样分布，并且不必担心任何假设。这种技术可以准确估计统计数据，这在使用数据做出决策时至关重要。
样本量
在机器学习中，通常使用与原始数据集相同的样本大小。
bootstrap 样本与原始数据集大小相同。结果，一些样本将在引导样本中多次表示，而其他样本则根本不会被选中。
— 第 72 页，应用预测建模，2013 年。
如果数据集庞大且计算效率存在问题，则可以使用较小的样本，例如数据集大小的 50% 或 80%。
重复
重复次数必须足够大，以确保可以对样本计算有意义的统计数据，例如平均值、标准差和标准误差。
最少可以重复 20 或 30 次。可以使用较小的值将进一步增加对估计值样本计算的统计数据的方差。
理想情况下，在给定时间资源的情况下，估计样本应尽可能大，重复数以百计或数千次。
必须注意的是，增加重新采样的数量_m_不会增加数据中的信息量。也就是说，将原始集合重采样 100,000 次并不比仅重采样 1,000 次有用。集合中的信息量取决于样本大小_n_，它将在每次重新采样期间保持不变。因此，更多重采样的好处是可以更好地估计采样分布。
了解自抽样方法的工作原理后，一个合乎逻辑的问题是“自抽样方法是否过于依赖观察到的数据？” 这是一个很好的问题，因为重采样来自初始样本。正因为如此，假设异常值会扭曲重新采样的估计值是合乎逻辑的。虽然这是真的，但如果考虑到传统的做法，可以看出，数据集中的异常值也会扭曲均值并夸大估计的标准误差。因此，虽然人们可能会认为异常值会在重新采样的数据中多次出现并扭曲结果，从而使传统方法更好，但自抽样方法与传统方法一样依赖数据。bootstrapping 的优点是它是一种直接得出标准误差和置信区间估计值的方法，而且它很方便，因为它避免了重复实验以获取其他组采样数据的成本。尽管不可能知道大多数问题的真实置信区间_。_
R语言基于Bootstrap的线性回归预测置信区间估计方法
我们知道参数的置信区间的计算，这些都服从一定的分布(t分布、正态分布），因此在标准误前乘以相应的t分值或Z分值。但如果我们找不到合适的分布时，就无法计算置信区间了吗？幸运的是，有一种方法几乎可以用于计算各种参数的置信区间，这就是Bootstrap 法。
本文使用BOOTSTRAP来获得预测的置信区间。我们将在线性回归基础上讨论。
> reg=lm(dist~speed,data=cars) > points(x,predict(reg,newdata= data.frame(speed=x)))
这是一个单点预测。当我们想给预测一个置信区间时，预测的置信区间取决于参数估计误差。
预测置信区间
让我们从预测的置信区间开始
> for(s in 1:500){ + indice=sample(1:n,size=n, + replace=TRUE) + points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")
蓝色值是通过在我们的观测数据库中重新取样获得的可能预测值。值得注意的是，在残差正态性假设下（回归线的斜率和常数估计值），置信区间（90%）如下所示：
predict(reg,interval ="confidence",
在这里，我们可以比较500个生成数据集上的值分布，并将经验分位数与正态假设下的分位数进行比较，
> hist(Yx,proba=TRUE > boxplot(Yx,horizontal=TRUE > polygon(c( x ,rev(x I]))))
可以看出，经验分位数与正态假设下的分位数是可以比较的。
> quantile(Yx,c(.05,.95)) 5% 95% 58.63689 70.31281 + level=.9,newdata=data.frame(speed=x)) fit lwr upr 1 65.00149 59.65934 70.34364
感兴趣变量的可能值
现在让我们看看另一种类型的置信区间，关于感兴趣变量的可能值。这一次，除了提取新样本和计算预测外，我们还将在每次绘制时添加噪声，以获得可能的值。
> for(s in 1:500){ + indice=sample(1:n,size=n, + base=cars[indice,] + erreur=residuals(reg) + predict(reg,newdata=data.frame(speed=x))+E
在这里，我们可以（首先以图形方式）比较通过重新取样获得的值和在正态假设下获得的值，
> hist(Yx,proba=TRUE) > boxplot(Yx) abline(v=U[2:3) > polygon(c(D$x[I,rev(D$x[I])
数值上给出了以下比较
> quantile(Yx,c(.05,.95)) 5% 95% 44.43468 96.01357 U=predict(reg,interval ="prediction" fit lwr upr 1 67.63136 45.16967 90.09305
这一次，右侧有轻微的不对称。显然，我们不能假设高斯残差，因为有更大的正值，而不是负值。考虑到数据的性质，这是有意义的（制动距离不能是负数）。
然后开始讨论在供应中使用回归模型。为了获得具有独立性，有人认为必须使用增量付款的数据，而不是累计付款。
可以创建一个数据库，解释变量是行和列。
> base=data.frame( + y > head(base,12) y ai bj 1 3209 2000 0 2 3367 2001 0 3 3871 2002 0 4 4239 2003 0 5 4929 2004 0 6 5217 2005 0 7 1163 2000 1 8 1292 2001 1 9 1474 2002 1 10 1678 2003 1 11 1865 2004 1 12 NA 2005 1
然后，我们可以从基于对数增量付款数据的回归模型开始，该模型基于对数正态模型
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.9471 0.1101 72.188 6.35e-15 *** as.factor(ai)2001 0.1604 0.1109 1.447 0.17849 as.factor(ai)2002 0.2718 0.1208 2.250 0.04819 * as.factor(ai)2003 0.5904 0.1342 4.399 0.00134 ** as.factor(ai)2004 0.5535 0.1562 3.543 0.00533 ** as.factor(ai)2005 0.6126 0.2070 2.959 0.01431 * as.factor(bj)1 -0.9674 0.1109 -8.726 5.46e-06 *** as.factor(bj)2 -4.2329 0.1208 -35.038 8.50e-12 *** as.factor(bj)3 -5.0571 0.1342 -37.684 4.13e-12 *** as.factor(bj)4 -5.9031 0.1562 -37.783 4.02e-12 *** as.factor(bj)5 -4.9026 0.2070 -23.685 4.08e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1753 on 10 degrees of freedom (15 observations deleted due to missingness) Multiple R-squared: 0.9975, Adjusted R-squared: 0.9949 F-statistic: 391.7 on 10 and 10 DF, p-value: 1.338e-11 > exp(predict(reg1, + newdata=base)+summary(reg1)$sigma^2/2) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 2871.2 1091.3 41.7 18.3 7.8 21.3 [2,] 3370.8 1281.2 48.9 21.5 9.2 25.0 [3,] 3768.0 1432.1 54.7 24.0 10.3 28.0 [4,] 5181.5 1969.4 75.2 33.0 14.2 38.5 [5,] 4994.1 1898.1 72.5 31.8 13.6 37.1 [6,] 5297.8 2013.6 76.9 33.7 14.5 39.3 > sum(py[is.na(y)]) [1] 2481.857
这与链式梯度法的结果略有不同，但仍然具有可比性。我们也可以尝试泊松回归（用对数链接）
glm(y~ + as.factor(ai)+ + as.factor(bj),data=base, + family=poisson) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 8.05697 0.01551 519.426 < 2e-16 *** as.factor(ai)2001 0.06440 0.02090 3.081 0.00206 ** as.factor(ai)2002 0.20242 0.02025 9.995 < 2e-16 *** as.factor(ai)2003 0.31175 0.01980 15.744 < 2e-16 *** as.factor(ai)2004 0.44407 0.01933 22.971 < 2e-16 *** as.factor(ai)2005 0.50271 0.02079 24.179 < 2e-16 *** as.factor(bj)1 -0.96513 0.01359 -70.994 < 2e-16 *** as.factor(bj)2 -4.14853 0.06613 -62.729 < 2e-16 *** as.factor(bj)3 -5.10499 0.12632 -40.413 < 2e-16 *** as.factor(bj)4 -5.94962 0.24279 -24.505 < 2e-16 *** as.factor(bj)5 -5.01244 0.21877 -22.912 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 46695.269 on 20 degrees of freedom Residual deviance: 30.214 on 10 degrees of freedom (15 observations deleted due to missingness) AIC: 209.52 Number of Fisher Scoring iterations: 4 > predict(reg2, newdata=base,type="response") > sum(py2[is.na(y)]) [1] 2426.985
预测结果与链式梯度法得到的估计值吻合。克劳斯·施密特（Klaus Schmidt）和安吉拉·温什（Angela Wünsche）于1998年在链式梯度法、边际和最大似然估计中建立了与最小偏差方法的联系。
