Generating Random Numbers in R

Introduction

Inverse Transform Method

Prosedur:

• Tentukan CDF dari distribusi yang dibangkitkan
• Diketahui u
• Cari inverse CDF atau quantile-nya dalam x
• Bangkitkan u berdasarkan distribusi Uniform(0,1)
• Diperoleh bilangan acak x dengan menghitung nilai quantile dalam u

Misal:

1. Diketahui bilangan acak distribusi Eksponensial(3) dengan amatan berjumlah 1000.
2. Menentukan inverse CDF x diperoleh F^-1(u) = -ln(1-u)/3
3. Bangkitkan u ~ U(0,1)

Bangkitan Manual

set.seed(1)
n <- 1000
u <- runif(n)
# memasukkan u ke fungsi inverse CDF
x <- -(log(1-u)/3)
head(x, 10)

##  [1] 0.10285903 0.15513747 0.28354264 0.79607595 0.07508273 0.76220341
##  [7] 0.96484508 0.36038629 0.33062022 0.02125917

Bangkitan dengan package R

set.seed(1)
y <- rexp(n, rate=3)
head(y, 10)

##  [1] 0.25172728 0.39388093 0.04856891 0.04659842 0.14535621 0.96498951
##  [7] 0.40985402 0.17989428 0.31885583 0.04901533

Plot Histogram

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

Acceptance Rejection Method

Prosedur:

• Menetapkan peubah acak Y sebagai target sebaran dan PDF-nya juga -> g(y)
• Membangkitkan peubah acak Y
• Membangkitkan u berdasarkan distribusi Uniform(0,1)
• Hitung c dengan mencari maksimum dari f(y)/g(y)
• Jika u <= f(y)/cg(y), jadikan Y sebagai X
• Jika tidak, ulangi langkah dari awal

Contoh:

1. Tetapkan Y sebagai target dan pdf nya
2. Bangkitkan Y

n <- 1000
set.seed(10)
y <- runif(n)

3. Bangkitkan u ~ Unif(0,1)

set.seed(10)
u <- runif(n, 0,1)

4. Menghitung nilai c dengan memaksimumkan f(y)/g(y)

h <- function(y) {6*y*(1-y)}
derivH <- optimize(f = h, interval = c(0,1), maximum = T)
derivH$objective

## [1] 1.5

f <- function(y) 6*y*(1-y)
g <- function(y) 1
nilaic <- derivH$objective
nilaic

## [1] 1.5

#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

x <- y[criteria]
head(x)

## [1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597 0.22543662

# memeriksa banyak amatan
length(x)

## [1] 739

Membangkitkan dengan fungsi rbeta

z <- rbeta(n, shape1 = 2, shape2 = 2)
par(mfrow=c(1,2)); hist(x, main="X ~ Beta(2,2) dari ARM"); hist(z, main="Z ~ Beta(2,2) dari rbeta")

Contoh lain

1. Diminta bangkitan X ~ Gamma(3/2, 1)
2. Dipilih Y ~ Exponential(3/2) karena memiliki domain sama

n <- 1500
set.seed(10)
y <- rexp(n, rate = 2/3)

3. Bangkitkan u ~ Unif(0,1)

set.seed(10)
u <- runif(n)

4. Hitung c dengan memaksimumkan f(y)/g(y)

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

Keterangan: optimize : untuk domain fungsi yang finite/terbatas optim : untuk domain fungsi infinite di salah satu batas atau keduanya Pada fungsi optim secara default dia menghasilkan nilai minimum, maka untuk menghasilkan nilai maksimum dalam selang, dikali -1

5. Terapkan kriteria jika u <= f(y)/cg(y)

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 untuk mengubah nilai minimum ke maksimum

#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)

## [1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE

6. Menjadikan Y yang memenuhi kriteria sebagai X

x <- y[criteria]
head(x)

## [1] 1.3803318 1.1282384 2.3625628 0.3474879 1.6300095 3.4914343

length(x)

## [1] 1398

#karena >1000 ambil 1000 saja
x <- x[1:1000]

Plot Data Bangkitan

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

Direct Transformation

Bangkitan dalam Regresi

Instruksi: Bangkitkan Y, X1, X2, X3 berdasarkan model linier berganda: Y = 10 + 3X1 + 5X3 + 7X3 + e dengan asumsi e ~ N(0,1) dengan n=1000.

1. Bangkitkan residual

n <- 1000
set.seed(10)
eps <- rnorm(n, 0, 1)

2. Bangkitkan Xi

set.seed(10)
X <- replicate(3, rnorm(n, 0, 1.5)) # ragam X harus sama
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

3. Generate Y sesuai model dan data bangkitan

beta <- matrix(c(10,3,5,7), nrow=4)
Xgab <- cbind(1, X)

Y <- Xgab%*%beta+eps

data_reg <- data.frame(Y, X)
colnames(data_reg) <- c("Y", "X1", "X2", "X3")
head(data_reg)

##           Y          X1         X2          X3
## 1 14.745624  0.02811926  1.5750205 -0.46179747
## 2 19.092204 -0.27637881  0.4291389  1.13712837
## 3 -1.763648 -2.05699582  0.3608472 -0.86079512
## 4  3.093050 -0.89875157  1.2490578 -1.40811669
## 5  9.656781  0.44181769 -0.3344748 -0.04154899
## 6  3.110839  0.58469145  0.4325163 -1.59937298

4. Menduga model

model_reg <- lm(Y~., data = data_reg)
summary(model_reg)

## 
## Call:
## lm(formula = Y ~ ., data = data_reg)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.549e-13 -2.410e-15 -1.070e-15  4.500e-16  4.296e-13 
## 
## Coefficients:
##              Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) 1.000e+01  6.940e-16 1.441e+16   <2e-16 ***
## X1          3.667e+00  4.672e-16 7.848e+15   <2e-16 ***
## X2          5.000e+00  4.437e-16 1.127e+16   <2e-16 ***
## X3          7.000e+00  4.551e-16 1.538e+16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.194e-14 on 996 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 1.474e+32 on 3 and 996 DF,  p-value: < 2.2e-16

Kesimpulan

• Kalau cara mengasumsikan modelnya benar, maka hasil menunjukkan bahwa simulasi yang dibuat sudah tepat.
• Simulasi yang dibuat sudah tepat mengindikasikan bahwa model regresi linier yang di R sudah akurat (hampir sama dengan koefisien di awal)
• Biasanya digunakan untuk memeriksa apakah modelnya sudah sesuai teori atau tidak.
• Dalam simulasi X nya nilainya fixed jadi tidak digabungkan dalam iterasi

Regresi dengan >1 iterasi pembangkitan data

1. Menentukan n dan tempat menaruh koefisien

n <- 1000
replic <- 100

replic.coeff <- NULL

2. Mendefinisikan nilai-nilai X yang fixed

set.seed(11)
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
Xgab <- cbind(1,X)

3. Membangkitkan data

set.seed(11)
for(i in seq(replic)){
  eps <- rnorm(n,mean = 0,sd = 1)
  Y <- Xgab%*%beta+eps
  data_reg <- data.frame(Y,X)
  colnames(data_reg) <- c("Y","X1","X2","X3")
  model_reg <- lm(Y~.,data = data_reg)
  # rbind digunakan untuk menggabungkan baris
  replic.coeff <- rbind(replic.coeff,coef(model_reg))
}

head(replic.coeff)

##      (Intercept)       X1       X2       X3
## [1,]   10.000000 3.666667 5.000000 7.000000
## [2,]   10.000000 3.000000 5.666667 7.000000
## [3,]   10.000000 3.000000 5.000000 7.666667
## [4,]   10.004285 2.979482 5.010359 6.980066
## [5,]   10.057015 2.971122 4.982448 6.968867
## [6,]    9.991595 3.009841 5.028527 7.017688

4. Menjadikan nilai tunggal dari koefisien-koefisien bangkitan

apply(replic.coeff, 2, mean) # 2 = berdasarkan kolom, 1 = berdasarkan baris

## (Intercept)          X1          X2          X3 
##   10.004456    3.007780    5.006718    7.005587

apply(replic.coeff, 2, sd)

## (Intercept)          X1          X2          X3 
##  0.03053835  0.06949131  0.07037874  0.06933019