Pendahuluan

Pembangkitan bilangan acak merupakan alat yang diperlukan dalam komputasi statistik umumnya untuk simulasi. Bilangan acak yang dibangkitkan merupakan pseudo random (acak semu). Pseudo random artinya jika dilakukan pengujian terhadap pengacakan, bilangan tersebut bersifat acak, tetapi sebenarnya dibangkitkan oleh sebuah fungsi. Bilangan acak yang dibangkitkan umumnya memenuhi sebaran statitik tertentu, sehingga kita memerlukan pdf/pmf atau cdf-nya. Pembangkitan bilang acak teorinya memerlukan bilangan acak uniform (0,1).

R telah menyiapkan banyak fungsi untuk membangkitkan bilangan acak berdasarkan sebaran yang sudah umum. Untuk membangkitkan bilangan acak dari suatu sebaran, nama fungsinya diawali dengan huruf r dikuti dengan nama sebarannya. Contoh data sebaran pseudo normal dibangkitkan dengan fungsi rnorm. Sewaktu membangkitkan bilangan acak, tidak benar-benar acak, karena kita membutuhkan suatu nilai awal, atau diistilahkan dengan nilai seed. Umumnya bilangan seed mengambil waktu di komputer sampai dengan milisecon. Tetapi, jika kita tetapkan nilai seed nya bilangan acak apapun nilainya akan tetap sama. Untuk menetapkan nilai seed menggunakan fungsi set seed.

Fungsi-fungsi peluang suatu sebaran

• Menghitung nilai density/ peluang fungsi kepekatan/ massa (pdf/pmf), dimulai dengan huruf d diikuti dengan nama sebarannya : dnorm dunif
  
• Menghitung fungsi peluang sebaran kumulatif (cdf), dimulai dengan huruf p diikuti dengan nama sebarannya: pnorm punif
  
• Fungsi invers/ quantile cdf, dimulai dengan huruf q diikuti dengan nama sebarannya: qnorm qunif

Setiap sebaran memiliki parameter yang mungkin berbeda. Walaupun secara umum setiap sebaran akan terkait tiga parameter yang berbeda. Paramater lokasi (\(\mu\)), dispersi (contoh ragam) dan shape (bentuk). Misal pada binom memerlukan parameter size dan prob. Beberapa fungsi sebaran memiliki default, misalkan sebaran normal dengan default \(\mu\)= 0 dan \(\sigma\)-nya adalah 1.

Contoh: membangkitkan bilangan acak sebanyak 10 dari distribusi normal menggunakan fungsi rnorm. rnorm mempunyai default normal baku, dengan \(\mu\) = 0 dan \(\sigma\) =1

rnorm(10)

 [1]  0.9999273  0.1350045 -1.7002207  0.1097219 -0.3376966 -0.3397540
 [7]  0.5986394 -0.3771589 -1.4077928 -0.8906223

Dengan mengatur set.seed sehingga menghasilkan bilangan acak yang diperoleh akan sama persis pada setiap pembangkitan

set.seed(1)
rnorm(10)

 [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684
 [7]  0.4874291  0.7383247  0.5757814 -0.3053884

# membangkita dengan mu = 5 dan standar deviasi =1
rnorm(10, mean = 5, sd = 1)

 [1] 6.511781 5.389843 4.378759 2.785300 6.124931 4.955066 4.983810 5.943836
 [9] 5.821221 5.593901

Membangkitkan bilangan acak dari distribusi uniform dengan fungsi runif.

# membangkitkan bilang acak sebanyak 10
runif(10)  #parameter default 0,1

 [1] 0.8209463 0.6470602 0.7829328 0.5530363 0.5297196 0.7893562 0.0233312
 [8] 0.4772301 0.7323137 0.6927316

# nilai minimum 1 maksimum 100
runif(10, 1, 100)

 [1] 48.284343 86.259738 44.371614 25.234930  7.997226 10.847150 32.310899
 [8] 52.344792 66.538503 41.276189

Mencari nilai peluang kumulatif peubah acak

p1 <- pnorm(1.645)
p1

[1] 0.9500151

Mencari nilai kuantil dari peluang peubah acak

q1 <- qnorm(0.975)
q1

[1] 1.959964

Mencari nilai density peubah acak

a <- seq(-4, 4, length = 1000)
da <- dnorm(a)
plot(a, da)

Teknik Pembangkitan Bilangan Acak

Jika fungsi untuk pembangkitan bilangan acak belum ada, karena sebarannya belum ada. Misalkan kita sudah memiliki suatu peluang sebaran acak, maka kita bisa melakukan pembangkitan bilangan acak. Teknik yang dapat digunakan adalah:

• Inverse-Transform method
  
• Acceptance-rejection method
  
• Direct Transformation –> istilah saja, metodenya tidak hanya satu

r <- runif(100)
x <- qnorm(r, mean = 10, sd = 1)  #invers dari fungsi sebaran kumulatif
hist(x, prob = TRUE)  # untuk melihat pola hasil datanya
sbx <- seq(7, 14, 0.01)  # karena datanya dari 8 - 13, peningkatan garis by 0.01 supaya smooth
lines(sbx, dnorm(sbx, mean = 10, sd = 1))  # gunakan denisty

n <- 1000
u <- runif(n)  #bangkitakan bilangan acak sebanyak 1000 dari distribusi uniform (0,1)
x <- u^(1/3)  # invers fungsi kumulatifnya
hist(x, prob = TRUE)  # periksa dengan histogram
sbx <- seq(0, 1, 0.01)  # sumbu x ibatasi 0 - 1
lines(sbx, 3 * sbx^2)  # buat garis sesuai fungsinya. hasil garis dan histogram sama.

u <- runif(10000)  # bangkitkan bilangan acak uniform
x <- qbinom(u, size = 5, prob = 0.5)  # menghitung invers cdf binom
barplot(table(x))  #cek dengan barplot karena peubah acak diskret

u <- runif(100)  #bangkitkan bilangan acak uniform
x <- ifelse(u <= 0.5, 0, ifelse(u <= 0.8, 1, 2))  # fungsi inverse
barplot(table(x))

set.seed(10)
# Eksponensial (lambda=3)

eks <- function(n, lambda) {
    U <- runif(n)
    X <- -log(1 - U)/lambda
    return(X)
}

# invers transform method
yy1 <- eks(1000, 3)

# fungsi bawaan R
yy2 <- rexp(1000, rate = 3)

par(mfrow = c(1, 2))

hist(yy1, main = "Exp dari Inverse Transform")
hist(yy2, main = "Exp dari Fungsi Rexp")

curve((3/2) * x^3 + (11/8) * x^2 + (1/6) * x + (1/12), 0, 1)

n <- 1000
j <- k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1)
    if (u < (((3/2) * (x^3)) + ((11/8) * (x^2)) + ((1/6) * x) + (1/12))/3.5 * (dunif(x, 
        0, 1))) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)
sbx <- seq(0, 1, 0.01)
lines(sbx, ((3/2) * (sbx^3)) + ((11/8) * (sbx^2)) + ((1/6) * sbx) + (1/12))

curve(6 * x * (1 - x), 0, 1)

x <- seq(0, 1, 0.01)
y <- 6 * x * (1 - x)
max(y)

[1] 1.5

n <- 1000
j <- k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1)
    if (3 * x * (1 - x) > u) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)
sbx <- seq(0, 1, 0.01)
lines(sbx, dbeta(sbx, 2, 2))

curve(3 * (x^2))

x <- seq(0, 1, 0.01)
y <- 3 * (x^2)
max(y)

[1] 3

n <- 1000
j <- 0
k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1)
    if (u < (3 * (x^2))/3.5 * dunif(x, 0, 1)) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)
sbx <- seq(0, 1, 0.01)
lines(sbx, (3 * (sbx^2)))

jumlahY <- length(y)
cat("Banyaknya y=", jumlahY)

Banyaknya y= 1000

curve(12 * (x^2) - 12 * (x^3), 0, 1)

x <- seq(0, 1, 0.01)
y <- 12 * (x^2) - 12 * (x^3)
max(y)

[1] 1.777644

set.seed(1501)

n <- 1000
j <- 0
k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1)
    if (u < (12 * (x^2) * (1 - x))/2 * dunif(x, 0, 1)) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)
sbx <- seq(0, 1, 0.01)
lines(sbx, (12 * (sbx^2 * (1 - sbx))))

jumlahY <- length(y)
cat("Banyaknya y=", jumlahY)

Banyaknya y= 1000

set.seed(1501)

n <- 1000
j <- 0
k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1)
    if ((6 * (x^2) * (1 - x)) > u) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)
sbx <- seq(0, 1, 0.01)
lines(sbx, (12 * (sbx^2 * (1 - sbx))))

cat("Jumlah iterasi= ", j)

Jumlah iterasi=  1991

cat("Jumlah bilangan acak= ", length(y))

Jumlah bilangan acak=  1000

curve((1/2) * x)

set.seed(1)
n <- 1000
j <- 0
k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1, 0, 2)
    if (u < ((1/2) * x)/1.1 * dunif(x, 0, 2)) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)

jumlahY <- length(y)
cat("Banyaknya y=", jumlahY)

Banyaknya y= 1000

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

Multivariate Normal Distribution

Terkadang kita ingin membangkitkan bilangan acak lebih dari satu yang di antara peubah bebas yang dibangkitkan tersebut memiliki hubungan. Pada pemrograman R untuk membangkitkan bilangan acak multivariate normal bisa dilakukan melalui library MASS. Jika korelasi 0 tinggal membangkitkan univariate sebanyak yang diinginkan, tapi jika membutuhkan korelasi dapat menggunakan multivariate normal.

Mendapatkan multivariate normal dengan teknik eigen

rmvn.eigen <- function(n, mu, Sigma) {
    d <- length(mu)
    ev <- eigen(Sigma, symmetric = TRUE)
    lambda <- ev$values
    V <- ev$vectors
    R <- V %*% diag(sqrt(lambda)) %*% t(V)
    Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
    X <- Z %*% R + matrix(mu, n, d, byrow = T)
    X
}
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.9, 0.9, 1), nrow = 2, ncol = 2)
X <- rmvn.eigen(1000, mu, Sigma)
cov(X)  #atauvar(X) > cor(X)

          [,1]      [,2]
[1,] 1.0981130 0.9892851
[2,] 0.9892851 1.0901014

head(X)

            [,1]        [,2]
[1,]  0.49105054 -0.14808204
[2,] -0.07872957 -1.03822367
[3,]  0.12238043  0.90993884
[4,]  0.45433461  0.55190300
[5,]  0.02913881 -0.01046732
[6,] -1.03893009 -0.29292777

Melalui library MASS untuk membangkitkan bilangan acak multivariate normal

library(MASS)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.9, 0.9, 1), nrow = 2, ncol = 2)
X1 <- mvrnorm(1000, mu, Sigma)
head(X1)

           [,1]       [,2]
[1,]  0.9741111  0.4666891
[2,]  0.2058594  0.5477797
[3,]  1.1359570  1.3911863
[4,] -0.4809639 -1.0854598
[5,] -1.1081501 -2.0159425
[6,]  0.7775855  1.0416556

X2 <- mvrnorm(1000, mu, Sigma, empirical = TRUE)  # empirical nilai korelasinya akan sama persis
var(X2)

     [,1] [,2]
[1,]  1.0  0.9
[2,]  0.9  1.0

Pembangkitan bilangan acak untuk model

Jika dalam multivariate antar X memiliki korelasi. Pada model, ada satu peubah yang merupakan fungsi dari peubah lainnya. Membangkitkan bilangan acak untuk model berarti membangkitkan bilangan acak dengan peubah lebih dari satu dan peubahnya memiliki keterkaitan fungsi. Dicontohkan dibangkitkan dari model regresi linear. Bangkitkan peubah Y, X1, X2 dan X3 berdasarkan model regresi linear berganda berikut ini:
\[\mathrm{Y} = 10+3X_1+5X_2+7X_3+e\] dengan mengasumsikan bahwa ϵ ~ N(0,1). Banyaknya amatan yang dibangkitkan adalah 1000

# bangkitkan residual berdistribusi normal

# banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n, mean = 0, sd = 1)

# bangkitkan peubah penjelas Xi 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

# bangkitkan peubah 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

# pengecekan dengan model regresi
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

Misalkan akan dibangkitkan model regresi linear sederhana, dengan:
\(Y =\) \(\beta_0\) + \(\beta_1\)\(X\) + \(ϵ\)
dalam hal ini \(ϵ\) ~ N(0,1), \(\beta_0\) = 1 dan \(\beta_1\)=1

b0 <- 1
b1 <- 1
b0hat <- NULL  # untuk menyimpan hasil dugaan
b1hat <- NULL
for (i in 1:100) {
    eps <- rnorm(10)
    X <- runif(10, 5, 10)
    Y <- b0 + b1 * X + eps
    obj <- lm(Y ~ X)
    obj
    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.250713 0.9661114
sd     1.955215 0.2639555

Misalkan akan dibangkitkan model regresi linear dengan X1 dan X2 sebesar 0.9
\(Y =\) \(\beta_0\) + \(\beta_1\)\(X_1\) + \(\beta_2\)\(X_2\) + \(ϵ\)
dalam hal ini \(ϵ\) ~ N(0,1), \(\beta_0\) = 1, \(\beta_1\)=1, \(\beta_2\)=1

library(MASS)
b0 <- 1
b1 <- 1
b2 <- 1
b0hat <- NULL
b1hat <- NULL
b2hat <- NULL
Sigma <- matrix(c(1, 0.9, 0.9, 1), nrow = 2, ncol = 2)
mu <- c(1, 1)
for (i in 1:100) {
    eps <- rnorm(10)
    X <- mvrnorm(10, mu, Sigma)
    Y <- b0 + b1 * X[, 1] * b2 * X[, 2] + eps
    obj <- lm(Y ~ X)
    b0hat <- c(b0hat, obj$coefficients[1])
    b1hat <- c(b1hat, obj$coefficients[2])
    b2hat <- c(b2hat, obj$coefficients[3])
}
hasil <- matrix(c(mean(b0hat), sd(b0hat), mean(b1hat), sd(b1hat), mean(b2hat), sd(b2hat)), 
    nrow = 2, ncol = 3)
rownames(hasil) <- c(" mean ", "sd")
colnames(hasil) <- c("b0", "b1", "b2")
hasil

              b0        b1        b2
 mean  0.7367504 0.9760073 0.9481515
sd     1.1343515 1.4791844 1.5394227

Referensi

Raharjo, Mulianto. (10 Maret 2021). STA561 Pemrograman Statistika: Pembangkitan Bilangan Acak. Retrieved from https://newlms.ipb.ac.id/

Soleh, A.M. (2021). STA561 Pemrograman Statistika: Pembangkitan Bilangan Acak. Retrieved from https://newlms.ipb.ac.id/

Rejection Sampling + R Demo. Retrieved from https://www.youtube.com/watch?v=kMb4JlvuGlw