Pembangkitan bilangan acak merupakan alat yang diperlukan dalam komputasi statistik, umumnya untuk simulasi. Bilangan acak yang dibangkitkan merupakan pseudorandom yang memenuhi sebaran statistik tertentu seperti fungsi massa peluang atau fungsi kepekatan peluang dan fungsi sebaran kumulatif.

Semua metode pembangkitan bilangan acak tergantung dari pembangkitan bilangan acak seragam.

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:

Prefix Deskripsi
r- pembangkitan bilangan acak dari suatu sebaran
d- pdf \(f(x)\) / pmf \(P(X=x)\)
p- cdf \(P(X \le x)\)
q- fungsi quantile/invers cdf

Fungsi Sebaran

R dilengkapi dengan sekumpulan fungsi untuk membangkitkan bilangan acak sebaran peluang yang umum digunakan seperti sebaran normal, poisson, binomial dan lain-lain. Beberapa contoh fungsi untuk membangkitkan data yang berasal dari sebaran tertentu antara lain (Peng 2020):

Distribution R name Arguments
beta beta shape1, shape2, ncp
binomial binom size, prob
Chauchy cauchy location, scale
chi-squared chisq df, ncp
exponential exp rate
F f df1, df2, ncp
gamma gamma shape, scale
geometric geom prob
hypergeometric hyper m, n, k
log-normal lnorm meanlog, sdlog
uniform unif min, max
neg binomial nbinom size, prob
normal norm mean, sd
Poisson pois lambda
Student’s t t df, ncp

Fungsi sebaran peluang diskret

Fungsi Deskripsi
-binom sebaran binomial \(Bi(n, \lambda)\)
-hyper sebaran hipergeometrik \(H(N, m, n)\)
-nbinom sebaran binomial negatif \(NB(k, pi)\)
-pois sebaran Poisson \(Po(\lambda, t)\)

Fungsi sebaran peluang kontinu

Fungsi Deskripsi
-unif sebaran seragam \(U(a, b)\)
-exp sebaran eksponensial \(Exp(\lambda)\)
-norm sebaran normal \(N(\pi, \sigma)\)

Teknik Pembangkitan Bilangan Acak

Dalam beberapa kondisi tertentu diperlukan juga teknik pembangkitan bilangan acak selain yang disediakan di dalam pustaka R. Teknik umum dalam pembangkitan bilangan acak, antara lain

  1. Inverse-transform method
  2. Acceptance-rejection method
  3. Direct Transformation

The Inverse Transform Method

Untuk membangkitkan contoh acak \(X\) dengan sebaran tertentu:

  • Tentukan fungsi sebaran kumulatif \(F(x)\) dari sebaran yang diinginkan
  • Hitung inverse dari \(F\) atau \(F^{-1}(x)\)
  • Bangkitkan contoh acak \(u_1, u_2, u_3, ..., u_n\) yang menyebar \(Seragam(0,1)\)
  • Hitung \(x_1 = F^{-1}(u_1)\), \(x_2 = F^{-1}(u_2)\), \(x_3 = F^{-1}(u_3)\), …, \(x_n = F^{-1}(u_n)\)

\(x_1, x_2, x_3\) merupakan contoh acak yang saling bebas dari peubah acak \(X\).

Contoh lain

U <- runif(100)
y_inverse <- -log(1-U)/3

y_default <- rexp(100, rate=3) # fungsi sebaran exp dengan package R
par(mfrow=(c(1,2)))
hist(y_default, main="Plot dari fungsi exp")
hist(y_inverse, main="Plot dari fungsi inverse")

Contoh 1

Contoh membangkitkan sebaran eksponensial.

Misal \(X \sim Eksponensial(\lambda)\)

  • Fungsi sebaran kumulatif

\[ F(X) = 1 - e^{-{\lambda}x}; \quad x \ge 0 \]

  • Fungsi inverse

\[ \begin{aligned} 1 - e^{-{\lambda}x} &= u \\ e^{-{\lambda}x} &= 1 - u \\ {-{\lambda}x} &= ln(1-u) \\ x &= -\frac{1}{\lambda}ln(1-u) \\ \end{aligned} \] Sehingga, \(F^{-1}(u) = -ln(1-u)/\lambda\)

  • Bangkitkan contoh acak \(Seragam(0,1)\)
  • Terapkan fungsi inverse untuk contoh acak yang telah dibangkitkan tersebut
inv.exp <- function(u, lambda){
  -log(1-u)/lambda
}
  
rand.exp <- function(n = 1, lambda = 1){
  u <- runif(n)
  x <- inv.exp(u, lambda)
  return(x)
}
set.seed(123)
rand.exp(n = 10, lambda = 2)
##  [1] 0.16954209 0.77630468 0.26295011 1.07286505 1.41061464 0.02331342
##  [7] 0.37549990 1.11475582 0.40085086 0.30496835
set.seed(123)
X <- rand.exp(n = 1000, lambda = 2)

hist(X, freq=F, breaks =20, xlab='X', main='Contoh Acak Eksponensial')
curve(dexp(x, rate=2) , 0, 3, lwd=2, xlab = "", ylab = "", add = T)

Contoh 2

Bangkitkan contoh acak dengan sebaran yang mempunyai fungsi kepekatan peluang

\[ f(x) = 3x^2, \quad 0 < x < 1 \]

Fungsi sebaran kumulatif

\[ F(x) = \int_0^x f_x(t)\,dt = x^3 \] Sehingga diperoleh fungsi inverse \(F^{-1}(u) = u^{1/3}\)

# http://www.di.fc.ul.pt/~jpn/r/ECS/index.html

inv.f <- function(u){
  u^(1/3)
}
  
inv.transform <- function(n = 1, inv.f){
  u <- runif(n)
  x <- inv.f(u)
  return(x)
}
set.seed(123)
X <- inv.transform(n = 10000, inv.f)

hist(X, freq=F, breaks =20, main="f(x)=" ~3*x^2~". Sample vs Density")
curve(3*x^2, 0, 1, lwd=2, xlab = "", ylab = "", add = T, col="red")

Referensi: Inverse Transform Sampling dan Simulation and Modeling to Understand Change

Acceptance-Rejection Method

Ilustrasi sederhana pada metode acceptance-rejection

a <- 5
b <- 12
target <- function(x) dbeta(x, a, b)
proposal <- dunif

mode = (a-1)/(a+b-2)
M = target(mode)

n = 1000
points <- runif(n)
uniforms <- runif(n)
accept <- uniforms < (target(points)/(M*proposal(points)))

curve(target, lwd=2)
curve(proposal, add=T, col="seagreen", lwd=2)
curve(M*proposal(x), add=T, col="seagreen", lty=2, lwd=2)
points(points, M*uniforms, pch=ifelse(accept, 1, 4),
  col=ifelse(accept, "blue", "red"), lwd=2       
)
legend("topright", c("target", "proposal", "accepted", "rejected"), 
       lwd = c(2, 2, NA, NA), col=c("black", "seagreen", "blue", "red"), 
       pch=c(NA, NA, 1, 4), bg="white")

Berikut ini adalah salah satu contoh algoritme untuk membangkitkan \(X \sim f(x); x_0 < x < x_1\) menggunakan metode Acceptance-rejection:

  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 3

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

ar <- function(n,x0,x1,f) {
  xx <- seq(x0,x1,length=10000)
  c <- max(f(xx))
  terima <- 0; iterasi <- 0
  hasil <- numeric(n)
  
  while(terima<n) {
    x <- runif(1,x0,x1)
    y1<- runif(1,0,c)
    y2<- f(x)
    if(y1<=y2) {
      terima <- terima+1
      hasil[terima]<-x 
    }
    iterasi <- iterasi+1
  }
  
  list(hasil=hasil,iterasi=iterasi)
}

set.seed(10)
f <- function(x) {3*x^2}
yyy <- ar(100,0,1,f)
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))
hist(yyy$hasil, main="f(x)="~3*x^2~" Sample vs Density", prob=T)
x <- seq(0, 1, 0.01)
lines(x, f(x), lwd=2, col="red")

Contoh 4

Algoritma

  1. Carilah peubah acak \(Y \sim g\) sehingga \(f(t)/g(t) \le c\) untuk semua \(t\) di mana \(f(t) > 0\)
  2. Bangkitkan bilangan acak \(y\) dari sebaran dengan fungsi kepekatan/masa peluang \(g\)
  3. Bangkitkan bilangan acak \(u\) dari sebaran \(Seragam(0,1)\)
  4. Jika \(u < f(y)/(c.g(y))\), terima \(y\) sebagai bilangan acak, selainnya tolak \(y\) dan ulangi langkah 2-4

Bangkitkan peubah acak yang menyebar \(Beta(2,2)\)

# http://www.di.fc.ul.pt/~jpn/r/ECS/index.html
accept.reject <- function(f, c, g, rg, n) { 
  n.accepts     <- 0
  result.sample <- rep(NA, n)
  
  while (n.accepts < n) {
    y <- rg(1)               # step 1
    u <- runif(1,0,1)        # step 2
    if (u < f(y)/(c*g(y))) { # step 3 (accept)
      n.accepts <- n.accepts+1
      result.sample[n.accepts] = y
    }
  }
  
  result.sample
}
f  <- function(x) 6*x*(1-x)     # pdf of Beta(2,2), maximum density is 1.5
g  <- function(x) x/x           # g(x) = 1 but in vectorized version
rg <- function(n) runif(n,0,1)  # uniform, in this case
c  <- 2                         # c=2 since f(x) <= 2 g(x)

vals <- accept.reject(f, c, g, rg, 10000) 

hist(vals, breaks=30, freq=FALSE, main="Beta(2,2). Sample vs true Density")
xs <- seq(0, 1, len=100)
lines(xs, dbeta(xs,2,2), col="red", lwd=2)

Contoh 5

Misalnya akan dibangkitkan bilangan acak dengan fungsi kepekatan peluang

\[ f(x) = \frac{3}{2}x^3 + \frac{11}{8}x^2 + \frac{1}{6}x + \frac{1}{12}; \quad 0 \le x \le 1 \]

f  <- function(x) (3/2)*(x^3)+(11/8)*(x^2)+(1/6)*(x)+(1/12) 
g  <- function(x) x/x             # g(x) = 1
rg <- function(n) runif(n, 0, 1)  # uniform, in this case
c  <- 4                           # 

vals <- accept.reject(f, c, g, rg, 10000) 

hist(vals, breaks=30, freq=FALSE, main="Beta(2,2). Sample vs true Density")
xs <- seq(0, 1, len=100)
lines(xs, f(xs), col="red", lwd=2)

Direct Transformation

Beberapa transformasi sebaran dari sebaran lain dapat digunakan untuk membangkitkan bilangan acak:
* Jika \(Z \sim N(0,1)\), maka \(V = Z^2 \sim \chi^2_{(1)}\)
* Jika \(U \sim \chi^2_{(m)}\), \(V \sim \chi^2_{(n)}\), \(U\) dan \(V\) saling bebas, maka \(F = (U/m) / (V/n) \sim F_{(m,n)}\)
* Jika \(Z \sim N(0,1)\), \(V \sim \chi^2_{(n)}\), \(Z\) dan \(V\) saling bebas, maka \(T = \frac{Z}{\sqrt{\frac{V}{n}}} \sim t_{(n)}\)
* dan sebagainya

Contoh 6

Jika \(Z \sim N(0,1)\), maka \(Y = Z^2 \sim \chi^2_{(1)}\)

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

Membangkitkan Bilangan Acak untuk Regresi

Contoh 7

Bangkitkan bilangan acak untuk membangun persamaan regresi linier sederhana antara \(X\) terhadap \(Y\), dengan \(b_0 = 1, b_1 = 1\).

b0 <- 1
b1 <- 1

b0hat <- NULL
b1hat <- NULL

for (i in 1:100) {
  eps <- rnorm (10)
  X <- runif (10 ,5 ,10)
  Y <- b0 + b1*X + eps
  
  obj <- lm(Y~X)
  b0hat <- c(b0hat, obj$coefficients[1])
  b1hat <- c(b1hat, obj$coefficients[2])
}

hasil <- matrix (c(mean(b0hat), sd(b0hat), mean(b1hat) , sd(b1hat)), nrow =2, ncol =2)
rownames(hasil) <- c("mean","sd")
colnames(hasil) <- c("b0", "b1")
hasil
##            b0        b1
## mean 1.416092 0.9471990
## sd   1.912003 0.2529918

Contoh 8

Bangkitkan bilangan acak untuk membangun persamaan regresi linier berganda antara \(X_1\) dan \(X_2\) terhadap \(Y\), sehingga diperoleh \(b_0 = 10, b_1 = 2.3, b_2 = 0.7\).

set.seed(123)
X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,1)

model <- lm(Y~X1+X2)
summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.64768 -0.51104 -0.05084  0.56224  2.28118 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) 8.819021   1.489143   5.922           0.00000584 ***
## X1          2.347724   0.045223  51.914 < 0.0000000000000002 ***
## X2          0.702732   0.006784 103.594 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9363 on 22 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.9978 
## F-statistic:  5509 on 2 and 22 DF,  p-value: < 0.00000000000000022

  1. Mahasiswa Pascasarjana IPB, ↩︎