2d. Temukan mean dan varians dari distribusi sampling t ̂_yr. Bagaimana nilai-nilai ini dibandingkan dengan mean dan varians dari Nȳ? Berapa bias dari t ̂_yr?

# Data
x <- c(13, 7, 11, 12, 4, 3, 11, 3, 5)
y <- c(10, 7, 13, 17, 8, 1, 15, 7, 4)

# Total t_x dan N (jumlah unit)
tx <- sum(x)
N <- length(y)

# Paket yang diperlukan
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
# Fungsi untuk menghitung B, t_yr, dan N * y_bar
calculate_B_tyr_Nybar <- function(sample_indices) {
  x_sample <- x[sample_indices]
  y_sample <- y[sample_indices]
  
  x_bar_sample <- mean(x_sample)
  y_bar_sample <- mean(y_sample)
  
  # Menghitung B
  B <- y_bar_sample / x_bar_sample
  
  # Menghitung t_hat_yr
  t_hat_yr <- B * tx
  
  # Menghitung N * y_bar
  N_y_bar <- N * y_bar_sample
  
  return(c(x_bar_sample, y_bar_sample, B, t_hat_yr, N_y_bar))
}

# Generate all combinations of samples of size 3
library(combinat)
## 
## Attaching package: 'combinat'
## The following object is masked from 'package:utils':
## 
##     combn
samples <- combn(1:N, 3)

# Tabel distribusi
results <- data.frame(Sample = character(), 
                      x_bar = numeric(), 
                      y_bar = numeric(), 
                      B = numeric(), 
                      t_hat_yr = numeric(), 
                      N_y_bar = numeric())

# Mengisi tabel dengan perhitungan dari setiap sampel
for (i in 1:ncol(samples)) {
  sample_indices <- samples[, i]
  calculation <- calculate_B_tyr_Nybar(sample_indices)
  
  results <- rbind(results, data.frame(Sample = paste(sample_indices, collapse = ","),
                                       x_bar = calculation[1],
                                       y_bar = calculation[2],
                                       B = calculation[3],
                                       t_hat_yr = calculation[4],
                                       N_y_bar = calculation[5]))
}

# Hitung mean dan varians dari distribusi t_hat_yr
mean_t_hat_yr <- mean(results$t_hat_yr)
var_t_hat_yr <- var(results$t_hat_yr)

# Hitung mean dan varians dari distribusi N_y_bar
mean_N_y_bar <- mean(results$N_y_bar)
var_N_y_bar <- var(results$N_y_bar)

# Tampilkan hasilnya
cat("Mean dari t_hat_yr: ", mean_t_hat_yr, "\n")
## Mean dari t_hat_yr:  82.79241
cat("Varians dari t_hat_yr: ", var_t_hat_yr, "\n")
## Varians dari t_hat_yr:  196.9583
cat("Mean dari N_y_bar: ", mean_N_y_bar, "\n")
## Mean dari N_y_bar:  82
cat("Varians dari N_y_bar: ", var_N_y_bar, "\n")
## Varians dari N_y_bar:  489.3253

bias tyr

# Total populasi y (t_y)
t_y <- sum(y)

# Hitung bias
bias_t_hat_yr <- mean_t_hat_yr - t_y

# Tampilkan bias
cat("Bias dari t_hat_yr: ", bias_t_hat_yr, "\n")
## Bias dari t_hat_yr:  0.7924067