# ==============================================================================
# SCRIPT UTAMA: SIMULASI DISTRIBUSI SAMPLING & VISUALISASI HISTOGRAM (ggplot2)
# ==============================================================================

# 1. Mengaktifkan Library yang Dibutuhkan
if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.5.3
if (!require("gridExtra")) install.packages("gridExtra")
## Loading required package: gridExtra
## Warning: package 'gridExtra' was built under R version 4.5.3
library(ggplot2)
library(gridExtra)

# Mengunci nilai acak agar hasil simulasi selalu konsisten
set.seed(42) 

# ==============================================================================
# LANGKAH 1: MEMBUAT POPULASI ASLI (DENGAN DISTRIBUSI SANGAT MENCENG/SKEWED)
# ==============================================================================
# Kita buat populasi besar berdistribusi Eksponensial (Rate = 0.1)
# Rata-rata populasi asli (mu) seharusnya adalah 1/rate = 10
N_populasi <- 100000
populasi <- rexp(N_populasi, rate = 0.1)

mu_populasi <- mean(populasi)
sd_populasi <- sd(populasi)

cat("--- PARAMETER POPULASI ASLI ---\n")
## --- PARAMETER POPULASI ASLI ---
cat(sprintf("Rata-rata Populasi (mu)   : %.4f\n", mu_populasi))
## Rata-rata Populasi (mu)   : 9.9818
cat(sprintf("Standar Deviasi (sigma)   : %.4f\n\n", sd_populasi))
## Standar Deviasi (sigma)   : 10.0059
# ==============================================================================
# LANGKAH 2: PROSES SIMULASI SAMPLING BERULANG (MONTE CARLO SAMPLING)
# ==============================================================================
S_replikasi <- 10000 # Jumlah pengulangan pengambilan sampel

# Inisialisasi vektor kosong untuk menampung rata-rata sampel
rata_rata_n2   <- numeric(S_replikasi)
rata_rata_n10  <- numeric(S_replikasi)
rata_rata_n30  <- numeric(S_replikasi)
rata_rata_n100 <- numeric(S_replikasi)

# Melakukan sampling berulang menggunakan perulangan (loop)
for (s in 1:S_replikasi) {
  rata_rata_n2[s]   <- mean(sample(populasi, size = 2, replace = TRUE))
  rata_rata_n10[s]  <- mean(sample(populasi, size = 10, replace = TRUE))
  rata_rata_n30[s]  <- mean(sample(populasi, size = 30, replace = TRUE))
  rata_rata_n100[s] <- mean(sample(populasi, size = 100, replace = TRUE))
}


# ==============================================================================
# LANGKAH 3: MERAPIKAN DATA KEDALAM SATU DATAFRAME UNTUK GGPLOT2
# ==============================================================================
df_sampling <- data.frame(
  Rata_Rata = c(rata_rata_n2, rata_rata_n10, rata_rata_n30, rata_rata_n100),
  Ukuran_Sampel = factor(rep(c("n = 2", "n = 10", "n = 30", "n = 100"), each = S_replikasi),
                         levels = c("n = 2", "n = 10", "n = 30", "n = 100"))
)


# ==============================================================================
# LANGKAH 4: VISUALISASI HISTOGRAM DENGAN GGPLOT2
# ==============================================================================

# Plot A: Histogram Populasi Asli (Sangat Menceng Kanan)
plot_populasi <- ggplot(data.frame(x = populasi), aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "#ff6b6b", color = "white", alpha = 0.8) +
  geom_vline(xintercept = mu_populasi, color = "black", linetype = "dashed", size = 1) +
  labs(
    title = "1. Distribusi Populasi Asli (Eksponensial)",
    subtitle = paste("Bentuk awal sangat menceng kanan | Rata-rata asli =", round(mu_populasi, 2)),
    x = "Nilai Data", y = "Densitas"
  ) +
  theme_minimal(base_size = 12)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Plot B: Grid Histogram Perbandingan Efek Ukuran Sampel (CLT)
plot_sampling <- ggplot(df_sampling, aes(x = Rata_Rata, fill = Ukuran_Sampel)) +
  geom_histogram(aes(y = after_stat(density)), bins = 50, color = "white", alpha = 0.75) +
  geom_density(color = "black", size = 0.5, linetype = "dashed") +
  geom_vline(xintercept = mu_populasi, color = "red", linetype = "solid", size = 0.8) +
  facet_wrap(~Ukuran_Sampel, scales = "free", ncol = 2) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "2. Distribusi Sampling dari Rata-rata Sampel (10.000 Replikasi)",
    subtitle = "Perhatikan bagaimana bentuk histogram berubah menjadi lonceng (Normal) seiring bertambahnya n",
    x = "Nilai Rata-rata Sampel (x-bar)", y = "Densitas"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

# Menggabungkan kedua plot menjadi satu gambar vertikal
grid.arrange(plot_populasi, plot_sampling, ncol = 1, heights = c(1, 2.2))