Praktikum Pemrograman Statistika

Randomness in Action: Simulasi Distribusi dengan R

library(tidyverse)

Warning: package 'ggplot2' was built under R version 4.4.3

Warning: package 'dplyr' was built under R version 4.4.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Pendahuluan: Dari Kebetulan ke Simulasi

Dalam dunia statistik modern, randomness bukan sekadar kebetulan, melainkan sebuah alat sains. Pembangkitan bilangan acak menjadi fondasi penting dalam komputasi statistik, terutama ketika kita ingin melakukan simulasi, eksperimen virtual, atau pengujian model tanpa harus selalu mengandalkan data nyata.

Bilangan acak yang dibangkitkan komputer pada dasarnya bukan “benar-benar acak”, melainkan pseudorandom: angka-angka yang dihasilkan oleh algoritma, tetapi cukup menyerupai sifat acak alami. Hebatnya, pseudorandom ini bisa dikendalikan sehingga eksperimen bisa direplikasi dengan hasil identik.

Lebih jauh, bilangan acak yang kita bangkitkan tidak berdiri sendiri, melainkan mengikuti distribusi probabilitas tertentu:

• Fungsi massa peluang (pmf) untuk sebaran diskrit,
• Fungsi kepekatan peluang (pdf) untuk sebaran kontinu,
• serta Fungsi sebaran kumulatif (cdf) untuk melihat peluang kumulatif dari suatu sebaran.

Dengan memahami dan mempraktikkan pembangkitan bilangan acak, kita akan mampu:

• Mensimulasikan data dari berbagai distribusi (normal, uniform, binomial, dsb.).
• Menerapkan metode matematis seperti inverse transform, acceptance–rejection, dan direct transformation.
• Menggunakan bilangan acak untuk membangun dan mengevaluasi model regresi sederhana maupun berganda ataupun model-model lainnya.

Praktikum ini akan mengajak kita untuk melihat aksi randomness secara nyata: bagaimana teori peluang berubah menjadi deretan angka, grafik distribusi, hingga model statistik yang bisa dianalisis.

Fungsi-fungsi Peluang dalam R

Sebagai salah satu bahasa pemrograman yang lahir dari kebutuhan akan komputasi statistik, R menyediakan sekumpulan fungsi bawaan untuk bekerja dengan berbagai sebaran peluang. Fungsi ini punya pola nama yang konsisten dengan prefix berikut:

Prefix	Deskripsi
r-	Pembangkitan bilangan acak dari suatu sebaran
d-	Probability density function (pdf) atau probability mass function (pmf)
p-	Cumulative distribution function (cdf)
q-	Fungsi quantile / invers CDF

Sebagai contoh, dnorm() menghitung kepadatan peluang dari distribusi normal, sedangkan rbinom() digunakan untuk membangkitkan angka acak dari distribusi binomial.

Tabel Fungsi Distribusi di R

Tabel berikut menyajikan daftar distribusi peluang yang paling umum digunakan dalam analisis dan simulasi statistik di R. Masing-masing distribusi disajikan dalam bentuk:

• Notasi matematis formal yang menunjukkan bentuk distribusinya.
• Nama fungsi dasar di R (misalnya norm, binom, pois) yang digunakan dengan awalan d-, p-, q-, atau r- sesuai kebutuhan.
• Parameter-parameter yang harus ditentukan agar fungsi bekerja sesuai distribusi yang dimaksud.
• Rumus probabilitas atau kepadatan peluang (PMF/PDF) sebagai representasi matematis distribusi tersebut.

Pemahaman terhadap struktur distribusi dan cara R mengimplementasikannya akan menjadi bekal penting sebelum kita mulai melakukan simulasi data acak dari distribusi-distribusi ini.

Distribusi	Notasi	R Name	Parameter	PDF / PMF
Beta	\(X \sim \text{Beta}(\alpha, \beta)\)	beta	shape1 = \(\alpha\), shape2 = \(\beta\)	\(f(x) = \dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1}, \; 0 < x < 1\)
Binomial	\(X \sim \text{Bin}(n, p)\)	binom	size = \(n\), prob = \(p\)	\(P(X = x) = \binom{n}{x} p^x (1 - p)^{n - x}, \; x = 0, 1, \dots, n\)
Cauchy	\(X \sim \text{Cauchy}(\theta, \gamma)\)	cauchy	location = \(\theta\), scale = \(\gamma\)	\(f(x) = \dfrac{1}{\pi \gamma \left[1 + \left( \frac{x - \theta}{\gamma} \right)^2 \right]}\)
Chi-squared	\(X \sim \chi^2_v\)	chisq	df = \(v\)	\(f(x) = \dfrac{1}{2^{v/2} \Gamma(v/2)} x^{v/2 - 1} e^{-x/2}, \; x > 0\)
Exponential	\(X \sim \text{Exp}(\lambda)\)	exp	rate = \(\lambda\)	\(f(x) = \lambda e^{-\lambda x}, \; x \geq 0\)
F	\(X \sim F(d_1, d_2)\)	df	df1 = \(d_1\), df2 = \(d_2\)	\(f(x) = \dfrac{\Gamma\left(\frac{d_1 + d_2}{2}\right)}{\Gamma(d_1/2)\Gamma(d_2/2)} \left( \dfrac{d_1}{d_2} \right)^{d_1/2} \dfrac{x^{d_1/2 - 1}}{ \left( 1 + \dfrac{d_1}{d_2} x \right)^{(d_1 + d_2)/2}}, \; x > 0\)
Gamma	\(X \sim \text{Gamma}(\alpha, \beta)\)	gamma	shape = \(\alpha\), scale = \(\beta\)	\(f(x) = \dfrac{1}{\Gamma(\alpha) \beta^\alpha} x^{\alpha - 1} e^{-x/\beta}, \; x \geq 0\)
Geometrik	\(X \sim \text{Geom}(p)\)	geom	prob = \(p\)	\(P(X = x) = (1 - p)^{x - 1} p, \; x = 1, 2, \dots\)
Hypergeometrik	\(X \sim \text{H}(N, m, n)\)	hyper	m = \(m\), n = \(n\), N = \(N\)	\(P(X = x) = \dfrac{\binom{m}{x} \binom{N - m}{n - x}}{\binom{N}{n}}, \; x = 0, 1, \dots, n\)
Log-normal	\(X \sim \text{Lognormal}(\mu, \sigma^2)\)	lnorm	meanlog = \(\mu\), sdlog = \(\sigma\)	\(f(x) = \dfrac{1}{x \sigma \sqrt{2\pi}} e^{- \frac{(\log x - \mu)^2}{2\sigma^2}}, \; x > 0\)
Uniform	\(X \sim U(a, b)\)	unif	min = \(a\), max = \(b\)	\(f(x) = \dfrac{1}{b - a}, \; a \leq x \leq b\)
Neg. Binomial	\(X \sim \text{NB}(r, p)\)	nbinom	size = \(r\), prob = \(p\)	\(P(X = x) = \binom{x - 1}{r - 1} p^r (1 - p)^{x - r}, \; x = r, r+1, \dots\)
Normal	\(X \sim N(\mu, \sigma^2)\)	norm	mean = \(\mu\), sd = \(\sigma\)	\(f(x) = \dfrac{1}{\sigma \sqrt{2\pi}} e^{- \frac{(x - \mu)^2}{2\sigma^2}}, -\infty< x < \infty\)
Poisson	\(X \sim \text{Poisson}(\lambda)\)	pois	lambda = \(\lambda\)	\(P(X = x) = \dfrac{e^{-\lambda} \lambda^x}{x!}, \; x = 0, 1, 2, \dots\)
Student’s t	\(X \sim t(v)\)	t	df = \(v\)	\(f(x) = \dfrac{\Gamma\left( \frac{v+1}{2} \right)}{\sqrt{v\pi} \Gamma(v/2)} \left( 1 + \dfrac{x^2}{v} \right)^{-(v+1)/2}\)

Note

Catatan: Distribusi Eksponensial di R

Pada literatur statistika, distribusi eksponensial dengan dua parameter yang berbeda:

1. Parameter rate (\(\lambda\)):
   \[ f(x) = \lambda e^{-\lambda x}, \quad x \geq 0 \]
   dengan \(E[X] = 1/\lambda\).

2. Parameter scale (\(\theta\)):
   \[ f(x) = \frac{1}{\theta} e^{-x/\theta}, \quad x \geq 0 \]
   dengan \(E[X] = \theta\).

Keduanya konsisten karena \(\lambda = 1/\theta\).

Pada R (rexp(), dexp(), pexp(), qexp()), parameter yang digunakan adalah rate (\(\lambda\)), bukan scale.
Jika notasi yang dipakai adalah \(X \sim \text{Exp}(\theta)\), maka implementasinya di R harus ditulis dengan rate = 1/theta.

Contoh Penggunaan Fungsi Distribusi di R

Mari kita lihat bagaimana R bekerja dengan distribusi peluang melalui contoh distribusi Normal dan Binomial.

Distribusi Normal (Kontinu)

Distribusi Normal merupakan salah satu distribusi paling penting dalam statistika. Pada contoh ini, kita akan:

• Membangkitkan 1000 data acak dari distribusi \(X \sim N(\mu = 10, \sigma = 2)\) menggunakan rnorm(),
• Menampilkan sebagian data dengan head(),
• Memvisualisasikan distribusi data bangkitan melalui histogram,
• Menambahkan kurva PDF teoritis menggunakan dnorm().

set.seed(123)
data_norm <- rnorm(1000, mean = 10, sd = 2)
data_norm <- data.frame(x=data_norm)

head(data_norm)

ggplot(data_norm, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), 
                 fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = 10, sd = 2), 
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = "Distribusi Normal: Histogram vs Kurva PDF",
    x = "Nilai",
    y = "Kepadatan (Density)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
  )

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Pada grafik di atas, histogram mewakili distribusi data hasil simulasi, sedangkan kurva Vermilion menunjukkan bentuk teoretis dari PDF distribusi Normal. Kemiripan bentuk antara keduanya menunjukkan bahwa rnorm() berhasil membangkitkan data acak yang mengikuti karakteristik distribusi Normal yang diharapkan.

Distribusi Binomial (Diskrit)

Distribusi Binomial digunakan untuk memodelkan banyaknya keberhasilan dari sejumlah percobaan Bernoulli. Kita akan:

• Membangkitkan 1000 data acak dari distribusi \(X \sim \text{Bin}(n = 10, p = 0.3)\) menggunakan rbinom(),
• Menampilkan sebagian data dengan head(),
• Menyajikan frekuensi hasil simulasi menggunakan barplot,
• Menambahkan garis nilai probabilitas teoretis dari dbinom().

set.seed(123)
data_binom <- rbinom(1000, size = 10, prob = 0.3)
data_binom <- data.frame(x = data_binom)

head(data_binom)

# Hitung PMF teoritis
pmf_binom <- data.frame(
  x = 0:10,
  prob = dbinom(0:10, size = 10, prob = 0.3)
)

# Visualisasi
ggplot(data_binom, aes(x = factor(x))) +
  geom_bar(fill = "#56B4E9", alpha = 0.7) +
  geom_point(
    data = pmf_binom,
    aes(x = factor(x), y = prob * 1000),  # dikali 1000 agar skala sama
    color = "#D55E00", size = 2
  ) +
  labs(
    title = "Distribusi Binomial: Barplot vs PMF Teoritis",
    x = "Jumlah Keberhasilan (x)",
    y = "Frekuensi"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
  )

Karena distribusi Binomial bersifat diskrit, grafiknya disajikan sebagai barplot berdasarkan frekuensi hasil simulasi. Titik-titik Vermilion menunjukkan nilai probabilitas dari fungsi massa peluang (PMF) yang diperoleh dari dbinom() dan dikalikan dengan 1000 agar berada pada skala yang sama dengan frekuensi observasi. Distribusi simulasi dan probabilitas teoretis terlihat sejalan, memperlihatkan bahwa rbinom() membangkitkan data dengan distribusi yang sesuai secara empiris.

Teknik Pembangkitan Bilangan Acak

Komputer pada dasarnya adalah mesin deterministik—semua operasinya bisa diprediksi. Maka untuk menghasilkan bilangan acak, kita membutuhkan algoritma khusus yang menghasilkan pseudorandom number: bilangan yang tampak acak, namun sebenarnya dibangkitkan dari proses deterministik yang bisa direplikasi.

Ada beberapa teknik utama untuk membangkitkan bilangan acak, terutama dari distribusi tertentu:

Inverse-Transform Method

Konsep dasar: Jika kita memiliki fungsi distribusi kumulatif (CDF) \(F(x)\) dari suatu distribusi, dan sebuah bilangan acak \(U \sim \text{Unif}(0,1)\), maka bilangan acak \(X = F^{-1}(U)\) akan mengikuti distribusi \(F\) tersebut.

Prinsip Inverse-Transform Method

Untuk membangkitkan contoh acak \(X\) dengan sebaran tertentu:

• Tentukan fungsi sebaran kumulatif \(F(x)\) dari sebaran yang diinginkan.
• Hitung invers dari \(F\) atau \(F^{-1}(x)\).
• Bangkitkan contoh acak \(u_1, u_2, u_3, \ldots, u_n\) yang menyebar \(\text{Seragam}(0,1)\).
• Hitung \(x_1 = F^{-1}(u_1),\ x_2 = F^{-1}(u_2),\ x_3 = F^{-1}(u_3),\ \ldots,\ x_n = F^{-1}(u_n)\).

\(x_1, x_2, x_3, \ldots, x_n\) merupakan contoh acak yang saling bebas dari peubah acak \(X\).

Dari PDF ke CDF: Langkah Mendapatkan Fungsi Sebaran Kumulatif

Sebelum menggunakan metode inverse transform, kita perlu mengetahui fungsi distribusi kumulatif (CDF) dari distribusi yang dimaksud. CDF diperoleh dengan mengintegralkan fungsi kepadatan peluang (PDF).

Jika suatu peubah acak kontinu \(X\) memiliki fungsi kepadatan peluang \(f(x)\), maka fungsi distribusi kumulatif \(F(x)\) didefinisikan sebagai:

\[ F(x) = P(X \leq x) = \int_{-\infty}^{x} f(t)\,dt \]

CDF menggambarkan probabilitas kumulatif bahwa \(X\) bernilai kurang dari atau sama dengan \(x\).

Setelah fungsi \(F(x)\) diperoleh, kita bisa mencari fungsi inversnya \(F^{-1}(u)\), yang diperlukan untuk teknik inverse transform. Bila \(U \sim \text{Unif}(0,1)\), maka:

\[ X = F^{-1}(U) \]

akan mengikuti distribusi peluang \(F\).

Contoh: Membangkitan dari Distribusi Eksponensial

Misalkan kita ingin membangkitkan bilangan acak dari distribusi eksponensial dengan parameter \(\theta = 2\), yaitu:

\[ f(x) = \frac{1}{2} e^{-x/2}, \quad x \geq 0 \]

Fungsi distribusi kumulatifnya adalah:

\[ F(x) = \int_{0}^{x} \frac{1}{2} e^{-t/2} \, dt = 1 - e^{-x/2} \]

Kemudian kita cari invers dari CDF tersebut:

\[ F^{-1}(u) = -2 \ln(1 - u) \]

Karena \(U \sim \text{Unif}(0,1)\), maka \(1 - U \sim \text{Unif}(0,1)\) juga. Jadi formulasi dapat disederhanakan menjadi:

\[ X = -2 \ln(U) \]

set.seed(123)
n <- 1000
theta <- 2

# Bangkitkan U ~ Unif(0,1)
u <- runif(n)

# Transformasi inverse
x_exp <- -theta * log(u)
data_exp <- data.frame(x = x_exp)

head(data_exp)

# Visualisasi
ggplot(data_exp, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 fill = "#56B4E9", color = "white", alpha = 0.8) +
  stat_function(fun = dexp, args = list(rate = 1/theta),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = expression("Pembangkitan Acak dari Distribusi Eksponensial"~(theta == 2)),
    x = "Nilai X",
    y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Pada grafik di atas, histogram biru merepresentasikan hasil simulasi pembangkitan bilangan acak dari distribusi eksponensial menggunakan metode inverse transform. Sementara itu, kurva Vermilion adalah fungsi kepadatan peluang (PDF) teoritis dari distribusi eksponensial dengan parameter \(\theta = 2\). Keselarasan bentuk histogram dengan kurva PDF menunjukkan bahwa metode inverse transform berhasil menghasilkan data acak yang mengikuti distribusi eksponensial sesuai teori.

Selain menggunakan rumus manual \(X = -\theta \ln(U)\), R juga menyediakan fungsi qexp() yang merupakan invers dari fungsi distribusi kumulatif eksponensial. Jika \(U \sim \text{Unif}(0,1)\), maka qexp(U, rate = 1/theta) akan menghasilkan nilai acak dari distribusi eksponensial dengan parameter \(\theta\). Dengan kata lain, qexp() mengimplementasikan metode inverse transform secara langsung.

set.seed(123)
n <- 1000
theta <- 2

# Membangkitan langsung dengan qexp()
x_exp_qexp <- qexp(runif(n), rate = 1/theta)

head(x_exp_qexp)

[1] 0.67816835 3.10521872 1.05180043 4.29146021 5.64245855 0.09325366

data_exp_qexp <- data.frame(x = x_exp_qexp)

ggplot(data_exp_qexp, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 fill = "#56B4E9", color = "white", alpha = 0.8) +
  stat_function(fun = dexp, args = list(rate = 1/theta),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = expression("Pembangkitan Acak dengan qexp()"~(theta == 2)),
    x = "Nilai X",
    y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Contoh: Distribusi dengan PDF Khusus

Sekarang, misalkan kita ingin membangkitkan contoh acak dari sebuah distribusi yang tidak tersedia secara langsung dalam fungsi bawaan R, yaitu dengan fungsi kepadatan peluang:

\[ f(x)=3x^2, \quad 0<x<1 \]

Distribusi ini tidak memiliki fungsi r- bawaan di R, sehingga kita perlu menggunakan metode inverse transform.

Langkah 1: Tentukan Fungsi Distribusi Kumulatif (CDF)

\[ F(x) = \int_{0}^{x} 3t^2 dt= x^3 \]

Langkah 2: Cari Invers dari CDF

Jika \(U \sim \text{Unif}(0,1)\), maka:

\[ F(x) = u \;\Rightarrow\; x^3 = u \;\Rightarrow\; x = u^{\tfrac{1}{3}} \]

Langkah 3: Rumus Pembangkitan

Dengan demikian, untuk membangkitkan bilangan acak yang mengikuti distribusi \(f(x) = 3x^2\), kita dapat menggunakan:

\[ X = U^{\tfrac{1}{3}}, \quad U \sim \text{Unif}(0,1) \]

set.seed(123)
n <- 1000

# 1. Bangkitkan U ~ Uniform(0,1)
u <- runif(n)

# 2. Transformasi inverse
x_custom <- u^(1/3)
data_custom <- data.frame(x = x_custom)
head(data_custom)

ggplot(data_custom, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = function(x) 3*x^2, 
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = expression("Pembangkitan Acak dari"~f(x) == 3*x^2~","~~"0 < x < 1"),
    x = "Nilai X",
    y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Pada grafik di atas, histogram biru adalah hasil simulasi pembangkitan acak dengan metode inverse transform, sedangkan kurva Vermilion menunjukkan bentuk teoritis dari \(f(x) = 3x^2\). Keselarasan keduanya menunjukkan bahwa metode inverse transform dapat digunakan untuk membangkitkan distribusi non-standar yang tidak tersedia dalam fungsi bawaan R.

Acceptance-Rejection Method

Konsep Dasar: Metode Acceptance–Rejection digunakan ketika distribusi target sulit atau tidak praktis untuk dibangkitkan langsung dengan inverse transform. Idenya adalah:

• Kita pilih distribusi pembangkit yang lebih sederhana (disebut proposal distribution, misalnya uniform atau eksponensial).
• Kita bangkitkan data dari distribusi proposal.
• Kita terima atau tolak contoh tersebut berdasarkan probabilitas tertentu sehingga contoh yang diterima mengikuti distribusi target.

Prinsip Acceptance–Rejection Method

Misalkan:

• \(f(x)\) = PDF target (yang ingin kita bangkitkan).
• \(g(x)\) = PDF proposal (mudah dibangkitkan).
• \(C>0\) adalah konstanta sehingga

\[ f(x)≤Cg(x), \quad ∀x \] Algoritma:

1. Bangkitkan \(Y∼g(x)\).
2. Bangkitkan \(U∼Unif(0,1)\).
3. Jika \(U≤\frac{f(Y)}{C\times g(Y)}\) maka terima \(Y\) sebagai sampel dari distribusi target. Jika tidak, tolak \(Y\) dan ulangi.

set.seed(123)
n_target <- 2000
alpha <- 5
beta <- 12

# Target distribusi Beta(5,12)
f <- function(x) dbeta(x, alpha, beta)

# Proposal Uniform(0,1)
g <- function(x) dunif(x, 0, 1)

# Cari konstanta C (approx pakai maksimum empiris)
x_seq <- seq(0, 1, length.out = 1000)
C <- max(f(x_seq) / g(x_seq))

# Acceptance-Rejection sampai accepted = n_target
accepted_x <- c()
rejected_x <- c()
rejected_y <- c()
accepted_y <- c()

while(length(accepted_x) < n_target){
  y <- runif(1)          # proposal
  u <- runif(1)          # uniform(0,1)
  
  if(u <= f(y)/(C*g(y))){
    accepted_x <- c(accepted_x, y)
    accepted_y <- c(accepted_y, u*C*g(y))
  } else {
    rejected_x <- c(rejected_x, y)
    rejected_y <- c(rejected_y, u*C*g(y))
  }
}

# Buat data frame
df_points <- data.frame(
  x = c(accepted_x, rejected_x),
  y = c(accepted_y, rejected_y),
  accepted = c(rep(TRUE, length(accepted_x)), rep(FALSE, length(rejected_x)))
)

df_curves <- data.frame(
  x = x_seq,
  target = f(x_seq),
  proposal = C*g(x_seq)
)

# Hitung ringkasan
total_iter <- nrow(df_points)
total_accept <- length(accepted_x)
total_reject <- length(rejected_x)

# Plot
ggplot() +
  geom_point(data = df_points, aes(x = x, y = y, color = accepted, shape = accepted), alpha=0.6) +
  scale_color_manual(values = c("red", "blue"), labels=c("Rejected", "Accepted")) +
  scale_shape_manual(values = c(4, 1), labels=c("Rejected", "Accepted")) +
  geom_line(data = df_curves, aes(x = x, y = target), color="black", linewidth=0.8) +
  geom_line(data = df_curves, aes(x = x, y = proposal), color="darkgreen", linetype="dashed", linewidth=0.8) +
  labs(
    title = "Acceptance-Rejection Method (Target: Beta(5,12))",
    subtitle = paste0("Total iterasi = ", total_iter,
                      " | Total diterima = ", total_accept,
                      " | Total ditolak = ", total_reject),
    x = "x", y = "Density",
    color = "Status", shape = "Status"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size=14, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=11, hjust=0.5))

Gambar di atas memperlihatkan ilustrasi proses Acceptance–Rejection Method untuk membangkitkan sampel dari distribusi Beta(5,12). Kurva hitam merepresentasikan PDF target \(f(x)\), sedangkan garis hijau putus-putus menggambarkan \(Cg(x)\), yaitu distribusi proposal Uniform(0,1) yang telah diskalakan dengan konstanta \(C\) sehingga selalu berada di atas kurva target.

Titik-titik biru (lingkaran) menunjukkan sampel proposal yang diterima karena jatuh di bawah kurva target, sementara titik-titik merah (silang) adalah sampel yang ditolak. Hanya titik-titik yang diterima yang membentuk distribusi sesuai dengan distribusi target Beta(5,12).

Dari visualisasi ini dapat dipahami bahwa meskipun sebagian besar titik proposal ditolak, distribusi akhir dari titik-titik yang diterima akan mengikuti pola distribusi target. Efisiensi metode ini bergantung pada nilai konstanta \(C\); semakin kecil \(C\), semakin banyak titik yang diterima sehingga proses lebih efisien. Namun, jika \(C\) dipilih terlalu kecil, maka ada kemungkinan kurva target \(f(x)\) tidak sepenuhnya berada di bawah \(Cg(x)\). Akibatnya, sebagian sampel valid bisa salah ditolak, dan distribusi yang dihasilkan tidak lagi sesuai dengan target. Sebaliknya, jika \(C\) dipilih terlalu besar, semua sampel memang valid (selalu \(f(x)≤Cg(x)\)), tetapi efisiensi metode menjadi rendah karena proporsi titik yang diterima sangat kecil.

Oleh karena itu, nilai \(C\) biasanya dipilih sebagai nilai maksimum dari rasio \(f(x)/g(x)\).

Contoh: Membangkitkan Data Berdistribusi Beta

Pada contoh ini, kita akan membangkitkan sampel acak dari distribusi Beta(2,2) menggunakan metode Acceptance–Rejection. Distribusi target Beta(2,2) memiliki fungsi kepadatan peluang (PDF) sebagai berikut:

\[ $f(x) = \dfrac{\Gamma(2+2)}{\Gamma(2)\Gamma(2)} x^{2-1}(1-x)^{2-1}=6x(1-x), \; 0 < x < 1$ \]

Fungsi ini bersifat simetris dan berbentuk seperti parabola terbalik, dengan puncaknya di tengah interval. Untuk membangkitkan data dengan metode Acceptance–Rejection, diperlukan distribusi proposal \(g(x)\) yang mudah dibangkitkan dan mendominasi \(f(x)\) di seluruh domainnya. Dalam hal ini, digunakan distribusi Uniform(0,1) sebagai proposal, dengan PDF:

\[ g(x)=1, \quad 0<x<1 \]

Selanjutnya, perlu dicari nilai konstanta \(C\) sedemikian hingga \(f(x)≤Cg(x\)) untuk semua \(x∈(0,1)\). Dengan kata lain, kita mencari nilai maksimum dari fungsi rasio berikut:

\[ h(x)=\frac{f(x)}{g(x)}=\frac{6x(1-x)}{1}=6x(1-x) \]

Untuk memperoleh nilai maksimum dari \(h(x)\), diturunkan fungsi tersebut terhadap \(x\), yaitu:

\[ \frac{dh}{dx}=\frac{6x(1-x)}{dx}=6(1-x)-6x=6-12x \]

Dengan menyelesaikan \(\frac{dh}{dx}=0\), diperoleh nilai \(x=\frac{1}{2}\) sebagai titik stasioner. Substitusi kembali ke fungsi \(h(x)\) menghasilkan:

\[ h(\frac{1}{2})=6(\frac{1}{2})(1-\frac{1}{2}))=\frac{6}{4}=1.5 \]

Dengan demikian, konstanta \(C=1.5\). Nilai ini menunjukkan batas atas dari rasio antara fungsi target dan proposal di seluruh domain.

Untuk memverifikasi nilai maksimum tersebut secara numerik, dapat digunakan fungsi optimize() dalam R sebagai berikut:

h <- function(x) 6 * x * (1 - x)
optimize(h, interval = c(0, 1), maximum = TRUE)$objective

[1] 1.5

Hal ini mengonfirmasi bahwa nilai maksimum dari rasio \(\frac{f(x)}{g(x)}\) adalah \(1.5\), yang digunakan sebagai konstanta skala dalam algoritma Acceptance–Rejection.

Note

Penentuan Konstanta Skala ( C )

Dalam praktik komputasi, pendekatan untuk mendapatkan nilai konstanta ( C ) dapat dilakukan dengan mengevaluasi rasio pada grid titik yang rapat (misalnya 1000 titik pada interval ([0,1])). Dengan pemilihan seperti ini, ( C ) cukup besar untuk menjamin kondisi ( f(x) Cg(x) ) selalu terpenuhi, tetapi tidak terlalu besar sehingga tetap menjaga efisiensi proses acceptance–rejection.

set.seed(123)
n_iter <- 3000       # jumlah titik yang dicek
alpha <- 2
beta <- 2

# Fungsi target: Beta(2,2) → 6x(1-x)
f <- function(x) 6 * x * (1 - x)

# Fungsi proposal: Uniform(0,1) → g(x) = 1
g <- function(x) 1


# Cari konstanta C secara numerik
x_seq <- seq(0, 1, length.out = 1000)
C <- 1.5

# Simulasi Acceptance–Rejection
accepted_x <- c()
rejected_x <- c()
accepted_y <- c()
rejected_y <- c()

for (i in 1:n_iter) {
  y <- runif(1)  # dari proposal
  u <- runif(1)  # uniform(0,1)
  
  if (u <= f(y) / (C * g(y))) {
    accepted_x <- c(accepted_x, y)
    accepted_y <- c(accepted_y, u * C * g(y))
  } else {
    rejected_x <- c(rejected_x, y)
    rejected_y <- c(rejected_y, u * C * g(y))
  }
}

# Data untuk titik dan kurva
df_points <- data.frame(
  x = c(accepted_x, rejected_x),
  y = c(accepted_y, rejected_y),
  status = factor(
    c(rep("Accepted", length(accepted_x)), rep("Rejected", length(rejected_x))),
    levels = c("Accepted", "Rejected"))
)

df_curve <- data.frame(
  x = x_seq,
  target = f(x_seq),
  proposal = C * g(x_seq)
)

# Plot
library(ggplot2)

ggplot() +
  geom_point(data = df_points, aes(x = x, y = y, color = status, shape = status), alpha = 0.6) +
  scale_color_manual(values = c("blue", "red")) +
  scale_shape_manual(values = c(1, 4)) +
  geom_line(data = df_curve, aes(x = x, y = target), color = "black", linewidth = 0.8) +
  geom_line(data = df_curve, aes(x = x, y = proposal), color = "darkgreen", linetype = "dashed", linewidth = 0.8) +
  labs(
    title = "Acceptance-Rejection Method (Target: Beta(2,2))",
    subtitle = paste0("Total dicek = ", n_iter,
                      " | Diterima = ", length(accepted_x),
                      " | Ditolak = ", length(rejected_x),
                      " | Rasio diterima ≈ ", round(length(accepted_x)/n_iter, 3)),
    x = "x", y = "Kepadatan",
    color = "Status", shape = "Status"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5)
  )

Jika diperhatikan, banyaknya amatan pada hasil akhir tidak berjumlah sebanyak yang kita inginkan. Hal ini disebabkan karena tidak semua amatan yang dibangkitkan dari distribusi proposal memenuhi kriteria penerimaan (acceptance rule) dalam metode Acceptance–Rejection. Oleh karena itu, jika kita ingin memperoleh tepat 1000 sampel dari distribusi target, maka:

• Kita perlu membangkitkan lebih dari 3000 amatan awal, atau
• Kita bisa menentukan terlebih dahulu jumlah amatan target (misalnya 3000), lalu melakukan proses pembangkitan dan pengecekan berulang hingga jumlah tersebut tercapai.

Untuk mengevaluasi hasil dari metode Acceptance–Rejection secara visual, kita dapat membandingkan histogram dari sampel yang diterima dengan kurva fungsi kepadatan peluang (PDF) dari distribusi Beta(2,2). Grafik ini memperlihatkan seberapa baik sampel yang dihasilkan mengikuti bentuk distribusi target yang diinginkan.

# Histogram dari sampel yang diterima vs kurva PDF Beta(2,2)
data_beta <- data.frame(x = accepted_x)

ggplot(data_beta, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 40, fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dbeta, args = list(shape1 = alpha, shape2 = beta),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = expression("Histogram Sampel Diterima vs PDF Beta(2,2)"),
    x = "x",
    y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

Contoh: Membangkitkan Data Berdistribusi Gamma

Misalkan kita ingin membangkitkan data sebanyak 100 dari distribusi \(X \sim \text{Gamma}(\alpha=3, \beta=2)\) dengan PDF:

\[ f(x) = \dfrac{1}{\Gamma(3) 2^3} x^{3 - 1} e^{-x/2}, \; x \geq 0 \]

Untuk proposal \(g(x)\), kita bisa pilih distribusi Eksponensial (karena domainnya sama-sama \(x>0\), dan bentuk kurvanya mirip).

• Proposal: \(g(x)=Exp(\lambda=1/2)\).

Lalu kita cari konstanta \(C\) sehingga \(f(x)≤Cg(x)\).

set.seed(123)
n_target <- 100
alpha <- 3
beta <- 2

# Target distribusi Gamma(alpha, beta)
f <- function(x) dgamma(x, shape = alpha, scale = beta)

# Proposal: Exponential dengan rate = 1/beta
g <- function(x) dexp(x, rate = 1/beta)

# Cari konstanta C (approx maksimum empiris)
x_seq <- seq(0, 30, length.out = 2000) # rentang cukup lebar
C <- max(f(x_seq) / g(x_seq))

# Acceptance-Rejection
accepted_x <- c()
rejected_x <- c()
accepted_y <- c()
rejected_y <- c()

while(length(accepted_x) < n_target){
  y <- rexp(1, rate = 1/beta) # proposal
  u <- runif(1)
  
  if(u <= f(y)/(C*g(y))){
    accepted_x <- c(accepted_x, y)
    accepted_y <- c(accepted_y, u*C*g(y))
  } else {
    rejected_x <- c(rejected_x, y)
    rejected_y <- c(rejected_y, u*C*g(y))
  }
}

# Data
df_points <- data.frame(
  x = c(accepted_x, rejected_x),
  y = c(accepted_y, rejected_y),
  accepted = c(rep(TRUE, length(accepted_x)), rep(FALSE, length(rejected_x)))
)

df_curves <- data.frame(
  x = x_seq,
  target = f(x_seq),
  proposal = C*g(x_seq)
)

head(df_curves)

# Hitung ringkasan
total_iter <- length(accepted_x) + length(rejected_x)
total_accept <- length(accepted_x)
total_reject <- length(rejected_x)

total_iter

[1] 12065

total_accept

[1] 100

total_reject

[1] 11965

# Plot
ggplot() +
  geom_point(data = df_points, aes(x = x, y = y, color = accepted, shape = accepted),
             alpha = 0.5, size = 1.2) +
  scale_color_manual(values = c("red", "blue"), labels = c("Rejected", "Accepted")) +
  scale_shape_manual(values = c(4, 1), labels = c("Rejected", "Accepted")) +
  geom_line(data = df_curves, aes(x = x, y = target), color = "black", linewidth = 0.8) +
  geom_line(data = df_curves, aes(x = x, y = proposal), color = "darkgreen", 
            linetype = "dashed", linewidth = 0.8) +
  labs(
    title = "Acceptance-Rejection Method (Target: Gamma(3,2))",
    subtitle = paste0("Total iterasi = ", total_iter,
                      " | Total diterima = ", total_accept,
                      " | Total ditolak = ", total_reject),
    x = "x", y = "Density",
    color = "Status", shape = "Status"
  ) +
  coord_cartesian(xlim = c(0, 20)) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 11, hjust = 0.5))

Pada visualisasi diatas, terlihat bahwa sebaran titik hasil proposal didominasi oleh titik-titik yang ditolak. Karena jumlah sampel ditolak sangat besar, sumbu-y secara otomatis menyesuaikan ke nilai maksimum dari proposal \(Cg(x)\). Akibatnya, kurva target Gamma(3,2) dan titik-titik accepted (biru) tampak kecil dan kurang jelas.

Untuk mengatasi hal ini, kita perlu membatasi rentang sumbu-y agar fokus pada area di sekitar kurva target. Dengan cara ini, bentuk distribusi target serta titik-titik accepted akan lebih mudah terlihat tanpa “tertutup” oleh skala yang terlalu besar akibat titik rejected.

# Plot
y_max <- max(df_curves$target) * 1.2

ggplot() +
  geom_point(data = df_points, aes(x = x, y = y, color = accepted, shape = accepted),
             alpha = 0.5, size = 1.2) +
  scale_color_manual(values = c("red", "blue"), labels = c("Rejected", "Accepted")) +
  scale_shape_manual(values = c(4, 1), labels = c("Rejected", "Accepted")) +
  geom_line(data = df_curves, aes(x = x, y = target), color = "black", linewidth = 0.8) +
  geom_line(data = df_curves, aes(x = x, y = proposal), color = "darkgreen", 
            linetype = "dashed", linewidth = 0.8) +
  labs(
    title = "Acceptance-Rejection Method (Target: Gamma(3,2))",
    subtitle = paste0("Total iterasi = ", total_iter,
                      " | Total diterima = ", total_accept,
                      " | Total ditolak = ", total_reject),
    x = "x", y = "Density",
    color = "Status", shape = "Status"
  ) +
  coord_cartesian(xlim = c(0, 20), ylim = c(0, y_max)) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 11, hjust = 0.5))

# Data hanya dari sampel accepted
data_gamma <- data.frame(x = accepted_x)

# Histogram vs PDF target
ggplot(data_gamma, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dgamma, args = list(shape = alpha, scale = beta),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = "Histogram Sampel Accepted vs PDF Gamma(3,2)",
    x = "x",
    y = "Density"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Histogram di atas memperlihatkan distribusi dari sampel yang diterima hasil metode Acceptance–Rejection untuk target distribusi Gamma(3,2). Histogram biru menunjukkan sebaran data simulasi, sedangkan kurva oranye merepresentasikan fungsi kepadatan peluang (PDF) teoretis Gamma(3,2).

Meskipun jumlah sampel yang diterima relatif sedikit (hanya 100), pola histogram sudah mulai mengikuti bentuk kurva target: naik di awal, mencapai puncak di sekitar nilai 4–5, kemudian menurun secara eksponensial. Jika jumlah sampel diterima diperbesar (misalnya 2000 atau lebih), bentuk histogram akan semakin halus dan mendekati kurva PDF teoretis.

Hal ini menunjukkan bahwa metode Acceptance–Rejection berhasil menghasilkan data acak yang distribusinya konsisten dengan distribusi target, meskipun efisiensinya sangat dipengaruhi oleh pemilihan konstanta \(C\) dan distribusi proposal \(g(x)\).

Contoh: Membangkitkan Data Berdistribusi Normal(0,1)

Distribusi target yang ingin dibangkitkan adalah distribusi Normal standar (\(X \sim N(\mu=0, \sigma^2=1)\)) , dengan fungsi kepadatan peluang \[ f(x) = \dfrac{1}{\sqrt{2\pi}} e^{- \frac{x^2}{2}}, \quad -\infty< x < \infty \]

Kita pilih distribusi proposal yang lebih sederhana. Dalam hal ini digunakan distribusi seragam pada interval \([−5,5]\), dengan fungsi kepadatan peluang

\[ f(x) = \dfrac{1}{10}, \quad -5 \leq x \leq 5 \]

Kita memerlukan suatu konstanta C>0 yang menjamin bahwa kurva target selalu berada di bawah kurva proposal yang telah diskalakan, atau dengan kata lain

\[ f(x)≤Cg(x), \quad -5 \leq x \leq 5 \]

set.seed(123)
n_target <- 1000
mu <- 0
sigma <- 1

# Target: Normal(0,1)
f <- function(x) dnorm(x, mean = mu, sd = sigma)

# Proposal: Uniform(-5,5)
g <- function(x) dunif(x, min = -5, max = 5)

# Cari konstanta C
x_seq <- seq(-5, 5, length.out = 2000)
C <- max(f(x_seq) / g(x_seq))

# Acceptance-Rejection
accepted_x <- c(); rejected_x <- c()
accepted_y <- c(); rejected_y <- c()

while(length(accepted_x) < n_target){
  y <- runif(1, min = -5, max = 5) # proposal
  u <- runif(1)
  
  if(u <= f(y)/(C*g(y))){
    accepted_x <- c(accepted_x, y)
    accepted_y <- c(accepted_y, u*C*g(y))
  } else {
    rejected_x <- c(rejected_x, y)
    rejected_y <- c(rejected_y, u*C*g(y))
  }
}

# Data
df_points <- data.frame(
  x = c(accepted_x, rejected_x),
  y = c(accepted_y, rejected_y),
  accepted = c(rep(TRUE, length(accepted_x)), rep(FALSE, length(rejected_x)))
)

df_curves <- data.frame(
  x = x_seq,
  target = f(x_seq),
  proposal = C*g(x_seq)
)

# Ringkasan
total_iter <- length(accepted_x) + length(rejected_x)
total_accept <- length(accepted_x)
total_reject <- length(rejected_x)

# Plot proses AR
ggplot() +
  geom_point(data = df_points, aes(x=x, y=y, color=accepted, shape=accepted), alpha=0.5) +
  scale_color_manual(values=c("red","blue"), labels=c("Rejected","Accepted")) +
  scale_shape_manual(values=c(4,1), labels=c("Rejected","Accepted")) +
  geom_line(data = df_curves, aes(x=x, y=target), color="black", linewidth=0.8) +
  geom_line(data = df_curves, aes(x=x, y=proposal), color="darkgreen", linetype="dashed", linewidth=0.8) +
  labs(
    title="Acceptance-Rejection Method (Target: Normal(0,1))",
    subtitle=paste0("Total iterasi = ", total_iter,
                    " | Total diterima = ", total_accept,
                    " | Total ditolak = ", total_reject),
    x="x", y="Density",
    color="Status", shape="Status"
  ) +
  coord_cartesian(xlim=c(-5,5)) +
  theme_minimal() +
  theme(plot.title=element_text(size=14, face="bold", hjust=0.5),
        plot.subtitle=element_text(size=11, hjust=0.5))

# Data hanya dari sampel accepted
data_norm_ar <- data.frame(x = accepted_x)

# Histogram vs PDF target Normal(0,1)
ggplot(data_norm_ar, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = "Histogram Sampel Accepted vs PDF Normal(0,1)",
    x = "x",
    y = "Density"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Histogram biru menunjukkan distribusi sampel yang diterima dari proses Acceptance–Rejection, sedangkan kurva vermilion adalah fungsi kepadatan peluang Normal(0,1). Dengan jumlah sampel yang cukup (misalnya 1000 atau lebih), histogram terlihat semakin mendekati bentuk kurva Normal standar. Hal ini membuktikan bahwa metode Acceptance–Rejection mampu menghasilkan data acak yang distribusinya sesuai dengan distribusi target meskipun menggunakan proposal sederhana seperti Uniform(-5,5).

Setelah melihat bagaimana metode Acceptance–Rejection dapat digunakan untuk membangkitkan sampel dari distribusi Normal(0,1), langkah selanjutnya adalah memahami pengaruh ukuran sampel terhadap hasil simulasi. Jumlah sampel yang diterima sangat menentukan seberapa halus histogram hasil simulasi dibandingkan dengan bentuk teoretis distribusi target. Dengan jumlah sampel yang kecil, histogram cenderung kasar dan tidak sepenuhnya mencerminkan pola distribusi target. Sebaliknya, semakin banyak sampel yang diterima, semakin baik pula histogram mendekati kurva kepadatan peluang (PDF) teoretis. Untuk mengilustrasikan hal ini, mari kita bandingkan hasil Acceptance–Rejection dengan jumlah sampel diterima sebanyak \(n=100\) dan \(n=2000\).

set.seed(123)
mu <- 0
sigma <- 1

# Target dan proposal
f <- function(x) dnorm(x, mean = mu, sd = sigma)
g <- function(x) dunif(x, min = -5, max = 5)

# Fungsi AR untuk menghasilkan n_target sampel
AR_sampler <- function(n_target, f, g){
  # Cari konstanta C
  x_seq <- seq(-5, 5, length.out = 2000)
  C <- max(f(x_seq) / g(x_seq))
  
  accepted_x <- c()
  while(length(accepted_x) < n_target){
    y <- runif(1, min = -5, max = 5) # proposal
    u <- runif(1)
    if(u <= f(y)/(C*g(y))){
      accepted_x <- c(accepted_x, y)
    }
  }
  return(accepted_x)
}

# Bangkitkan data
accepted_100 <- AR_sampler(100, f, g)
accepted_2000 <- AR_sampler(2000, f, g)

# Buat data frame gabungan
df_hist <- data.frame(
  x = c(accepted_100, accepted_2000),
  Sampel = c(rep("n = 100", length(accepted_100)),
             rep("n = 2000", length(accepted_2000)))
)

# Plot histogram berdampingan
ggplot(df_hist, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma),
                color = "#D55E00", linewidth = 0.8) +
  facet_wrap(~ Sampel, ncol = 2) +
  labs(
    title = "Perbandingan jumlah sampel n=100 dan n=2000",
    x = "x",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5)
  )

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Pada gambar di atas ditampilkan hasil simulasi distribusi Normal(0,1) dengan metode Acceptance–Rejection untuk dua ukuran sampel yang berbeda. Panel kiri menunjukkan histogram dari \(n=100\) sampel yang diterima, sedangkan panel kanan memperlihatkan histogram dari \(n=2000\) sampel.

Dapat dilihat bahwa ketika jumlah sampel relatif kecil (\(n=100\)), histogram masih terlihat kasar, bergelombang, dan belum sepenuhnya mengikuti bentuk kurva teoretis (garis merah). Namun, ketika jumlah sampel diperbesar (\(n=2000\)), histogram semakin halus dan bentuknya jauh lebih menyerupai kurva kepadatan peluang Normal standar.

Hal ini menegaskan prinsip dasar simulasi Monte Carlo: semakin banyak sampel yang dibangkitkan, semakin baik distribusi hasil simulasi mendekati distribusi target yang sebenarnya. Materi Monte Carlo akan dibahas pada pertemuan lainnya.

Acceptance–Rejection untuk Distribusi dengan PDF Khusus

Pada praktiknya, tidak semua distribusi target memiliki fungsi r- bawaan di R atau atau merupakan distribusi standar yang sudah dikenal namanya. Ada kalanya kita hanya mengetahui bentuk fungsi kepadatan peluang (PDF), tetapi distribusi tersebut tidak memiliki nama khusus atau tidak tersedia fungsi pembangkitnya di perangkat lunak statistik. Dalam situasi seperti ini, metode Acceptance–Rejection sangat berguna karena tetap memungkinkan kita membangkitkan sampel acak hanya berdasarkan bentuk PDF target yang diketahui.

Sebagai ilustrasi, Misalkan target distribusi acak kita memiliki fungsi kepadatan peluang:

\[ f(x)=\frac{3}{2}x^3+\frac{11}{8}x^2+\frac{1}{6}x+\frac{1}{12},\quad 0≤x≤1 \]

Distribusi ini tidak memiliki nama khusus dalam literatur statistik, tetapi merupakan PDF valid karena memenuhi dua syarat: 1. \(f(x)≥0\) untuk semua \(x∈[0,1]\), 2. \(\int_{0}^{1} f(x)\,dx=1\)

Proposal yang sesuai adalah distribusi Uniform(0,1):

\[ g(x)=1, \quad 0≤x≤1 \] Dengan proposal sederhana ini, konstanta \(C\) diperoleh dari:

\[ C = \max_{x \in (0,1]} \frac{f(x)}{g(x)} = \max_{x \in (0,1]} f(x) \]

Karena \(g(x)=1\), maka \(C\) cukup dicari dengan menghitung maksimum fungsi \(f(x)\) pada interval [0,1]. Setelah itu, algoritma Acceptance–Rejection dapat digunakan untuk menghasilkan sampel dari distribusi target polinomial ini.

# Target: f(x) = 3/2 x^3 + 11/8 x^2 + 1/6 x + 1/12, 0<=x<=1
f <- function(x) {
  (3/2)*x^3 + (11/8)*x^2 + (1/6)*x + (1/12)
}

# Proposal: Uniform(0,1), g(x)=1
g <- function(x) dunif(x, min=0, max=1)

# Cari konstanta C
x_seq <- seq(0, 1, length.out = 1000)
C <- max(f(x_seq) / g(x_seq))  # karena g(x)=1, cukup max f(x)

# Acceptance-Rejection
set.seed(123)
n_target <- 2000
accepted_x <- c(); rejected_x <- c()
accepted_y <- c(); rejected_y <- c()

while(length(accepted_x) < n_target){
  y <- runif(1, min=0, max=1) # proposal
  u <- runif(1)
  
  if(u <= f(y)/(C*g(y))){
    accepted_x <- c(accepted_x, y)
    accepted_y <- c(accepted_y, u*C*g(y))
  } else {
    rejected_x <- c(rejected_x, y)
    rejected_y <- c(rejected_y, u*C*g(y))
  }
}

# Data untuk plotting
df_points <- data.frame(
  x = c(accepted_x, rejected_x),
  y = c(accepted_y, rejected_y),
  accepted = c(rep(TRUE, length(accepted_x)), rep(FALSE, length(rejected_x)))
)

df_curves <- data.frame(
  x = x_seq,
  target = f(x_seq),
  proposal = C*g(x_seq)
)

# Ringkasan
total_iter <- length(accepted_x) + length(rejected_x)
total_accept <- length(accepted_x)
total_reject <- length(rejected_x)

# Plot proses AR
ggplot() +
  geom_point(data = df_points, aes(x=x, y=y, color=accepted, shape=accepted), alpha=0.5) +
  scale_color_manual(values=c("red","blue"), labels=c("Rejected","Accepted")) +
  scale_shape_manual(values=c(4,1), labels=c("Rejected","Accepted")) +
  geom_line(data = df_curves, aes(x=x, y=target), color="black", linewidth=0.8) +
  geom_line(data = df_curves, aes(x=x, y=proposal), color="darkgreen", linetype="dashed", linewidth=0.8) +
  labs(
    title = expression("Acceptance–Rejection Method:"~~
                   f(x) == frac(3, 2)*x^3 + frac(11, 8)*x^2 + frac(1, 6)*x + frac(1, 12)),
    subtitle=paste0("Total iterasi = ", total_iter,
                    " | Total diterima = ", total_accept,
                    " | Total ditolak = ", total_reject),
    x="x", y="Density",
    color="Status", shape="Status"
  ) +
  coord_cartesian(xlim=c(0,1), ylim=c(0,C*1.1)) +
  theme_minimal() +
  theme(plot.title=element_text(size=14, face="bold", hjust=0.5),
        plot.subtitle=element_text(size=11, hjust=0.5))

# Histogram sampel accepted vs PDF target
df_accept <- data.frame(x = accepted_x)

ggplot(df_accept, aes(x)) +
  geom_histogram(aes(y=after_stat(density)),
                 bins = 30, fill="#56B4E9", color="white", alpha=0.7) +
  stat_function(fun = f, color="#D55E00", linewidth=0.8) +
  labs(
    title = expression("Histogram Sampel Accepted vs"~~
                   f(x) == frac(3, 2)*x^3 + frac(11, 8)*x^2 + frac(1, 6)*x + frac(1, 12)),
    x = "x", y = "Density"
  ) +
  theme_minimal() +
  theme(plot.title=element_text(size=14, face="bold", hjust=0.5))

Hasil visualisasi pada histogram menunjukkan bahwa sampel yang dihasilkan melalui metode Acceptance–Rejection mengikuti bentuk distribusi target yang telah ditentukan. Histogram biru yang merepresentasikan sebaran sampel acak memperlihatkan kesesuaian dengan kurva fungsi kepadatan peluang (PDF) target yang digambarkan oleh garis oranye. Hal ini menegaskan bahwa metode Acceptance–Rejection mampu menghasilkan sampel yang mendekati distribusi target, meskipun distribusi tersebut tidak memiliki nama standar atau fungsi pembangkit bawaan dalam perangkat lunak statistik. Seiring dengan bertambahnya jumlah sampel, bentuk histogram cenderung semakin halus dan semakin mendekati kurva teoritis, sehingga akurasi pendekatan distribusi semakin baik. Perbandingan antara histogram dengan jumlah sampel \(n=2000\) dan \(n=100000\) ditampilkan pada gambar di bawah untuk memperlihatkan pengaruh ukuran sampel terhadap kualitas pendekatan distribusi.

# Target PDF: polinomial
f <- function(x) {
  (3/2)*x^3 + (11/8)*x^2 + (1/6)*x + (1/12)
}

# Proposal: Uniform(0,1)
g <- function(x) dunif(x, min=0, max=1)

# Fungsi AR untuk menghasilkan n_target sampel
AR_sampler <- function(n_target, f, g){
  x_seq <- seq(0, 1, length.out = 2000)
  C <- max(f(x_seq) / g(x_seq))
  
  accepted_x <- c()
  while(length(accepted_x) < n_target){
    y <- runif(1, min=0, max=1) # proposal
    u <- runif(1)
    if(u <= f(y)/(C*g(y))){
      accepted_x <- c(accepted_x, y)
    }
  }
  return(accepted_x)
}

# Bangkitkan data
set.seed(123)
accepted_2000 <- AR_sampler(2000, f, g)
accepted_100000 <- AR_sampler(100000, f, g)

# Buat data frame gabungan
df_hist <- data.frame(
  x = c(accepted_2000, accepted_100000),
  Sampel = c(rep("n = 2000", length(accepted_2000)),
             rep("n = 100000", length(accepted_100000)))
)

# Plot histogram berdampingan
ggplot(df_hist, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 30, fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = f, color = "#D55E00", linewidth = 0.8) +
  facet_wrap(~ Sampel, ncol = 2) +
  labs(
    title = "Perbandingan jumlah sampel n=2000 dan n=100000",
    x = "x",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5)
  )

Direct Transform Method

Metode Direct Transform digunakan untuk membangkitkan sampel acak dari distribusi target dengan cara mentransformasi peubah acak lain yang lebih sederhana. Teknik ini memanfaatkan sifat-sifat hubungan antar distribusi dalam teori probabilitas. Dengan demikian, pembangkitan tidak lagi melalui inversi fungsi distribusi kumulatif atau prosedur acceptance–rejection, melainkan melalui transformasi langsung dari distribusi yang sudah diketahui algoritma pembangkitannya. Diagram berikut memperlihatkan hubungan matematis antara berbagai distribusi peluang.

Contoh: Normal dari Dua Uniform (Metode Box-Muller)

Jika \(U_1, U_2 \sim \text{Unif}(0,1)\), maka kita bisa membangkitkan dua variabel normal independen \(Z_1, Z_2 \sim \mathcal{N}(0,1)\) menggunakan Transformasi Box-Muller:

\[ Z_1 = \sqrt{-2 \ln(U_1)} \cos(2\pi U_2) \\ Z_2 = \sqrt{-2 \ln(U_1)} \sin(2\pi U_2) \]

set.seed(123)
n <- 10000
u1 <- runif(n)
u2 <- runif(n)

z1 <- sqrt(-2 * log(u1)) * cos(2 * pi * u2)
z2 <- sqrt(-2 * log(u1)) * sin(2 * pi * u2)

data_normal <- data.frame(z1=z1,z2=z2)

data_normal_long <- data.frame(
  z = c(z1, z2),
  Variabel = rep(c("Z1", "Z2"), each = n)
)

# Facet
ggplot(data_normal_long, aes(x = z)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 50, fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1),
                color = "#D55E00", linewidth = 0.8) +
  facet_wrap(~ Variabel, ncol = 2) +
  labs(
    title = expression("Distribusi"~Normal(0,1)~"dari Transformasi Box–Muller"),
    subtitle = expression("Membandingkan"~Z[1]~"dan"~Z[2]~"dari dua"~Uniform(0,1)),
    x = "Nilai", y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    plot.subtitle = element_text(hjust = 0.5, size = 11)
  )

Contoh: Normal ke Chi-Square

Jika \(Z \sim \mathcal{N}(0,1)\), maka \(X = Z^2\) mengikuti distribusi Chi-square dengan 1 derajat bebas: \(X \sim \chi^2(1)\).

set.seed(123)
z <- rnorm(10000)
x_chisq <- z^2
data_chisq <- data.frame(x = x_chisq)

ggplot(data_chisq, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), 
                 fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dchisq, args = list(df = 1),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = expression("Distribusi"~chi^2*"(1), dengan"~Z^2 %~% N(0,1)),
    x = "x", y = "Kepadatan"
  ) +
  coord_cartesian(xlim = c(0, 10), ylim = c(0, 1)) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Contoh: t-Student dari Normal dan Chi-Square

Jika: \(Z \sim \mathcal{N}(0,1)\) dan \(V \sim \chi^2(\nu)\), dimana \(Z, V\) saling bebas maka:

\[ T = \frac{Z}{\sqrt{V/v}}\sim t(v) \]

set.seed(123)
n <- 10000
df_t <- 5

z <- rnorm(n)
v <- rchisq(n, df = df_t)
t_stat <- z / sqrt(v / df_t)
data_t <- data.frame(x = t_stat)

ggplot(data_t, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 50, fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = dt, args = list(df = df_t),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = bquote("Distribusi t("~.(df_t)~") dari Z dan Chi"^2),
    x = "t", y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Contoh: F-Distribution dari Dua Chi-Square

Jika \(U \sim \chi^2(d_1)\) dan \(V \sim \chi^2(d_2)\), dimana \(U\) dan \(V\) bebas, maka:

\[ F = \frac{U/d_1}{V/d_2}\sim F(d_1,d_2) \]

set.seed(123)
n <- 10000
df1 <- 5
df2 <- 10

u <- rchisq(n, df = df1)
v <- rchisq(n, df = df2)

f_stat <- (u / df1) / (v / df2)
data_f <- data.frame(x = f_stat)

ggplot(data_f, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 60, fill = "#56B4E9", color = "white", alpha = 0.7) +
  stat_function(fun = df, args = list(df1 = df1, df2 = df2),
                color = "#D55E00", linewidth = 0.8) +
  labs(
    title = bquote("Distribusi"~F(.(df1),~.(df2))~"dari Dua"~chi^2),
    x = "F", y = "Kepadatan"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Multivariate Normal Distribution

Distribusi Normal Multivariat (Multivariate Normal Distribution, MND) adalah generalisasi dari distribusi Normal univariat ke dimensi lebih tinggi. Sebuah vektor acak berdimensi \(p\) dikatakan berdistribusi normal multivariat jika setiap kombinasi linear dari komponennya mengikuti distribusi Normal univariat.

Misalkan \(\mathbf{X}\) adalah vektor acak berdimensi \(p\):

\[ \mathbf{X} = \begin{bmatrix} X_1 \\ X_2 \\ \vdots \\ X_p \end{bmatrix} \sim \mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \]

dengan:

• \(\boldsymbol{\mu} = \mathbb{E}[\mathbf{X}]\): vektor mean berdimensi \(p\),
• \(\boldsymbol{\Sigma} = \text{Cov}(\mathbf{X})\): matriks kovarians berdimensi \(p \times p\), yang simetris dan semi-definit positif.
• Fungsi ini terdefinisi untuk semua \(\mathbf{x} \in \mathbb{R}^p\).

Fungsi kepadatan dari distribusi Normal multivariat diberikan oleh:

\[ f(\mathbf{x}) = \frac{1}{(2\pi)^{p/2} |\boldsymbol{\Sigma}|^{1/2}} \exp \left\{ -\frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right\} \]

di mana:

• \(|\boldsymbol{\Sigma}|\) adalah determinan dari matriks kovarians \(\boldsymbol{\Sigma}\),
• \((\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\) disebut kuadrat Mahalanobis,

Beberapa properti penting dari distribusi Normal multivariat:

• Marginal: Subvektor dari \(\mathbf{X}\) juga berdistribusi normal.
• Kondisional: Distribusi kondisional dari subset \(\mathbf{X}\), diberikan subset lainnya, juga merupakan distribusi normal.
• Simetri Eliptik: Kontur kepadatan berbentuk elips di ruang multidimensi.
• Independensi dan Nol Kovarians: Untuk MND, independensi ekuivalen dengan nol kovarians (berlaku dua arah).

Misalkan \(\mathbf{X} = [X_1, X_2]^\top\) mengikuti distribusi normal bivariat dengan parameter:

\[ \boldsymbol{\mu} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \quad \boldsymbol{\Sigma} = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} \]

Maka fungsi kepadatannya menjadi:

\[ f(x_1, x_2) = \frac{1}{2\pi \sqrt{1 - \rho^2}} \exp\left\{ -\frac{1}{2(1 - \rho^2)} \left( x_1^2 - 2\rho x_1 x_2 + x_2^2 \right) \right\} \]

Pembangkitan bilangan acak dari distribusi \(\mathcal{N}_d(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) dilakukan dalam dua tahap utama

Pertama, kita membangkitkan vektor acak standar. Bangkitkan vektor \(\mathbf{Z} = (Z_1, Z_2, \dots, Z_d)\) dari distribusi Normal standar multivariat \(\mathcal{N}_d(\mathbf{0}, \mathbf{I})\).

Kedua, kita melakukan transformasi linear Transformasi vektor \(\mathbf{Z}\) menggunakan rata-rata \(\boldsymbol{\mu}_d\) dan matriks kovarians \(\boldsymbol{\Sigma}\) yang diinginkan:

\[ \mathbf{X} = C \mathbf{Z} + \boldsymbol{\mu} \sim \mathcal{N}_d(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \]

dengan \(C\) merupakan matriks faktor kovarians sedemikian hingga:

\[ \boldsymbol{\Sigma} = C C^\top \]

Pemfaktoran \(\boldsymbol{\Sigma}\) dapat dilakukan melalui beberapa metode dekomposisi matriks:

• Dekomposisi Eigen: eigen() di R
• Dekomposisi Cholesky: chol() di R (paling umum jika \(\Sigma\) positif-definit)
• Dekomposisi Singular Value (SVD): svd() di R (jika matriks singular atau numerik tidak stabil)

Dekomposisi Eigen

Dekomposisi eigen digunakan untuk membentuk transformasi linear terhadap vektor acak \(\mathbf{Z} \sim \mathcal{N}_d(\mathbf{0}, \mathbf{I})\) sehingga menghasilkan vektor \(\mathbf{X} \sim \mathcal{N}_d(\boldsymbol{\mu}, \boldsymbol{\Sigma})\).

Jika \(\boldsymbol{\Sigma}\) adalah matriks kovarians simetris dan semi-definit positif, maka ia dapat didekomposisi sebagai:

\[ \boldsymbol{\Sigma} = V \Lambda V^\top \]

dengan:

• \(V\): matriks eigenvektor (\(d \times d\)), bersifat ortonormal
  
• \(\Lambda\): matriks diagonal yang berisi nilai eigen positif (akar variansi)

Dari sini, kita dapat membentuk matriks akar kovarians \(\mathbf{C}\) sebagai:

\[ C = V \Lambda^{1/2} \]

Transformasi ke vektor target dilakukan dengan:

\[ \mathbf{X} = C \mathbf{Z} + \boldsymbol{\mu} = V \Lambda^{1/2} \mathbf{Z} + \boldsymbol{\mu} \]

# Fungsi untuk membangkitkan multivariate normal dengan dekomposisi eigen
rmvn.eigen <- function(n, mu, Sigma) {
  d <- length(mu)  # Dimensi
  ev <- eigen(Sigma, symmetric = TRUE)  # Dekomposisi eigen
  lambda <- ev$values     # Nilai eigen (varians)
  V <- ev$vectors         # Matriks eigenvektor
  R <- V %*% diag(sqrt(lambda)) %*% t(V)  # Matriks akar kovarians

  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)  # Sampel acak standar
  X <- Z %*% R + matrix(mu, n, d, byrow = TRUE)  # Transformasi ke distribusi target
  return(X)
}

# Parameter distribusi
mu <- c(0, 0)  # Vektor rata-rata
Sigma <- matrix(c(1, 0.9,
                  0.9, 1), nrow = 2, byrow = TRUE)  # Matriks kovarians

# Pembangkitan data
set.seed(123)
X <- rmvn.eigen(n = 10000, mu = mu, Sigma = Sigma)

colMeans(X)      # Harus mendekati mu = c(0, 0)

[1] -0.006845915 -0.008975630

Nilai mean dari data hasil simulasi sangat dekat dengan nilai teoritis (0, 0). Penyimpangan kecil ini wajar terjadi karena sifat acak simulasi (sample error).

print(cov(X))           # Harus mendekati Sigma

          [,1]      [,2]
[1,] 1.0043492 0.9062085
[2,] 0.9062085 1.0069044

Kovarians empiris mendekati matriks \(\Sigma\) yang diberikan, dengan korelasi antara dua variabel mendekati 0.9.

print(cor(X))           # Harus mendekati korelasi 0.9

          [,1]      [,2]
[1,] 1.0000000 0.9011387
[2,] 0.9011387 1.0000000

Nilai korelasi empiris antara dua variabel X1 dan X2 adalah 0.915, mendekati nilai target 0.9.

Dekomposisi Cholesky

Dekomposisi Cholesky adalah metode untuk memfaktorkan matriks kovarians yang simetris dan positif-definit \(\boldsymbol{\Sigma}\) menjadi hasil perkalian dua matriks segitiga:

\[ \boldsymbol{\Sigma} = LL^\top \]

di mana:

• \(L\) adalah matriks segitiga bawah berdimensi \(d\times d\),
• \(L^\top\) adalah transpose dari \(L\)
• Seluruh elemen diagonal dari \(L\) bernilai positif.

Agar dekomposisi ini valid, matriks \(\boldsymbol{\Sigma}\) harus memenuhi dua syarat:

• Simetris: \(\boldsymbol{\Sigma\top}=\boldsymbol{\Sigma}\)
• Positif-definit: \(\forall \mathbf{x} \in \mathbb{R}^d \setminus \{0\}, \quad\mathbf{x}^\top \boldsymbol{\Sigma} \mathbf{x} > 0\)

Dalam konteks transformasi distribusi normal multivariat, kita ingin memperoleh akar kovarians \(\boldsymbol{\C}\) sehingga:

\[ \boldsymbol{\Sigma} = \mathbf{C} \mathbf{C}^\top \]

Dengan dekomposisi Cholesky, akar kovarians tersebut diperoleh dari \(\boldsymbol{C} = L\) sehingga \(\boldsymbol{\Sigma} = CC^\top=LL^\top\)

# Fungsi untuk membangkitkan multivariate normal dengan dekomposisi Cholesky
rmvn.chol <- function(n, mu, Sigma) {
  d <- length(mu)
  L <- chol(Sigma)  # Transpos dari hasil chol() agar berbentuk segitiga bawah
  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
  X <- Z %*% L + matrix(mu, n, d, byrow = TRUE)
  return(X)
}

# Parameter distribusi
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.9,
                  0.9, 1), nrow = 2, byrow = TRUE)

# Pembangkitan data
set.seed(123)
X <- rmvn.chol(n = 10000, mu = mu, Sigma = Sigma)

colMeans(X)  # Harus mendekati mu

[1] -0.002371702 -0.006103943

cov(X)       # Harus mendekati Sigma

          [,1]     [,2]
[1,] 0.9972751 0.900173
[2,] 0.9001730 1.003115

cor(X)       # Harus mendekati 0.9

          [,1]      [,2]
[1,] 1.0000000 0.9000014
[2,] 0.9000014 1.0000000

Singular Value Decomposition

Dekomposisi Singular Value (SVD) adalah metode numerik yang sangat stabil dan umum digunakan dalam berbagai aplikasi aljabar linier, termasuk pembangkitan bilangan acak dari distribusi normal multivariat.

Jika \(\boldsymbol{\Sigma}\) adalah matriks kovarians simetris dan semi-definit positif berdimensi \(d \times d\), maka dapat dilakukan dekomposisi SVD sebagai berikut:

\[ \boldsymbol{\Sigma} = UDU^\top \]

dengan:

• \(U\): matriks ortogonal berisi vektor singular kiri (untuk matriks simetris positif-definit, \(U\) sama dengan matriks eigenvektor),
• \(D\): matriks diagonal berisi nilai singular positif (sama dengan nilai eigen jika \(\boldsymbol{\Sigma}\) simetris positif-definit).

Dari dekomposisi tersebut, kita bentuk akar kovarians sebagai:

\[ \mathbf{C} = \mathbf{U} \mathbf{D}^{1/2} \mathbf{U}^\top \] Sehingga transformasi dari vektor acak standar \(Z \sim \mathcal{N}_d(\boldsymbol{\mu}, \boldsymbol{I})\) menjadi:

\[ \mathbf{X} = C \mathbf{Z} + \boldsymbol{\mu} = UD^{1/2}U \mathbf{Z} + \boldsymbol{\mu} \] dengan hasil akhir \(X \sim \mathcal{N}_d(\boldsymbol{\mu}, \boldsymbol{\Sigma})\).

# Fungsi untuk membangkitkan multivariate normal dengan dekomposisi SVD
rmvn.svd <- function(n, mu, Sigma) {
  d <- length(mu)
  sv <- svd(Sigma)  # Dekomposisi SVD
  U <- sv$u
  D <- sv$d
  C <- U %*% diag(sqrt(D)) %*% t(U)  # Matriks akar kovarians simetris

  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)  # Sampel acak standar
  X <- Z %*% C + matrix(mu, n, d, byrow = TRUE)  # Transformasi
  return(X)
}

# Parameter distribusi
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.9,
                  0.9, 1), nrow = 2, byrow = TRUE)

# Pembangkitan data
set.seed(123)
X <- rmvn.svd(n = 10000, mu = mu, Sigma = Sigma)

colMeans(X)  # Harus mendekati mu

[1] -0.006845915 -0.008975630

cov(X)       # Harus mendekati Sigma

          [,1]      [,2]
[1,] 1.0043492 0.9062085
[2,] 0.9062085 1.0069044

cor(X)       # Harus mendekati 0.9

          [,1]      [,2]
[1,] 1.0000000 0.9011387
[2,] 0.9011387 1.0000000

MASS untuk Distribusi Normal Multivariat

Pada R, telah tersedia library MASS yang menyediakan berbagai fungsi statistik klasik, salah satunya adalah mvrnorm() yang digunakan untuk membangkitkan data dari distribusi Normal multivariat. Untuk menggunakan fungsi ini, kita perlu terlebih dahulu memuat library MASS berikut:

library(MASS)

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

Fungsi mvrnorm(n, mu, Sigma) digunakan untuk menghasilkan n observasi dari distribusi \(\mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), dengan:

• n: jumlah observasi yang ingin disimulasikan,
• mu: vektor mean berdimensi \(p\),
• Sigma: matriks kovarians berdimensi \(p \times p\).

set.seed(123)

# Parameter distribusi
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.8,
                  0.8, 1), nrow = 2, byrow = TRUE)

# Simulasi 100 observasi dari distribusi Normal multivariat
sim_data1 <- mvrnorm(n = 100000, mu = mu, Sigma = Sigma)

colMeans(sim_data1)  # Harus mendekati mu

[1] -0.0007226531  0.0025759037

cov(sim_data1)       # Harus mendekati Sigma

          [,1]      [,2]
[1,] 1.0019027 0.7987391
[2,] 0.7987391 0.9986999

cor(sim_data1)       # Harus mendekati 0.9

          [,1]      [,2]
[1,] 1.0000000 0.7984995
[2,] 0.7984995 1.0000000

Parameter tambahan empirical = TRUE memaksa hasil simulasi memiliki rata-rata dan kovarians sampel yang persis sama dengan nilai mu dan Sigma. Hal ini bermanfaat jika ingin memastikan struktur statistik eksak dari data simulasi.

set.seed(123)

# Parameter distribusi
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.8,
                  0.8, 1), nrow = 2, byrow = TRUE)

# Simulasi 100 observasi dari distribusi Normal multivariat
sim_data2 <- mvrnorm(n = 100000, mu = mu, Sigma = Sigma, empirical = TRUE)

colMeans(sim_data2)  # Harus mendekati mu

[1] -2.304122e-17 -2.770374e-17

cov(sim_data2)       # Harus mendekati Sigma

     [,1] [,2]
[1,]  1.0  0.8
[2,]  0.8  1.0

cor(sim_data2)       # Harus mendekati 0.9

     [,1] [,2]
[1,]  1.0  0.8
[2,]  0.8  1.0

Pembangkitan Bilangan Acak untuk Regresi Linier

Dalam analisis regresi, sering kali kita memerlukan data simulasi untuk menguji performa metode estimasi, validasi model, atau tujuan pembelajaran. Data ini dapat dibangkitkan secara acak berdasarkan model regresi linear yang kita tentukan.

Model Regresi Linier Sederhana

Model regresi linear sederhana:

\[ Y=\beta_0+\beta_1X_1+\mu \] dengan:

• \(Y\): peubah respon,
• \(X\): peubah penjelas,
• \(\beta_0, \beta_1\): parameter regresi,
• \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\): error acak berdistribusi Normal.

Semisal, kita ingin membangkitkan data simulasi yang mengikuti model regresi linear dengan:

• Intersep \(\beta_0 = 5\),
• Koefisien slope \(\beta_1 = 2\),
• Error acak mengikuti distribusi Normal dengan rata-rata 0 dan simpangan baku \(\sigma = 1\).

Untuk membangkitkan data simulasi sebanyak \(n = 1000\) observasi, kita dapat memilih nilai prediktor \(X\) dari distribusi tertentu (misalnya seragam/Uniform), lalu menghitung nilai \(Y\) berdasarkan rumus model dengan menambahkan error acak.

set.seed(123)

# Tentukan parameter
n <- 1000             # Jumlah observasi
beta0 <- 5           # Intersep
beta1 <- 2           # Koefisien slope
sigma <- 1           # Simpangan baku error

# Bangkitkan X ~ Uniform(0, 10)
x <- runif(n, min = 0, max = 10)

# Bangkitkan error ~ N(0, sigma^2)
epsilon <- rnorm(n, mean = 0, sd = sigma)

# Hitung Y
y <- beta0 + beta1 * x + epsilon

# Gabungkan ke dalam data frame
data_regresi <- data.frame(x = x, y = y)

Kita dapat memvisualisasikan data hasil simulasi dalam bentuk scatter plot:

ggplot(data_regresi, aes(x = x, y = y)) +
  geom_point(color = "#56B4E9", alpha = 0.7) +
  labs(
    title = "Data Simulasi Regresi Linier Sederhana",
    x = "X",
    y = "Y"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

Setelah memperoleh data simulasi, kita bisa mengestimasi parameter \(\beta_0\) dan \(\beta_1\) (diestimasi dengan \(\hat{\beta}_0\) dan \(\hat{\beta}_1\)) menggunakan fungsi lm():

model_lm <- lm(y ~ x, data = data_regresi)
summary(model_lm)

Call:
lm(formula = y ~ x, data = data_regresi)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8431 -0.7007  0.0136  0.6558  3.4072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.04094    0.06333    79.6   <2e-16 ***
x            1.99417    0.01103   180.9   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.002 on 998 degrees of freedom
Multiple R-squared:  0.9704,    Adjusted R-squared:  0.9704 
F-statistic: 3.271e+04 on 1 and 998 DF,  p-value: < 2.2e-16

Perhatikan bahwa hasil estimasi Intercept dan x dalam output summary() seharusnya mendekati nilai parameter sebenarnya, yaitu \(\beta_0 = 5\) dan \(\beta_1 = 2\), dengan sedikit variasi akibat pengaruh error acak.

# Visualisasi 2: Scatter plot + garis regresi dari hasil lm()
ggplot(data_regresi, aes(x = x, y = y)) +
  geom_point(aes(color = "Data Simulasi"), alpha = 0.7) +
  geom_abline(aes(intercept = coef(model_lm)[1],
                  slope = coef(model_lm)[2],
                  color = "Garis Regresi (Hasil Estimasi)"),
              linewidth = 0.8) +
  scale_color_manual(
    name = "Keterangan",
    values = c("Data Simulasi" = "#56B4E9", "Garis Regresi (Hasil Estimasi)" = "#D55E00")
  ) +
  labs(
    title = "Garis Regresi dari Hasil Estimasi Model Linier",
    x = "X",
    y = "Y"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Terkadang, kita ingin mengetahui bagaimana hasil estimasi parameter regresi berubah-ubah ketika data dibangkitkan ulang berkali-kali dalam kondisi yang sama. Pendekatan ini membantu mengevaluasi rata-rata (bias) dan keragaman (variabilitas) dari estimator parameter.

Misalnya, kita ingin mengevaluasi kestabilan hasil estimasi dari model regresi sederhana melalui simulasi sebanyak 100 kali, dengan skenario sebagai berikut:

• Ukuran sampel kecil (\(n = 10\)),
• Intersep sebenarnya \(\beta_0 = 5\), slope \(\beta_1 = 2\),
• Error acak mengikuti distribusi Normal: \(\varepsilon \sim \mathcal{N}(0, 1)\),
• Prediktor \(X\) diambil dari distribusi Uniform(5, 10).

set.seed(123)

# Parameter simulasi
n_simulasi <- 100
n_obs <- 10
beta0 <- 5
beta1 <- 2
sigma <- 1

# Vektor untuk menyimpan hasil estimasi
b0_hat <- numeric(n_simulasi)
b1_hat <- numeric(n_simulasi)

# Simulasi berulang 100 kali
for (i in 1:n_simulasi) {
  x <- runif(n_obs, min = 5, max = 10)
  eps <- rnorm(n_obs, mean = 0, sd = sigma)
  y <- beta0 + beta1 * x + eps
  
  model <- lm(y ~ x)
  b0_hat[i] <- coef(model)[1]
  b1_hat[i] <- coef(model)[2]
}

# Ringkasan hasil
hasil_simulasi <- data.frame(
  Estimator = c("Intercept (b0)", "Slope (b1)"),
  Mean = c(mean(b0_hat), mean(b1_hat)),
  SD = c(sd(b0_hat), sd(b1_hat))
)

hasil_simulasi

Setelah melakukan simulasi berulang sebanyak 100 kali, kita memperoleh 100 nilai estimasi untuk masing-masing parameter regresi, yaitu intersep (\(\beta_0\)) dan slope (\(\beta_1\)). Nilai-nilai estimasi ini dapat bervariasi karena adanya komponen error acak (\(\varepsilon\)) dalam setiap replikasi data.

Untuk memahami sebaran, bias, dan stabilitas estimasi parameter yang dihasilkan dari proses simulasi tersebut, kita dapat menggunakan boxplot. Visualisasi ini memberikan gambaran ringkas mengenai:

• Letak pusat estimasi (median dan mean secara kasar),
• Sebaran estimasi (varians),
• Potensi outlier,
• Apakah hasil estimasi cenderung bias terhadap parameter sebenarnya.

Dengan boxplot, kita bisa dengan cepat melihat apakah estimasi parameter regresi konsisten dan mendekati nilai sebenarnya (\(\beta_0 = 5\), \(\beta_1 = 2\)) dalam konteks simulasi.

# Buat data boxplot
df_boxplot <- data.frame(
  Estimator = factor(rep(c("Intercept (b0)", "Slope (b1)"), each = n_simulasi)),
  Nilai = c(b0_hat, b1_hat)
)

ggplot(df_boxplot, aes(x = Estimator, y = Nilai, fill = Estimator)) +
  geom_boxplot(alpha = 0.8, width = 0.5, outlier.shape = 16, outlier.size = 1.5) +
  
  # Garis horizontal untuk nilai sebenarnya
  geom_hline(aes(yintercept = 5, color = "beta0"),
             linetype = "dashed", linewidth = 0.8) +
  geom_hline(aes(yintercept = 2, color = "beta1"),
             linetype = "dashed", linewidth = 0.8) +
  
  scale_color_manual(
    name = NULL,
    values = c("beta0" = "#FF6666", "beta1" = "#66CC66"),
    labels = c(
      "beta0" = expression("True value of " * beta[0]),
      "beta1" = expression("True value of " * beta[1])
    )
  ) +
  
  # Warna boxplot tanpa legend
  scale_fill_manual(
    values = c("Intercept (b0)" = "#56B4E9", "Slope (b1)" = "#D55E00"),
    guide = "none"  # Hapus legend fill
  ) +
  
  labs(
    title = "Variabilitas Estimator Parameter Regresi dari 100 Simulasi",
    x = NULL,
    y = "Nilai Estimasi"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Model Regresi Linier Berganda

Model regresi linier berganda merupakan perluasan dari regresi linier sederhana dengan melibatkan lebih dari satu peubah penjelas:

\[ Y=\beta_0+\beta_1X_1+\beta_2X_2+\dots+\beta_pX_p+\mu \]

dengan:

• \(Y\): peubah respon,
• \(X_1, X_2, \dots, X_p\): peubah penjelas,
• \(\beta_0, \beta_1, \dots, \beta_p\): parameter model,
• \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\): error acak berdistribusi Normal.

Kita akan mensimulasikan data sebanyak \(n = 1000\) observasi dengan

• \(\beta_0 = 3\), \(\beta_1 = 1.5\), \(\beta_2 = -2\),
• Peubah penjelas \(X_1 \sim \text{Uniform}(0, 10)\) dan \(X_2 \sim \text{Normal}(5, 2)\),
• Error acak \(\varepsilon \sim \mathcal{N}(0, 1)\).

set.seed(123)

# Parameter
n <- 1000
beta0 <- 3
beta1 <- 1.5
beta2 <- -2
sigma <- 1

# Bangkitkan peubah penjelas
x1 <- runif(n, min = 0, max = 10)      # X1 ~ Uniform
x2 <- rnorm(n, mean = 5, sd = 2)       # X2 ~ Normal

# Error
epsilon <- rnorm(n, mean = 0, sd = sigma)

# Respon Y
y <- beta0 + beta1 * x1 + beta2 * x2 + epsilon

# Gabungkan ke dalam data frame
data_regresi_ganda <- data.frame(x1 = x1, x2 = x2, y = y)

# Estimasi model
model_ganda <- lm(y ~ x1 + x2, data = data_regresi_ganda)

# Output model
summary(model_ganda)

Call:
lm(formula = y ~ x1 + x2, data = data_regresi_ganda)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9879 -0.6392 -0.0180  0.6932  3.2486 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.92260    0.10044    29.1   <2e-16 ***
x1           1.48509    0.01082   137.3   <2e-16 ***
x2          -1.97037    0.01553  -126.9   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.983 on 997 degrees of freedom
Multiple R-squared:  0.9727,    Adjusted R-squared:  0.9727 
F-statistic: 1.776e+04 on 2 and 997 DF,  p-value: < 2.2e-16

Seperti sebelumnya, kita ingin mengetahui bagaimana estimasi \(\hat{\beta}_0\), \(\hat{\beta}_1\), dan \(\hat{\beta}_2\) berubah saat model dibangkitkan ulang sebanyak 100 kali dengan kondisi yang sama. Pendekatan ini membantu mengevaluasi rata-rata, sebaran, dan potensi bias dari setiap estimator dalam regresi linier berganda.

set.seed(123)

# Parameter
n <- 10
beta0 <- 3
beta1 <- 1.5
beta2 <- -2
sigma <- 1

n_simulasi <- 100
b0_hat <- numeric(n_simulasi)
b1_hat <- numeric(n_simulasi)
b2_hat <- numeric(n_simulasi)

for (i in 1:n_simulasi) {
  x1 <- runif(n, 0, 10)
  x2 <- rnorm(n, 5, 2)
  eps <- rnorm(n, 0, sigma)
  y <- beta0 + beta1 * x1 + beta2 * x2 + eps

  model <- lm(y ~ x1 + x2)
  b0_hat[i] <- coef(model)[1]
  b1_hat[i] <- coef(model)[2]
  b2_hat[i] <- coef(model)[3]
}

# Ringkasan hasil
hasil_regresi_ganda <- data.frame(
  Estimator = c("Intercept (b0)", "Slope X1 (b1)", "Slope X2 (b2)"),
  Mean = c(mean(b0_hat), mean(b1_hat), mean(b2_hat)),
  SD = c(sd(b0_hat), sd(b1_hat), sd(b2_hat))
)

hasil_regresi_ganda

df_box_ganda <- data.frame(
  Estimator = factor(rep(c("Intercept (b0)", "Slope X1 (b1)", "Slope X2 (b2)"), each = n_simulasi)),
  Nilai = c(b0_hat, b1_hat, b2_hat)
)

ggplot(df_box_ganda, aes(x = Estimator, y = Nilai, fill = Estimator)) +
  geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 1.5) +
  
  geom_hline(aes(yintercept = 3, color = "b0"), linetype = "dashed", linewidth = 0.8) +
  geom_hline(aes(yintercept = 1.5, color = "b1"), linetype = "dashed", linewidth = 0.8) +
  geom_hline(aes(yintercept = -2, color = "b2"), linetype = "dashed", linewidth = 0.8) +
  
  scale_color_manual(
    name = NULL,
    values = c(
      "b0" = "#009E73",
      "b1" = "#D55E00",
      "b2" = "#0072B2"
    ),
    labels = c(
      "b0" = expression("True value of " * beta[0]),
      "b1" = expression("True value of " * beta[1]),
      "b2" = expression("True value of " * beta[2])
    )
  ) +
  
  # Warna boxplot (tanpa legend)
  scale_fill_manual(
    values = c("#E69F00", "#56B4E9", "#009E73"),
    guide = "none"
  ) +
  
  labs(
    title = "Variabilitas Estimator Parameter pada Regresi Linier Berganda",
    x = NULL,
    y = "Nilai Estimasi"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Pada contoh sebelumnya, peubah penjelas \(X_1\) dan \(X_2\) dibangkitkan secara independen. Dalam kenyataannya, peubah penjelas sering kali saling berkorelasi, misalnya pendapatan dan pengeluaran rumah tangga. Untuk merepresentasikan kondisi tersebut, kita dapat membangkitkan \(X_1\) dan \(X_2\) dari distribusi Normal multivariat dengan kovarians tidak nol.

Model yang digunakan:

\[ Y=\beta_0+\beta_1X_1+\beta_2X_2+\mu \]

dengan:

• \(\beta_0 = 3, \beta_1 = 1.5, \beta_2 = -2\),
• \((X_1, X_2) \sim \mathcal{N}(\mu, \Sigma)\) dengan

\[ \mu = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \quad \Sigma = \begin{bmatrix} 1 & 0.7 \\ 0.7 & 1 \end{bmatrix} \]

• \(\varepsilon \sim \mathcal{N}(0,1)\).

set.seed(123)
library(MASS)

# Parameter
n <- 100
beta0 <- 3
beta1 <- 1.5
beta2 <- -2
sigma <- 1

# Parameter distribusi multivariat normal untuk X1 dan X2
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.7,
                  0.7, 1), nrow = 2, byrow = TRUE)

# Bangkitkan X1 dan X2 yang berkorelasi
X <- mvrnorm(n = n, mu = mu, Sigma = Sigma)
x1 <- X[,1]
x2 <- X[,2]

# Bangkitkan error
epsilon <- rnorm(n, mean = 0, sd = sigma)

# Hitung Y
y <- beta0 + beta1 * x1 + beta2 * x2 + epsilon

# Gabungkan ke data frame
data_corr <- data.frame(x1 = x1, x2 = x2, y = y)

# Estimasi model regresi
model_corr <- lm(y ~ x1 + x2, data = data_corr)
summary(model_corr)

Call:
lm(formula = y ~ x1 + x2, data = data_corr)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8730 -0.6607 -0.1245  0.6214  2.0798 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.13507    0.09614   32.61   <2e-16 ***
x1           1.39704    0.13729   10.18   <2e-16 ***
x2          -2.04148    0.14244  -14.33   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9513 on 97 degrees of freedom
Multiple R-squared:  0.6799,    Adjusted R-squared:  0.6733 
F-statistic:   103 on 2 and 97 DF,  p-value: < 2.2e-16

Kita bisa ulangi proses simulasi sebanyak 100 kali untuk melihat stabilitas estimasi parameter dalam kondisi multikolinearitas moderat.

set.seed(123)

n_simulasi <- 100
b0_hat <- numeric(n_simulasi)
b1_hat <- numeric(n_simulasi)
b2_hat <- numeric(n_simulasi)

n <- 10
beta0 <- 3
beta1 <- 1.5
beta2 <- -2
sigma <- 1

# Parameter distribusi multivariat normal untuk X1 dan X2
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.7,
                  0.7, 1), nrow = 2, byrow = TRUE)

for (i in 1:n_simulasi) {
  X <- mvrnorm(n = n, mu = mu, Sigma = Sigma)
  x1 <- X[,1]
  x2 <- X[,2]
  eps <- rnorm(n, mean = 0, sd = sigma)
  y <- beta0 + beta1*x1 + beta2*x2 + eps
  
  model <- lm(y ~ x1 + x2)
  b0_hat[i] <- coef(model)[1]
  b1_hat[i] <- coef(model)[2]
  b2_hat[i] <- coef(model)[3]
}

# Ringkasan hasil estimasi
hasil_corr <- data.frame(
  Estimator = c("Intercept (b0)", "Slope X1 (b1)", "Slope X2 (b2)"),
  Mean = c(mean(b0_hat), mean(b1_hat), mean(b2_hat)),
  SD = c(sd(b0_hat), sd(b1_hat), sd(b2_hat))
)

hasil_corr

df_corr <- data.frame(
  Estimator = factor(rep(c("Intercept (b0)", "Slope X1 (b1)", "Slope X2 (b2)"), each = n_simulasi)),
  Nilai = c(b0_hat, b1_hat, b2_hat)
)

ggplot(df_corr, aes(x = Estimator, y = Nilai, fill = Estimator)) +
  geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 1.5) +
  geom_hline(aes(yintercept = 3, color = "b0"), linetype = "dashed", linewidth = 0.8) +
  geom_hline(aes(yintercept = 1.5, color = "b1"), linetype = "dashed", linewidth = 0.8) +
  geom_hline(aes(yintercept = -2, color = "b2"), linetype = "dashed", linewidth = 0.8) +
  scale_color_manual(
    name = NULL,
    values = c("b0" = "#009E73", "b1" = "#D55E00", "b2" = "#0072B2"),
    labels = c(
      "b0" = expression("True value of " * beta[0]),
      "b1" = expression("True value of " * beta[1]),
      "b2" = expression("True value of " * beta[2])
    )
  ) +
  scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73"), guide = "none") +
  labs(
    title = "Variabilitas Estimator dengan Prediktor Berkorelasi",
    x = NULL,
    y = "Nilai Estimasi"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Muhammad Yusran
Program Studi Magister Statistika dan Sains Data
Sekolah Sains Data, Matematika, dan Informatika
IPB University

Email: muhammadyusran@apps.ipb.ac.id
LinkedIn: Muhammad Yusran
Instagram: mhammad.yusran_