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.

# banyaknya amatan
n <- 1000
# membangkitkan U(0,1)
u <- runif(n)
#fungsi inverse CDF
set.seed(10) #membuat seed
x <- -(log(1-u)) #Log merupakan Ln di R
head(x)
## [1] 2.933906694 1.060397922 0.655152269 1.437412851 0.004792286 1.284341819
#memeriksa banyaknya amatan
length(x)
## [1] 1000

#membangkitkan 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

#Menampilkan histogram di R

par(mfrow=c(1,2)) 
hist(x,main="Exp dari Inverse Transform", col = c("pink"))
hist(y,main="Exp dari fungsi rexp", col = c("yellow"))

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!

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 \({f(y)}/{g(y)}\)

misalkan saja \(h(y)= f(y)/g(y)\), sehingga \(h(y)=6y(1-y)/1=6y(1-y)\) untuk mencari nilai maksimum dari \(h(y)\) maka kita harus menyelesaikan persamaan \(dh(y)/dy=0\) Dengan menggunakan aturan perkalian pada turunan maka didapatkan: \(6(1-y)-6y=0\) \(6-6y-6y=0\) \(6-12y=0\) \(y=1/2\)

Subtitusi \(y=1/2\) ke \(h(y)\) sehingga didapatkan: \(h(1/2)=6(1/2)(1-1/2)=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<=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.

#menampilkan Histogram

z <- rbeta(n,shape1 = 2,shape2 = 2)
par(mfrow=c(1,2)) 
hist(x,main="beta(2,2) dari AR", col=c("cyan"))
hist(z,main="beta(2,2) dari rbeta", col =c("pink"))

Latihan 3

Bangkitkan peubah acak $X~Gamma(3/2,1) yang memiliki fungsi PDF

$f(x)=1/Γ(3/2)x1/2e-x,0<x<∞

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

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

Dalam Kasus ini kita pilih \(Y~exponential(3/2)\) karena memiliki domain yang sama dengan peubah acak \(X\) yaitu \(0<x<∞\) dan karena distribusi exponensial merupakan kasus khusus dari distribusi gamma. Pdf dari Y adalah

\(g(y)=!/(3/2)e ^-y/(3/2), 0<=y<=∞\)

atau \(g(y) = 2/3e^-2/3y, 0<=y<=∞\)

  1. Bangkitkan peubah acak \(Y\) tersebut.
# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)
  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 \(f(x)/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

Beberapa catatan yang perlu diperhatikan untuk menggunakan fungsi optim adalah secara default fungsi ini hanya bisa mencari nilai minimum dari suatu fungsi. Oleh karena itu, fungsi \(h\) perlu dikalikan dengan \(1\). Argumen par digunakan untuk menentukan nilai awal.

  1. Jika \(u<=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=c("plum"))
hist(z,main="beta(3/2,1) dari rgamma", col=c("coral"))

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.

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

#Menampilkan Histogram di R

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

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

  1. Bangkitkan residual berdasarkan distribusi normal \(ϵ∼N(0,σ)\)
#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 1)
  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
  1. 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

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

Latihan 6

Berdasarkan Latihan 5 lakukan pengulangan sebanyak 100 ulangan dan kemudian hitung rata-rata dan simpangan baku dari masing-masing koefisien regresi.

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