srcmini本文概述
Bootstrap是一种使用样本数据推断总体的方法。布拉德利·埃夫隆(Bradley Efron)于1979年在论文中首次介绍了它。Bootstrap依赖于采样, 并从样本数据中进行替换。该技术可用于估计任何统计信息的标准误差, 并为其获取置信区间(CI)。如果CI没有封闭的形式, 或者具有非常复杂的形式, 则Bootstrap尤其有用。
假设我们有一个n个元素的样本：X = {x1, x2, …, xn}, 我们对CI的一些统计量T = t(X)感兴趣。 Bootstrap框架很简单。我们只重复R次以下方案：对于第i次重复, 从可用样本中替换n个元素的样本(其中一些元素将被多次提取)。将此新样本称为第i个Bootstrap样本Xi, 并计算所需的统计量Ti = t(Xi)。
结果, 我们将获得统计数据的R值：T1, T2, …, TR。我们称它们为T的Bootstrap实现或T的Bootstrap分布。基于此, 我们可以计算T的CI。有多种实现方法。采取百分位数似乎是最简单的方法。
让我们(再次使用)著名的虹膜数据集。看一下前几行：
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
假设我们要找到中位数Sepal.Length, 中位数Sepal.Width和Spearman的秩相关系数的CI。我们将使用R的启动程序包和一个名为… boot的函数。要使用其功能, 我们必须创建一个函数, 该函数可以根据重采样数据计算统计信息。它至少应具有两个参数：一个数据集和一个向量, 其中包含一个数据集中元素的索引, 这些元素被挑选来创建引导程序样本。
如果我们希望一次计算多个统计量的CI, 则我们的函数必须将它们作为单个向量返回。
对于我们的示例, 它可能看起来像这样：
library(boot)
foo <- function(data, indices){
dt<-data[indices, ]
c(
cor(dt[, 1], dt[, 2], method='s'), median(dt[, 1]), median(dt[, 2])
)
}
foo从数据中选择所需的元素(将数字存储在索引中)并计算前两列的相关系数(方法=’s’用于选择Spearman的秩系数, method =’p’将产生Pearson系数)及其中位数。
我们还可以添加一些额外的参数, 例如。让用户选择od类型的相关系数：
foo <- function(data, indices, cor.type){
dt<-data[indices, ]
c(
cor(dt[, 1], dt[, 2], method=cor.type), median(dt[, 1]), median(dt[, 2])
)
}
为了使本教程更加通用, 我将使用foo的最新版本。
现在, 我们可以使用启动功能了。我们必须为其提供名称od数据集, 刚刚创建的函数, 重复次数(R)以及函数的其他任何参数(例如cor.type)。下面, 我使用set.seed来重现此示例。
set.seed(12345)
myBootstrap <- boot(iris, foo, R=1000, cor.type='s')
引导功能返回一个名为…(是的, 你是对的！)引导的类的对象。它有两个有趣的元素。 $ t包含由引导程序(T的引导实现)生成的统计值的R值：
head(myBootstrap$t)
## [, 1] [, 2] [, 3]
## [1, ] -0.26405188 5.70 3
## [2, ] -0.12973299 5.80 3
## [3, ] -0.07972066 5.75 3
## [4, ] -0.16122705 6.00 3
## [5, ] -0.20664808 6.00 3
## [6, ] -0.12221170 5.80 3
$ t0包含原始, 完整数据集中的统计值：
myBootstrap$t0
## [1] -0.1667777 5.8000000 3.0000000
将启动对象打印到控制台可提供一些进一步的信息：
myBootstrap
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = foo, R = 1000, cor.type = "s")
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.1667777 0.002546391 0.07573983
## t2* 5.8000000 -0.013350000 0.10295571
## t3* 3.0000000 0.007900000 0.02726414
original与$ t0相同。偏差是引导实现的平均值(来自$ t的引导实现)之间的差, 称为T的引导估计和原始数据集中的值(来自$ t0的引导估计)。
colMeans(myBootstrap$t)-myBootstrap$t0
## [1] 0.002546391 -0.013350000 0.007900000
标准误差是引导程序估计的标准误差, 等于引导程序实现的标准偏差。
apply(myBootstrap$t, 2, sd)
## [1] 0.07573983 0.10295571 0.02726414
在开始使用CI之前, 始终值得一看引导实现的分布。我们可以使用plot函数, 通过索引告诉我们希望查看foo中计算的统计量。此处index = 1是萼片长度和宽度之间的Spearman相关系数, index = 2是中位数od sepal长度, index = 3是中位数od sepal宽度。
plot(myBootstrap, index=1)
Bootstrap相关系数的分布似乎很正常。让我们为其找到CI。我们可以使用boot.ci。默认值为95％的配置项, 但可以使用conf参数进行更改。
boot.ci(myBootstrap, index=1)
## Warning in boot.ci(myBootstrap, index = 1): bootstrap variances needed for
## studentized intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 1)
##
## Intervals :
## Level Normal Basic
## 95% (-0.3178, -0.0209 ) (-0.3212, -0.0329 )
##
## Level Percentile BCa
## 95% (-0.3007, -0.0124 ) (-0.3005, -0.0075 )
## Calculations and Intervals on Original Scale
boot.ci提供5种类型的引导CI。其中之一, 学生间隔时间是唯一的。它需要估计Bootstrap方差。我们没有提供它, 因此R会显示警告：学生化间隔所需的引导程序差异。方差估计可以通过二级引导程序获得, 也可以通过折刀技术轻松获得。这在某种程度上超出了本教程的范围, 因此让我们集中讨论其余四种类型的引导CI。
如果我们不想全部看到它们, 可以在type参数中选择相关的参数。可能的值为norm, basic, stud, perc, bca或这些的向量。
boot.ci(myBootstrap, index=1, type=c('basic', 'perc'))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, type = c("basic", "perc"), index = 1)
##
## Intervals :
## Level Basic Percentile
## 95% (-0.3212, -0.0329 ) (-0.3007, -0.0124 )
## Calculations and Intervals on Original Scale
boot.ci函数创建类的对象…(是的, 你猜对了！)bootci。就像在类型实参中使用的CI类型一样, 调用其元素。 $ norm是包含置信度水平和CI边界的3元素向量。
boot.ci(myBootstrap, index=1, type='norm')$norm
## conf
## [1, ] 0.95 -0.3177714 -0.02087672
$ basic, $ stud, $ perc和$ bca是由5个元素组成的向量, 还包含用于计算CI的百分位数(我们稍后会再介绍)：
boot.ci(myBootstrap, index=1, type='basic')$basic
## conf
## [1, ] 0.95 975.98 25.03 -0.3211981 -0.03285178
一点符号(对不起！)
要了解不同类型的CI是什么, 我们需要引入一些符号。所以让：
- t⋆是引导估计(引导实现的平均值),
- t0是原始数据集中的统计值,
- 这是引导估计的标准误差,
- b是Bootstrap估计的偏差b =t⋆-t0
- α为置信度, 通常为α= 0.95,
- zα是标准正态分布的$ 1- \ frac \ alpha 2 $-分位数,
- θα是引导实现的分布的α百分数。
百分位数CI
使用上述符号, 百分位数CI为：
(θ(1-α)/ 2, θ1-(1-α)/ 2)
所以这只需要相关的百分位数。而已。
典型的Wald CI为：
t0±zα⋅se⋆
但是在Bootstrap情况下, 我们应该对其进行校正以消除偏差。这样就变成：
$$ t_0-b \ pm z_ alpha \ cdot se ^ \ star \\ 2t_0-t ^ \ star \ pm alpha_ cdot se ^ \ star $$
通常不建议使用百分位数CI, 因为它在处理怪异的分布时表现不佳。基本配置项(也称为关键配置项或经验配置项)对此更为健壮。其背后的原理是计算每个引导复制和t0之间的差异, 并使用其分布的百分位数。可以找到完整的详细信息, 例如在L. Wasserman的所有统计资料中
基本CI的最终公式为：
(2t0-θ1-(1-α)/ 2, 2t0-θ(1-α)/ 2)
BCα来自偏差校正, 加速。它的公式不是很复杂, 但是有点不直观, 因此我将跳过它。如果你对详细信息感兴趣, 请参阅Thomas J. DiCiccio和Bradley Efron的文章。
方法名称中提到的加速要求使用引导实现的特定百分比。有时可能会出现一些极端的百分位数, 可能是异常值。在这种情况下, BCα可能不稳定。
让我们来看看BCαCI的花瓣宽度中位数。在原始数据集中, 该中值正好等于3。
boot.ci(myBootstrap, index=3)
## Warning in boot.ci(myBootstrap, index = 3): bootstrap variances needed for
## studentized intervals
## Warning in norm.inter(t, adj.alpha): extreme order statistics used as
## endpoints
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 3)
##
## Intervals :
## Level Normal Basic
## 95% ( 2.939, 3.046 ) ( 2.900, 3.000 )
##
## Level Percentile BCa
## 95% ( 3.0, 3.1 ) ( 2.9, 2.9 )
## Calculations and Intervals on Original Scale
## Warning : BCa Intervals used Extreme Quantiles
## Some BCa intervals may be unstable
我们得到BCαCI(2.9, 2.9)。这很奇怪, 但是幸运的是R告诫我们极端顺序统计用作端点。让我们看看这里到底发生了什么：
plot(myBootstrap, index=3)
引导实现的分布并不常见。其中绝大多数(超过90％)为3秒。
table(myBootstrap$t[, 3])
##
## 2.9 2.95 3 3.05 3.1 3.15 3.2
## 1 1 908 24 63 1 2
在这种情况下, BCα方法中的”加速度”很容易”跳”至极值。在此, 构成BCαCI基础的百分位数CI为(3, 3.1)。现在, 当我们将其”向左”扩展时, 必须使用0.002％。 “向右”的扩展也达到极限：0.997％。
有时, 我们需要重新创建引导复制。如果我们可以使用R, 这不是问题。 set.seed函数解决了这个问题。在其他软件中重新创建引导程序复制将更加复杂。而且, 对于大的R, R中的重新计算也不是一种选择(例如, 由于时间不足)。
我们可以解决这个问题, 保存形成每个引导程序样本的原始数据集元素的索引。这就是boot.array函数(带有index = T参数)的作用。
tableOfIndices<-boot.array(myBootstrap, indices=T)
每行是一个引导程序样本。例如, 我们的第一个示例包含以下元素：
tableOfIndices[1, ]
## [1] 109 12 143 65 41 28 105 62 23 102 37 34 55 53 38 3 11
## [18] 146 31 122 85 141 118 116 117 78 100 13 98 49 59 88 24 56
## [35] 136 4 59 67 126 118 104 101 70 13 12 19 21 149 73 133 67
## [52] 86 6 88 131 105 121 145 121 83 70 68 71 111 76 73 122 116
## [69] 85 144 114 139 89 124 64 27 81 78 58 39 17 50 37 117 76
## [86] 97 114 64 17 93 141 65 124 80 137 54 37 57 70 128 10 55
## [103] 13 53 59 45 116 14 77 118 108 138 50 78 49 104 49 1 85
## [120] 28 43 82 74 64 33 55 32 59 62 112 8 11 96 30 14 30
## [137] 38 85 66 85 97 107 4 18 76 35 31 133 27 69
设置index = F(默认值), 我们将得到以下问题的答案：”每个引导样本中原始数据集的每个元素出现了多少次？”。例如, 在第一个样本中：数据集的第一个元素出现一次, 根本没有第二个元素出现, 第三个元素出现一次, 第四个元素出现两次, 依此类推。
tableOfAppearances<-boot.array(myBootstrap)
tableOfAppearances[1, ]
## [1] 1 0 1 2 0 1 0 1 0 1 2 2 3 2 0 0 2 1 1 0 1 0 1 1 0 0 2 2 0 2 2 1 1 1 1
## [36] 0 3 2 1 0 1 0 1 0 1 0 0 0 3 2 0 0 2 1 3 1 1 1 4 0 0 2 0 3 2 1 2 1 1 3
## [71] 1 0 2 1 0 3 1 3 0 1 1 1 1 0 5 1 0 2 1 0 0 0 1 0 0 1 2 1 0 1 1 1 0 2 2
## [106] 0 1 1 1 0 1 1 0 2 0 3 2 3 0 0 2 2 0 2 0 1 0 1 0 0 1 0 2 0 0 1 1 1 1 0
## [141] 2 0 1 1 1 1 0 0 1 0
这样的表使我们可以在R外部重新创建引导实现。或者, 当我们不想使用set.seed并再次执行所有计算时, 可以在R自身中重新创建引导实现。
onceAgain<-apply(tableOfIndices, 1, foo, data=iris, cor.type='s')
让我们检查结果是否相同：
head(t(onceAgain))
## [, 1] [, 2] [, 3]
## [1, ] -0.26405188 5.70 3
## [2, ] -0.12973299 5.80 3
## [3, ] -0.07972066 5.75 3
## [4, ] -0.16122705 6.00 3
## [5, ] -0.20664808 6.00 3
## [6, ] -0.12221170 5.80 3
head(myBootstrap$t)
## [, 1] [, 2] [, 3]
## [1, ] -0.26405188 5.70 3
## [2, ] -0.12973299 5.80 3
## [3, ] -0.07972066 5.75 3
## [4, ] -0.16122705 6.00 3
## [5, ] -0.20664808 6.00 3
## [6, ] -0.12221170 5.80 3
all(t(onceAgain)==myBootstrap$t)
## [1] TRUE
对, 他们是！
如果你想了解有关R中机器学习的更多信息, 请参加srcmini的机器学习工具箱课程。
评论前必须登录！
注册