经验首页 前端设计 程序设计 Java相关 移动开发 数据库/运维 软件/图像 大数据/云计算 其他经验
当前位置:技术经验 » 程序设计 » R语言 » 查看文章
R语言学习笔记缺失数据的Bootstrap与Jackknife方法
来源:jb51  时间:2021/11/8 19:51:41  对本文有异议

一、题目

请添加图片描述

下面再加入缺失的情况来继续深入探讨,同样还是如习题1.6的构造方式来加入缺失值,其中a=2, b = 0

请添加图片描述

我们将进行如下几种操作:

请添加图片描述

二、解答

a)Bootstrap与Jackknife进行估计

首先构建生成数据函数。

  1. # 生成数据
  2. # 生成数据
  3. GenerateData <- function(a = 0, b = 0) {
  4. y <- matrix(nrow = 3, ncol = 100)
  5. z <- matrix(rnorm(300), nrow = 3)
  6. y[1, ] <- 1 + z[1, ]
  7. y[2, ] <- 5 + 2 * z[1, ] + z[2, ]
  8. u <- a * (y[1, ] - 1) + b * (y[2, ] - 5) + z[3, ]
  9. # m2 <- 1 * (u < 0)
  10. y[3, ] <- y[2, ]
  11. y[3, u < 0] <- NA
  12. dat_comp <- data.frame(y1 = y[1, ], y2 = y[2, ])
  13. dat_incomp <- data.frame(y1 = y[1, ], y2 = y[3, ])
  14. # dat_incomp <- na.omit(dat_incomp)
  15. return(list(dat_comp = dat_comp, dat_incomp = dat_incomp))
  16. }

Bootstrap与Jackknife的函数:

  1. Bootstrap1 <- function(Y, B = 200, fun) {
  2. Y_len <- length(Y)
  3. mat_boots <- matrix(sample(Y, Y_len * B, replace = T), nrow = B, ncol = Y_len)
  4. statis_boots <- apply(mat_boots, 1, fun)
  5. boots_mean <- mean(statis_boots)
  6. boots_sd <- sd(statis_boots)
  7. return(list(mean = boots_mean, sd = boots_sd))
  8. }
  9.  
  10. Jackknife1 <- function(Y, fun) {
  11. Y_len <- length(Y)
  12. mat_jack <- sapply(1:Y_len, function(i) Y[-i])
  13. redu_samp <- apply(mat_jack, 2, fun)
  14. jack_mean <- mean(redu_samp)
  15. jack_sd <- sqrt(((Y_len - 1) ^ 2 / Y_len) * var(redu_samp))
  16. return(list(mean = jack_mean, sd = jack_sd))
  17. }

进行重复试验所需的函数:

  1. RepSimulation <- function(seed = 2018, fun) {
  2. set.seed(seed)
  3. dat <- GenerateData()
  4. dat_comp_y2 <- dat$dat_comp$y2
  5. boots_sd <- Bootstrap1(dat_comp_y2, B = 200, fun)$sd
  6. jack_sd <- Jackknife1(dat_comp_y2, fun)$sd
  7. return(c(boots_sd = boots_sd, jack_sd = jack_sd))
  8. }

下面重复100次实验进行 Y2​的均值与变异系数标准差的估计:

  1. nrep <- 100
  2. ## 均值
  3. fun = mean
  4. mat_boots_jack <- sapply(1:nrep, RepSimulation, fun)
  5. apply(mat_boots_jack, 1, function(x) paste(round(mean(x), 3), '±', round(sd(x), 3)))
  1. ## 变异系数
  2. fun = function(x) sd(x) / mean(x)
  3. mat_boots_jack <- sapply(1:nrep, RepSimulation, fun)
  4. apply(mat_boots_jack, 1, function(x) paste(round(mean(x), 3), '±', round(sd(x), 3)))

从上面可以发现,Bootstrap与Jackknife两者估计结果较为相近,其中对均值标准差的估计,Jackknife的方差更小。这其实较为符合常识:Jackknife估计每次只取出一个样本,用剩下的样本来作为样本整体;而Bootstrap每次都会比较随机地重抽样,随机性相对较高,所以重复100次模拟实验,导致其方差相对较大。

下面我们用计算公式来进行推导。

b)均值与变异系数(大样本)的标准差解析式推导与计算

均值

请添加图片描述

变异系数(大样本近似)

  1. ## 变异系数
  2. sd(sapply(1:10000, function(x) {
  3. set.seed(x)
  4. dat <- GenerateData(a = 0, b = 0)
  5. sd(dat$dat_comp$y2) / mean(dat$dat_comp$y2)
  6. }))

变异系数大样本近似值为:0.03717648,说明前面的Bootstrap与Jackknife两种方法估计的都较为准确。

c)缺失插补后的Bootstrap与Jackknife

构造线性填补的函数,并进行线性填补。

  1. DatImputation <- function(dat_incomp) {
  2. dat_imp <- dat_incomp
  3. lm_model = lm(y2 ~ y1, data = na.omit(dat_incomp))
  4. # 找出y2缺失对应的那部分data
  5. na_ind = is.na(dat_incomp$y2)
  6. na_dat = dat_incomp[na_ind, ]
  7. # 将缺失数据进行填补
  8. dat_imp[na_ind, 'y2'] = predict(lm_model, na_dat)
  9. return(dat_imp)
  10. }
  11.  
  12. dat <- GenerateData(a = 2, b = 0)
  13. dat_imp <- DatImputation(dat$dat_incomp)
  1. fun = mean
  2. Bootstrap1(dat_imp$y2, B = 200, fun)$sd
  1. Jackknife1(dat_imp$y2, fun)$sd
  1. fun = function(x) sd(x) / mean(x)
  2. Bootstrap1(dat_imp$y2, B = 200, fun)$sd
  1. Jackknife1(dat_imp$y2, fun)$sd

Bootstrap与Jackknife的填补结果,很大一部分是由于数据的缺失会造成距离真实值较远。但单从两种方法估计出来的值比较接近。

c)缺失插补前的Bootstrap与Jackknife

先构建相关的函数:

  1. Array2meancv <- function(j, myarray) {
  2. dat_incomp <- as.data.frame(myarray[, j, ])
  3. names(dat_incomp) <- c('y1', 'y2')
  4. dat_imp <- DatImputation(dat_incomp)
  5. y2_mean <- mean(dat_imp$y2)
  6. y2_cv <- sd(dat_imp$y2) / y2_mean
  7. return(c(mean = y2_mean, cv = y2_cv))
  8. }
  9.  
  10. Bootstrap_imp <- function(dat_incomp, B = 200) {
  11. n <- nrow(dat_incomp)
  12. array_boots <- array(dim = c(n, B, 2))
  13. mat_boots_ind <- matrix(sample(1:n, n * B, replace = T), nrow = B, ncol = n)
  14.  
  15. array_boots[, , 1] <- sapply(1:B, function(i) dat_incomp$y1[mat_boots_ind[i, ]])
  16. array_boots[, , 2] <- sapply(1:B, function(i) dat_incomp$y2[mat_boots_ind[i, ]])
  17. mean_cv_imp <- sapply(1:B, Array2meancv, array_boots)
  18. boots_imp_mean <- apply(mean_cv_imp, 1, mean)
  19. boots_imp_sd <- apply(mean_cv_imp, 1, sd)
  20. return(list(mean = boots_imp_mean, sd = boots_imp_sd))
  21. }
  22.  
  23. Jackknife_imp <- function(dat_incomp) {
  24. n <- nrow(dat_incomp)
  25. array_jack <- array(dim = c(n - 1, n, 2))
  26. array_jack[, , 1] <- sapply(1:n, function(i) dat_incomp[-i, 'y1'])
  27. array_jack[, , 2] <- sapply(1:n, function(i) dat_incomp[-i, 'y2'])
  28. mean_cv_imp <- sapply(1:n, Array2meancv, array_jack)
  29. jack_imp_mean <- apply(mean_cv_imp, 1, mean)
  30. jack_imp_sd <- apply(mean_cv_imp, 1, function(x) sqrt(((n - 1) ^ 2 / n) * var(x)))
  31. return(list(mean = jack_imp_mean, sd = jack_imp_sd))
  32. }

然后看看两种方式估计出来的结果:

  1. Bootstrap_imp(dat$dat_incomp)$sd
  1. Jackknife_imp(dat$dat_incomp)$sd

缺失插补前进行Bootstrap与Jackknife也还是有一定的误差,标准差都相对更大,表示波动会比较大。具体表现情况下面我们多次重复模拟实验,通过90%置信区间来看各个方法的优劣。

d)比较各种方式的90%置信区间情况(重复100次实验)

  1. RepSimulationCI <- function(seed = 2018, stats = 'mean') {
  2. mean_true <- 5
  3. cv_true <- sqrt(5) / 5
  4. myjudge <- function(x, value) {
  5. return(ifelse((x$mean - qnorm(0.95) * x$sd < value) & (x$mean + qnorm(0.95) * x$sd > value), 1, 0))
  6. }
  7. if(stats == 'mean') {
  8. fun = mean
  9. value = mean_true
  10. } else if(stats == 'cv') {
  11. fun = function(x) sd(x) / mean(x)
  12. value = cv_true
  13. }
  14. set.seed(seed)
  15. boots_after_ind <- boots_before_ind <- jack_after_ind <- jack_before_ind <- 0
  16. dat <- GenerateData(a = 2, b = 0)
  17. dat_incomp <- dat$dat_incomp
  18. # after imputation
  19. dat_imp <- DatImputation(dat_incomp)
  20.  
  21. boots_after <- Bootstrap1(dat_imp$y2, B = 200, fun)
  22. boots_after_ind <- myjudge(boots_after, value)
  23. jack_after <- Jackknife1(dat_imp$y2, fun)
  24. jack_after_ind <- myjudge(jack_after, value)
  25. # before imputation
  26. boots_before <- Bootstrap_imp(dat_incomp)
  27. jack_before <- Jackknife_imp(dat_incomp)
  28. if(stats == 'mean') {
  29. boots_before$mean <- boots_before$mean[1]
  30. boots_before$sd <- boots_before$sd[1]
  31. jack_before$mean <- jack_before$mean[1]
  32. jack_before$sd <- jack_before$sd[1]
  33. } else if(stats == 'cv') {
  34. boots_before$mean <- boots_before$mean[2]
  35. boots_before$sd <- boots_before$sd[2]
  36. jack_before$mean <- jack_before$mean[2]
  37. jack_before$sd <- jack_before$sd[2]
  38. }
  39. boots_before_ind <- myjudge(boots_before, value)
  40. jack_before_ind <- myjudge(jack_before, value)
  41. return(c(boots_after = boots_after_ind,
  42. boots_before = boots_before_ind,
  43. jack_after = jack_after_ind,
  44. jack_before = jack_before_ind))
  45. }
  46.  

重复100次实验,均值情况:

  1. nrep <- 100
  2. result_mean <- apply(sapply(1:nrep, RepSimulationCI, 'mean'), 1, sum)
  3. names(result_mean) <- c('boots_after', 'boots_before', 'jack_after', 'jack_before')
  4. result_mean

变异系数情况:

  1. result_cv <- apply(sapply(1:nrep, RepSimulationCI, 'cv'), 1, sum)
  2. names(result_cv) <- c('boots_after', 'boots_before', 'jack_after', 'jack_before')
  3. result_cv

上面的数字越表示90%置信区间覆盖真实值的个数,数字越大表示覆盖的次数越多,也就说明该方法会相对更好。

填补之前进行Bootstrap或Jackknife

无论是均值还是变异系数,通过模拟实验都能看出,在填补之前进行Bootstrap或Jackknife,其估计均会远优于在填补之后进行Bootstrap或Jackknife。而具体到Bootstrap或Jackknife,这两种方法相差无几。

填补之后进行Bootstrap或Jackknife

在填补之后进行Bootstrap或Jackknife,效果都会很差,其实仔细思考后也能够理解,本身缺失了近一半的数据,然后填补会带来很大的偏差,此时我们再从中抽样,有很大可能抽出来的绝大多数都是原本填补的有很大偏差的样本,这样估计就会更为不准了。

当然,从理论上说,填补之前进行Bootstrap或Jackknife是较为合理的,这样对每个Bootstrap或Jackknife样本,都可以用当前的观测值去填补当前的缺失值,这样每次填补可能花费的时间将对较长,但实际却更有效。

以上就是R语言学习笔记缺失数据的Bootstrap与Jackknife方法的详细内容,更多关于R语言学习笔记的资料请关注w3xue其它相关文章!

 友情链接:直通硅谷  点职佳  北美留学生论坛

本站QQ群:前端 618073944 | Java 606181507 | Python 626812652 | C/C++ 612253063 | 微信 634508462 | 苹果 692586424 | C#/.net 182808419 | PHP 305140648 | 运维 608723728

W3xue 的所有内容仅供测试,对任何法律问题及风险不承担任何责任。通过使用本站内容随之而来的风险与本站无关。
关于我们  |  意见建议  |  捐助我们  |  报错有奖  |  广告合作、友情链接(目前9元/月)请联系QQ:27243702 沸活量
皖ICP备17017327号-2 皖公网安备34020702000426号