5 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.

6 Jawaban 1

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

  1. Tentukan inverse CDF \(F(X)\)

inverse CDF dari \(exponential(3) F^{-1}(u) = \frac{-ln(1-u)}{3} , x>=0\)

  1. Bangkitkan \(u\) ~ \(U(0,1)\)
# banyaknya amatan
n <- 1000
# membangkitan U(0,1)
u <- runif(n)
  1. Dapatkan bilangan acak x dengan menghitung F^{-1}(u)
# fungsi inverse CDF
set.seed(10) # membuat seed
x <- -(log(1-u)/3) #log merupakan ln di R
head(x)
## [1] 0.26696179 0.23752991 0.05479975 0.24991206 1.04225958 0.32507957
# 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 =c("lightpink"))
hist(y,main="Exp dari fungsi rexp", col =c("lightblue"))

fungsi par 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:

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

  2. Bangkitkan peubah acak dari \(Y\) tersebut.

  3. Bangkitkan \(u\) berdasarkan distribusi \(uniform(0,1)\) atau \(u\) ~ \(U(0,1)\)

  4. Hitung nilai \(c\) dengan mencari nilai maksimum dari \(\frac{f(y}{g(y)}\)

  5. Jika \(u\) <= \(\frac{f(y}{g(y)}\), jadikan \(Y\) sebagai \(X\).

  6. Jika kondisi ke-4 tidak terpenuhi kembali ke-1.

7 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!

8 Jawaban 2

Untuk menjawab soal Latihan 3, kita akan terapkan algoritme dari Acceptance-Rejection Method sebagai berikut:

  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

\(g(y)=1,0<=y<=1\)

  1. Bangkitkan peubah acak \(Y\) tersebut.
# banyaknya amatan
n=1000
set.seed(10)
y = runif(n,0,1)
  1. Bangkitkan \(u\) berdasarkan distribusi \(uniform(0,1)\) atau \(u~U(0,1)\)
#set.seed(30)
u = runif(n,0,1)
  1. Hitung nilai \(c\) dengan mencari nilai maksimum dari \(\frac{f(y}{g(y)}\)

Misalkan saja \(h(y) =\frac{6y(1-y)}{1}=6y(1-y)\) Untuk mencari nilai maksimum dari \(h(y)\) maka kita harus menyelesaikan persamaan

=0$

Dengan menggunakan aturan perkalian pada turunan maka didapatkan:

\(6y(1-y)-6y=0\)

\(6-6y-6y=0\)

\(6-12y=0\)

\(y= \frac{1}{2}\)

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

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)\).

  1. Jika \(u\) <= , 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 FALSE  TRUE
#jadikan Y sebagai X
x <- y[criteria]
head(x)
## [1] 0.5074782 0.3067685 0.4269077 0.6931021 0.2254366 0.2745305
# banyaknya amatan
length(x)
## [1] 673

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= c("lightblue"))
hist(z,main="beta(2,2) dari rbeta", col= c("lightgreen"))

10 Jawaban 3

# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)
set.seed(30)
u = runif(n,0,1)
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
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]
z <- rgamma(n,shape = 3/2,scale = 1)
par(mfrow=c(1,2)) 
hist(x,main="gamma(3/2,1) dari AR", col= c("lightblue"))
hist(z,main="beta(3/2,1) dari rgamma", col= c("lightgreen"))

12 Jawaban 4

#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
set.seed(5)
v <- rt(n,df=2)
par(mfrow=c(1,2)) 
hist(x,main="t(2) dari Direct", col= c("lightpink"))
hist(v,main="t(2) dari rt", col= c("lightblue"))

14 Jawaban 5

#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 1)
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
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
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