Pembangkitan Bilangan Acak

Author

Muhammad Syafiq

Rangkuman Materi

Random Generators in R

Sebagai salah satu bahasa pemrograman yang hadir karena kebutuhan akan komputasi statistik, R telah menyiapkan banyak fungsi untuk membangkitkan data berdasarkan sebaran. Fungsi peluang suatu sebaran dalam R ditandai dengan prefix sebagai berikut:

Inverse Transform Method

Algoritme Inverse Transform Method adalah sebagai berikut:

  1. Tentukan CDF \(F(x)\) dari distribusi yang ingin dibangkitkan
  2. Carilah inverse CDF (quantile) \(F^{-1}(x)\)
  3. Bangkitkan \(u\) berdasarkan distribusi \(Uniform(0,1)\) atau \(u \sim U(0,1)\)
  4. Dapatkan bilangan acak \(x\) dengan menghitung \(F^{-1}(u)\)

Contoh

Bangkitkan bilangan acak yang berdistribusi \(exponential(3)\) dengan Inverse Transform Method, yang amatannya berjumlah \(1000\). Bandingkan hasilnya dengan fungsi bawaan R rexp dengan menggunakan histogram.

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

PDF dari \(\text{Exponential}(3)\)

\[ f(x) = 3e^{-3x}, x \geq 0 \]

CDF dari \(\text{Exponential}(3)\)

\[ F(x) = 1-e^{-3x}, x \geq 0 \]

2.Tentukan inverse CDF \(F(x)\)

inverse CDF dari \(\text{Exponential}(3)\)

\[ F^{-1}(u) = -\frac{ln(1-u)}{3}, x \geq 0 \]

3.Bangkitkan \(u \sim U(0,1)\)

# banyaknya amatan
n <- 1000
# membangkitan U(0,1)
u <- runif(n)

4.Dapatkan bilangan acak \(x\) dengan menghitung \(F^{-1}(u)\)

# fungsi inverse CDF
set.seed(10) # membuat seed
x <- -(log(1-u)/3) #log merupakan ln di R
head(x)
[1] 0.032350283 0.033035477 0.063608827 0.006058708 0.270181804 0.364711443
# memeriksa banyaknya amatan
length(x)
[1] 1000

Membangkitan dengan fungsi bawaan R

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
# memeriksa banyaknya amatan
length(y)
[1] 1000

Menamplikan histogram di R

par(mfrow=c(1,2)) 
hist(x,main="Exp dari Inverse Transform")
hist(y,main="Exp dari fungsi rexp")

fungsi par digunakan untuk menampilkan beberapa plot dalam satu layout, argumen mfrow digunakan untuk menyatakan bagaimana pembagian layout tersebut.

Argumen main dalam fungsi hist digunakan untuk memberi judul pada plot.

Acceptance-Rejection Method

Algoritme Acceptance-Rejection Method adalah sebagai berikut:

  1. Tetapkan peubah acak dari \(Y\) sebagai target sebaran beserta PDFnya \(g(y)\)
  2. Bangkitkan peubah acak \(Y\) tersebut.
  3. Bangkitkan \(u\) berdasarkan distribusi \(Uniform(0,1)\) atau \(u ~ U(0,1)\)
  4. Hitung nilai \(c\) dengan mencari nilai maksimum dari \(\frac{f(y)}{g(y)}\)
  5. Jika \(u\leq \frac{f(y)}{c g(y)}\), jadikan \(Y\) sebagai \(X\).
  6. Jika kondisi ke-4 tidak terpenuhi kembali ke-1.

Contoh

Bangkitkan peubah acak \(X \sim Beta(2,2)\) yang memiliki fungsi PDF

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

dengan menggunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histrogramnya!

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)
head(u)
[1] 0.09878282 0.48823179 0.36403673 0.42061913 0.30096439 0.14763513
length(u)
[1] 1000

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:

\[ \begin{aligned} 6(1-y)-6y &= 0 \\ 6-6y-6y &= 0 \\ 6-12y &= 0 \\ y &=\frac{1}{2} \\ \end{aligned} \]

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)\).

5.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.

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

Direct Transformation

y <- rnorm(1000)
x <- y^2
hist(x, prob=T, col="lightgreen")
sbx <- seq(0,13,0.01)
lines(sbx, dchisq(sbx,df=1), col="red")

Referensi: rpubs.com/alfanugraha/random-var