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")