1 Latar Belakang & Tujuan

1.1 Ringkasan

Bootstrap adalah teknik resampling non-parametrik yang diperkenalkan oleh Efron (1979). Metode ini menjawab pertanyaan mendasar dalam statistika inferensial:

“Bagaimana kita dapat mengukur ketidakpastian sebuah estimasi statistik ketika kita hanya memiliki satu sampel dari populasi?”

Pada minggu ini, kita mempelajari bootstrap melalui dua latihan:

Latihan Tugas
1 Bangkitkan data normal (\(n=1000\), \(\mu=30\), \(\sigma=2.5\)), lalu hitung rata-rata dari 50 sampel bootstrap
2 Visualisasikan distribusi rata-rata bootstrap vs distribusi data asli dalam satu panel gabungan yang mewah

1.2 Landasan Teoritis

Prinsip Bootstrap dapat diringkas sebagai berikut:

\[\hat{\theta}^* = T(X^*_1, X^*_2, \ldots, X^*_B)\]

di mana \(X^*_b\) adalah sampel bootstrap ke-\(b\) (diambil dengan pengembalian dari sampel asli \(X\)), dan \(\hat{\theta}^*\) adalah distribusi estimator yang dihasilkan.

Standard Error Bootstrap (estimasi ketidakpastian):

\[\widehat{\text{SE}}_B(\hat{\theta}) = \sqrt{\frac{1}{B-1} \sum_{b=1}^{B} \left(\hat{\theta}^*_b - \bar{\theta}^*\right)^2}\]

Perbandingan dengan SE Analitik untuk mean:

\[\text{SE}_{\text{analitik}}(\bar{X}) = \frac{\sigma}{\sqrt{n}} = \frac{2.5}{\sqrt{1000}} \approx 0.0791\]

Jika bootstrap berjalan dengan benar, \(\widehat{\text{SE}}_B\) akan mendekati nilai analitik ini.

Teorema Limit Pusat (CLT) memastikan bahwa distribusi bootstrap dari \(\bar{X}^*\) akan mendekati distribusi normal:

\[\bar{X}^* \xrightarrow{d} \mathcal{N}\!\left(\bar{X},\, \frac{s^2}{n}\right) \quad \text{ketika } B \to \infty\]

2 Persiapan: Memuat Paket

Bagian ini memuat semua paket yang diperlukan. ggplot2 untuk visualisasi berkualitas tinggi, dplyr & tidyr untuk manipulasi data, viridis & scales untuk palet warna yang aksesibel, patchwork untuk menyusun multi-panel plot, dan knitr untuk tabel yang rapi. Paket-paket ini adalah fondasi ekosistem tidyverse yang menjadi standar industri analisis data modern.

library(ggplot2)
library(dplyr)
library(tidyr)
library(viridis)
library(patchwork)
library(scales)
library(knitr)

3 Bagian 1 — Pembangkitan Data Asli (myData)

Sebelum melakukan bootstrap, kita membutuhkan data populasi awal. Di sini dibangkitkan 1000 observasi dari distribusi normal dengan \(\mu = 30\) dan \(\sigma = 2.5\) menggunakan rnorm(). Nilai set.seed(150) memastikan reproduktibilitas penuh — siapapun yang menjalankan kode ini akan mendapatkan hasil yang identik.

Penting untuk dipahami bahwa rnorm() menghasilkan nilai acak di sekitar parameter yang diberikan — nilai mean dan SD sampel tidak akan persis sama dengan \(\mu\) dan \(\sigma\) populasi, melainkan mendekatinya seiring bertambahnya \(n\) (Hukum Bilangan Besar).

# ── Parameter ────────────────────────────────────────────────
set.seed(150)          # Seed untuk reproduktibilitas
N_OBS  <- 1000        # Jumlah observasi
MU     <- 30          # Rata-rata populasi
SIGMA  <- 2.5         # Standar deviasi populasi

# ── Bangkitkan data ──────────────────────────────────────────
myData <- rnorm(N_OBS, mean = MU, sd = SIGMA)

3.1 Sanity Check Data Asli

Tiga pemeriksaan awal (sanity checks) dilakukan untuk memvalidasi data yang baru dibangkitkan. Kita berharap: jumlah observasi = 1000, mean ≈ 30, dan SD ≈ 2.5. Penyimpangan kecil adalah hal normal dan mencerminkan variasi sampling alami.

# ── Sanity checks ────────────────────────────────────────────
cat("Jumlah observasi :", length(myData), "\n")
## Jumlah observasi : 1000
cat("Mean sampel      :", round(mean(myData), 5), "\n")
## Mean sampel      : 29.92068
cat("Std. Dev. sampel :", round(sd(myData),   5), "\n")
## Std. Dev. sampel : 2.47517
cat("Min              :", round(min(myData),   4), "\n")
## Min              : 22.2416
cat("Max              :", round(max(myData),   4), "\n")
## Max              : 36.7845
# ── Tabel ringkasan ──────────────────────────────────────────
data.frame(
  Statistik = c("Jumlah Observasi", "Mean", "Std. Deviasi",
                "Min", "Q1", "Median", "Q3", "Max"),
  `Nilai Sampel` = round(c(
    length(myData), mean(myData), sd(myData),
    quantile(myData, 0)[[1]], quantile(myData, 0.25)[[1]],
    median(myData), quantile(myData, 0.75)[[1]],
    quantile(myData, 1)[[1]]
  ), 4),
  `Nilai Teoritis` = c(1000, 30, 2.5, NA, NA, 30, NA, NA),
  check.names = FALSE
) %>%
  kable(align = "c",
        caption = "Ringkasan Statistik myData vs Parameter Populasi Teoritis")
Ringkasan Statistik myData vs Parameter Populasi Teoritis
Statistik Nilai Sampel Nilai Teoritis
Jumlah Observasi 1000.0000 1000.0
Mean 29.9207 30.0
Std. Deviasi 2.4752 2.5
Min 22.2416 NA
Q1 28.2559 NA
Median 29.8430 30.0
Q3 31.6312 NA
Max 36.7845 NA

3.2 Visualisasi Data Asli

Histogram berikut memvisualisasikan distribusi myData. Tiga elemen kunci ditampilkan: (1) kurva densitas empiris (biru), (2) kurva distribusi normal teoritis \(\mathcal{N}(30, 2.5^2)\) (oranye putus-putus) sebagai pembanding, dan (3) garis vertikal mean sampel. Kedekatan antara kurva empiris dan teoritis mengkonfirmasi data terbangkit dengan benar.

ggplot(data.frame(x = myData), aes(x = x)) +
  geom_histogram(
    aes(y = ..density..),
    bins  = 40,
    fill  = "#3949ab", color = "white", alpha = 0.75
  ) +
  geom_density(color = "#1565c0", size = 1.2, linetype = "solid") +
  stat_function(
    fun  = dnorm, args = list(mean = MU, sd = SIGMA),
    color = "#e65100", size = 1.1, linetype = "dashed"
  ) +
  geom_vline(
    xintercept = mean(myData),
    color = "#c62828", size = 1.2, linetype = "solid"
  ) +
  annotate(
    "text",
    x     = mean(myData) + 0.25,
    y     = Inf, vjust = 2, hjust = 0,
    label = paste0("Mean = ", round(mean(myData), 4)),
    color = "#c62828", fontface = "bold", size = 4
  ) +
  scale_x_continuous(breaks = pretty_breaks(8)) +
  labs(
    title    = "Distribusi Data Asli — myData",
    subtitle = paste0(
      "N = ", N_OBS, "  |  μ = ", MU, "  |  σ = ", SIGMA,
      "  |  Garis merah = Mean sampel  |  Kurva oranye = N(30, 2.5²) teoritis"
    ),
    x = "Nilai Observasi",
    y = "Densitas"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1a237e", size = 14),
    plot.subtitle = element_text(color = "#5c6bc0", size = 9.5)
  )

4 Bagian 2 — Latihan 1: Simulasi Bootstrap 50 Sampel

Latihan 1 meminta kita melakukan bootstrap dengan 50 ulangan (\(B = 50\)), di mana setiap ulangan mengambil sampel berukuran 1000 dari myData dengan pengembalian (replace = TRUE), lalu menghitung mean-nya. Hasilnya disimpan dalam vektor bootstrap.results.

Kunci konseptual: pengambilan dengan pengembalian berarti satu observasi bisa terpilih berkali-kali dalam satu sampel bootstrap, yang menciptakan variasi antar sampel meski semua berasal dari data yang sama. Inilah yang memungkinkan kita mensimulasikan variabilitas sampling tanpa perlu data baru.

4.1 Menjalankan Bootstrap

# ── Parameter bootstrap ──────────────────────────────────────
set.seed(150)                      # Seed identik dengan soal
sample.size       <- N_OBS         # Ukuran tiap sampel bootstrap = 1000
n.samples         <- 50            # Jumlah ulangan bootstrap (B = 50)
bootstrap.results <- numeric(n.samples)  # Vektor penyimpan hasil

# ── Loop bootstrap ───────────────────────────────────────────
for (i in 1:n.samples) {
  # Ambil indeks acak DENGAN PENGEMBALIAN
  obs <- sample(seq_len(sample.size), size = sample.size, replace = TRUE)

  # Hitung mean dari sub-sampel bootstrap ke-i
  bootstrap.results[i] <- mean(myData[obs])
}

4.2 Sanity Check & Ringkasan Hasil Bootstrap

Pemeriksaan panjang vektor memastikan 50 nilai mean tersimpan. Standar deviasi dari bootstrap.results adalah Standard Error Bootstrap — estimasi empiris dari ketidakpastian mean. Nilai ini kemudian kita bandingkan dengan SE analitik \(\sigma/\sqrt{n}\).

# ── Sanity check ─────────────────────────────────────────────
cat("Jumlah sampel bootstrap  :", length(bootstrap.results), "\n")
## Jumlah sampel bootstrap  : 50
cat("Mean dari bootstrap means:", round(mean(bootstrap.results), 5), "\n")
## Mean dari bootstrap means: 29.92612
cat("SE Bootstrap (SD means)  :", round(sd(bootstrap.results),   5), "\n")
## SE Bootstrap (SD means)  : 0.07028
cat("SE Analitik (σ/√n)       :", round(SIGMA / sqrt(N_OBS),     5), "\n")
## SE Analitik (σ/√n)       : 0.07906
# ── Tabel perbandingan ────────────────────────────────────────
data.frame(
  Statistik = c("B (jumlah ulangan)", "Mean bootstrap",
                "SD bootstrap (≈ SE)", "SE analitik σ/√n",
                "Selisih absolut SE", "Min mean", "Max mean", "Rentang mean"),
  Nilai = round(c(
    length(bootstrap.results),
    mean(bootstrap.results),
    sd(bootstrap.results),
    SIGMA / sqrt(N_OBS),
    abs(sd(bootstrap.results) - SIGMA / sqrt(N_OBS)),
    min(bootstrap.results),
    max(bootstrap.results),
    max(bootstrap.results) - min(bootstrap.results)
  ), 5),
  check.names = FALSE
) %>%
  kable(align = "c",
        caption = paste0(
          "Ringkasan Distribusi Bootstrap Mean (B = ", n.samples,
          ", n = ", sample.size, ")"
        ))
Ringkasan Distribusi Bootstrap Mean (B = 50, n = 1000)
Statistik Nilai
B (jumlah ulangan) 50.00000
Mean bootstrap 29.92612
SD bootstrap (≈ SE) 0.07028
SE analitik σ/√n 0.07906
Selisih absolut SE 0.00878
Min mean 29.73110
Max mean 30.07557
Rentang mean 0.34447

4.3 Seluruh Nilai Mean Bootstrap

Tabel berikut menampilkan seluruh 50 nilai mean yang diperoleh dari 50 sampel bootstrap. Fluktuasi kecil antar nilai mencerminkan variasi sampling alami — inilah inti dari apa yang diestimasi oleh bootstrap.

data.frame(
  `Sampel ke-` = 1:n.samples,
  `Mean Bootstrap` = round(bootstrap.results, 5),
  `Deviasi dari Grand Mean` = round(bootstrap.results - mean(bootstrap.results), 5),
  check.names = FALSE
) %>%
  kable(align = "c",
        caption = "Seluruh 50 Nilai Mean Bootstrap")
Seluruh 50 Nilai Mean Bootstrap
Sampel ke- Mean Bootstrap Deviasi dari Grand Mean
1 29.88877 -0.03735
2 29.95125 0.02514
3 29.91789 -0.00823
4 29.91696 -0.00916
5 29.90099 -0.02513
6 29.83304 -0.09308
7 29.85498 -0.07114
8 29.95973 0.03361
9 29.86083 -0.06529
10 30.03870 0.11258
11 29.84256 -0.08356
12 29.91217 -0.01395
13 29.87484 -0.05128
14 29.90375 -0.02237
15 29.94881 0.02269
16 29.87566 -0.05045
17 29.89902 -0.02710
18 29.87050 -0.05562
19 29.91005 -0.01607
20 29.88778 -0.03834
21 30.00650 0.08038
22 29.92184 -0.00428
23 29.98550 0.05939
24 29.73110 -0.19502
25 29.96483 0.03871
26 29.91879 -0.00733
27 29.85054 -0.07558
28 30.02939 0.10327
29 29.95772 0.03160
30 30.06538 0.13926
31 29.98891 0.06280
32 29.95251 0.02639
33 29.79974 -0.12638
34 29.93682 0.01071
35 29.91790 -0.00822
36 29.98175 0.05563
37 29.98773 0.06161
38 29.91802 -0.00809
39 29.83295 -0.09316
40 29.99206 0.06594
41 29.91268 -0.01344
42 30.00392 0.07780
43 29.90862 -0.01750
44 30.00265 0.07653
45 30.07557 0.14945
46 29.92908 0.00296
47 29.88484 -0.04128
48 30.02723 0.10112
49 29.96257 0.03646
50 29.81047 -0.11564

5 Bagian 3 — Latihan 2: Visualisasi Gabungan

Latihan 2 meminta dua histogram: satu untuk distribusi mean bootstrap, satu untuk data asli, digabung dalam satu panel. Di sini kita melampaui permintaan minimal dengan membangun visualisasi multi-layer yang kaya informasi menggunakan ggplot2 dan patchwork.

5.1 Plot A — Distribusi Mean Bootstrap

Histogram ini menampilkan distribusi dari 50 mean bootstrap. Sesuai CLT, distribusi ini seharusnya mendekati normal meskipun \(B\) hanya 50. Garis merah adalah grand mean bootstrap, sementara area berbayang menunjukkan rentang ±1 SE.

df_boot <- data.frame(mean_bs = bootstrap.results)
grand_mean_bs <- mean(bootstrap.results)
se_bs         <- sd(bootstrap.results)

pA <- ggplot(df_boot, aes(x = mean_bs)) +
  # Rentang ±1 SE
  annotate(
    "rect",
    xmin  = grand_mean_bs - se_bs,
    xmax  = grand_mean_bs + se_bs,
    ymin  = -Inf, ymax = Inf,
    alpha = 0.12, fill = "#c62828"
  ) +
  geom_histogram(
    aes(y = ..density..),
    bins  = 15,
    fill  = "#c62828", color = "white", alpha = 0.80
  ) +
  geom_density(color = "#7b1fa2", size = 1.2) +
  geom_vline(
    xintercept = grand_mean_bs,
    color = "#880e4f", size = 1.2, linetype = "dashed"
  ) +
  annotate(
    "text",
    x = grand_mean_bs + 0.008, y = Inf,
    vjust = 2, hjust = 0,
    label = paste0("Grand Mean = ", round(grand_mean_bs, 4),
                   "\nSE = ", round(se_bs, 5)),
    color = "#880e4f", fontface = "bold", size = 3.8
  ) +
  scale_x_continuous(breaks = pretty_breaks(6)) +
  labs(
    title    = "Plot A — Distribusi Mean Bootstrap",
    subtitle = paste0("B = ", n.samples, " ulangan  |  Area merah = ±1 SE  |  ",
                      "SE analitik σ/√n = ", round(SIGMA/sqrt(N_OBS), 5)),
    x        = "Nilai Mean Bootstrap",
    y        = "Densitas"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#880e4f", size = 13),
    plot.subtitle = element_text(color = "#ad1457", size = 9)
  )

5.2 Plot B — Distribusi Data Asli

Untuk perbandingan, histogram data asli myData ditampilkan dengan skala densitas yang sama. Perhatikan perbedaan rentang sumbu-X yang dramatis: data asli tersebar dari ~20 hingga ~40, sedangkan distribusi mean bootstrap sangat terkonsentrasi di sekitar 30. Inilah visualisasi langsung dari Standard Error vs Standard Deviation.

df_raw <- data.frame(nilai = myData)

pB <- ggplot(df_raw, aes(x = nilai)) +
  geom_histogram(
    aes(y = ..density..),
    bins  = 40,
    fill  = "#1565c0", color = "white", alpha = 0.80
  ) +
  geom_density(color = "#0d47a1", size = 1.2) +
  stat_function(
    fun  = dnorm, args = list(mean = MU, sd = SIGMA),
    color = "#e65100", size = 1.1, linetype = "dashed"
  ) +
  geom_vline(
    xintercept = mean(myData),
    color = "#1a237e", size = 1.2, linetype = "dashed"
  ) +
  annotate(
    "text",
    x = mean(myData) + 0.15, y = Inf,
    vjust = 2, hjust = 0,
    label = paste0("Mean = ", round(mean(myData), 4),
                   "\nSD = ",  round(sd(myData), 4)),
    color = "#1a237e", fontface = "bold", size = 3.8
  ) +
  scale_x_continuous(breaks = pretty_breaks(8)) +
  labs(
    title    = "Plot B — Distribusi Data Asli (myData)",
    subtitle = paste0("N = ", N_OBS, "  |  Kurva oranye = N(",
                      MU, ", ", SIGMA, "²) teoritis"),
    x        = "Nilai Observasi",
    y        = "Densitas"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1565c0", size = 13),
    plot.subtitle = element_text(color = "#1976d2", size = 9)
  )

5.3 Panel Gabungan (Jawaban Latihan 2)

Ini adalah panel gabungan yang diminta Latihan 2, dibangun dengan patchwork. Susunan vertikal (dua baris) memudahkan perbandingan langsung antara distribusi mean bootstrap (atas) dan distribusi data asli (bawah). Perhatikan perbedaan skala sumbu-X yang sangat kontras — ini secara visual mendemonstrasikan bahwa SE \((\approx 0.08)\) jauh lebih kecil dari SD \((= 2.5)\).

panel_latihan2 <- pA / pB +
  plot_annotation(
    title    = "LATIHAN 2 — Bootstrap vs Distribusi Asli",
    subtitle = paste0(
      "Atas: Distribusi mean dari ", n.samples,
      " sampel bootstrap  ·  Bawah: Distribusi data asli (N = ", N_OBS, ")"
    ),
    caption  = paste0(
      "Muhammad Dzaky Pratama — 2404010021  |  set.seed(150)  |  ",
      "μ = ", MU, ", σ = ", SIGMA
    ),
    theme = theme(
      plot.title    = element_text(size = 16, face = "bold",
                                   color = "#0d0d2b", hjust = 0.5),
      plot.subtitle = element_text(size = 10, color = "#3e3e5e", hjust = 0.5),
      plot.caption  = element_text(color = "grey55", size = 8.5, hjust = 0.5)
    )
  )

panel_latihan2

6 Bagian 4 — Analisis Mendalam & Tambahan

6.1 Konvergensi SE Bootstrap terhadap SE Analitik

Salah satu pertanyaan teoritis penting adalah: bagaimana akurasi SE Bootstrap meningkat seiring bertambahnya \(B\)? Blok ini mensimulasikan SE bootstrap untuk berbagai nilai \(B\) (dari 10 hingga 500) dan membandingkannya dengan nilai analitik. Ini adalah demonstrasi visual dari konvergensi metode bootstrap.

set.seed(42)
B_grid <- c(10, 25, 50, 100, 200, 300, 500)
se_hasil <- numeric(length(B_grid))

for (j in seq_along(B_grid)) {
  B_j <- B_grid[j]
  boot_j <- numeric(B_j)
  for (i in 1:B_j) {
    obs      <- sample(seq_len(N_OBS), size = N_OBS, replace = TRUE)
    boot_j[i] <- mean(myData[obs])
  }
  se_hasil[j] <- sd(boot_j)
}

se_analitik <- SIGMA / sqrt(N_OBS)

df_conv <- data.frame(B = B_grid, SE_bootstrap = se_hasil)

ggplot(df_conv, aes(x = B, y = SE_bootstrap)) +
  geom_hline(
    yintercept = se_analitik,
    color = "#e65100", size = 1.1, linetype = "dashed"
  ) +
  geom_ribbon(
    aes(ymin = se_analitik * 0.85, ymax = se_analitik * 1.15),
    fill = "#e65100", alpha = 0.08
  ) +
  geom_line(color = "#1a237e", size = 1.3) +
  geom_point(
    fill = "#3949ab", color = "white",
    size = 4, shape = 21, stroke = 2
  ) +
  geom_text(
    aes(label = round(SE_bootstrap, 4)),
    vjust = -1.1, size = 3.2, color = "#1a237e", fontface = "bold"
  ) +
  annotate(
    "text",
    x = max(B_grid) * 0.75, y = se_analitik + 0.0008,
    label = paste0("SE Analitik σ/√n = ", round(se_analitik, 5)),
    color = "#e65100", fontface = "bold", size = 4
  ) +
  scale_x_continuous(breaks = B_grid) +
  labs(
    title    = "Konvergensi SE Bootstrap terhadap SE Analitik",
    subtitle = paste0("Semakin besar B → SE Bootstrap mendekati σ/√n = ",
                      round(se_analitik, 5),
                      "  |  Area oranye = toleransi ±15%"),
    x        = "Jumlah Ulangan Bootstrap (B)",
    y        = "Standard Error Bootstrap"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1a237e", size = 13),
    plot.subtitle = element_text(color = "#5c6bc0", size = 9.5)
  )

6.2 Interval Kepercayaan Bootstrap (Metode Persentil)

Interval Kepercayaan Bootstrap (metode persentil) diperoleh langsung dari distribusi empiris estimator bootstrap — tidak memerlukan asumsi distribusi apapun. Untuk CI 95%, kita ambil persentil ke-2.5 dan ke-97.5 dari bootstrap.results. Ini merupakan keunggulan utama bootstrap dibanding metode parametrik.

Secara formal: \(\text{CI}_{95\%} = [\hat{\theta}^*_{(0.025)},\; \hat{\theta}^*_{(0.975)}]\)

# ── CI Bootstrap: Metode Persentil ───────────────────────────
ci_lower <- quantile(bootstrap.results, 0.025)
ci_upper <- quantile(bootstrap.results, 0.975)
ci_width <- ci_upper - ci_lower

# ── CI Analitik (z) untuk perbandingan ───────────────────────
se_anal  <- SIGMA / sqrt(N_OBS)
ci_anal_lo <- mean(myData) - qnorm(0.975) * se_anal
ci_anal_hi <- mean(myData) + qnorm(0.975) * se_anal

cat("── Interval Kepercayaan 95% ──────────────────────────────\n")
## ── Interval Kepercayaan 95% ──────────────────────────────
cat(sprintf("Bootstrap (persentil): [%.5f,  %.5f]  | lebar = %.5f\n",
            ci_lower, ci_upper, ci_width))
## Bootstrap (persentil): [29.80215,  30.05937]  | lebar = 0.25722
cat(sprintf("Analitik (z)         : [%.5f,  %.5f]  | lebar = %.5f\n",
            ci_anal_lo, ci_anal_hi, ci_anal_hi - ci_anal_lo))
## Analitik (z)         : [29.76573,  30.07563]  | lebar = 0.30990
# ── Tabel perbandingan ────────────────────────────────────────
data.frame(
  Metode           = c("Bootstrap (Persentil)", "Analitik (z-interval)"),
  `Batas Bawah`    = round(c(ci_lower, ci_anal_lo), 5),
  `Batas Atas`     = round(c(ci_upper, ci_anal_hi), 5),
  `Lebar CI`       = round(c(ci_width, ci_anal_hi - ci_anal_lo), 5),
  check.names = FALSE
) %>%
  kable(align = "c",
        caption = "Perbandingan CI 95% — Bootstrap Persentil vs Analitik (z)")
Perbandingan CI 95% — Bootstrap Persentil vs Analitik (z)
Metode Batas Bawah Batas Atas Lebar CI
2.5% Bootstrap (Persentil) 29.80215 30.05937 0.25722
Analitik (z-interval) 29.76573 30.07563 0.30990

6.3 Visualisasi CI Bootstrap

Grafik berikut menampilkan secara visual posisi dan lebar CI Bootstrap (merah) dan CI Analitik (biru) pada distribusi mean bootstrap. Kedekatan keduanya mengkonfirmasi keakuratan metode bootstrap meskipun hanya dengan \(B = 50\) ulangan.

df_boot2 <- data.frame(mean_bs = bootstrap.results)

ggplot(df_boot2, aes(x = mean_bs)) +
  # CI Bootstrap
  annotate("rect",
           xmin = ci_lower, xmax = ci_upper,
           ymin = -Inf, ymax = Inf,
           fill = "#c62828", alpha = 0.15) +
  # CI Analitik
  annotate("rect",
           xmin = ci_anal_lo, xmax = ci_anal_hi,
           ymin = -Inf, ymax = Inf,
           fill = "#1565c0", alpha = 0.15) +
  geom_histogram(
    aes(y = ..density..),
    bins = 15, fill = "#ff8f00", color = "white", alpha = 0.75
  ) +
  geom_density(color = "#e65100", size = 1.2) +
  geom_vline(xintercept = ci_lower,   color = "#c62828", linetype = "dashed", size = 0.9) +
  geom_vline(xintercept = ci_upper,   color = "#c62828", linetype = "dashed", size = 0.9) +
  geom_vline(xintercept = ci_anal_lo, color = "#1565c0", linetype = "dotted", size = 0.9) +
  geom_vline(xintercept = ci_anal_hi, color = "#1565c0", linetype = "dotted", size = 0.9) +
  geom_vline(xintercept = mean(bootstrap.results),
             color = "#880e4f", size = 1.3) +
  annotate("text", x = ci_lower - 0.003, y = Inf, vjust = 2, hjust = 1,
           label = paste0("Boot\n", round(ci_lower, 3)),
           color = "#c62828", fontface = "bold", size = 3.2) +
  annotate("text", x = ci_upper + 0.003, y = Inf, vjust = 2, hjust = 0,
           label = paste0("Boot\n", round(ci_upper, 3)),
           color = "#c62828", fontface = "bold", size = 3.2) +
  labs(
    title    = "CI 95%: Bootstrap Persentil (merah) vs Analitik-z (biru)",
    subtitle = "Garis merah putus = batas CI Bootstrap  |  Garis biru titik = batas CI Analitik",
    x        = "Mean Bootstrap", y = "Densitas"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1a237e", size = 13),
    plot.subtitle = element_text(color = "#5c6bc0", size = 9.5)
  )

6.4 Dotplot Seluruh 50 Mean Bootstrap

Dotplot di bawah menampilkan posisi setiap mean bootstrap secara individual, diurutkan, sehingga mudah melihat pola sebaran. Titik berwarna di zona hijau berada dalam rentang ±1 SE dari grand mean, sedangkan titik di zona merah adalah outlier relatif.

df_dot <- data.frame(
  idx     = 1:n.samples,
  mean_bs = sort(bootstrap.results),
  zona    = ifelse(
    abs(sort(bootstrap.results) - grand_mean_bs) <= se_bs,
    "Dalam ±1 SE", "Di luar ±1 SE"
  )
)

ggplot(df_dot, aes(x = idx, y = mean_bs, color = zona)) +
  geom_hline(yintercept = grand_mean_bs,
             color = "#1a237e", linetype = "dashed", size = 0.8) +
  geom_hline(yintercept = grand_mean_bs + se_bs,
             color = "#c62828", linetype = "dotted", size = 0.7) +
  geom_hline(yintercept = grand_mean_bs - se_bs,
             color = "#c62828", linetype = "dotted", size = 0.7) +
  geom_segment(aes(xend = idx, yend = grand_mean_bs),
               color = "grey80", size = 0.5) +
  geom_point(size = 3.5, shape = 21, fill = "white", stroke = 2.5) +
  scale_color_manual(
    values = c("Dalam ±1 SE" = "#2e7d32", "Di luar ±1 SE" = "#c62828"),
    name   = "Zona"
  ) +
  annotate("text", x = 1, y = grand_mean_bs + 0.001, vjust = -0.5,
           label = paste0("Grand Mean = ", round(grand_mean_bs, 4)),
           color = "#1a237e", fontface = "bold", size = 3.5, hjust = 0) +
  labs(
    title    = "Dotplot — 50 Mean Bootstrap (Diurutkan)",
    subtitle = "Garis putus biru = grand mean  |  Garis titik merah = ±1 SE",
    x        = "Peringkat (1 = terkecil)",
    y        = "Nilai Mean Bootstrap"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1a237e", size = 13),
    plot.subtitle = element_text(color = "#5c6bc0", size = 9.5),
    legend.position = "top"
  )

6.5 Panel Akhir — Ringkasan Visual Menyeluruh

Panel komprehensif ini menggabungkan empat perspektif dalam satu tampilan: distribusi data asli, distribusi mean bootstrap, plot konvergensi SE, dan dotplot individual. Panel ini dirancang sebagai ringkasan satu halaman yang dapat digunakan dalam laporan atau presentasi.

# ── Rebuild konvergensi plot untuk panel ─────────────────────
p_conv <- ggplot(df_conv, aes(x = B, y = SE_bootstrap)) +
  geom_hline(yintercept = se_analitik, color = "#e65100",
             size = 1.1, linetype = "dashed") +
  geom_ribbon(aes(ymin = se_analitik * 0.85, ymax = se_analitik * 1.15),
              fill = "#e65100", alpha = 0.08) +
  geom_line(color = "#1a237e", size = 1.2) +
  geom_point(fill = "#3949ab", color = "white", size = 4, shape = 21, stroke = 2) +
  scale_x_continuous(breaks = B_grid) +
  labs(title = "Konvergensi SE Bootstrap",
       subtitle = paste0("Garis oranye = SE analitik σ/√n = ",
                         round(se_analitik, 5)),
       x = "B (jumlah ulangan)", y = "SE Bootstrap") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold", color = "#1a237e"))

# ── Panel gabungan ────────────────────────────────────────────
panel_besar <- (pB | pA) /
               (p_conv | ggplot(df_dot, aes(x = idx, y = mean_bs, color = zona)) +
                  geom_hline(yintercept = grand_mean_bs, color = "#1a237e",
                             linetype = "dashed", size = 0.8) +
                  geom_segment(aes(xend = idx, yend = grand_mean_bs),
                               color = "grey80") +
                  geom_point(size = 3, shape = 21, fill = "white", stroke = 2) +
                  scale_color_manual(
                    values = c("Dalam ±1 SE" = "#2e7d32", "Di luar ±1 SE" = "#c62828"),
                    name = "Zona") +
                  labs(title = "Dotplot 50 Mean Bootstrap",
                       x = "Peringkat", y = "Mean Bootstrap") +
                  theme_minimal(base_size = 11) +
                  theme(plot.title = element_text(face = "bold", color = "#1a237e"),
                        legend.position = "top")
               ) +
  plot_annotation(
    title    = "ANALISIS BOOTSTRAP LENGKAP — WEEK 6",
    subtitle = paste0(
      "Data: N(μ=30, σ=2.5), n=1000  ·  ",
      "Bootstrap: B=50 sampel dengan pengembalian  ·  ",
      "set.seed(150)"
    ),
    caption  = paste0(
      "Muhammad Dzaky Pratama — 2404010021  |  ",
      "SE Bootstrap = ", round(sd(bootstrap.results), 5),
      "  |  SE Analitik = ", round(se_analitik, 5)
    ),
    theme = theme(
      plot.title    = element_text(size = 18, face = "bold",
                                   color = "#0d0d2b", hjust = 0.5),
      plot.subtitle = element_text(size = 10, color = "#3e3e5e", hjust = 0.5),
      plot.caption  = element_text(color = "grey55", size = 9, hjust = 0.5)
    )
  )

panel_besar

7 Bagian 5 — Interpretasi Komprehensif

7.1 Latihan 1: Hasil Numerik

Rangkuman Seluruh Hasil Numerik — Latihan 1 & 2
Item Nilai
Ukuran data asli (N) 1000.00000
Parameter populasi μ 30.00000
Parameter populasi σ 2.50000
Mean sampel myData 29.92068
SD sampel myData 2.47517
Jumlah sampel bootstrap (B) 50.00000
Ukuran tiap sampel bootstrap 1000.00000
Grand Mean bootstrap 29.92612
SE Bootstrap (SD of means) 0.07028
SE Analitik (σ/√n) 0.07906
Selisih SE absolut 0.00878
CI 95% batas bawah (Bootstrap) 29.80215
CI 95% batas atas (Bootstrap) 30.05937

7.2 Diskusi Teoritis

1. Mengapa menggunakan pengembalian (replacement)?

Tanpa pengembalian, setiap sampel bootstrap akan identik dengan data asli dan tidak ada variasi yang dapat diukur. Dengan pengembalian, setiap sampel ke-\(b\) merupakan realisasi berbeda dari “distribusi empiris” \(\hat{F}\), sehingga menghasilkan variasi antar \(\hat{\theta}^*_b\) yang mencerminkan ketidakpastian estimasi.

2. Mengapa \(B = 50\) cukup untuk estimasi mean?

Karena mean adalah estimator smooth dengan varians yang kecil, 50 ulangan sudah cukup stabil. Untuk estimator yang lebih tidak stabil (misal kuantil ekstrem, koefisien regresi), \(B \geq 1000\) lebih disarankan.

3. Interpretasi SE Bootstrap:

\[\widehat{\text{SE}}_{\text{bootstrap}} = 0.07028 \approx \frac{\sigma}{\sqrt{n}} = \frac{2.5}{\sqrt{1000}} = 0.07906\]

Perbedaan sebesar 0.00878 adalah akibat variasi sampling dalam bootstrap itu sendiri — yang akan semakin kecil seiring bertambahnya \(B\).

4. Perbandingan SD vs SE:

Statistik Nilai Makna
SD data asli 2.4752 Variabilitas antar individu dalam populasi
SE bootstrap 0.07028 Variabilitas estimasi mean antar sampel
Rasio SD/SE 35.2 Sesuai teori: ≈ √1000 ≈ 31.6

Rasio ini membuktikan secara empiris bahwa \(\text{SE}(\bar{X}) = \text{SD}(X) / \sqrt{n}\).

7.3 Kesimpulan Akhir

# Temuan Implikasi
1 Bootstrap mean ≈ mean populasi Estimator mean tidak bias
2 SE Bootstrap ≈ SE Analitik Bootstrap valid untuk estimasi ketidakpastian
3 Distribusi mean bootstrap → Normal CLT berlaku meskipun \(B\) kecil
4 SE meningkat akurasinya seiring \(B\) Untuk presisi tinggi, gunakan \(B \geq 500\)
5 CI Bootstrap ≈ CI Analitik Bootstrap dapat menggantikan asumsi distribusi

Simpulan: Bootstrap adalah metode yang sangat fleksibel karena tidak memerlukan asumsi distribusi parametrik. Dengan hanya menggunakan data yang ada dan prinsip resampling sederhana, kita dapat membangun estimasi ketidakpastian yang valid untuk hampir semua jenis estimator statistik.