# R code for bootstrap resampling to approximate 
# sampling distribution of X-bar and s^2

# Defining the sample data:

my.data <- c(7, 9, 13, 12, 4, 6, 8, 10, 10, 7)
my.n <- length(my.data)

#Defining the number of resamples:

my.m <- 1000

# Setting up the matrix to hold bootstrap-sample values

setup.data.matrix <- matrix(my.data, nrow=my.m, ncol=my.n, byrow=T)

# carrying out the sampling (with replacement):

bootstrap.data.matrix <- apply(setup.data.matrix, 1, sample, size=my.n, replace=TRUE)

# Transposing to get back to same dimensions as setup.data.matrix

bootstrap.data.matrix <- t(bootstrap.data.matrix)

# Calculating the sample mean for each of the bootstrap samples

my.means <- apply(bootstrap.data.matrix, 1, mean)

# Calculating the sample variance for each of the bootstrap samples

my.vars <- apply(bootstrap.data.matrix, 1, var)


par(mfrow=c(2,2))
hist(my.means)
plot(density(my.means))

hist(my.vars)
plot(density(my.vars))

# sample estimates for mu and sigma^2:

print(mean(my.data)); print(var(my.data))

# Bootstrap point estimates for mu and sigma^2:

print(median(my.means)); print(median(my.vars))

# Bootstrap interval estimates for mu and sigma^2:

lower.upper.mean <- quantile(my.means, probs=c(0.025, 0.975))

print(paste("95% bootstrap interval for mu: (", round(lower.upper.mean[1],3), ", ", round(lower.upper.mean[2],3), ")"))

lower.upper.variance <- quantile(my.vars, probs=c(0.025, 0.975))

print(paste("95% bootstrap interval for variance: (", round(lower.upper.variance[1],3), ", ", round(lower.upper.variance[2],3), ")"))