10 Jawaban 3

Untuk menjawab soal Latihan 4, kita akan terapkan algoritme dari Acceptance-Rejection Method sebagai berikut:

1. Tetapkan peubah acak dari \(Y\) sebagai target sebaran beserta PDFnya \(g(y)\)

Dalam Kasus ini kita pilih \(Y \sim exponential(\frac{3}{2})\) karena memiliki domain yang sama dengan peubah acak \(X\) yaitu \(0<x< \infty\) dan karena distribusi exponensial merupakan kasus khusus dari distribus gamma. Pdf dari Y adalah

\(g(y) = \frac{1}{\frac{3}{2}} e^{-\frac{y}{\frac{3}{2}}}, 0\leq y \leq \infty\)

atau

\(g(y) = \frac{2}{3} e^{-\frac{2}{3}y}, 0\leq y \leq \infty\) 2. Bangkitkan peubah acak \(Y\) tersebut.

# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)

Note : Distribusi \(exponential(\theta)\) di R memiliki pdf

\(f(x) = \theta e^{-\theta y}, 0\leq y \leq \infty\)

sedangkan yang kita gunakan \(exponential(\theta)\) memiliki pdf

\(f(x) = \frac{1}{\theta} e^{-\frac{y}{\theta}}, 0\leq y \leq \infty\)

sehingga \(exponential(\frac{3}{2})\) yang dimaksud sama dengan \(exponential(\frac{2}{3})\) di R.

3. Bangkitkan \(u\) berdasarkan distribusi \(Uniform(0,1)\) atau \(u ~ U(0,1)\)

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

4. Hitung nilai \(c\) dengan mencari nilai maksimum dari \(\frac{f(y)}{g(y)}\)

Untuk menghitung manual nilai \(c\) kita lewati agar bisa dikerjakan secara mandiri. Kita langsung saja menggunakan fungsi optim di R. Fungsi optimize tidak bisa digunakan disini karena domain dari pdf \(f(y)\) dan \(g(y)\) terdapat unsur ketakhinggaan \((\infty)\).

fungsi \(h(x)\) yang diperoleh adalah

\(h(x) = \frac{3}{2 \Gamma \left(\frac{3}{2}\right)} y^{\frac{1}{2}} e^{-\frac{y}{3}}\)

h <- function(y) (-1) * (3/(2*gamma(3/2))) * (y)^(1/2) * exp(-y/3)
derivH <- optim(par = 0,fn=h,lower = 0,method = "L-BFGS-B")
derivH$value

## [1] -1.257317

Beberapa catatan yang perlu diperhatikan untuk menggunakan fungsi optim adalah secara default fungsi ini hanya bisa mencari nilai minimum dari suatu fungsi. Oleh karena itu, fungsi \(h\) perlu dikalikan dengan \(-1\). Argumen par digunakan untuk menentukan nilai awal.

4. Jika \(u\leq \frac{f(y)}{c g(y)}\), jadikan \(Y\) sebagai \(X\).

f <- function(y) (3/(2*gamma(3/2))) * (y)^(1/2) * exp(-y)
g <- function(y) (2/3) * exp(-(2/3) * y)
nilaic <- -derivH$value #jangan lupa tambahkan negatif
#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.02243461 1.38033181 1.12823841 2.36256278 0.34748792 1.63000950

# banyaknya amatan
length(x)

## [1] 1393

#karena >1000 ambil 1000 saja
x <- x[1:1000]

• Menamplikan histogram di R

z <- rgamma(n,shape = 3/2,scale = 1)
par(mfrow=c(1,2)) 
hist(x,main="gamma(3/2,1) dari AR")
hist(z,main="beta(3/2,1) dari rgamma")