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

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
  2. Tentukan inverse CDF F(x) F^-1 (u) = -ln(1-u)/3 , x>=0
  3. 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) {r echo=FALSE} # fungsi inverse CDF set.seed(10) # membuat seed x <- -(log(1-u)/3) #log merupakan ln di R head(x)

  2. Memeriksa banyak amatan {r echo=FALSE} # memeriksa banyaknya amatan length(x)

  3. Membangkitkan dengan fungsi bawaan R {r echo=FALSE} set.seed(10) # membuat seed y <- rexp(n,rate = 3) head(y)

{r echo=FALSE} # memeriksa banyaknya amatan length(y)

  1. Menampilkan histogram di R

{r echo=FALSE} par(mfrow=c(1,2)) hist(x,main=“Exp dari Inverse Transform”) hist(y,main=“Exp dari fungsi rexp”)

Latihan 2

Bangkitkan peubah acak Y ~ 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 histogramnya!

Jawaban 2

  1. Tetapkan peubah acak dari X sebagai target sebaran beserta PDFnya g(y)=1,0<=y<=1

  2. Bangkitkan peubah acak Y {r echo=FALSE} # banyaknya amatan n=1000 set.seed(10) y = runif(n,0,1)

  3. Bangkitkan u berdasarkan distribusi Uniform(0,1) {r echo=FALSE} set.seed(30) u = runif(n,0,1)

  4. Hitung nilai c dengan mencari maksimum dari f(y)/g(y) dan didapatkan nilai C = 1.5 {r echo=FALSE} h <- function(y) 6y(1-y) derivH <- optimize(f=h,interval=c(0,1),maximum=T) derivH$objective

  5. Jika u <= f(y)/g(y), jadikan Y sebagai X {r echo=FALSE} f <- function(y) 6y(1-y) g <- function(y) 1 nilaic <- derivH$objective #kriteria penerimaan criteria <- u < f(y)/(nilaic*g(y)) head(criteria)

{r echo=FALSE} #jadikan Y sebagai X x <- y[criteria] head(x)

banyaknya amatan {r echo=FALSE} length(x)

  1. Menampilkan histogram di R

{r echo=FALSE} z <- rbeta(n,shape1 = 2,shape2 = 2) par(mfrow=c(1,2)) hist(x,main=“beta(2,2) dari AR”) hist(z,main=“beta(2,2) dari rbeta”)

Latihan 3

Bangkitkan peubah acak X ~ Gamma(3/2, 1) yang memiliki fungsi PDF f(x) = 1/r(3/2) x^1/2 e^-x, 0<x<~ dengan menggunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histogramnya!

Jawaban 3

  1. Tetapkan peubah acak dari Y sebagai target sebaran beserta PDFnya g(y). Didapatkan g(y) = 2/3e^(-2/3y), 0<=y<=~

  2. Bangkitkan peubah acak Y tersebut {r echo=FALSE} # banyaknya amatan n=1500 set.seed(10) y = rexp(n,rate = 2/3)

  3. Bangkitkan u berdasarkan distribusi Uniform(0,1) {r echo=FALSE} set.seed(30) u = runif(n,0,1)

  4. Hitung nilai c dengan mencari nilai maksimum dari f(y)/g(y) {r echo=FALSE} h <- function(y) (-1) * (3/(2gamma(3/2))) (y)^(1/2) * exp(-y/3) derivH <- optim(par = 0,fn=h,lower = 0,method = “L-BFGS-B”) derivH$value

  5. Jika u <= f(y)/cg(y), jadikan Y sebagai X {r echo=FALSE} f <- function(y) (3/(2gamma(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)

{r echo=FALSE} x <- y[criteria] head(x)

  1. Banyaknya amatan {r echo=FALSE} length(x)

{r echo=FALSE} #karena >1000 ambil 1000 saja x <- x[1:1000]

  1. Menampilkan histogram di R

{r echo=FALSE} 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”)

Latihan 4

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

Jawaban 4

  1. Transformasi distribusi x^2 dan N(0,1) melalui persamaan X=z/sqrt(Y/2) {r echo=FALSE} #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. Menampilkan histogram di R {r echo=FALSE} set.seed(5) v <- rt(n,df=2) par(mfrow=c(1,2)) hist(x,main=“t(2) dari Direct”) hist(v,main=“t(2) dari rt”)

Latihan 5

Bangkitkan peubah Y, X1, X2, dan 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

Jawaban 5

  1. Bangkitkan residual berdasarkan distribusi normal E~N(0,sigma) {r echo=FALSE} #banyaknya amatan n <- 1000 set.seed(12) epsilon <- rnorm(n,mean = 0,sd = 1)

  2. Bangkitkan peubah penjelas Xi yang saling bebas {r echo=FALSE} set.seed(10) X <- replicate(3,rnorm(n,mean=0,sd=1.5)) head(X)

  3. Bangkitkan peubah respon Y dengan menggunakan model regresi {r echo=FALSE} 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)

  1. Memodelkan data dengan fungsi lm yang ada di R {r echo=FALSE} modelRegresi <- lm(Y~.,data = dataRegresi) summary(modelRegresi)