4 Pembangkitan Bilangan Acak

4.1 Teknik Pembangkitan Bilangan Acak

Secara umum terdapat tiga metode untuk membangkitkan bilangan acak, yaitu:

  1. Inverse Transform Method

  2. Acceptance-Rejection Method

  3. Direct Transformation

4.2 Algoritme Inverse Transform Method

Algoritme Inverse Transform Method adalah sebagai berikut:

  1. Tentukan CDF dari distribusi yang ingin dibangkitkan

Latihan 1

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.

Jawaban 1

  1. Tentukan CDF \(F(x)\) PDF dari \(exponential(3)\) \(f(x)=3e^{-3x},x \geq 0\) CDF dari \(exponential(3)\) \(f(x)=1-e^{-3x},x \geq 0\)

  2. Tentukan inverse CDF \(F(x)\) inverse CDF dari \(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)
  1. 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.11591076 1.10580059 0.63204987 0.89916898 0.45320418 0.07437934
# 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.

Algoritme 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)}{cg(y)}\)

  6. Jika kondisi ke-4 tidak terpenuhi kembali ke-1.

Latihan 2

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 histogramnya!

Jawaban 2

Untuk menjawab soal latihan 2, Kita akan terapkan algoritma dari Acceptance-Rejection Method sebagai berikut:

  1. Tetapkan peubah acak dari \(Y\) sebagai sebaran dari PDFnya \(g(y)\) Dalam kasus ini kita pilih \(Y \sim U(0,1)\) karena memiliki domainyang 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)
  1. Bangkitkan \(u\) berdasarkan distribusi \(Uniform(0,1)\) atau \(u U(0,1)\)
set.seed(30)
u = runif(n,0,1)
  1. 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{dh(u)}{dy}=0\)

Dengan menggunakan aturan perkalian pada turunan maka didapatkan:

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

\(6-6y-6y=0\)

\(6-12y=0\)

\(y={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\).

  1. Jika \(u \leq \frac{f(y)}{cg(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")

Latihan 3

Bangkitkan peubah acak \(X \sim Gamma(\frac{3}{2},1)\) yang memiliki fungsi PDF

\(f(x)= \frac{1}{\Gamma \left(\frac{3}{2} \right)} x^{\frac{1}{2}} e^{-x}, 0<x< \infty\)

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

Jawaban 3

Untuk menjawab soal latihan 3, kita akan terapkan algoritma 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 khusu dari distribusi gamma PDF dari Y adalah

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

  1. Bangkitkan peubah acak \(Y\) tersebut
# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)

Note : Distribusi\(exponential(θ)\) di R memiliki PDF \(f(x)=θe^{-θy},0\leq y \leq \infty\) sedangkan yang kita gunakan \(exponential(θ)\) memiliki PDF \(f(x)=\frac{1}{θ}e^-\frac{y}{0},0\leq y \leq \infty\)

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

  1. Bangkitkan \(u\) berdasarkan distribusi \(Uniform(0,1)\) atau \(u U(0,1)\)
set.seed(30)
u = runif(n,0,1)
  1. 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 (\frac{3}{2})}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.

  1. Jika \(u \leq \frac{f(y)}{cg(y)}\), jadika \(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")

Latihan 4

Bangkitkan peubah acak \(X \sim t(2)\) sebanyak 1000 amatan yang didapatkan dari hasil transformasi langsung peubah acak \(Y \sim X^2 (2)\) dan \(Z \sim N(0,1)\).

Buatlah histogramnya dan bandingkan hasilnya dengan hasik bangkitan yang diperoleh dari fungsi rt

Jawaban 4

Distribusi \(t\) bisa didapatkan dari hasil transformasi distribusi \(\chi^2\) dan \(N(0,1)\)melalui persamaan dibawah ini.

\(X = \frac{Z}{\sqrt{\frac{Y}{2}}}\)

Sehingga berdasarkan persamaan tersebut kita bisa menuliskan program R untuk membangkitkan distribusi \(t(2)\) sebagai berikut

#banyaknya amatan
n <- 1000
set.seed(15)
z <- rnorm(n,mean = 0,sd = 1)
y <- rchisq(n,df = 2)

x <- z/sqrt(y/2)
head(x)
## [1]  0.4249926  2.2981978 -0.3956900  0.8938340  0.9930945 -1.4325087

Menampilkan histogram di R

set.seed(5)
v <- rt(n,df=2)
par(mfrow=c(1,2)) 
hist(x,main="t(2) dari Direct")
hist(v,main="t(2) dari rt")

Algoritme pembangkitan bilangan acak untuk Regresi Linear

Algoritme pembangkitan bilangan acak untuk Regresi Linear adalah sebagai berikut:

  1. Bangkitkan residual berdasarkan distribusi normal \(\epsilon \sim N(0,\sigma)\)
  2. Bangkitkan peubah penjelas \(X_i\) yang saling bebas.
  3. Bangkitkan peubah respon \(Y\) dengan menggunakan model regresi

Latihan 5

Bangkitkan peubah \(Y,X_1,X_2\) dan \(X_3\) berdasarkan model regresi linear berganda berikut ini: \(Y = 10+3X_1+5X_2+7X_3+ \epsilon\) dengan mengasumsikan bahwa \(\epsilon N(0,1)\). Banyaknya amatan yang dibangkitkan adalah 1000

Jawaban 5

  1. Bangkitkan residual berdasarkan distribusi normal \(\epsilon \sim N(0,\sigma)\)
#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 1)
  1. Bangkitkan peubah penjelas \(X_i\) yang saling bebas.
set.seed(10)
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
head(X)
##             [,1]       [,2]        [,3]
## [1,]  0.02811926  1.5750205 -0.46179747
## [2,] -0.27637881  0.4291389  1.13712837
## [3,] -2.05699582  0.3608472 -0.86079512
## [4,] -0.89875157  1.2490578 -1.40811669
## [5,]  0.44181769 -0.3344748 -0.04154899
## [6,]  0.58469145  0.4325163 -1.59937298

Pembangkitan peubah penjelas \(X_i\) dapat dilakukan juga dengan fungsi mvnorm seperti pada slide kuliah.

  1. Bangkitkan peubah respon \(Y\) dengan menggunakan model regresi
betas <- matrix(c(10,3,5,7),nrow = 4,ncol = 1)
Xgab <- cbind(1,X)

Y <- Xgab%*%betas+epsilon

dataRegresi <- data.frame(Y,X)
colnames(dataRegresi) <- c("Y","X1","X2","X3")
head(dataRegresi)
##           Y          X1         X2          X3
## 1 13.246310  0.02811926  1.5750205 -0.46179747
## 2 20.853626 -0.27637881  0.4291389  1.13712837
## 3 -1.349062 -2.05699582  0.3608472 -0.86079512
## 4  2.772212 -0.89875157  1.2490578 -1.40811669
## 5  7.364594  0.44181769 -0.3344748 -0.04154899
## 6  2.448749  0.58469145  0.4325163 -1.59937298

Note : sintaks Xgab <- cbind(1,X) digunakan untuk menggabungkan peubah penjelas yang ada di matrik \(X\)dengan vektor yang isinya 1. Vektor ini digunakan sebagai pengkali konstanta model (dalam hal ini bernilai 10).

Untuk membuktikan bahwa bangkitan kita benar, maka kita akan memodelkan data tersebut dengan fungsi lm yang ada id R.

modelRegresi <- lm(Y~.,data = dataRegresi)
summary(modelRegresi)
## 
## Call:
## lm(formula = Y ~ ., data = dataRegresi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0364 -0.6052 -0.0207  0.6087  3.1774 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.97484    0.03033   328.8   <2e-16 ***
## X1           2.97225    0.02042   145.6   <2e-16 ***
## X2           4.97282    0.01939   256.5   <2e-16 ***
## X3           7.00493    0.01989   352.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9589 on 996 degrees of freedom
## Multiple R-squared:  0.9955, Adjusted R-squared:  0.9954 
## F-statistic: 7.275e+04 on 3 and 996 DF,  p-value: < 2.2e-16

Berdasarkan output regresi diatas, nilai koefisien regresi yang dihasilkan mendekati nilai koefisien regresi yang kita bangkitkan. Hal ini membuktikan bahwa data sudah sesuai dengan model regresi yang kita inginkan.

Latihan 6

Berdasarkan Latihan 5 lakukan pengulangan sebanyak 100 ulangan dan kemudian hitung rata-rata dan simpangan baku dari masing-masing koefisien regresi.

Jawaban 6

#banyaknya amatan
n <- 1000
ulangan <- 100

# membuat tempat penyimpanan
# koefisien tiap ulangan
koefisien_ulangan <- NULL

set.seed(13)
for(i in seq(ulangan)){
epsilon <- rnorm(n,mean = 0,sd = 1)
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
Xgab <- cbind(1,X)
Y <- Xgab%*%betas+epsilon
dataRegresi <- data.frame(Y,X)
colnames(dataRegresi) <- c("Y","X1","X2","X3")
modelRegresi <- lm(Y~.,data = dataRegresi)
# rbind digunakan untuk menggabungkan baris
koefisien_ulangan <- rbind(koefisien_ulangan,coef(modelRegresi))
}
head(koefisien_ulangan)
##      (Intercept)       X1       X2       X3
## [1,]    9.998451 2.997131 4.961536 7.019022
## [2,]    9.999899 2.981410 5.012151 6.977783
## [3,]    9.971732 2.986735 4.992031 7.019945
## [4,]    9.940730 3.013790 5.006416 7.003125
## [5,]   10.029303 3.033154 5.009336 6.972909
## [6,]   10.008416 3.019092 5.010800 6.999371
#menghitung mean dari setiap koefisien
colMeans(koefisien_ulangan)
## (Intercept)          X1          X2          X3 
##    9.997236    3.002950    5.000115    7.001103
#menghitung sd dari setiap koefisien
apply(koefisien_ulangan,2,sd)
## (Intercept)          X1          X2          X3 
##  0.02972359  0.01960207  0.01923794  0.02379381

Pengulangan seperti ini biasanya digunakan ketika kita ingin melakukan kajian pada model-model statistik. Menghitung mean dari semua koefisien dimaksudkan untuk melihat sebaik mana model yang kita gunakan dalam mengestimasi parameter. Parameter dalam regresi pada data bangkitan ini adalah koefisien yang kita tentukan di awal. Semakin dekat nilai rata-rata semua koefisien tersebut dengan parameter maka bisa disimpulkan modelnya sudah baik. Sedangkan simpangan baku dari semua koefisien dapat diinterpretasikan sebagai sebagai galat pendugaan dari masing-masing ulangan.

Note

fungsi apply digunakan untuk menerapkan fungsi di R berdasarkan baris dan kolom. Coba gunakan fungsi mean, min, max, dan as.character dan lihat bagaimana hasilnya. Angka 2 pada fungsi apply berarti operasi fungsi dilakukan berdasarkan kolom. Jika diubah menjadi angka 1 berarti operasi fungsi dilakukan berdasarkan baris.