8 Jawaban 2

Untuk menjawab soal Latihan 3, 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 U(0,1)$ karena memiliki domain yang sama dengan peubah acak $X$ yaitu $0<x<1$. Pdf dari Y adalah

$g(y) = 1, 0\leq y \leq 1$

2. Bangkitkan peubah acak $Y$ tersebut.

# banyaknya amatan
n=1000
set.seed(10)
y = runif(n,0,1)

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)}$

Misalkan saja $h(y) = \frac{f(y)}{g(y)}$, sehingga

$h(y) = \frac{6y(1-y)}{1} = 6y(1-y)$ Untuk mencari nilai maksimum dari $h(y)$ maka kita harus menyelesaikan persamaan

$\frac{d h(y)}{dy}=0$

Dengan menggunakan aturan perkalian pada turunan maka didapatkan:

$6(1-y)-6y=0$

$6-6y-6y=0$

$6-12y=0$

$y=\frac{1}{2}$

Subtitusi $y=\frac{1}{2}$ ke $h(y)$ sehingga didapatkan:

$h(\frac{1}{2}) = 6(\frac{1}{2})(1-\frac{1}{2})=\frac{3}{2}$

sehingga nilai $c=1.5$

Di R terdapat fungsi untuk mendapatkan nilai $c$ yaitu dengan menggunakan fungsi optimize

h <- function(y) 6*y*(1-y)
derivH <- optimize(f=h,interval=c(0,1),maximum=T)
derivH$objective

## [1] 1.5

argumen interval adalah domain dari $f(y)$ dan $g(y)$, maximum=TRUE menujukkan bahwa fungsi optimze mencari nilai maksimum dari $h(y)$.

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

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

# banyaknya amatan
length(x)

## [1] 668

Jika diperhatikan banyaknya amatan pada $x$ tidak berjumlah 1000, dikarenakan ada beberapa amatan dalam $y$ yang tidak diloloskan oleh kriteria penerimaan. Oleh karena itu, jika jumlah amatan yang harus dibangkitkan sebanyak 1000 maka banyaknya amatan awal yang ditentukan harus lebih dari 1000. Selamat mencoba!!!

Note: Program di atas adalah alternatif lain dari program yang sudah ada di kuliah. Silahkan gunakan program yang anda mengerti

• Menamplikan histogram di R

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