Teknik umum dalam pembangkitan bilangan acak ada 3, yaitu:

  1. Inverse-Transform Method

  2. Acceptance-rejection method

  3. Direct Transformation

1. Inverse-Transform Method

Jika X adalah peubah acak kontinu dengan cdf \(F(x)\), maka \(U = F(x) \sim Uniform (0,1)\). Jika \(X = F^{-1}(u)\) maka \(X \sim f(x)\). Jadi, algoritma untuk membangkitkan data menggunakan metode ITM:

  1. Tentukan fsk \(F(x)\) dari sebaran yang diinginkan

  2. Tentukan \(F^{-1}(x)\)

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

  4. Transformasi \(X = F^{-1}(u)\)

Contoh 1

Membangkitkan bilangan acak \(Eksponensial(\lambda)\)

  • \(X \sim eks (\lambda)\)

  • \(f(x) = \lambda e ^{- \lambda x}\) , untuk \(x \ge 0\)

  • \(F(x) = 1-e ^{- \lambda x}\) , untuk \(x \ge 0\)

  • \(U=1-e ^{- \lambda x}\) , untuk \(x \ge 0\)

  • \(X=-ln(1-U)/\lambda\)

Algoritma :

  1. Bangkitkan U, bilangan acak Seragam(0, 1)

  2. Hitung \(X=-ln(1-U)/\lambda\)

  3. Ulangi berkali-kali sesuai dengan banyaknya bilangan yang diinginkan

#Eksponensial(lambda=3)
eks <- function(n,lambda){
  U <- runif(n) #membangkitkan bilangan acak seragam sebanyak n
  X <- -log(1-U)/lambda #menentukan inversnya dengan fungsi log
  return(X) #
}
yy1 <- eks(1000,3) #inverse transform method #membangkitkan 1000 bilangan dengan lambda=3
yy2 <- rexp(1000,rate=3) #membangkitkan 1000 bilangan dengan rate=3 menggunakan fungsi bawaan R
par(mfrow=c(1,2)) #membagi layout menjadi 1 baris dan 2 kolom
hist(yy1,main="Exp dari Inverse Transform") #membuat histogram dari invers transform method
hist(yy2,main="Exp dari fungsi rexp") #membuat histogram dari fungsi bawaab R

Kesimpulan:

Baik mmenggunakan invers transform method maupun fungsi bawaan R, hasilnya mirip.

Contoh 2

Bangkitkan \(X \sim f (x) = 3x^2\), \(0 < x < 1\), atau dengan kata lain, \(F_X(x) = x^3\), sehingga \(F^{-1}(u) = u^{1/3}\), \(0 < u < 1\)

n <- 1000 #jumlah bilangan ada 1000
u <- runif(n) #membangkitkan bilang acak sebaran seragam sebanyak n=1000
x <- u^(1/3) # inversnya dari f(x)
hist(x, prob=TRUE) ##membuat histogram dari invers transform method pada f(x)
sbx <- seq(0,1,.01) #sumbu x dari 0 sampai 1 dengan jarak 0.01
lines(sbx, 3*sbx^2) #menambah garis

2. Acceptance-rejection method

Algoritme untuk membangkitkan \(X \sim f(x) ; x_0 < x < x_1\) :

  1. Bangkitkan \(x \sim U(x_0,x_1)\)

  2. Tentukan c sehingga \(f(x) \le c\) untuk \(x_0 < x < x_1\)

  3. Bangkitkan \(y_1 \sim U(0,c)\)

  4. Tentukan \(y_2 = f(x)\)

  5. Jika \(y_1 \le y_2\), terima x

Contoh 1

Bangkitkan bilangan acak yang memiliki fkp \(f_Y(y)=3y^2 ; 0<y<1\) menggunakan acceptance-rejection method !

ar <- function(n,x0,x1,f) {
  xx <- seq(x0,x1,length=10000) #baris bilangan dari x0 sampai x1 dengan 10000 titik
  c <- max(f(xx)) #nilai maksimum dari fungsi f(xx)
  terima <- 0; iterasi <- 0 #menyimpan berapa banyak iterasi yang dibutuhkan dan yg termasuk ke dalam kriteria
  hasil <- numeric(n)
  while(terima<n) { #menjalankan fungsi selama terima < n
    x <- runif(1,x0,x1) #membangkitkan bilangan seragam dari x0 sampai x1
    y1<- runif(1,0,c) #membangkitkan bilangan seragam dari 0 sampai c
    y2<- f(x) #nilai f(x) untuk setiap x yang dibangkitkan 
    if(y1<=y2) { #jika y1<=y2, maka
      terima <- terima+1
      hasil[terima]<-x } #hasil x disimpan pada vektor hasil
    iterasi <- iterasi+1 }
  list(hasil=hasil,iterasi=iterasi) #hasil dari iterasi disimpan dalam bentuk list
}
set.seed(10)
f <- function(x) {3*x^2} #f(x)
yyy <- ar(100,0,1,f) #membangkitkan 100 bilangan dari 0 sampai 1 pada f(x)
yyy
## $hasil
##   [1] 0.8647212 0.7751099 0.8382877 0.7707715 0.5355970 0.8613824 0.2036477
##   [8] 0.7979930 0.7438394 0.3443435 0.9837322 0.6935082 0.6331153 0.8880315
##  [15] 0.7690405 0.6483695 0.8795432 0.9360689 0.7233519 0.7620444 0.9868082
##  [22] 0.8760261 0.7240640 0.8140516 0.5588949 0.8900940 0.7456896 0.8480646
##  [29] 0.8703302 0.8223331 0.8508123 0.7709219 0.8953595 0.5803863 0.5982260
##  [36] 0.9235285 0.7367755 0.6898170 0.8301572 0.9293209 0.9095163 0.5347576
##  [43] 0.3478601 0.8759762 0.7286815 0.8749293 0.6988356 0.8312562 0.5572238
##  [50] 0.6647687 0.7400502 0.9806898 0.3800746 0.7553169 0.5184889 0.8879149
##  [57] 0.9177773 0.8084086 0.8537441 0.4232184 0.7604306 0.3405763 0.3886568
##  [64] 0.4774175 0.5387605 0.9485434 0.7124685 0.9081691 0.9457656 0.7716899
##  [71] 0.6946655 0.5368832 0.8481593 0.8242752 0.5123742 0.3152032 0.9924487
##  [78] 0.9327120 0.9892809 0.6283590 0.5254605 0.8810815 0.5291748 0.5765517
##  [85] 0.7231807 0.8761180 0.3995670 0.8986123 0.9335217 0.7859216 0.7784128
##  [92] 0.6955333 0.9060413 0.9916424 0.4729846 0.9770567 0.9386110 0.9959093
##  [99] 0.8543663 0.8309547
## 
## $iterasi
## [1] 322
par(mfrow=c(1,1)) #membuat layout 1 kolom 1 baris
hist(yyy$hasil,
     main="f(x)=3*x^2", prob=T) #membuat histogram dari f(x) untuk 0 sampai 1
x <- seq(0, 1, 0.01)
lines(x, f(x), lwd=2, col=4) #menambah garis

3. Direct Transformation

Memanfaatkan beberapa fungsi transformasi dari berbagai sebaran yang ada.

Contoh 1

Jika \(X \sim U(0,1)\) menggunakan \(Y=b-aX+a\) maka \(Y \sim U(a,b)\). Membangkitkan bilangan acak \(Y \sim U(3,5) ; Y =2X+3\):

x <- runif(1000)
y <- 2*x+3

Contoh 2

Bangkitkan bilangan acak \(\sim X^2_{(1)}\) !

y <- rnorm(1000) #membangkitkan bilangan acak sebaran normal sebanyak 1000 
x <- y^2  #mengkuadratkan bilangan acak normal yang sudah diperoleh agar berdistribusi chi-square
hist(x,prob=T) #membuat histogram dari x
sbx <- seq(0,13,0.01) ##sumbu x dari 0 sampai 13 dengan jarak 0.01 
lines(sbx,dchisq(sbx,df =1),col="red") #menambah garis