library(boot)
y <- c(2,2,1,4,1,0,5,3,1,6,0,0,3,1,3,0,3,0,2,20,0,2,3,1,25)
y
n <- length(y)
n
thetahat = mean(y)
thetahat
Brep = 10000

# Bootstrap the sample mean
sampmean <- function(y,i) mean(y[i])

bootmean <- boot(data=y,statistic=sampmean,R=Brep)
bootmean
boot.ci(bootmean,conf=.95,type=c("norm"))
boot.ci(bootmean,conf=.95,type=c("perc"))
boot.ci(bootmean,conf=.95,type=c("bca"))
par(mfrow=c(2,1))

hist(bootmean$t,main="Bootstrap Sample Means")
plot(ecdf(bootmean$t),main="Empirical CDF of Bootstrap Means")