Homework 1 (Boostrapping Methods)

Tải Thư Viện

Code

library(tidyverse)
library(dplyr)
library(MASS)
library(boot)
library(ggplot2)
library(gridExtra)
# suppressPackageStartupMessages(library(tidyverse)) # Suppress warnings from Tidyverse

Các Hàm Phụ Trợ

Code

mytheme = list(
  theme_classic()+
    theme(panel.background = element_blank(),
          strip.background = element_rect(colour=NA, fill=NA),
          panel.border = element_rect(fill = NA, color = "black"),
          legend.title = element_blank(),legend.position="bottom", 
          strip.text = element_text(face="bold", size=9),
          axis.text=element_text(face="bold"),
          axis.title = element_text(face="bold"),
          plot.title = element_text(face = "bold", hjust = 0.5,size=13)))

Bài Tập 5

Cho biết \(X\) là một biến ngẫu nhiên tuân theo phân phối chuẩn với trung bình bằng \(36\) và độ lệch chuẩn bằng \(8\). Ta có phân phối mẫu cho trung bình mẫu với cỡ mẫu \(n\) như sau:

\[ \bar{X} \sim \mathcal{N}(\mu, \sigma^2/n) \] với \(\mu_{\bar{X}} = \mu\) và \(\sigma_{\bar{X}} = \sigma\). Như vậy, với \(n = 200\), phân phối mẫu cho \(\bar{X}\) là \(\mathcal{N}(36, 8^2 / \sqrt{200})\).

Code

bootstrap_mean <- function(data, n_bootstrap = 1000) {
  n <- length(data)
  bootstrap_means <- numeric(n_bootstrap)
  
  for (i in 1:n_bootstrap) {
    bootstrap_sample <- sample(data, n, replace = TRUE)
    bootstrap_means[i] <- mean(bootstrap_sample)
  }
  
  return(bootstrap_means)
}

# Theoretical mean and standard deviation of the sample mean
mu <- 36
sigma <- 8
theoretical_mean <- mu
theoretical_sd <- sigma / sqrt(200)

sample_200 <- rnorm(200, mean = mu, sd = sigma)

# Boostrapping mean and standard deviation of the sample mean
bootstrap_means_200 <- bootstrap_mean(sample_200)
boostrap_mean <- mean(bootstrap_means_200)
boostrap_sd <- sd(bootstrap_means_200)

print(paste("Theoretical Mean:", theoretical_mean))

[1] "Theoretical Mean: 36"

Code

print(paste("Theoretical SD:", round(theoretical_sd, 3)))

[1] "Theoretical SD: 0.566"

Code

print(paste("Sample Mean:", round(mean(sample_200), 3)))

[1] "Sample Mean: 34.736"

Code

print(paste("Sample SD:", round(sd(sample_200), 3)))

[1] "Sample SD: 7.924"

Code

print(paste("Bootstrap Mean:", round(boostrap_mean, 3)))

[1] "Bootstrap Mean: 34.714"

Code

print(paste("Bootstrap SD:", round(boostrap_sd, 3)))

[1] "Bootstrap SD: 0.547"

Code

# List of sample sizes
sample_sizes <- c(200, 50, 10)
plots <- list()
sampling_dist <- list()
# For each sample size, generate data, perform bootstrap, and create plot
for (n in sample_sizes) {
  set.seed(123)
  sample_data <- rnorm(n, mean = 36, sd = 8)
  
  sample_mean <- mean(sample_data)
  sample_sd <- sd(sample_data)
  
  bootstrap_means <- bootstrap_mean(sample_data)
  
  bootstrap_df <- data.frame(bootstrap_means = bootstrap_means)
  
  plot <- ggplot(bootstrap_df, aes(x = bootstrap_means)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "lightblue", color = "blue") +
    stat_function(fun = dnorm, args = list(mean = 36, sd = 8 / sqrt(n)), color = "red", linewidth = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = mean(bootstrap_means)), color = "green", linewidth = 1) +
    labs(title = paste("n=", n, sep=""), x = "Sample Mean", y = "Density") +
    mytheme +
    scale_color_manual(values = c("Theoretical Density" = "red", "Bootstrap Mean" = "green")) +
    theme(legend.position = "top") +
    guides(color = guide_legend(override.aes = list(size = 1.5))) +
    labs(color = "")
  
  plots[[as.character(n)]] <- plot
}

# Arrange plots
grid.arrange(plots[[1]], plots[[2]], plots[[3]], ncol = 3)

Ta thấy phân phối mẫu của trung bình mẫu được khi dùng boostrap gần như tương đồng với phân phối mẫu lý thuyết với sai biệt trị tuyệt đối chỉ khoảng \(1\) cho trung bình và \(0.02\) cho độ lệch chuẩn. Tuy nhiên khi dữ liệu của ta không đủ lớn (trong trường hợp \(n=50\) và \(n=10\)), phân phối mẫu do boostrap thiết lập cho sự biến thiên lớn và tính chính xác của xấp xỉ này là không cao. Điều này có thể lý giải thông qua cách xây dựng của boostrap phụ thuộc và số lượng dữ liệu. Khi dữ liệu càng lớn thì xấp xỉ của hàm phân phối thực nghiệm càng tiến về hàm phân phối thực sự của dữ liệu (chiếu theo góc nhìn từ định lý Glivenko Cantelli).

Bài Tập 6

Giả sử ta không có dữ liệu thực tế, ta sẽ tạo dữ liệu giả bằng đoạn mã sau:

Code

# Simulate data 
set.seed(123) 
scores <- floor(rnorm(50, mean = 60, sd = 3))
group <- sample(c(0, 1), 50, replace = TRUE)  # 0: not infected, 1: infected

df <- data.frame(scores, group)

Các điểm số sẽ được tạo ngẫu nhiên từ phân phối chuẩn với trung bình \(60\) và độ lệch chuẩn \(3\). Vì điểm số chỉ có thể nhận giá trị nguyên dương nên ta sẽ lấy phân nguyên nhỏ hơn của các giá trị khởi tạo. Các nhóm được lấy ngẫu nhiên theo phân phối Bernoulli với xác suất là bằng nhau. Trong trường hợp biến nhị phân thì xác xuât \(p(X = 0) = p(X = 1) = 1/2\) và ta có thể thực hiện với lệnh sample. Như vậy, ta đã có dữ liệu giả để thực hiện trong phần bài tập này.

Code

ggplot(df, aes(x = factor(group), y = scores)) +
  geom_boxplot() +
  labs(x = "Group", y = "Score", title = "Compare scores between groups") + mytheme

Ta thấy rằng trung vị phần trăm tế bào CD3-positive lymphocytes của nhóm bị nhiễm vi khuẩn cao hơn nhóm không nhiễm. Điều này ủng hộ cho giả thuyết khả năng việc sử dung phân biệt giữa nhóm bị nhiễm và nhóm không bị nhiễm.

Với \(X, Y\) lần lượt là giá trị của điểm trong nhóm \(0\) và \(1\), ta cần ước lượng cho \(\theta = P(X < Y)\) như là một thước đo độ chính xác trong chuẩn đoán bệnh nhân brucellois bằng CD3. Ta có ước lượng này được tính như sau \[ \hat{\theta} = \frac{1}{nm} \sum_{i=1}^{n} \sum_{j=1}^{n} \mathcal{I}_{\left(X_i < Y_j\right)} \] với \(\mathcal{I}(.)\) là hàm chỉ mục. \(\hat{\theta}\) có thể được tính như sau trong R:

Code

X <- df$scores[df$group == 0]
Y <- df$scores[df$group == 1]

theta_hat <- mean(outer(X, Y, "<"))

estimate_theta <- function(df) {
  X <- df$scores[df$group == 0]
  Y <- df$scores[df$group == 1]
  
  n <- length(X)
  m <- length(Y)
  
  count <- 0
  for (i in 1:n) {
    for (j in 1:m) {
      if (X[i] < Y[j]) {
        count <- count + 1
      }
    }
  }
  
  theta_hat <- count / (n * m)
  return(theta_hat)
}

theta_hat_loop <- estimate_theta(df)
theta_hat_loop

[1] 0.6095076

Code

abs(theta_hat - theta_hat_loop)

[1] 0

Ta có thể tính nhanh hai tổng bằng cách tạo ma trận với số lượng \(X\) là số hàng và số lượng \(Y\) là số cột sau đó thực hiện so sánh \(X_i\) và \(Y_j\). Như vậy ma trận của chúng ta chỉ gồm các giá trị nhị nhân TRUE/FALSE. Khi này, ta chỉ cần đếm số lượng TRUE và chia cho \(nm\). So sánh với cách tiếp cận dùng vòng lặp, ta thấy kết quả chó ra đều như nhau.

Ta thực hiện boostrap để xác định phân phối mẫu của \(\hat{\theta}\) như sau:

1. Với số lần boostrap \(B\), tạo một mẫu ngẫu nhiên với kích cỡ là \(n\) với \(n\) là số lượng mẫu trong dữ liệu gốc có lặp lại.
2. Trích xuất \(X\) và \(Y\) tương ứng với \(X\) là nhóm người không bị nhiễm và \(Y\) là nhóm người bị nhiễm.
3. Tính \(\hat{\theta} = \frac{1}{nm} \sum_{i=1}^{n} \sum_{j=1}^{n} \mathcal{I}_{\left(X_i < Y_j\right)}\).
4. Lặp lại bước 2 và 3 \(B\) lần.

Đoạn mã R sau đây mô phỏng quá trình boostrap trên cũng như phân phối mẫu cho \(\hat{\theta}\).

Code

# Boostrap 
B <- 10000
theta_bootstrap <- numeric(B)
for (i in 1:B) {
  idx <- sample(nrow(df), replace = TRUE)
  X_boot <- df$scores[idx][df$group[idx] == 0]
  Y_boot <- df$scores[idx][df$group[idx] == 1]
  theta_bootstrap[i] <- mean(outer(X_boot, Y_boot, "<"))
}

# Mean Boostrap
mean_bootstrap <- mean(theta_bootstrap)

bootstrap_df <- data.frame(theta = theta_bootstrap)
ci <- quantile(theta_bootstrap, c(0.025, 0.975))

# Plotting 
ggplot(bootstrap_df, aes(x = theta)) +
  geom_histogram(binwidth = 0.01, fill = "lightblue", color = "white") +
  geom_vline(xintercept = mean_bootstrap, color = "red", linetype = "dashed") +
  labs(x = expression(theta), y = "Frequency", 
       title = "Boostrap distribution") +
  annotate("text", x = mean_bootstrap, y = Inf,
           label = paste("Bootstrap mean:", 
                         round(mean_bootstrap, 3)), 
           vjust = 1, hjust = -0.1, color = "red") + mytheme