1 Tugas Optimasi

1.1 Soal Transformasi inverse

Membangkitkan 1000 bilangan acak dari distribusi Exponential(3) menggunakan Inverse Transform Method dan membandingkannya dengan fungsi bawaan rexp() menggunakan histogram.

1.2 Teori

Distribusi exponential dengan laju \(\lambda\) memiliki fungsi densitas dan fungsi distribusi sebagai berikut:

\[ f(x) = \lambda e^{-\lambda x}, \quad x \ge 0 \]

\[ F(x) = 1 - e^{-\lambda x}, \quad x \ge 0 \]

Untuk membangkitkan bilangan acak dari distribusi exponential menggunakan Inverse Transform Method, kita gunakan hubungan: \[ F^{-1}(u) = -\frac{1}{\lambda}\ln(1-u) \] dengan \(u \sim U(0,1)\).

1.3 Jawaban 1

  1. Tentukan CDF \(F(x)\)

    PDF dari \(exponential(3)\): \[ f(x) = 3e^{-3x}, \quad x \ge 0 \]

    Maka CDF-nya adalah: \[ F(x) = 1 - e^{-3x}, \quad x \ge 0 \]

  2. Tentukan inverse CDF \(F^{-1}(u)\)

    Inverse CDF dari \(exponential(3)\): \[ F^{-1}(u) = -\frac{\ln(1-u)}{3}, \quad x \ge 0 \]

Bangkitkan u∼U(0,1)

# banyaknya amatan
n <- 1000
u <- runif(n)
set.seed(10) # membuat seed
x <- -(log(1-u)/3) #log merupakan ln di R
head(x)
## [1] 0.05771407 0.57252996 0.22396868 0.21951921 0.19273004 0.26469711
set.seed(10) # membuat seed
y <- rexp(n,rate = 3)
head(y)
## [1] 0.004985469 0.306740402 0.250719646 0.525013950 0.077219539 0.362224334
# 4. Tampilkan dua grafik berdampingan
par(mfrow = c(1, 2))  # 1 baris, 2 kolom

# Grafik 1: Inverse Transform
hist(x, breaks = 20, probability = TRUE, 
     main = "Exp : Metode Invers-Transform",
     xlab = "exp1", col = "lightcoral", border = "white")

# Grafik 2: Fungsi rexp()
hist(y, breaks = 20, probability = TRUE, 
     main = "Exp : Fungsi Rexp",
     xlab = "exp2", col = "lightcoral", border = "white")

1.4 Soal Acceptance-Rejection

Membangkitkan 1000 bilangan acak dari distribusi Beta(2,2) menggunakan Metode Acceptance–Rejection dan membandingkannya dengan fungsi bawaan rbeta() melalui histogram.

##Teori Distribusi Beta(α, β) dengan parameter \(\alpha = 2\) dan \(\beta = 2\) memiliki fungsi kepadatan probabilitas (PDF):

\[ f(x) = 6x(1-x), \quad 0 < x < 1 \]

Kita menggunakan distribusi Uniform(0,1) sebagai fungsi proposal:

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

Rasio antara fungsi target dan proposal:

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

Nilai maksimumnya: \[ M = \max_{x \in (0,1)} 6x(1-x) = 1.5 \]

Syarat acceptance–rejection:

\[ u \leq \frac{f(x)}{M g(x)} = \frac{6x(1-x)}{1.5} = 4x(1-x) \]

dengan \(u \sim U(0,1)\).

1.5 Langkah 1: Menetapkan Peubah Acak dan Fungsi PDF-nya

Tetapkan peubah acak dari \(Y\) sebagai target sebaran beserta PDF-nya \(g(y)\).

Dalam kasus ini, kita pilih: \[ Y \sim U(0,1) \] karena memiliki domain yang sama dengan peubah acak \(X\), yaitu \(0 < x < 1\).

PDF dari \(Y\) adalah: \[ g(y) = 1, \quad 0 \le y \le 1 \]

1.6 Langkah 2: Membangkitan Peubah Acak \(Y\)

n=1000
set.seed(10)
y = runif(n,0,1)
set.seed(30)
u = runif(n,0,1)

Hitung nilai \(c\) dengan mencari nilai maksimum dari rasio: \[ \frac{f(y)}{g(y)} \]

Misalkan: \[ h(y) = \frac{f(y)}{g(y)} \]

Untuk distribusi \(X \sim \text{Beta}(2,2)\) dengan \[ f(y) = 6y(1 - y), \quad 0 < y < 1 \] dan karena \(Y \sim U(0,1)\) maka \(g(y) = 1\), diperoleh: \[ h(y) = 6y(1 - y) \]

Untuk mencari nilai maksimum dari \(h(y)\), kita cari titik kritis dengan menurunkan \(h(y)\) terhadap \(y\): \[ \frac{dh(y)}{dy} = 0 \]

Turunan dari \(h(y) = 6y(1 - y)\) adalah: \[ \frac{dh(y)}{dy} = 6(1 - y) - 6y = 6 - 12y \]

Dengan menyetarakan turunan ke nol: \[ 6 - 12y = 0 \Rightarrow y = \frac{1}{2} \]

Substitusi nilai \(y = \frac{1}{2}\) ke dalam fungsi \(h(y)\): \[ h\left(\frac{1}{2}\right) = 6\left(\frac{1}{2}\right)\left(1 - \frac{1}{2}\right) = 6 \times \frac{1}{2} \times \frac{1}{2} = \frac{3}{2} \]

Sehingga diperoleh: \[ c = 1.5 \]

h <- function(y) 6*y*(1-y)
derivH <- optimize(f=h,interval=c(0,1),maximum=T)
derivH$objective
## [1] 1.5
f <- function(y) 6*y*(1-y)
g <- function(y) 1
nilaic <- derivH$objective
#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
#jadikan Y sebagai X
x <- y[criteria]
head(x)
## [1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597 0.22543662
z <- rbeta(n,shape1 = 2,shape2 = 2)
par(mfrow=c(1,2)) 
hist(x,main="beta(2,2) dari AR")
hist(z,main="beta(2,2) dari rbeta")