######***###### UNIVARIATE DATA ######***######
### Example 1: univariate statistic (median)
# generate 100 standard normal observations
set.seed(1)
n <- 100
x <- rnorm(n)
# nonparametric bootstrap
npbs <- np.boot(x = x, statistic = median)
npbs
### Example 2: multivariate statistic (quartiles)
# generate 100 standard normal observations
set.seed(1)
n <- 100
x <- rnorm(n)
# nonparametric bootstrap
npbs <- np.boot(x = x, statistic = quantile,
probs = c(0.25, 0.5, 0.75))
npbs
######***###### MULTIVARIATE DATA ######***######
### Example 1: univariate statistic (correlation)
# correlation matrix square root (with rho = 0.5)
rho <- 0.5
val <- c(sqrt(1 + rho), sqrt(1 - rho))
corsqrt <- matrix(c(val[1], -val[2], val), 2, 2) / sqrt(2)
# generate 100 bivariate observations (with rho = 0.5)
n <- 100
set.seed(1)
data <- cbind(rnorm(n), rnorm(n)) %*% corsqrt
# define statistic function
statfun <- function(x, data) cor(data[x,1], data[x,2])
# nonparametric bootstrap
npbs <- np.boot(x = 1:n, statistic = statfun, data = data)
npbs
### Example 2: multivariate statistic (variances and covariance)
# correlation matrix square root (with rho = 0.5)
rho <- 0.5
val <- c(sqrt(1 + rho), sqrt(1 - rho))
corsqrt <- matrix(c(val[1], -val[2], val), 2, 2) / sqrt(2)
# generate 100 bivariate observations (with rho = 0.5)
n <- 100
set.seed(1)
data <- cbind(rnorm(n), rnorm(n)) %*% corsqrt
# define statistic function
statfun <- function(x, data) {
cmat <- cov(data[x,])
ltri <- lower.tri(cmat, diag = TRUE)
cvec <- cmat[ltri]
names(cvec) <- c("var(x1)", "cov(x1,x2)", "var(x2)")
cvec
}
# nonparametric bootstrap
npbs <- np.boot(x = 1:n, statistic = statfun, data = data)
npbs
if (FALSE) {
######***###### REGRESSION ######***######
### Example 1: bootstrap cases
# generate 100 observations
n <- 100
set.seed(1)
x <- seq(0, 1, length.out = n)
y <- 1 + 2 * x + rnorm(n)
data <- data.frame(x = x, y = y)
# define statistic function
statfun <- function(x, data) {
xmat <- cbind(1, data$x[x])
xinv <- solve(crossprod(xmat)) %*% t(xmat)
coef <- as.numeric(xinv %*% data$y[x])
names(coef) <- c("(Intercept)", "x")
coef
}
# nonparametric bootstrap
npbs <- np.boot(x = 1:n, statistic = statfun, data = data)
npbs
### Example 2: bootstrap residuals
# generate 100 observations
n <- 100
set.seed(1)
x <- seq(0, 1, length.out = n)
y <- 1 + 2 * x + rnorm(n)
# prepare data
xmat <- cbind(1, x)
xinv <- solve(crossprod(xmat)) %*% t(xmat)
fit <- xmat %*% xinv %*% y
data <- list(fit = fit, resid = y - fit, xinv = xinv, x = x)
# define statistic function
statfun <- function(x, data) {
ynew <- data$fit + data$resid[x]
coef <- as.numeric(data$xinv %*% ynew)
names(coef) <- c("(Intercept)", "x")
coef
}
# define jackknife function
jackfun <- function(x, data){
ynew <- data$fit[x] + data$resid[x]
xmat <- cbind(1, data$x[x])
xinv <- solve(crossprod(xmat)) %*% t(xmat)
coef <- as.numeric(xinv %*% ynew)
names(coef) <- c("(Intercept)", "x")
coef
}
# nonparametric bootstrap
npbs <- np.boot(x = 1:n, statistic = statfun, data = data,
jackknife = jackfun)
npbs
}
Run the code above in your browser using DataLab