自助法示例

心灵

1. n <- 1000 #样本数据集数量
   
2. B <- 1000 # boot.samples
   
3. alpha <- 0.05 #置信概率
   
4. df <- 6
   
5. true.mu <- 0 # 真是均值
   
6. true.sd <- 1 # 真实标准差
   
7. 
   
8. # 模拟数据生成
   
9. x <- rnorm(n,mean = true.mu, sd = true.sd)
   
10. 
    
11. sample.mean.x <- mean(x) #样本均值
    
12. sample.sd.x <- sd(x) #样本标准差
    
13. sample.se.x <- sample.sd.x/sqrt(n) #样本标准误
    
14. sample.conf.x <- sample.mean.x + c(-1,1)*qnorm(alpha/2,lower.tail=F)*sample.se.x #样本95%置信区间
    
15. 
    
16. # 为自助法预留缓存向量
    
17. boot.x <- vector(mode="numeric",length=B)
    
18. boot.mean <- vector(mode="numeric",length=B)
    
19. boot.sd <- vector(mode="numeric",length=B)
    
20. boot.t <- vector(mode="numeric",length=B)
    
21. 
    
22. for(i in 1:B)
    
23. {
    
24.   boot.x <- sample(x, size=n, replace=T)
    
25.   boot.mean[i] <- mean(boot.x)
    
26.   boot.sd[i] <- sd(boot.x)
    
27.   boot.t[i] <- (sample.mean.x-boot.mean[i])/(boot.sd[i]/sqrt(n))
    
28. }
    
29. 
    
30. # bootstrap 统计量
    
31. boot.mean.x <- mean(boot.mean)
    
32. boot.sd.x <- mean(boot.sd)
    
33. boot.se.x <- boot.sd.x/sqrt(B)
    
34. boot.t.quantile <- quantile(boot.t,probs = c(alpha/2,1-alpha/2))
    
35. boot.mean.quantile <- quantile(boot.mean,probs = c(alpha/2,1-alpha/2))
    
36. 
    
37. boot.conf.x <- boot.mean.x + c(-1,1)*qt(alpha/2,df=B-1,lower.tail=F)*boot.se.x
    
38. boot.conf.x2 <- boot.mean.x + boot.t.quantile*boot.se.x
    
39. 
    
40. true.conf.x <- 2*sample.mean.x - c(boot.mean.quantile[2],boot.mean.quantile[1])
    
41. 
    
42. #boot.mean.x
    
43. #boot.sd.x
    
44. #boot.se.x
    
45. #boot.conf.x
    
46. #boot.conf.x2
    
47. 
    
48. #true.conf.x
    
49. 
    
50. #sample.mean.x
    
51. #sample.sd.x
    
52. #sample.se.x
    
53. #sample.conf.x
    
54. 
    
55. 
    
56. hist(boot.mean, col = "lightblue", border = "pink",freq = F, main="Normal distribution")
    
57. lines(density(boot.mean),col="blue", lwd=2)
    
58. abline(v=0, col="red",lwd=2)
    
59. 
    
60. message("boot.mean.x = ", format(boot.mean.x,digits = 4),"           sample.mean.x = ", format(sample.mean.x,digits = 4),"\n",
    
61.         "boot.sd.x = ", format(boot.sd.x,digits = 4),"               sample.sd.x = ", format(sample.sd.x,digits = 4),"\n",
    
62.         "boot.conf.x = ", format(boot.conf.x,digits = 4),"      sample.conf.x = ", format(sample.conf.x,digits = 4),"\n",
    
63.         "boot.conf.x2 = ", format(boot.conf.x2,digits = 4),"      true.conf.x = ", format(true.conf.x,digits = 4),"\n")
    
64. 
    

复制代码

这段代码主要简单得阐述了bootstrap技术，并对比了在bootstrap统计量和样本统计量以及（真实）统计量的差异。

ｙōｙō琴緣

不错哦 ，么么哒

Freedzhang

{:soso_e179:}{:soso_e179:}

txyw

bootstrap和蒙特卡洛间的区别大不大？

wjyalx

这个的确不错啊，把统计学知识和机器学习结合起来，值得学习哈。

suppppperrr

您好~~~这个用于对比的truemean应该是sample的mean还是总体population的mean呀？用不用先rnorm一个总体，在从总体里抽样出来作为样本（sample）?正在做一项相关的作业，这里感到很困惑。