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