# Menetapkan seed spesifik untuk reproducibility
set.seed(84291)
n_sample <- 100
# BAGIAN A: DISTRIBUSI GUMBEL
# Parameter unik: mu = 1.5, beta = 0.8
mu_val <- 1.5
beta_val <- 0.8
# Fungsi Kepekatan Peluang (PDF) Gumbel
dgumbel_custom <- function(x, mu, beta) {
z <- (x - mu) / beta
(1 / beta) * exp(-z - exp(-z))
}
# 1. Inverse-Transform Method
# Formula: X = mu - beta * ln(-ln(U))
sim_gumbel_inversion <- function(n, mu, beta) {
u <- runif(n)
x <- mu - beta * log(-log(u))
return(x)
}
# 2. Acceptance-Rejection Method
# Proposal: Distribusi Logistik Logistic(mu, beta)
# g(x) = exp(-(x-mu)/beta) / (beta * (1 + exp(-(x-mu)/beta))^2)
# Rasio f(x)/g(x) maksimum pada x = mu, c = 4 * exp(-1) = 4/e ≈ 1.4715
sim_gumbel_ar <- function(n, mu, beta) {
c_bound <- 4 * exp(-1) # Nilai optimal supremum
samples <- numeric(n)
accepted <- 0
total_iter <- 0
while (accepted < n) {
total_iter <- total_iter + 1
y_candidate <- rlogis(1, location = mu, scale = beta)
u <- runif(1)
# Rasio penerimaan: f(y) / (c * g(y))
fy <- dgumbel_custom(y_candidate, mu, beta)
gy <- dlogis(y_candidate, location = mu, scale = beta)
if (u <= fy / (c_bound * gy)) {
accepted <- accepted + 1
samples[accepted] <- y_candidate
}
}
attr(samples, "acceptance_rate") <- n / total_iter
return(samples)
}
# Eksekusi Pembangkitan Data Soal A
data_gumbel_inv <- sim_gumbel_inversion(n_sample, mu_val, beta_val)
data_gumbel_ar <- sim_gumbel_ar(n_sample, mu_val, beta_val)
# BAGIAN B: DISTRIBUSI KUMARASWAMY
# Parameter unik: a = 2.5, b = 3.5
a_val <- 2.5
b_val <- 3.5
# Fungsi Kepekatan Peluang (PDF) Kumaraswamy
dkuma_custom <- function(y, a, b) {
a * b * (y^(a - 1)) * ((1 - y^a)^(b - 1))
}
# 1. Inverse-Transform Method
# Formula: Y = (1 - (1 - U)^(1/b))^(1/a)
sim_kuma_inversion <- function(n, a, b) {
u <- runif(n)
y <- (1 - (1 - u)^(1 / b))^(1 / a)
return(y)
}
# 2. Acceptance-Rejection Method
# Proposal: Uniform(0, 1) -> g(y) = 1
# Puncak maksimum f(y) tercapai pada y* = ((a - 1)/(a*b - 1))^(1/a)
sim_kuma_ar <- function(n, a, b) {
y_mode <- ((a - 1) / (a * b - 1))^(1 / a)
c_bound <- dkuma_custom(y_mode, a, b)
samples <- numeric(n)
accepted <- 0
total_iter <- 0
while (accepted < n) {
total_iter <- total_iter + 1
y_candidate <- runif(1, min = 0, max = 1)
u <- runif(1)
# Rasio penerimaan: f(y) / (c * 1)
if (u <= dkuma_custom(y_candidate, a, b) / c_bound) {
accepted <- accepted + 1
samples[accepted] <- y_candidate
}
}
attr(samples, "acceptance_rate") <- n / total_iter
return(samples)
}
# Eksekusi Pembangkitan Data Soal B
data_kuma_inv <- sim_kuma_inversion(n_sample, a_val, b_val)
data_kuma_ar <- sim_kuma_ar(n_sample, a_val, b_val)Assignment 1 - Random Number Generation in R
Penurunan Matematis
Bagian A: Distribusi Gumbel
Fungsi kepekatan peluang (PDF) dari distribusi Gumbel didefinisikan sebagai:
\[f(x) = \frac{1}{\beta} \exp\left(-\frac{x - \mu}{\beta}\right) \exp\left(-e^{-(x - \mu)/\beta}\right), \quad -\infty < x < \infty\]
Parameter yang ditetapkan: \(\mu = 1.5\) dan \(\beta = 0.8\).
Metode Transformasi Invers (Inverse-Transform Method)
Fungsi distribusi kumulatif (CDF):
\[F(x) = \int_{-\infty}^{x} f(t) \, dt = \exp\left(-e^{-(x - \mu)/\beta}\right)\]
Misalkan \(U \sim \text{Uniform}(0, 1)\) dan tetapkan \(u = F(x)\):
\[\ln(u) = -e^{-(x - \mu)/\beta} \implies -\ln(u) = e^{-(x - \mu)/\beta}\]
\[\ln(-\ln(u)) = -\frac{x - \mu}{\beta} \implies x = \mu - \beta \ln(-\ln(u))\]
Untuk \(\mu = 1.5\) dan \(\beta = 0.8\):
\[X = 1.5 - 0.8 \ln(-\ln(U))\]
Metode Penerimaan-Penolakan (Acceptance-Rejection Method)
Distribusi proposal yang digunakan adalah distribusi Logistik dengan lokasi 𝜇 = 1.5 dan skala 𝛽= 0.8:
Rasio densitas target terhadap densitas proposal:
\[\frac{f(x)}{g(x)} = \left(1 + e^{-(x - \mu)/\beta}\right)^2 \exp\left(-e^{-(x - \mu)/\beta}\right)\]
Rasio mencapai maksimum saat \(x = \mu = 1.5\), sehingga diperoleh konstanta pembatas optimum \(c\):
\[c = \sup_{x} \frac{f(x)}{g(x)} = (1 + 1)^2 e^{-1} = \frac{4}{e} \approx 1.4715\]
Kandidat \(Y \sim \text{Logistic}(1.5, 0.8)\) diterima jika bilangan acak \(U \sim \text{Uniform}(0, 1)\) memenuhi:
\[U \le \frac{f(Y)}{c \cdot g(Y)}\]
Bagian B: Distribusi Kumaraswamy
Fungsi kepekatan peluang (PDF) dari distribusi Kumaraswamy didefinisikan sebagai:
\[f(y) = a b y^{a-1} (1 - y^a)^{b-1}, \quad 0 < y < 1\]Parameter yang ditetapkan: \(a = 2.5\) dan \(b = 3.5\).
Metode Transformasi Invers (Inverse-Transform Method)
Fungsi distribusi kumulatif (CDF):
\[F(y) = \int_0^y a b t^{a-1} (1 - t^a)^{b-1} \, dt = 1 - (1 - y^a)^b\]
Misalkan \(U \sim \text{Uniform}(0, 1)\) dan tetapkan \(u = F(y)\):
\[1 - u = (1 - y^a)^b \implies (1 - u)^{1/b} = 1 - y^a\]
\[y = \left(1 - (1 - u)^{1/b}\right)^{1/a}\]
Untuk \(a = 2.5\) dan \(b = 3.5\):
\[Y = \left(1 - (1 - U)^{1/3.5}\right)^{1/2.5}\]
Metode Penerimaan-Penolakan (Acceptance-Rejection Method)
Distribusi proposal yang digunakan adalah \(\text{Uniform}(0, 1)\) dengan PDF konstan \(g(y) = 1\). Titik puncak modus kurva dicari melalui turunan pertama:
\[\frac{d}{dy}\ln f(y) = \frac{a - 1}{y} - \frac{a(b - 1)y^{a-1}}{1 - y^a} = 0 \implies y^* = \left(\frac{a - 1}{ab - 1}\right)^{1/a}\]
Untuk \(a = 2.5\) dan \(b = 3.5\):
\[y^* = \left(\frac{2.5 - 1}{(2.5)(3.5) - 1}\right)^{1/2.5} = \left(\frac{1.5}{7.75}\right)^{0.4} \approx 0.5181\]
Konstanta pembatas optimum \(c\) bernilai:
\[c = f(y^*) = (2.5)(3.5)(0.5181)^{1.5} \left(1 - 0.5181^{2.5}\right)^{2.5} \approx 1.8310\]
Kandidat \(Y \sim \text{Uniform}(0, 1)\) diterima jika \(U \sim \text{Uniform}(0, 1)\) memenuhi:
\[U \le \frac{f(Y)}{c}\]
Answer #1
Visualisasi
par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1.5))
# 1. Gumbel: Inverse Transform
hist(data_gumbel_inv, probability = TRUE, breaks = 12, col = "#A2D2FF",
border = "white", xlab = "x", main = "Gumbel: Inversion (n = 100)",
xlim = c(-1, 6), ylim = c(0, 0.6))
curve(dgumbel_custom(x, mu_val, beta_val), col = "#003049", lwd = 2.5, add = TRUE)
legend("topright", legend = c("Empirical", "True Density"), fill = c("#A2D2FF", NA),
border = c("white", NA), lty = c(NA, 1), col = c(NA, "#003049"), lwd = c(NA, 2.5), bty = "n")
# 2. Gumbel: Acceptance-Rejection
hist(data_gumbel_ar, probability = TRUE, breaks = 12, col = "#FFAFCC",
border = "white", xlab = "x", main = "Gumbel: Accept-Reject (n = 100)",
xlim = c(-1, 6), ylim = c(0, 0.6))
curve(dgumbel_custom(x, mu_val, beta_val), col = "#003049", lwd = 2.5, add = TRUE)
legend("topright", legend = c("Empirical", "True Density"), fill = c("#FFAFCC", NA),
border = c("white", NA), lty = c(NA, 1), col = c(NA, "#003049"), lwd = c(NA, 2.5), bty = "n")
# 3. Kumaraswamy: Inverse Transform
hist(data_kuma_inv, probability = TRUE, breaks = 10, col = "#C8E6C9",
border = "white", xlab = "y", main = "Kumaraswamy: Inversion (n = 100)",
xlim = c(0, 1), ylim = c(0, 2.5))
curve(dkuma_custom(x, a_val, b_val), col = "#2E7D32", lwd = 2.5, add = TRUE)
legend("topright", legend = c("Empirical", "True Density"), fill = c("#C8E6C9", NA),
border = c("white", NA), lty = c(NA, 1), col = c(NA, "#2E7D32"), lwd = c(NA, 2.5), bty = "n")
# 4. Kumaraswamy: Acceptance-Rejection
hist(data_kuma_ar, probability = TRUE, breaks = 10, col = "#FFE082",
border = "white", xlab = "y", main = "Kumaraswamy: Accept-Reject (n = 100)",
xlim = c(0, 1), ylim = c(0, 2.5))
curve(dkuma_custom(x, a_val, b_val), col = "#2E7D32", lwd = 2.5, add = TRUE)
legend("topright", legend = c("Empirical", "True Density"), fill = c("#FFE082", NA),
border = c("white", NA), lty = c(NA, 1), col = c(NA, "#2E7D32"), lwd = c(NA, 2.5), bty = "n")
Penjelasan Visualisasi
Distribusi Gumbel (\(\mu = 1.5, \beta = 0.8\)):
Kedua metode pembangkitan berhasil mereplikasi pola asimetris khas distribusi Gumbel yang condong ke kanan (right-skewed). Konsentrasi data empiris terpusat di sekitar \(x \approx 1.5\) dan distribusi histogram berhimpit secara konsisten dengan kurva kepadatan teoritis (True Density).
Distribusi Kumaraswamy (\(a = 2.5, b = 3.5\)):
Seluruh 100 nilai tergenerasi berada strictly di dalam interval terbuka \((0, 1)\). Kurva empiris membentuk distribusi unimodal dengan puncak di sekitar \(y \approx 0.52\), selaras dengan nilai modus analitis \(y^* \approx 0.5181\).
Perbandingan Metode:
Metode Inverse-Transform bekerja secara deterministik satu-ke-satu tanpa adanya penolakan sampel, sedangkan metode Acceptance-Rejection terbukti menghasilkan distribusi target yang identik melalui proses seleksi titik acak di bawah selubung pembatas
Penjelasan Perbedaan Grafik: Jika kedua histogram metode tersebut terlihat agak berbeda satu sama lain, itu bukanlah suatu kesalahan. Perbedaan bentuk ini sangat wajar terjadi karena kedua metode menggunakan cara kerja acak yang berbeda, dan karena ukuran sampelnya sangat kecil (hanya 100 observasi). Ukuran sampel yang kecil selalu memunculkan variasi acak (sampling variation).