# Efron and Tibshirani, 1994, An Introduction to the Bootstrap, chapter 11 # ======================================== # jackknife() # compute theta.hat for each jackknife data set jackknife <- function(x, theta, xdata, ...){ n <- length(x) u <- numeric(n) for(i in 1:n) { u[i] <- theta(x[ - i], xdata) } thetahat <- theta(x, xdata) jack.bias <- (n - 1) * (mean(u) - thetahat) jack.se <- sqrt(((n - 1)/n) * sum((u - mean(u))^2)) return(c(u, jack.bias=jack.bias, jack.se=jack.se)) } # ======================================== # test score data score <- read.table("scor.data") names(score) <- c("mechanics", "vectors", "algebra", "analysis", "statistics") head(score) str(score) summary(score) xdata <- as.matrix(score) n <- nrow(xdata) # theta calculates ratio of largest eigenvalue to sum of eigenvalues of covariance matrix library(panel) theta <- function(x, xdata){ cov.star <- cov(xdata[x, ]) eigensys.star <- eddcmp(cov.star) lambdas.star <- eigensys.star\$evalues return(lambdas.star[1] / sum(lambdas.star)) } theta.hat <- theta(1:n, xdata) # => 0.619115 # jackknife values # compare with figure 11.1, p.144 (results <- jackknife(1:n, theta, xdata)) results.theta.hat <- results[1:n] results.bias <- results[n+1] # => 0.001069139 results.se <- results[n+2] # => 0.04955231 hist(results.theta.hat, col="springgreen2", xlim=c(.45, .75), xlab="jackknife") theta.hat.mean <- mean(results.theta.hat) # => 0.6191273 scaled.results <- mu + sqrt(87) * (results.theta.hat - theta.hat.mean) hist(scaled.results, col="thistle4", xlim=c(.45, .75), xlab="inflated jackknife") # bootstrap values of theta.hat bootstrap <- function(x, nboot, theta, ...){ data <- matrix(sample(x, size=length(x) * nboot, replace=TRUE), nrow=nboot) return(apply(data, 1, theta, xdata)) } (bootstrap.results <- bootstrap(1:n, nboot=88, theta, xdata)) hist(bootstrap.results, col="palevioletred2", xlim=c(.45, .75), xlab="bootstrap") # ========================================