R语⾔与抽样技术学习笔记（bootstrap）
R语⾔与抽样技术学习笔记（Randomize，Jackknife，bootstrap）
Bootstrap⽅法
  Bootstrap⼀词来源于西⽅神话故事“The adventures of Baron Munchausen”归结出的短语“to pull oneself up by one's bootstrap"，意味着不靠外界⼒量，依靠⾃⾝提升性能。
  Bootstrap的基本思想是：因为观测样本包含了潜在样本的全部的信息，那么我们不妨就把这个样本看做“总体”。那么相关的统计⼯作（估计或者检验）的统计量的分布可以从“总体”中利⽤Monte Carlo模拟得到。其做法可以简单地概括为：既然样本是抽出来的，那我何不从样本中再抽样。bootstrap检验方法
bootstrap基本⽅法
  1、采⽤重抽样技术从原始样本中抽取⼀定数量（⾃⼰给定）的样本，此过程允许重复抽样。
  2、根据抽出的样本计算给定的统计量T。
  3、重复上述N次（⼀般⼤于1000），得到N个统计量T。其均值可以视作统计量T的估计。
  4、计算上述N个统计量T的样本⽅差，得到统计量的⽅差。
  上述的估计我们可以看成是Bootstrap的⾮参数估计形式，它基本的思想是⽤频率分布直⽅图来估计概率分布。当然Bootstrap也有参数形式，在已知分布下，我们可以先利⽤总体样本估计出对应参数，再利⽤估计出的分布做Monte Carlo模拟，得到统计量分布的推断。  值得⼀提的是，参数化的Bootstrap⽅法虽然不够稳健，但是对于不平滑的函数，参数⽅法往往要⽐⾮参数办法好，当然这是基于你对样本的分布有⼀个初步了解的基础上的。
  例如：我们要考虑均匀分布\( U(\theta) \)的参数\( \theta \)的估计。我们采⽤似然估计。
data.sim <- runif(10)
theta.hat <- max(data.sim)
theta.boot1 <- replicate(1000, expr = {
y <- sample(data.sim, size = 10, replace = TRUE)
max(y)
})
theta.boot1.estimate <- mean(theta.boot1)
cat("the original estimate is ", theta.hat, "after bootstrap is ", theta.boot1.estimate)
## the original estimate is  0.7398 after bootstrap is  0.7138
hist(theta.boot1)
  从结果来看，倒不是说估计有多不好，只是说⽅差⽐较⼤，⽽且它的经验分布真的不太像真正的分布，这个近似很糟糕，导致的直接结果是⽅差也很⼤。
  如果采⽤参数⽅法，我们再来看看：
theta.boot2 <- replicate(1000, expr = {
y <- runif(1000, 0, theta.hat)
max(y)
})
theta.boot2.estimate <- mean(theta.boot2)
cat("the original estimate is ", theta.hat, "after bootstrap is ", theta.boot2.estimate)
## the original estimate is  0.7398 after bootstrap is  0.7391
hist(theta.boot2)
  结果从直⽅图来看是更优秀了，估计也更好⼀些，关键是⽅差变⼩了，从⾮参数的0.0402减少到了7.3944 × 10-4。bootstrap推断与bootstrap置信区间
  既然我们已经得到了Bootstrap估计量的经验分布函数，那么⼀个⾃然的结果就是我们可以利⽤这个分布对统计量做出⼀些统计推断。例如可以推测估计量的⽅差，估计量的偏差，估计量的置信区间等。
现在，我们就来考虑如何做Bootstrap的统计推断。
利⽤Bootstrap估计偏差
  既然Bootstrap将获得的样本样本看成了”总体“，那么估计量T⾃然是⼀个⽆偏的估计，Bootstrap数据集构造的”样本“的统计量T 与原始估计量T的偏差⾃然就是估计量偏差的⼀个很好的估计。
  具体做法是：
1. 从原始样本\( x_{1},\cdots,x_{n} \)中有放回的抽取n个样本构成⼀个Bootstrap数据集，重复这个过程m次，得到m个数据集。
2. 对于每个Bootstrap数据集，计算估计量T的值，记为\( T^{*j} \)。
3. \( T^{*j} \)的均值是T的⽆偏估计，⽽其与T的差是偏差的估计。
利⽤Bootstrap估计⽅差
  估计量T的⽅差的估计可以看做每个Bootstrap数据集的统计量T的值的⽅差。
  以我们遗留的问题，求1到100中随机抽取10个数的中位数的⽅差为例来说明。
n <- 10
x <- sample(1:100, size = n)
Mboot <- replicate(1000, expr = {
y <- sample(x, size = n, replace = TRUE)
median(y)
})
print(var(Mboot))
## [1] 334.2
  这个应该是⼀个正确的估计了。Efron指出要得到标准差的估计并不需要⾮常多的Bootstrap数据集（m不需要过分的⼤），通常50已经不错了，m>200是⽐较少见的（区间估计可能需要多⼀些）
  在R中，bootstrap包的函数bootstrap可以帮助你完成这⼀过程。bootstrap函数的调⽤格式如下：
bootstrap(x,nboot,theta,…, func=NULL) 
  参数说明：
x：原始抽样数据 
theta：统计量T 
nboot：构造Bootstrap数据集个数
library(bootstrap)
theta <- function(x) {
median(x)
}
results <- bootstrap(x, 100, theta)
print(var(results$thetastar))
## [1] 393.2
  可以看到两个的结果是相近的，所以，利⽤这个函数还是不错的选择。类似的还有boot包的boot函数。我们在相关数据的Bootstrap 推断中会⽤到。
相关数据的Bootstrap推断
回归数据的Bootstrap推断
  我们之所以可以采⽤Bootstrap去做这些估计，蕴含了⼀个很重要的假设，这些样本是近似iid的，然⽽我们不能保证需要推断的数据都是近似独⽴同分布的，对于相关数据的Bootstrap推断，我们常⽤的⽅法有配对的Bootstrap（paired Bootstrap）与残差法。
  先说paired Bootstrap，它的基本想法是，对于观测构成的数据框，虽然观测的每⼀⾏数据是相关的，但是每⾏是独⽴的，我们Bootstrap抽样，每次抽取⼀⾏，⽽不是单独的抽⼀个数即可。
例如，数据集women列出了美国⼥性的平均⾝⾼与体重，我们以体重为响应变量，⾝⾼为协变量进⾏回归，得到回归系数的估计。
  使⽤paired Bootstrap：
m <- 200
n <- nrow(women)
beta <- numeric(m)
for (b in 1:m) {
i <- sample(1:n, size = n, replace = TRUE)
weight <- women$weight[i]
height <- women$height[i]
beta[b] <- lm(weight ~ height)$coef[2]
}
cat("the estimate of beta is", lm(weight ~ height, data = women)$coef[2], "paired bootstrap estimate is",
mean(beta))
## the estimate of beta is 3.45 paired bootstrap estimate is 3.452
cat("the bias is", lm(weight ~ height, data = women)$coef[2] - mean(beta), "the stand error is",
sd(beta))
## the bias is -0.002468 the stand error is 0.126
  我们可以看到，估计量是⽆偏的，但是这个办法估计的⽅差变化较⼩，可能导致区间估计是不够稳健。我们可以利⽤boot包的boot函数来解决。
beta <- function(x, i) {
xi <- x[i, ]
coef(lm(xi[, 2] ~ xi[, 1]))[2]
}
library(boot)
obj <- boot(data = women, statistic = beta, R = 2000)
obj
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = women, statistic = beta, R = 2000)
##
##
## Bootstrap Statistics :
##    original    bias    std. error
## t1*    3.45 -0.002534      0.1208
  接下来我们说说残差法：
1. 由观测数据拟合模型.
2. 获得响应\( \hat{y}_{i} \)与残差\( \epsilon_{i} \)
3. 从残差数据集中有放回的抽取残差，构成Bootstrap残差数据集\( \hat{\epsilon}_{i} \)（这是近似独⽴的）
4. 构造伪响应\( Y_{i}^{*}=y_{i}+\hat{\epsilon}_{i} \)
5. 对x回归伪响应Y，得到希望得到的统计量，重复多次，得到希望的统计量的经验分布，利⽤它做统计推断
  我们将women数据集的例⼦利⽤残差法在做⼀次，R代码如下：
< <- lm(weight ~ height, data = women)
y.fit <- )
y.bootstrap <- rep(1, 15)
beta <- NULL
for (p in 1:100) {
for (i in 1:15) {
res <- s, 1, replace = TRUE)
y.bootstrap[i] <- y.fit[i] + res
}
beta[p] <- lm(y.bootstrap ~ women$height)$coef[2]
}
cat("the estimate of beta is", lm(weight ~ height, data = women)$coef[2], "paired bootstrap estimate is",
mean(beta))
## the estimate of beta is 3.45 paired bootstrap estimate is 3.436
cat("the bias is", lm(weight ~ height, data = women)$coef[2] - mean(beta), "the stand error is",
sd(beta))
## the bias is 0.01436 the stand error is 0.08561
  可以看到，利⽤残差法得到的⽅差更为稳健，做出的估计也更为的合理。
  这⾥需要指出⼀点，Bootstrap虽然可以处理相关数据，但是在变量筛选⽅⾯，其效果远不如Cross Validation准则好。
时间序列数据中的Bootstrap⽅法
  还有⼀类数据的相关性是上述假定也不满⾜的，那就是时间序列数列。那么如何利⽤Bootstrap来推断时间序列呢？我们以1947年–1991年美国GNP季度增长率数据为例进⾏说明。这个数据来⾃《⾦融时间序列分析》⼀书，数据可以在下载。
我们先来了解⼀下这个数据：

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。