Latihan 1

Bangkitkan bilangan acak yang berdistribusi exponential dengan Inverse Transform Method, yang amatannya berjumlah 1000. Bandingkan hasilnya dengan fungsi bawaan R rexp dengan menggunakan histogram.

n <- 1000
u <- runif(n)
set.seed(10)
x <- -(log(1-u)/3)
head(x)
## [1] 0.1587204 0.2264414 0.6661987 0.5897032 0.3685879 0.4080832
length(x)
## [1] 1000
set.seed(10) 
y <- rexp(n,rate = 3)
head(y)
## [1] 0.004985469 0.306740402 0.250719646 0.525013950 0.077219539 0.362224334
length(y)
## [1] 1000
par(mfrow=c(1,2)) 
hist(x,main="Exp dari Inverse Transform", col= "skyblue")
hist(y,main="Exp dari fungsi rexp", col= "lightgreen")

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!

n=1000
set.seed(10)
y = runif(n,0,1)
set.seed(30)
u = runif(n,0,1)
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
#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
length(x)
## [1] 668
z <- rbeta(n,shape1 = 2,shape2 = 2)
par(mfrow=c(1,2)) 
hist(x,main="beta(2,2) dari AR", col= "lavender")
hist(z,main="beta(2,2) dari rbeta", col= "skyblue")

Latihan 3

Bangkitkan peubah acak X ~ Gamma(3/2,1) yang memiliki fungsi PDF dengan mengunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histrogramnya!

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
x <- y[criteria]
head(x)
## [1] 0.02243461 1.38033181 1.12823841 2.36256278 0.34748792 1.63000950
length(x)
## [1] 1393
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 = "cadetblue")
hist(z,main="beta(3/2,1) dari rgamma", col = "magenta")

Latihan 4

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

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 = "slateblue")
hist(v,main="t(2) dari rt", col = "goldenrod")

Latihan 5

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

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