library(MASS) data(Boston) y <- Boston[,14] x <- Boston[,-14] n <- nrow(Boston) g <- data.frame(crim=exp(seq(log(min(x[,1])),log(max(x[,1])),length=500))) p <- matrix(0,1,500) p2 <- matrix(0,1,500) B <- 100 plot(log(x$crim),y) for (b in 1:B) { bss <- sample(1:n, n , replace=TRUE) rp <- rpart(y[bss] ~ log(crim), data=x[bss,], control=rpart.control(maxdepth=3)) pred <- predict(rp,g) p <- p + pred p2 <- p2 + pred^2 lines(log(g$crim),pred,col=b,lwd=.5) } p <- p/B sd <- sqrt((p2/B - p^2)*B/(B-1)) lines(log(g$crim),p,col=1,lwd=10) lines(log(g$crim),p+sd,col=1,lwd=5) lines(log(g$crim),p-sd,col=1,lwd=5)