Fungsi Pembangkit Random

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

Tentukan CDF $F(x)$ PDF dari $exponential(3)$ $f(x)=3e−3x,x≥ 0$ CDF dari $exponential(3)$ $F(x)=1−e−3x,x≥0$
Tentukan inverse CDF $F(x)$ inverse CDF dari $exponential(3) F^{−1}(u)=−\frac{ln(1−u)}{3},x≥0$
Bangkitkan $u∼U(0,1)$

# banyaknya amatan
n <- 1000
# membangkitan U(0,1)
u <- runif(n)

# fungsi inverse CDF
set.seed(10) # membuat seed
x <- -(log(1-u)/3) #log merupakan ln di R
head(x)

## [1] 0.25729764 0.11191569 0.76838511 0.25356351 0.04066605 0.04837109

Dapatkan bilangan acak $x$ dengan menghitung $F^−1(u)$

# 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",col="red")
hist(y,main="Exp dari fungsi rexp",col="blue")

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.

6.1 Algoritme Acceptance-Rejection Method

Algoritme Acceptance-Rejection Method adalah sebagai berikut:

Tetapkan peubah acak dari $Y$ sebagai target sebaran beserta PDFnya $g(y)$
Bangkitkan peubah acak $Y$ tersebut.
Bangkitkan $u$ berdasarkan distribusi $Uniform(0,1)$ atau $u$ $U(0,1)$
Hitung nilai $c$ dengan mencari nilai maksimum dari $\frac{f(y)}{g(y)}$
Jika $u ≤ \frac{f(y)}{cg(y)}$ , jadikan $Y$ sebagai $X$.
Jika kondisi ke-4 tidak terpenuhi kembali ke-1.

LATIHAN 2

Bangkitkan peubah acak $X∼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 histrogramnya!

JAWABAN 2

1.Tetapkan peubah acak dari $Y$ sebagai target sebaran beserta PDFnya $g(y)$

Dalam Kasus ini kita pilih $Y∼U(0,1)$ karena memiliki domain yang sama dengan peubah acak $X$ yaitu $0<x<1$ . Pdf dari $Y$ adalah

Bangkitkan peubah acak $Y$ tersebut.

# banyaknya amatan
n=1000
set.seed(10)
y = runif(n,0,1)

Bangkitkan u berdasarkan distribusi $Uniform(0,1)$ atau $u$ $U(0,1)$

set.seed(30)
u = runif(n,0,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(y)}{dy}=0$ $6(1−y)−6y=0$ $6−6y−6y=0$ $6−12y=0$ $y=12$

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

Dengan menggunakan aturan perkalian pada turunan maka didapatkan: 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)$ .

Jika $u ≤ 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",col="red")
hist(z,main="beta(2,2) dari rbeta",col="blue")

LATIHAN 3

Bangkitkan peubah acak $X∼Gamma(\frac{3}{2},1)$ yang memiliki fungsi PDF $f(x)=\frac{1}{(Γ(\frac{3}{2})} )x^\frac{1}{2}e^{−x},0<x<∞$ dengan menggunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histrogramnya!

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

JAWABAN 3

Tetapkan peubah acak dari $Y$ sebagai target sebaran beserta PDFnya $g(y)$

# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)

Note : Distribusi exponential(θ) di R memiliki pdf

Bangkitkan u berdasarkan distribusi $Uniform(0,1)$ atau $u$ $U(0,1)$

set.seed(30)
u = runif(n,0,1)

Hitung nilai c dengan mencari nilai maksimum dari $\frac{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

Jika $u≤\frac{f(y)}{cg(y)}$ , jadikan $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",col="blue")
hist(z,main="beta(3/2,1) dari rgamma",col="red")

LATIHAN 4

Bangkitkan peubah acak $X∼t(2)$ sebanyak 1000 amatan yang didapatkan dari hasil transformasi langsung peubah acak $Y∼χ^2(2)$ dan $Z∼N(0,1)$ . Buatlah histogramnya dan bandingkan hasilnya dengan hasil bangkitan yang diperoleh dari fungsi rt.

JAWABAN 4

Distribusi t bisa didapatkan dari hasil transformasi distribusi $χ^2$ dan $N(0,1)$ melalui persamaan dibawah ini: $X=\frac{Z}{√Y^{2}}$

Sehingga berdasarkan persamaa 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

Menamplikan histogram di R

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

Algoritme pembangkitan bilangan acak untuk Regresi Linear

Algoritme pembangkitan bilangan acak untuk Regresi Linear adalah sebagai berikut: 1. Bangkitkan residual berdasarkan distribusi normal $ϵ∼N(0,σ)$ 2. Bangkitkan peubah penjelas $Xi$ yang saling bebas. 3. Bangkitkan peubah respon $Y$ dengan menggunakan model regresi

LATIHAN 5

Bangkitkan peubah $Y,X1,X2$ dan $X3$ berdasarkan model regresi linear berganda berikut ini: $Y=10+3X1+5X2+7X3+ϵ$ dengan mengasumsikan bahwa $ ϵN(0,1)$ . Banyaknya amatan yang dibangkitkan adalah 1000

Jawaban 5

Bangkitkan residual berdasarkan distribusi normal $ϵ∼N(0,σ)$

#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

Pembangkitan peubah penjelas $Xi$ dapat dilakukan juga dengan fungsi $mvnorm$ seperti pada slide kuliah.

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.

Fungsi Pembangkit Random

Moh. Aditya Sirojul Hilmi

2022-10-19

LATIHAN 1

JAWABAN 1

Membangkitan dengan fungsi bawaan R

Menamplikan histogram di R

6.1 Algoritme Acceptance-Rejection Method

LATIHAN 2

JAWABAN 2

Menamplikan histogram di R

LATIHAN 3

JAWABAN 3

Menamplikan histogram di R

LATIHAN 4

JAWABAN 4

Menamplikan histogram di R

Algoritme pembangkitan bilangan acak untuk Regresi Linear

LATIHAN 5

Jawaban 5

LATIHAN 6

Jawaban 6