-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBOOTSTRAP JACKKNIFE.R
52 lines (48 loc) · 1.52 KB
/
BOOTSTRAP JACKKNIFE.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
set.seed(123)
mean_x <- 0
mean_y <- 0
var_x <- 1
var_y <- 1
correlation <- 0.3
# Generate bivariate normal data
n <- 10
data <- MASS::mvrnorm(n, mu = c(mean_x, mean_y), Sigma = matrix(c(var_x, correlation, correlation, var_y), ncol = 2))
# Function to calculate the sample statistic of interest
calculate_statistic <- function(data) {
return(var(data[, 1]))
}
# Bootstrap function
bootstrap <- function(data, B) {
statistics <- numeric(B)
for (i in 1:B) {
resampled_data <- data[sample(1:n, replace = TRUE), ]
statistics[i] <- calculate_statistic(resampled_data)
}
return(statistics)
}
# Number of bootstrap samples
B <- 1000
bootstrap_results <- bootstrap(data, B)
original_statistic <- calculate_statistic(data)
bias <- mean(bootstrap_results) - original_statistic
variance <- var(bootstrap_results)
median_estimate <- median(bootstrap_results)
#######################
# Function to calculate the sample statistic of interest
calculate_statistic <- function(data) {
return(var(data[, 1]))
}
# Jackknife function
jackknife <- function(data) {
n <- nrow(data)
statistics <- numeric(n)
for (i in 1:n) {
subset_data <- data[-i, ]
statistics[i] <- calculate_statistic(subset_data)
}
return(statistics)
}
jackknife_results <- jackknife(data)
original_statistic <- calculate_statistic(data)
bias_jackknife <- (n - 1) * (mean(jackknife_results) - original_statistic)
variance_jackknife <- ((n - 1) / n) * sum((jackknife_results - mean(jackknife_results))^2)