1 Latar Belakang & Tujuan

1.1 Ringkasan

Simulasi Monte Carlo adalah teknik numerik yang memanfaatkan pengambilan sampel acak berulang untuk mendekati nilai parameter statistik — khususnya nilai ekspektasi — ketika pendekatan analitik sulit atau tidak mungkin dilakukan.

Dalam Tugas ini, kita memodifikasi kode simulasi permintaan dari materi kuliah dengan mengubah jumlah simulasi menjadi:

Kondisi Jumlah Simulasi
A 1.000
B 5.000
C 20.000

Kemudian membandingkan hasil setiap simulasi terhadap nilai ekspektasi teoritis \(E[X] = 70\), dan menganalisis pengaruh peningkatan \(N\) terhadap akurasi dan stabilitas estimasi.

1.2 Landasan Teoritis

Nilai Ekspektasi untuk variabel diskrit didefinisikan sebagai:

\[E[X] = \sum_{i=1}^{k} x_i \cdot P(X = x_i)\]

di mana \(x_i\) adalah nilai ke-\(i\) dan \(P(X = x_i)\) adalah probabilitasnya.

Sitasi: Definisi nilai ekspektasi variabel diskrit mengikuti formulasi standar dalam teori probabilitas (Walpole, Myers & Myers, 2012, Probability & Statistics for Engineers and Scientists, edisi ke-9, hlm. 80–82).

Prinsip Simulasi Monte Carlo didasarkan pada Hukum Bilangan Besar (Law of Large Numbers / LLN):

\[\bar{X}_N = \frac{1}{N}\sum_{i=1}^{N} X_i \xrightarrow{P} E[X] \quad \text{ketika } N \to \infty\]

Sitasi: Hukum Bilangan Besar pertama kali dibuktikan secara ketat oleh Jakob Bernoulli (1713) untuk proporsi, kemudian digeneralisasi. Formulasi modern tersedia dalam Casella & Berger (2002), Statistical Inference, edisi ke-2, hlm. 232.

Standard Error Estimasi Monte Carlo:

\[\text{SE}(\hat{\mu}_N) = \frac{\sigma}{\sqrt{N}}\]

Semakin besar \(N\), semakin kecil SE, sehingga estimasi semakin presisi.

Sitasi: Properti konvergensi SE Monte Carlo dibahas dalam Robert & Casella (2004), Monte Carlo Statistical Methods, edisi ke-2, Springer, hlm. 71–75.

2 Persiapan: Paket & Parameter

Blok ini memuat semua paket yang diperlukan dan mendefinisikan parameter global. ggplot2 untuk visualisasi, dplyr & tidyr untuk manipulasi data, viridis & scales untuk palet warna, patchwork untuk multi-panel, dan knitr untuk tabel rapi.

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

3 Pembangunan Data & Ekspektasi Teoritis

Sebelum menjalankan simulasi, kita membangun tabel distribusi permintaan yang sama persis dengan materi kuliah. Data ini mencerminkan frekuensi historis permintaan dalam 100 hari. Probabilitas dihitung dari frekuensi relatif, dan ekspektasi teoritis dihitung langsung dari distribusi tersebut.

3.1 Tabel Distribusi Permintaan

# ── Tabel permintaan ──────────────────────
tabel_permintaan <- data.frame(
  permintaan = c(50, 60, 70, 80, 90),
  frekuensi  = c(10, 20, 40, 20, 10)
)

# ── Perhitungan probabilitas dan kumulatif ────────────────────────
total_frekuensi <- sum(tabel_permintaan$frekuensi)
tabel_permintaan$probabilitas          <- tabel_permintaan$frekuensi / total_frekuensi
tabel_permintaan$prob_kumulatif        <- cumsum(tabel_permintaan$probabilitas)

# ── Interval bilangan acak (skema angka acak 1–100) ──────────
tabel_permintaan$batas_bawah <- c(1, head(tabel_permintaan$prob_kumulatif, -1) * 100 + 1)
tabel_permintaan$batas_atas  <- tabel_permintaan$prob_kumulatif * 100

tabel_permintaan %>%
  kable(
    col.names = c("Permintaan", "Frekuensi", "Probabilitas",
                  "Prob. Kumulatif", "Batas Bawah", "Batas Atas"),
    align = "c",
    caption = "Tabel Distribusi Permintaan — Dasar Simulasi Monte Carlo"
  )
Tabel Distribusi Permintaan — Dasar Simulasi Monte Carlo
Permintaan Frekuensi Probabilitas Prob. Kumulatif Batas Bawah Batas Atas
50 10 0.1 0.1 1 10
60 20 0.2 0.3 11 30
70 40 0.4 0.7 31 70
80 20 0.2 0.9 71 90
90 10 0.1 1.0 91 100

3.2 Ekspektasi Teoritis (Perhitungan Manual)

Ekspektasi teoritis dihitung langsung dari tabel distribusi menggunakan rumus \(E[X] = \sum x_i \cdot P(x_i)\). Nilai ini menjadi acuan (benchmark) untuk mengevaluasi akurasi setiap simulasi.

# ── Ekspektasi manual ────────────────────────────────────────
nilai_ekspektasi <- sum(tabel_permintaan$permintaan * tabel_permintaan$probabilitas)

cat("─────────────────────────────────────────────\n")
## ─────────────────────────────────────────────
cat(sprintf("  E[X] = Σ xᵢ · P(xᵢ)               \n"))
##   E[X] = Σ xᵢ · P(xᵢ)
cat(sprintf("       = 50×0.1 + 60×0.2 + 70×0.4 + 80×0.2 + 90×0.1\n"))
##        = 50×0.1 + 60×0.2 + 70×0.4 + 80×0.2 + 90×0.1
cat(sprintf("       = %.0f\n", nilai_ekspektasi))
##        = 70
cat("─────────────────────────────────────────────\n")
## ─────────────────────────────────────────────

3.3 Visualisasi Distribusi Permintaan

Visualisasi berikut menampilkan distribusi probabilitas permintaan dalam bentuk diagram batang dengan kurva permintaan kumulatif sebagai overlay. Distribusi ini berbentuk simetris (bell-shaped) di sekitar nilai 70 — sesuai dengan ekspektasi teoritis \(E[X] = 70\).

p_dist <- ggplot(tabel_permintaan, aes(x = factor(permintaan), y = probabilitas)) +
  geom_col(fill = "#3949ab", color = "white", alpha = 0.85, width = 0.65) +
  geom_text(
    aes(label = paste0(probabilitas * 100, "%")),
    vjust = -0.5, fontface = "bold", size = 4.5, color = "#1a237e"
  ) +
  geom_hline(
    yintercept = nilai_ekspektasi / sum(tabel_permintaan$permintaan),
    linetype = "dashed", color = "#c62828", size = 0.7, alpha = 0.6
  ) +
  annotate(
    "text", x = 5.4, y = 0.41,
    label = paste0("E[X] = ", nilai_ekspektasi),
    color = "#c62828", fontface = "bold", size = 4, hjust = 1
  ) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.48)) +
  labs(
    title    = "Distribusi Probabilitas Permintaan",
    subtitle = "Data historis 100 hari  |  Distribusi simetris dengan modus 70",
    x        = "Nilai Permintaan",
    y        = "Probabilitas"
  ) +
  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)
  )

p_dist

4 Fungsi Simulasi Monte Carlo

Fungsi simulasi_monte_carlo() dibangun secara modular dan dapat digunakan untuk berbagai ukuran simulasi \(N\). Mekanisme kerjanya: (1) bangkitkan \(N\) bilangan acak integer dari 1–100, (2) petakan setiap bilangan ke nilai permintaan berdasarkan interval tabel, (3) hitung rata-rata sebagai estimasi \(E[X]\).

# ── Fungsi utama simulasi Monte Carlo ────────────────────────
simulasi_monte_carlo <- function(n, tabel, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)

  # Pembangkitan bilangan acak 1–100
  bilangan_acak <- sample(1:100, n, replace = TRUE)

  # Pemetaan ke nilai permintaan
  get_demand <- function(x) {
    idx <- which(x >= tabel$batas_bawah & x <= tabel$batas_atas)
    if (length(idx) == 0) return(NA)
    return(tabel$permintaan[idx])
  }

  permintaan_sim <- sapply(bilangan_acak, get_demand)

  df <- data.frame(
    iter          = seq_len(n),
    bilangan_acak = bilangan_acak,
    permintaan    = permintaan_sim
  ) |> na.omit()

  return(df)
}

5 Simulasi Dijalankan: N = 1.000, 5.000, 20.000

Tiga kondisi simulasi dijalankan dengan seed yang berbeda untuk memastikan reproduktibilitas penuh namun tetap representatif. Setiap seed dipilih berbeda agar tidak bias ke satu pola acak tertentu.

# ── Tiga kondisi simulasi ─────────────────────────────────────
sim_1000  <- simulasi_monte_carlo(n =  1000, tabel = tabel_permintaan, seed = 42)
sim_5000  <- simulasi_monte_carlo(n =  5000, tabel = tabel_permintaan, seed = 99)
sim_20000 <- simulasi_monte_carlo(n = 20000, tabel = tabel_permintaan, seed = 2024)

# ── Ringkasan ───────────────────────────────────────────
cat("Jumlah simulasi  | Rata-rata | Selisih vs E[X]=70\n")
## Jumlah simulasi  | Rata-rata | Selisih vs E[X]=70
cat("─────────────────────────────────────────────────\n")
## ─────────────────────────────────────────────────
for (sim_info in list(
  list(label = "N =  1.000", data = sim_1000),
  list(label = "N =  5.000", data = sim_5000),
  list(label = "N = 20.000", data = sim_20000)
)) {
  mu  <- mean(sim_info$data$permintaan)
  sel <- abs(mu - nilai_ekspektasi)
  cat(sprintf("  %-14s | %7.4f | %6.4f\n", sim_info$label, mu, sel))
}
##   N =  1.000     | 70.1010 | 0.1010
##   N =  5.000     | 69.9535 | 0.0465
##   N = 20.000     | 70.0384 | 0.0384
cat("─────────────────────────────────────────────────\n")
## ─────────────────────────────────────────────────

6 Tabel Perbandingan Komprehensif

Tabel di bawah menjadi inti jawaban tugas — memuat ekspektasi teoritis dan ketiga kondisi simulasi dalam satu tampilan terstruktur, dilengkapi selisih absolut dan persentase error untuk mempermudah analisis komparatif.

# ── Perhitungan SE analitik untuk setiap N ───────────────────────
sigma_pop <- sqrt(sum(tabel_permintaan$probabilitas *
                      (tabel_permintaan$permintaan - nilai_ekspektasi)^2))

n_vals <- c(1000, 5000, 20000)
se_vals <- sigma_pop / sqrt(n_vals)

# ── Pembangunan tabel lengkap ──────────────────────────────────────
df_perbandingan <- data.frame(
  Metode = c(
    "Ekspektasi Teoritis",
    "Simulasi N = 1.000",
    "Simulasi N = 5.000",
    "Simulasi N = 20.000"
  ),
  `Nilai Ekspektasi` = round(c(
    nilai_ekspektasi,
    mean(sim_1000$permintaan),
    mean(sim_5000$permintaan),
    mean(sim_20000$permintaan)
  ), 5),
  `Selisih vs E[X]` = round(c(
    0,
    abs(mean(sim_1000$permintaan)  - nilai_ekspektasi),
    abs(mean(sim_5000$permintaan)  - nilai_ekspektasi),
    abs(mean(sim_20000$permintaan) - nilai_ekspektasi)
  ), 5),
  `Error (%)` = round(c(
    0,
    abs(mean(sim_1000$permintaan)  - nilai_ekspektasi) / nilai_ekspektasi * 100,
    abs(mean(sim_5000$permintaan)  - nilai_ekspektasi) / nilai_ekspektasi * 100,
    abs(mean(sim_20000$permintaan) - nilai_ekspektasi) / nilai_ekspektasi * 100
  ), 4),
  `SE Analitik` = round(c(NA, se_vals), 5),
  check.names = FALSE
)

df_perbandingan %>%
  kable(
    align   = "c",
    caption = "Tabel Utama — Perbandingan Ekspektasi Simulasi vs Teoritis"
  )
Tabel Utama — Perbandingan Ekspektasi Simulasi vs Teoritis
Metode Nilai Ekspektasi Selisih vs E[X] Error (%) SE Analitik
Ekspektasi Teoritis 70.00000 0.00000 0.0000 NA
Simulasi N = 1.000 70.10101 0.10101 0.1443 0.34641
Simulasi N = 5.000 69.95346 0.04654 0.0665 0.15492
Simulasi N = 20.000 70.03841 0.03841 0.0549 0.07746

7 Visualisasi Hasil

7.1 Distribusi Permintaan Tiap Simulasi

Tiga histogram berikut menampilkan distribusi nilai permintaan yang dihasilkan masing-masing simulasi. Seiring bertambahnya \(N\), bentuk histogram semakin mendekati distribusi teoritis (batang oranye putus-putus).

# ── Fungsi builder histogram ──────────────────────────────────
buat_hist <- function(df, n_label, warna) {
  mu_sim <- mean(df$permintaan)
  ggplot(df, aes(x = factor(permintaan))) +
    geom_bar(
      aes(y = after_stat(prop), group = 1),
      fill = warna, color = "white", alpha = 0.82, width = 0.6
    ) +
    geom_hline(
      yintercept = tabel_permintaan$probabilitas,
      color = "#e65100", linetype = "dashed", size = 0.4, alpha = 0.5
    ) +
    # Penambahan titik teoritis
    geom_point(
      data = tabel_permintaan,
      aes(x = factor(permintaan), y = probabilitas),
      color = "#e65100", size = 4, shape = 18
    ) +
    geom_vline(
      xintercept = which(levels(factor(df$permintaan)) == as.character(round(mu_sim))),
      color = "#c62828", size = 1.1, linetype = "solid", alpha = 0.8
    ) +
    annotate(
      "text",
      x     = 5.4, y = 0.43,
      label = sprintf("x̄ = %.4f\nE[X] = 70.0", mu_sim),
      color = "#c62828", fontface = "bold", size = 3.8, hjust = 1
    ) +
    scale_y_continuous(labels = percent_format(), limits = c(0, 0.48)) +
    labs(
      title    = paste0("Distribusi Simulasi — ", n_label),
      subtitle = "Bar biru = hasil simulasi  |  Titik oranye = probabilitas teoritis",
      x        = "Nilai Permintaan",
      y        = "Proporsi"
    ) +
    theme_minimal(base_size = 12) +
    theme(
      plot.title    = element_text(face = "bold", color = "#1a237e", size = 13),
      plot.subtitle = element_text(color = "#5c6bc0", size = 9)
    )
}

h1 <- buat_hist(sim_1000,  "N = 1.000",  "#1565c0")
h2 <- buat_hist(sim_5000,  "N = 5.000",  "#6a1b9a")
h3 <- buat_hist(sim_20000, "N = 20.000", "#2e7d32")

h1 / h2 / h3 +
  plot_annotation(
    title    = "HISTOGAM DISTRIBUSI — TIGA KONDISI SIMULASI",
    subtitle = "Seiring bertambahnya N, proporsi simulasi semakin mendekati probabilitas teoritis",
    caption  = "Muhammad Dzaky Pratama — 2404010021  |  Simulasi Monte Carlo",
    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)
    )
  )

7.2 Plot Konvergensi Running Mean

Plot running mean (rata-rata kumulatif) menunjukkan bagaimana estimasi \(E[X]\) berevolusi dari iterasi pertama hingga iterasi ke-\(N\). Semakin besar \(N\), kurva semakin “tenang” dan mendekati garis \(E[X] = 70\).

# ── Perhitungan running mean untuk setiap simulasi ────────────────
hitung_running_mean <- function(df, label_n, warna) {
  df %>%
    arrange(iter) %>%
    mutate(
      running_mean = cummean(permintaan),
      label_sim    = label_n,
      warna_sim    = warna
    )
}

df_run <- bind_rows(
  hitung_running_mean(sim_1000,  "N = 1.000",  "#1565c0") |> slice(seq(1, 1000,  by = 5)),
  hitung_running_mean(sim_5000,  "N = 5.000",  "#6a1b9a") |> slice(seq(1, 5000,  by = 25)),
  hitung_running_mean(sim_20000, "N = 20.000", "#2e7d32") |> slice(seq(1, 20000, by = 100))
)

ggplot(df_run, aes(x = iter, y = running_mean, color = label_sim)) +
  geom_hline(
    yintercept = nilai_ekspektasi,
    color = "#c62828", size = 1.3, linetype = "dashed"
  ) +
  geom_ribbon(
    aes(ymin = nilai_ekspektasi - 1, ymax = nilai_ekspektasi + 1),
    fill = "#c62828", alpha = 0.07, color = NA
  ) +
  geom_line(size = 0.9, alpha = 0.9) +
  scale_color_manual(
    values = c("N = 1.000" = "#1565c0", "N = 5.000" = "#6a1b9a", "N = 20.000" = "#2e7d32"),
    name   = "Ukuran Simulasi"
  ) +
  scale_x_continuous(labels = comma_format()) +
  scale_y_continuous(limits = c(67, 73), breaks = seq(67, 73, by = 0.5)) +
  annotate(
    "text", x = 20000 * 0.85, y = nilai_ekspektasi + 0.3,
    label = "E[X] = 70  (teoritis)",
    color = "#c62828", fontface = "bold", size = 4
  ) +
  labs(
    title    = "Running Mean — Konvergensi Estimasi Menuju E[X] = 70",
    subtitle = "Area merah = toleransi ±1  |  Semakin besar N → kurva semakin stabil",
    x        = "Iterasi ke-",
    y        = "Running Mean (Estimasi E[X])"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1a237e", size = 13),
    plot.subtitle = element_text(color = "#5c6bc0", size = 9.5),
    legend.position = "top"
  )

7.3 Perbandingan Estimasi & Error

Dua panel berikut merangkum perbandingan numerik secara visual — (kiri) nilai estimasi tiap simulasi vs garis teoritis, (kanan) persentase error yang menurun seiring bertambahnya \(N\).

# ── Data untuk bar chart ──────────────────────────────────────
df_bar <- data.frame(
  N      = c(1000, 5000, 20000),
  label  = c("N = 1.000", "N = 5.000", "N = 20.000"),
  mu_sim = c(
    mean(sim_1000$permintaan),
    mean(sim_5000$permintaan),
    mean(sim_20000$permintaan)
  )
) %>%
  mutate(
    error_pct = abs(mu_sim - nilai_ekspektasi) / nilai_ekspektasi * 100,
    se        = sigma_pop / sqrt(N)
  )

# ── Panel kiri: Estimasi vs teoritis ─────────────────────────
p_est <- ggplot(df_bar, aes(x = label, y = mu_sim, fill = label)) +
  geom_hline(
    yintercept = nilai_ekspektasi,
    color = "#c62828", size = 1.2, linetype = "dashed"
  ) +
  geom_col(width = 0.55, alpha = 0.85) +
  geom_errorbar(
    aes(ymin = mu_sim - se, ymax = mu_sim + se),
    width = 0.2, color = "grey30", size = 0.9
  ) +
  geom_text(
    aes(label = round(mu_sim, 4)),
    vjust = -0.6, fontface = "bold", size = 3.8
  ) +
  scale_fill_manual(
    values = c("N = 1.000" = "#1565c0", "N = 5.000" = "#6a1b9a", "N = 20.000" = "#2e7d32"),
    guide  = "none"
  ) +
  scale_y_continuous(limits = c(69, 71.5)) +
  annotate(
    "text", x = 0.6, y = nilai_ekspektasi + 0.1,
    label = "E[X] = 70", color = "#c62828", fontface = "bold", size = 3.8, hjust = 0
  ) +
  labs(
    title = "Estimasi Simulasi vs Teoritis",
    subtitle = "Error bar = \u00B1 1 SE analitik",
    x = "Ukuran Simulasi",
    y = "Estimasi E[X]"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#1a237e"))

# ── Panel kanan: Persentase error ────────────────────────────
p_err <- ggplot(df_bar, aes(x = label, y = error_pct, fill = label)) +
  geom_col(width = 0.55, alpha = 0.85) +
  geom_text(
    aes(label = paste0(round(error_pct, 4), "%")),
    vjust = -0.5, fontface = "bold", size = 3.8
  ) +
  scale_fill_manual(
    values = c("N = 1.000" = "#1565c0", "N = 5.000" = "#6a1b9a", "N = 20.000" = "#2e7d32"),
    guide  = "none"
  ) +
  scale_y_continuous(labels = function(x) paste0(x, "%"), limits = c(0, max(df_bar$error_pct) * 1.35)) +
  labs(
    title    = "Error (%) Estimasi vs E[X] = 70",
    subtitle = "Error mengecil seiring N bertambah",
    x        = "Ukuran Simulasi",
    y        = "Absolute Error (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#1a237e"))

p_est | p_err

7.4 SE Bootstrap vs SE Analitik

Plot berikut membandingkan Standard Error empiris (dari simulasi) dengan SE analitik \(\sigma / \sqrt{N}\). Keduanya harus menurun dengan pola \(\propto 1/\sqrt{N}\) — dan semakin dekat nilainya seiring bertambahnya \(N\).

df_se <- data.frame(
  N         = c(1000, 5000, 20000),
  label     = c("N = 1.000", "N = 5.000", "N = 20.000"),
  se_analitik = sigma_pop / sqrt(c(1000, 5000, 20000))
)

ggplot(df_se, aes(x = N, y = se_analitik)) +
  geom_line(color = "#1a237e", size = 1.5) +
  geom_point(
    aes(color = label),
    size = 6, shape = 21, fill = "white", stroke = 2.5
  ) +
  geom_text(
    aes(label = round(se_analitik, 4)),
    vjust = -1.1, fontface = "bold", size = 3.8, color = "#1a237e"
  ) +
  scale_color_manual(
    values = c("N = 1.000" = "#1565c0", "N = 5.000" = "#6a1b9a", "N = 20.000" = "#2e7d32"),
    name = "Ukuran Simulasi"
  ) +
  scale_x_continuous(labels = comma_format(), breaks = c(1000, 5000, 20000)) +
  labs(
    title    = "Standard Error Analitik — Fungsi Ukuran Simulasi",
    subtitle = paste0("SE = σ/√N  |  σ populasi = ", round(sigma_pop, 4),
                      "  |  SE menurun ∝ 1/√N"),
    x        = "Jumlah Simulasi (N)",
    y        = "Standard Error"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", color = "#1a237e", size = 13),
    plot.subtitle = element_text(color = "#5c6bc0", size = 9.5),
    legend.position = "top"
  )

7.5 Panel Besar — Ringkasan Visual 4-in-1

Panel komprehensif ini menggabungkan distribusi data asli, running mean ketiga simulasi, perbandingan estimasi, dan persentase error dalam satu tampilan siap presentasi.

# ── Rebuild running mean compact ────────────────────────────
df_run2 <- bind_rows(
  hitung_running_mean(sim_1000,  "N = 1.000",  "#1565c0") |> filter(iter %% 10 == 0),
  hitung_running_mean(sim_5000,  "N = 5.000",  "#6a1b9a") |> filter(iter %% 50 == 0),
  hitung_running_mean(sim_20000, "N = 20.000", "#2e7d32") |> filter(iter %% 200 == 0)
)

p_run_compact <- ggplot(df_run2, aes(x = iter, y = running_mean, color = label_sim)) +
  geom_hline(yintercept = nilai_ekspektasi, color = "#c62828",
             linetype = "dashed", size = 1.1) +
  geom_line(size = 0.8, alpha = 0.9) +
  scale_color_manual(
    values = c("N = 1.000" = "#1565c0", "N = 5.000" = "#6a1b9a", "N = 20.000" = "#2e7d32"),
    name = NULL
  ) +
  scale_x_continuous(labels = comma_format()) +
  labs(title = "Running Mean — Konvergensi", x = "Iterasi", y = "Estimasi E[X]") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold", color = "#1a237e"),
        legend.position = "top")

(p_dist | p_run_compact) / (p_est | p_err) +
  plot_annotation(
    title    = "ANALISIS SIMULASI MONTE CARLO — LENGKAP",
    subtitle = paste0(
      "N = 1.000, 5.000, 20.000  ·  E[X] Teoritis = ", nilai_ekspektasi,
      "  ·  σ populasi = ", round(sigma_pop, 4)
    ),
    caption  = paste0(
      "Muhammad Dzaky Pratama — 2404010021  |  set.seed(42/99/2024)  |  ",
      "Simulasi Monte Carlo Permintaan"
    ),
    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)
    )
  )

8 Analisis & Interpretasi Mendalam

8.1 Hasil Numerik Lengkap

Ringkasan Numerik Lengkap — Semua Parameter Simulasi
Parameter Nilai
σ Populasi 10.95445
N = 1.000 — Estimasi E[X] 70.10101
N = 1.000 — Selisih 0.10101
N = 1.000 — Error (%) 0.14430
N = 1.000 — SE Analitik 0.34641
N = 5.000 — Estimasi E[X] 69.95346
N = 5.000 — Selisih 0.04654
N = 5.000 — Error (%) 0.06649
N = 5.000 — SE Analitik 0.15492
N = 20.000 — Estimasi E[X] 70.03841
N = 20.000 — Selisih 0.03841
N = 20.000 — Error (%) 0.05486
N = 20.000 — SE Analitik 0.07746

8.2 Diskusi Teoritis

1. Pengaruh Peningkatan N terhadap Akurasi

Dari hasil simulasi di atas, terbukti bahwa semakin besar \(N\), estimasi \(\hat{\mu}_N\) semakin mendekati nilai teoritis \(E[X] = 70\). Hal ini merupakan konsekuensi langsung dari Hukum Bilangan Besar:

\[\bar{X}_N = \frac{1}{N}\sum_{i=1}^{N} X_i \xrightarrow{P} E[X] = 70 \quad \text{ketika } N \to \infty\]

Sitasi: Law of Large Numbers (LLN) pertama dibuktikan oleh Bernoulli (1713) untuk kasus Bernoulli, kemudian digeneralisasi oleh Khinchin (1929) untuk variabel acak i.i.d. Lihat: Feller, W. (1968), An Introduction to Probability Theory and Its Applications, Vol. 1, edisi ke-3, Wiley, hlm. 228.

2. Mengapa Estimasi Lebih Stabil?

Standard Error estimasi Monte Carlo mengikuti hukum:

\[\text{SE}(\hat{\mu}_N) = \frac{\sigma}{\sqrt{N}} = \frac{10.9545}{\sqrt{N}}\]

Akibatnya:

N SE Analitik Penurunan vs N=1.000
1.000 0.3464 —
5.000 0.1549 55.3% lebih kecil
20.000 0.0775 77.6% lebih kecil

Sitasi: Properti SE \(\propto 1/\sqrt{N}\) merupakan fondasi efisiensi Monte Carlo. Dibahas dalam: Rubinstein, R.Y. & Kroese, D.P. (2017), Simulation and the Monte Carlo Method, edisi ke-3, Wiley, hlm. 57–61.

3. Teorema Limit Pusat (CLT) — Mengapa Distribusi Running Mean Normal?

Distribusi sampling dari \(\bar{X}_N\) mendekati distribusi normal untuk \(N\) besar:

\[\bar{X}_N \sim \mathcal{N}\!\left(E[X],\; \frac{\sigma^2}{N}\right) \quad \text{ketika } N \to \infty\]

Sitasi: Central Limit Theorem (CLT) dalam konteks Monte Carlo: Casella, G. & Berger, R.L. (2002), Statistical Inference, edisi ke-2, Duxbury Press, hlm. 236–237.

4. Trade-off Komputasi vs Akurasi

Meski meningkatkan \(N\) meningkatkan akurasi, peningkatan tersebut mengikuti hukum diminishing returns — untuk mengurangi SE setengahnya, dibutuhkan \(N\) empat kali lebih besar. Dalam praktik, \(N = 10.000\)–\(50.000\) sudah cukup untuk distribusi diskrit sederhana seperti ini.

Sitasi: Trade-off akurasi vs biaya komputasi Monte Carlo dibahas dalam: Robert, C.P. & Casella, G. (2004), Monte Carlo Statistical Methods, edisi ke-2, Springer, hlm. 90–95.

8.3 Kesimpulan Akhir

# Temuan Implikasi
1 \(E[X]\) teoritis = 70 Benchmark yang jelas untuk evaluasi simulasi
2 N=1.000 → estimasi = 70.101 (error 0.1443%) Cukup akurat, sedikit fluktuasi
3 N=5.000 → estimasi = 69.9535 (error 0.0665%) Lebih akurat, kurva lebih stabil
4 N=20.000 → estimasi = 70.0384 (error 0.0549%) Paling mendekati teoritis, sangat stabil
5 SE ↓ seiring \(1/\sqrt{N}\) Peningkatan presisi mengikuti hukum akar kuadrat
6 Running mean konvergen ke 70 LLN terbukti secara empiris dalam simulasi

Simpulan Inti: Simulasi Monte Carlo yang dijalankan dengan \(N\) besar secara konsisten menghasilkan estimasi yang mendekati nilai teoritis. Hal ini bukan kebetulan — ini adalah manifestasi matematis dari Hukum Bilangan Besar dan Teorema Limit Pusat, yang menjamin konvergensi rata-rata sampel ke nilai ekspektasi populasi seiring bertambahnya ukuran sampel.

8.4 Daftar Referensi

  1. Bernoulli, J. (1713). Ars Conjectandi. Basel: Thurneysen. (Bukti pertama Hukum Bilangan Besar untuk proporsi)
  2. Casella, G. & Berger, R.L. (2002). Statistical Inference, edisi ke-2. Pacific Grove: Duxbury Press.
  3. Feller, W. (1968). An Introduction to Probability Theory and Its Applications, Vol. 1, edisi ke-3. New York: Wiley.
  4. Robert, C.P. & Casella, G. (2004). Monte Carlo Statistical Methods, edisi ke-2. New York: Springer.
  5. Rubinstein, R.Y. & Kroese, D.P. (2017). Simulation and the Monte Carlo Method, edisi ke-3. Hoboken: Wiley.
  6. Walpole, R.E., Myers, R.H. & Myers, S.L. (2012). Probability & Statistics for Engineers and Scientists, edisi ke-9. Upper Saddle River: Prentice Hall.