Laporan 5 Praktikum Pemrograman Statistik - Tugas Satria June Adwendi
Link Rpubs : klik disini
Praktikum Pertemuan 11 - Komputasi Untuk Regresi
Prosedur Sweep Operator (Eliminasi)
Secara numerik pendugaan koefisien regresi dilakukan dengan suatu algoritma yang disebut Sweep Operator. Algoritma ini berdasarkan eliminasi (Gauss elimination) terhadap matriks simetrik, A, yang berukuran m×m. Unsur diagonal utama matriks A, yaitu akk≠0 , k=1,2,…,m
Prosedur ini digunakan untuk menentukan :
- matriks kebalikan (inverse) -> A-1
- penduga koefisien regresi ->𝜷duga
- JKG (Jumlah Kuadrat Galat)
Latihan
library(fastmatrix)
y = c(1, 3, 3, 2, 2, 1) # vector y
x1 = c(1, 2, 3, 1, 2, 3) # vector x1
x2 = c(1, 1, 1, -1, -1, -1) # vector x2
s = c(1,1,1,1,1,1) # vector ‘satu’
X=cbind(s,x1,x2) # menggabungkan vector s, x1, dan x2
XX=t(X)%*%X # X’X
Xy=t(X)%*%y # X’y
yX=t(y)%*%X # y’X
yy=t(y)%*%y # y’y
b1=cbind(XX,Xy) # menggabungkan X’X dan X’y
b2=cbind(yX,yy) # menggabungkan y’X dan y’y
b=rbind(b1,b2) # menggabungkan X’X , X’y , y’X , dan y’y
sweep.operator(b, 1:3)## s x1 x2
## s 1.166667 -0.50 0.0000000 1.5000000
## x1 -0.500000 0.25 0.0000000 0.2500000
## x2 0.000000 0.00 0.1666667 0.3333333
## -1.500000 -0.25 -0.3333333 3.0833333round(sweep.operator(b, 1:3),3)## s x1 x2
## s 1.167 -0.50 0.000 1.500
## x1 -0.500 0.25 0.000 0.250
## x2 0.000 0.00 0.167 0.333
## -1.500 -0.25 -0.333 3.083k1=sweep.operator(b, k = 1)
k1## s x1 x2
## s 0.1666667 2 0 2
## x1 -2.0000000 4 0 1
## x2 0.0000000 0 6 2
## -2.0000000 1 2 4k2=sweep.operator(b, k = 1:2)
k2## s x1 x2
## s 1.166667 -0.50 0 1.50
## x1 -0.500000 0.25 0 0.25
## x2 0.000000 0.00 6 2.00
## -1.500000 -0.25 2 3.75k3=sweep.operator(b, k = 1:3)
k3## s x1 x2
## s 1.166667 -0.50 0.0000000 1.5000000
## x1 -0.500000 0.25 0.0000000 0.2500000
## x2 0.000000 0.00 0.1666667 0.3333333
## -1.500000 -0.25 -0.3333333 3.0833333k4=sweep.operator(b, k=c(1,3))
k4## s x1 x2
## s 0.1666667 2 0.0000000 2.0000000
## x1 -2.0000000 4 0.0000000 1.0000000
## x2 0.0000000 0 0.1666667 0.3333333
## -2.0000000 1 -0.3333333 3.3333333Dekomposisi Cholesky
Menghasilkan matriks segitiga bawah (akar kuadrat) dari matriks simetrik.
Penyelesaian persamaan linier berdasarkan matriks segitiga bawah (sederhana).
Untuk regresi linier sederhana dan berganda.
Faktorisasi matriks.
Latihan
#Solve the equation X’X β = X’y
#Dekomposisi Cholesky terhadap X’X -> LL’, Ll = X’y -> l , dan L’β = l -> 𝜷duga.
cho=t(chol(XX)) # LL’
el=solve(cho,Xy) # Lc = X’y (el adalah c ) -> L’β = c
el## [,1]
## s 4.8989795
## x1 0.5000000
## x2 0.8164966beta=solve(t(cho),el) # L’β = c (beta adalah 𝜷 ̂)
beta## [,1]
## s 1.5000000
## x1 0.2500000
## x2 0.3333333t(y)%*%y - t(el)%*%el # d2 = y’y-l’l = y’y -y’X(X’X)-1X’y = JKG## [,1]
## [1,] 3.083333Praktikum Pertemuan 12 - EM Algorithm
Suatu prosedur iteratif untuk pendugaan parameter menggunakan metode kemungkinan maksimum (maximum likehood). Konsep dasar: menggabungkan data pengamatan dengan data tidak teramati (hilang/tersembunyi) sehingga bentuk fungsi kemungkinannya lebih mudah untuk dianalisis.
Algoritma (EM):
begin
initialize θ0, T, i=0
do i <- i + 1
E step: compute Q(θ; θi)
M step: θi+1 <- arg max Q(θ; θi)
until [Q(θi+1; θi) - Q(θi; θi-1)] ≤ T
return 𝜃cap <- θi+1
end
Contoh 1
Pendugaan parameter pada data tabulasi dua arah dengan satu buah data hilang (Pawitan, 2001). Pendugaan berdasarkan model rancangan berikut: yij = Miyu + Alphai + Betaj + εij, i =1,2 dan j=1,2,3 Data pengamatan berdasarkan model itu adalah yij:
y11 10
y12 15
y13 17
y21 22
y22 23
y23 ==> data hilang
# Rumus EM dalam rancob
emrancob <- function(ftn, x,yij,i, j, tol = 1e-9, max.iter = 100) {
yold <- yij
ynew <- ftn(x,yold,i,j)
iter <- 1
cat("At iteration 1 value of x is:", ynew, "\n")
while ((abs(ynew-yold) > tol) && (iter < max.iter)) {
yold <- ynew;
ynew <- ftn(x,yold,i,j);
iter <- iter + 1
cat("At iteration", iter, "value of x is:", ynew, "\n")
}
if (abs(ynew-yold) > tol) {
cat("Algorithm failed to converge\n")
return(NULL)
} else {
cat("Algorithm converged\n")
return(ynew)
}
}
# Masukkan urutan matrix yang akan dicari (y23)
ftn <- function(x,yij,i,j){
x <- matrix(c(10,15,17,22,23,yij),ncol=3,byrow=TRUE)
mu <- mean(x)
alfai <- mean(x[i,])- mu
betaj <- mean(x[,j])- mu
yij <- mu + alfai + betaj
return(yij)
}
# Cari y23, dengan nilai awal 0, i=2, j=3
emrancob(ftn, x,0,2,3, tol = 1e-4, max.iter = 100) ## At iteration 1 value of x is: 9
## At iteration 2 value of x is: 15
## At iteration 3 value of x is: 19
## At iteration 4 value of x is: 21.66667
## At iteration 5 value of x is: 23.44444
## At iteration 6 value of x is: 24.62963
## At iteration 7 value of x is: 25.41975
## At iteration 8 value of x is: 25.9465
## At iteration 9 value of x is: 26.29767
## At iteration 10 value of x is: 26.53178
## At iteration 11 value of x is: 26.68785
## At iteration 12 value of x is: 26.7919
## At iteration 13 value of x is: 26.86127
## At iteration 14 value of x is: 26.90751
## At iteration 15 value of x is: 26.93834
## At iteration 16 value of x is: 26.95889
## At iteration 17 value of x is: 26.9726
## At iteration 18 value of x is: 26.98173
## At iteration 19 value of x is: 26.98782
## At iteration 20 value of x is: 26.99188
## At iteration 21 value of x is: 26.99459
## At iteration 22 value of x is: 26.99639
## At iteration 23 value of x is: 26.99759
## At iteration 24 value of x is: 26.9984
## At iteration 25 value of x is: 26.99893
## At iteration 26 value of x is: 26.99929
## At iteration 27 value of x is: 26.99952
## At iteration 28 value of x is: 26.99968
## At iteration 29 value of x is: 26.99979
## At iteration 30 value of x is: 26.99986
## Algorithm converged## [1] 26.99986Teladan 1
Pendugaan parameter pada data genetika (Rao 1973, hal 369), diasumsikan data fenotipe y=(y1,y2,y3,y4)=(125,18,20,34) memiliki distribusi multinomial dengan peluang
{1/2+θ/4,1−θ/4,1−θ/4,θ/4} sehingga fungsi peluang bagi y adalah
Dugalah data hilang x1 dan x2 dengan EM algorithm!
expectation1 <- function(theta) {125*( 1/2 / ((1/2)+(1/4)*theta) )}
expectation2 <- function(theta) {125*( ((1/4)*theta) / ((1/2)+(1/4)*theta) )}
maximization <- function(x2) {(x2 + x5) / (x2 + x3 + x4 + x5)}
x1 <- 10 #initial value
x2 <- 20 #initial value
x3 <- 18
x4 <- 20
x5 <- 34
niter <- 100
theta0 <- 100 #initial value
save_iter <- data.frame("iter"=0,"theta"=theta0,"x1"=x1,"x2"=x2)
for (i in 1:niter){
x1 <- expectation1(theta0)
x2 <- expectation2(theta0)
theta <- maximization(x2)
criteria <- abs(theta-theta0)
theta0 <- theta
save_iter <- rbind(save_iter,c(i,theta0,x1,x2))
if (criteria<10^-7) break
}
save_iter## iter theta x1 x2
## 1 0 100.0000000 10.00000 20.00000
## 2 1 0.8046765 2.45098 122.54902
## 3 2 0.6477018 89.13684 35.86316
## 4 3 0.6295520 94.42151 30.57849
## 5 4 0.6271833 95.07323 29.92677
## 6 5 0.6268695 95.15895 29.84105
## 7 6 0.6268279 95.17031 29.82969
## 8 7 0.6268223 95.17182 29.82818
## 9 8 0.6268216 95.17202 29.82798
## 10 9 0.6268215 95.17205 29.82795Praktikum Pertemuan 13 - Resampling
Proses pendugaan parameter selalu didasari dengan asumsi sebaran dan ukuran contoh yang memadai. Asumsi sebaran ini menjadi dasar penentuan selang kepercayaan bagi parameter (misal 𝜇). Selang kepercayaan bagi 𝜇 juga bisa diperoleh meskipun contoh acak berukuran kecil tetapi kita mengasumsikan populasinya memiliki sebaran Normal. Seandainya contoh acak berukuran kecil dan sebaran populasinya tidak diketahui, selang kepercayaan bagi parameter tidak bisa diperoleh tanpa sebaran.
Contoh acak berukuran kecil tidak jarang diperoleh dan bahkan tidak ada informasi tentang sebarannya.Pendugaan selang bagi parameter harus dengan metode yang didasarkan pada asumsi sebaran. Sebaran populasi data tidak diketahui sehingga kita tidak mempunyai dasar yang kuat untuk menduga selang kepecayaan tanpa sebaran. Bootstrap dan Jackknife dapat digunakan sebagai cara pendugaan selang dalam kondisi data berukuran kecil dan tidak memenuhi asumsi sebaran.
library(boot)Boostrap
Jacknife
Cross-Validation
Salah satu metode evaluasi pemodelan statistik
Kasus khusus dari validasi
Menggunakan konsep dari jackknife
Validasi vs validasi silang (CV)
Pada proses validasi, setiap observasi hanya berperan sebagai anggota data latih/training saja atau sebagai anggota data validasi saja. Jika dilakukan pengulangan, pengacakan pembagian data latih dan validasi dilakukan pada setiap ulangan.
Pada proses validasi silang, setiap observasi berperan menjadi data latih dan validasi, yaitu satu kali menjadi data validasi dan sisanya menjadi data latih. Pengacakan pembagian data menjadi k bagian dilakukan satu kali di awal, kemudian setiap bagian secara bergantian menjadi data latih dan data validasi.
Jenis CV yg banyak digunakan: LOO-CV (leave-one-out CV) dan k-fold CV
LOO-CV
Sisihkan satu buah amatan sebagai gugus data validasi, sedangkan (n-1) amatan lainnya sebagai gugus data training
Susun model analisis menggunakan data training, dan lakukan prediksi terhadap data validasi yang hanya berisi satu amatan tadi
3.Hitung galat prediksi. Misal jika respon berupa numerik: cv(−k) = 1/nk √(Σ(y_cap−y)^2)
4.Ulangi langkah 1-3 dengan menggunakan satu amatan lainnya sebagai gugus validasi
K-fold CV
Bagi gugus data menjadi k buah bagian yang berisi amatan dengan banyaknya yang sama (hampir sama). Misal k=10 dan n=250 maka tiap bagian berisi 25.
Gunakan salah satu bagian sebagai gugus data validasi dan (k-1) bagian lain sebagai gugus data training.
3.Susun model analisis menggunakan data training, dan lakukan prediksi respon pada data validasi. 4. Hitung galat prediksi. Misal jika respon berupa numerik: cv(−k) = 1/nk √(Σ(y_cap−y)^2)
5.Ulangi langkah 1-4 dengan menggunakan satu bagian lainnya sebagai gugus validasi dan sisanya sebagai data training. Ulangi hingga seluruh k buah bagian pernah menjadi data validasi.
Latihan 1
Tabel di bawah ini melaporkan jumlah populasi dari n=49 kota di USA pada tahun 1920 dan 1930, yang dinotasikan dengan u dan x.
Berdasarkan data ini, peniliti tertarik untuk menghitung rasio dari mean yang nantinya dapat digunakan untuk menduga populasi total USA pada tahun 1930 berdasarkan jumlah populasi pada tahun 1920. Jika kota-kota ini merupakan sampel acak dengan (U,X) menotasikan pasangan nilai populasi dari kota yang terpilih secara acak, maka total populasi 1930 merupakan hasil perkalian antara total populasi 1920 dengan rasio dari ekspetasi θ=E(X)/E(U), yang merupakan parameter rasio yang akan diduga. Dengan menggunakan bootstrap dan jackknife, hitunglah galat baku (standard error) dari θ_cap, yang merupakan √V(θ_cap), jika diketahui θ_cap=(x/u)!
Jawab:
Bootstrap
# Import data ke R
dataCity <- read.table(header = TRUE,text=
'u x
138 143
93 104
61 69
179 260
48 75
37 63
29 50
23 48
30 111
2 50
38 52
46 53
71 79
25 57
298 317
74 93
50 58
76 80
381 464
387 459
78 106
60 57
507 634
50 64
77 89
64 77
40 60
136 139
243 291
256 288
94 85
36 46
45 53
67 67
120 115
172 183
66 86
46 65
121 113
44 58
64 63
56 142
40 64
116 130
87 105
43 61
43 50
161 232
36 54
')- Cara Pertama : Boostraping dengan program R sendiri
- Bangkitan sampel baru dengan menggunakan metode penarikan sampel/contoh dengan pengembalian (with replacement) sebanyak B kali.
n <- nrow(dataCity)
B <- 1000
# sampel diambil berdasarkan urutan amatan/index
set.seed(123)
index <- replicate(B,sample.int(n,size =n,replace = TRUE))- Hitung penduga untuk setiap sampel baru
#penduga theta
thetaHat <- function(data) {mean(data$x) / mean(data$u)}
thetaBoost <- NULL
for (i in 1:B){
thetaBoost <- c(thetaBoost, thetaHat(dataCity[index[ ,i], ]))
}
head(thetaBoost,n=20)## [1] 1.332714 1.189381 1.248475 1.251276 1.323737 1.228723 1.228209 1.231559
## [9] 1.279538 1.263938 1.245061 1.191353 1.271465 1.215520 1.289922 1.309773
## [17] 1.266583 1.206254 1.194978 1.208488mean(thetaBoost)## [1] 1.23997- Gunakan hasil pada (b) untuk menghitung ragam V(θ_cap)
Ragam V(θ_cap) bisa didekati menggunakan ragam bootstrap V(θ^(B)).
# ragam
ragamTheta <- var(thetaBoost)
ragamTheta## [1] 0.001208174# galat baku
sqrt(ragamTheta)## [1] 0.03475879- Cara Kedua : Boostraping dengan menggunakan package
boot
thetaBoost2 <- function(data,index){
dataBoost <- data[index, ]
mean(dataBoost$x) / mean(dataBoost$u)
}
bootCity <- boot(data = dataCity,R = 1000,statistic = thetaBoost2)
bootCity##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = dataCity, statistic = thetaBoost2, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 1.239019 0.003905339 0.036673Jacknife
- Bangkitan sampel baru dengan menghilangkan satu amatan pada setiap ulangan.
n <- nrow(dataCity)
# menghilangkan satu amatan
removeOne <- function(i,data) seq(nrow(data))[-i]
set.seed(123)
index <- sapply(1:n,removeOne,data=dataCity)- Hitung penduga untuk setiap sampel baru
#penduga theta
thetaHat <- function(data) {mean(data$x) / mean(data$u)}
thetaJack <- NULL
for (i in 1:n){
thetaJack <- c(thetaJack, thetaHat(dataCity[index[ ,i], ]))
}
thetaJack## [1] 1.244711 1.241282 1.240336 1.231179 1.235917 1.235599 1.236219 1.235142
## [9] 1.224323 1.229612 1.238038 1.239816 1.240819 1.233844 1.250000 1.238755
## [17] 1.239808 1.241864 1.240745 1.243411 1.237138 1.242491 1.237739 1.238609
## [25] 1.240305 1.239479 1.236937 1.245018 1.241114 1.245102 1.245363 1.238741
## [33] 1.239569 1.242230 1.245845 1.245186 1.238172 1.237420 1.246503 1.238323
## [41] 1.242285 1.224490 1.236139 1.241798 1.239581 1.237478 1.239673 1.232373
## [49] 1.237146mean(thetaJack)## [1] 1.239054- Gunakan hasil pada b untuk menghitung ragam V(θ_cap)
Ragam V(θ_cap) bisa didekati menggunakan ragam jacknife V(θ_cap^(Jack)).
# ragam
ragamTheta2 <- (((n-1)^2)/n)*var(thetaJack)
ragamTheta2## [1] 0.001192619# standard error
sqrt(ragamTheta2)## [1] 0.03453431matrix(c(mean(thetaBoost),sqrt(ragamTheta),
mean(thetaJack),sqrt(ragamTheta2)),ncol=2,byrow=T,
dimnames=list(c("Bootstrap","Jackknife"),c("estimator","std error")))## estimator std error
## Bootstrap 1.239970 0.03475879
## Jackknife 1.239054 0.03453431Latihan 2
Gunakanlah data swiss yang ada di R, lalu hitunglah 5-folds CV dan LOO-CV (Leave-One-Out Cross Validation) berdasarkan regresi yang dilakukan pada data tersebut, dengan peubah respon adalah Fertility dan peubah lainnya merupakan peubah prediktor!
# Import data
data(swiss)
swiss## Fertility Agriculture Examination Education Catholic
## Courtelary 80.2 17.0 15 12 9.96
## Delemont 83.1 45.1 6 9 84.84
## Franches-Mnt 92.5 39.7 5 5 93.40
## Moutier 85.8 36.5 12 7 33.77
## Neuveville 76.9 43.5 17 15 5.16
## Porrentruy 76.1 35.3 9 7 90.57
## Broye 83.8 70.2 16 7 92.85
## Glane 92.4 67.8 14 8 97.16
## Gruyere 82.4 53.3 12 7 97.67
## Sarine 82.9 45.2 16 13 91.38
## Veveyse 87.1 64.5 14 6 98.61
## Aigle 64.1 62.0 21 12 8.52
## Aubonne 66.9 67.5 14 7 2.27
## Avenches 68.9 60.7 19 12 4.43
## Cossonay 61.7 69.3 22 5 2.82
## Echallens 68.3 72.6 18 2 24.20
## Grandson 71.7 34.0 17 8 3.30
## Lausanne 55.7 19.4 26 28 12.11
## La Vallee 54.3 15.2 31 20 2.15
## Lavaux 65.1 73.0 19 9 2.84
## Morges 65.5 59.8 22 10 5.23
## Moudon 65.0 55.1 14 3 4.52
## Nyone 56.6 50.9 22 12 15.14
## Orbe 57.4 54.1 20 6 4.20
## Oron 72.5 71.2 12 1 2.40
## Payerne 74.2 58.1 14 8 5.23
## Paysd'enhaut 72.0 63.5 6 3 2.56
## Rolle 60.5 60.8 16 10 7.72
## Vevey 58.3 26.8 25 19 18.46
## Yverdon 65.4 49.5 15 8 6.10
## Conthey 75.5 85.9 3 2 99.71
## Entremont 69.3 84.9 7 6 99.68
## Herens 77.3 89.7 5 2 100.00
## Martigwy 70.5 78.2 12 6 98.96
## Monthey 79.4 64.9 7 3 98.22
## St Maurice 65.0 75.9 9 9 99.06
## Sierre 92.2 84.6 3 3 99.46
## Sion 79.3 63.1 13 13 96.83
## Boudry 70.4 38.4 26 12 5.62
## La Chauxdfnd 65.7 7.7 29 11 13.79
## Le Locle 72.7 16.7 22 13 11.22
## Neuchatel 64.4 17.6 35 32 16.92
## Val de Ruz 77.6 37.6 15 7 4.97
## ValdeTravers 67.6 18.7 25 7 8.65
## V. De Geneve 35.0 1.2 37 53 42.34
## Rive Droite 44.7 46.6 16 29 50.43
## Rive Gauche 42.8 27.7 22 29 58.33
## Infant.Mortality
## Courtelary 22.2
## Delemont 22.2
## Franches-Mnt 20.2
## Moutier 20.3
## Neuveville 20.6
## Porrentruy 26.6
## Broye 23.6
## Glane 24.9
## Gruyere 21.0
## Sarine 24.4
## Veveyse 24.5
## Aigle 16.5
## Aubonne 19.1
## Avenches 22.7
## Cossonay 18.7
## Echallens 21.2
## Grandson 20.0
## Lausanne 20.2
## La Vallee 10.8
## Lavaux 20.0
## Morges 18.0
## Moudon 22.4
## Nyone 16.7
## Orbe 15.3
## Oron 21.0
## Payerne 23.8
## Paysd'enhaut 18.0
## Rolle 16.3
## Vevey 20.9
## Yverdon 22.5
## Conthey 15.1
## Entremont 19.8
## Herens 18.3
## Martigwy 19.4
## Monthey 20.2
## St Maurice 17.8
## Sierre 16.3
## Sion 18.1
## Boudry 20.3
## La Chauxdfnd 20.5
## Le Locle 18.9
## Neuchatel 23.0
## Val de Ruz 20.0
## ValdeTravers 19.5
## V. De Geneve 18.0
## Rive Droite 18.2
## Rive Gauche 19.35-folds CV
- Bagi amatan data menjadi k bagian.
#banyaknya amatan
n <- nrow(swiss)
# mengacak urutan amatan
set.seed(123)
randomIndex <- sample.int(n,n)
k <- 5
# menentukan folds/bagian untuk tiap amatan
folds <- rep(sample.int(k,k),length.out=n)
# membagi amatan berdasarkan folds yang ditentukan
indexfolds <- split(randomIndex,folds)
indexfolds## $`1`
## [1] 14 25 28 39 4 12 36 43 1
##
## $`2`
## [1] 3 26 9 7 41 46 30 20 2
##
## $`3`
## [1] 15 45 38 8 19 33 21 22 6 40
##
## $`4`
## [1] 31 37 5 44 34 11 13 16 32 24
##
## $`5`
## [1] 42 27 29 10 17 35 18 23 47Bagian ke-k digunakan untuk menduga respons, sedangkan bagian lainnya digabungkan untuk menduga model.
Hitung galat prediksi: cv(−k) = 1/nk √(Σ(y_cap − y)^2)
Ulangi sampai setiap bagian pernah digunakan untuk menduga galat prediksi
cv1 <- NULL
for(j in 1:k){
#menduga model
indexModel <- unlist(indexfolds[-j])
dataModel <- swiss[indexModel, ]
regresiSwiss <- lm(Fertility~.,data = dataModel)
# menduga respon
indexRespon <- unlist(indexfolds[j])
dataRespon <- swiss[indexRespon, ]
responseSwiss <- predict(regresiSwiss,newdata=dataRespon)
# menghitung galat prediksi
galat <- dataRespon$Fertility - responseSwiss
cv1 <- c(cv1,1/(n*k)*sqrt(sum(galat^2)))
}
cv1## [1] 0.08892357 0.08021139 0.13716277 0.11542136 0.09583984- Hitung RMSEP
RMSEP1 <- mean(cv1)
RMSEP1## [1] 0.1035118LOO-CV
Untuk LOO-CV bisa menggunakan program untuk k-fold dengan menetapkan nilai k=n. Atau bisa juga hanya mengambil program inti agar lebih sederhana.
cv2 <- NULL
for(j in 1:n){
#menduga model
regresiSwiss <- lm(Fertility~.,data = swiss[-j,])
# menduga respon
responseSwiss <- predict(regresiSwiss,newdata=swiss[j,])
# menghitung galat prediksi
galat <- swiss$Fertility[j] - responseSwiss
cv2 <- c(cv2,1/(n*k)*sqrt(sum(galat^2)))
}
cv2## [1] 0.0281845164 0.0028616702 0.0338942605 0.0415164841 0.0559431792
## [6] 0.0764369943 0.0203921731 0.0535268370 0.0049220519 0.0162970146
## [11] 0.0171026663 0.0235181736 0.0023376615 0.0137463865 0.0179616749
## [16] 0.0242204121 0.0004833863 0.0009529467 0.0233628193 0.0079161546
## [21] 0.0153134215 0.0501196972 0.0226890718 0.0321692177 0.0050210089
## [26] 0.0082059576 0.0027083520 0.0068050767 0.0244371035 0.0323969395
## [31] 0.0039094563 0.0359695286 0.0060425348 0.0275140671 0.0180155765
## [36] 0.0383084947 0.0760281997 0.0386535480 0.0221452050 0.0348836887
## [41] 0.0195380342 0.0582634572 0.0232653423 0.0280304009 0.0015825468
## [46] 0.0520888356 0.0735220543RMSEP2 <- mean(cv2)
RMSEP2## [1] 0.02602562Praktikum Pertemuan 14 - Metode Monte Carlo
library(ggplot2)
library(mvtnorm)
library(MBESS)Simulasi Monte Carlo yang memanfaatkan informasi mengenai sebaran data yang diketahui (dihipotesiskan, dianggap tahu) dengan pasti.
Pendugaan Peluang:
P(Y>a)atauP(Y<a)
dengan Y=g(x),X f(x)
Latihan 1
Seorang peneliti ingin melihat bagaimana kesalahan tipe 1 (α) dari uji t (t test) dan uji peringkat bertanda Wilcoxon (Wilcoxoan signed-rank test) [https://sixsigmastudyguide.com/1-sample-wilcoxon-non-parametric-hypothesis-test/] untuk satu populasi pada berbagai ukuran sampel. Hipotesis yang akan diujikan adalah sebagai berikut:
H0:μ=1500
H1:μ≠1500
Bangkitan data berdasarkan distribusi N(μ=1500,σ=200), α=0.05 dan ulangan sebanyak 10.000. Ukuran sampel yang diteliti adalah 5 sampai 50. Berilah kesimpulan atas hasilnya!
Jawaban:
Langkah 1 ada mendefinisikan parameter-parameter yang digunakan untuk simulasi
# mendefinisikan ukuran sampel dari 5 saampai 50
n <- seq(5,50)
# mendefinisikan miu distribusi normal
mu0 <- 1500
# mendefinisikan sigma distribusi normal
sigma0 <- 200
# mendefinisikan banyaknya ulangan
num_rep <- 10000
# mendefinisikan alpha
alpha <- 0.05Langkah 2 medefinisikan fungsi untuk mengextract p-value dari masing-masing uji (test).
t_test_pvalue <- function(x,alt="two.sided",mu=mu0){
t.test(x, alternative = alt, mu = mu)$p.value
}
wilcox_test_pvalue <- function(x,alt="two.sided",mu=mu0){
wilcox.test(x, alternative = alt, mu = mu)$p.value
}Langkah 3 Menjalankan Simulasi
# membuat objek data.frame dummy untuk menyimpan hasil simulasi
result_sim1 <- data.frame(n=0,alpha_empiris=0,test_name="")
set.seed(2020)
# iterasi simulasi berdasarkan ukuran sampel
for(i in seq_along(n)){
# membangkitan data berdistribusi normal dengan ulangan 10rb
sim_mc_norm <- replicate(num_rep, rnorm(n[i],mu0,sigma0))
# menghitung p value uji t dari 10rb data hasil bangkitan
t_test_res <- apply(sim_mc_norm,2,t_test_pvalue)
# menghitung p value uji wilcoxon dari 10rb data hasil bangkitan
wilcox_test_res <- apply(sim_mc_norm,2,wilcox_test_pvalue)
# menghitung alpha empiris dari uji t
t_p_value_sim <- mean(t_test_res<alpha)
# menghitung alpha empiris dari uji wilcoxon
wilcox_p_value_sim <- mean(wilcox_test_res<alpha)
# menyimpan hasil simulasi
result_sim1 <- rbind(result_sim1,list(n[i],t_p_value_sim,"t_test"),
list(n[i],wilcox_p_value_sim,"wilcox_test")
)
}
# membuang baris dummy
result_sim1 <- result_sim1[-1,]Langkah 4 Membuat Grafik dari hasil simulasi
ggplot(data = result_sim1,aes(n,alpha_empiris))+
geom_point(aes(color=test_name))+geom_hline(yintercept = 0.05)+scale_y_continuous(limits = c(0,0.1))
Latihan 2
Seorang peneliti ingin melihat bagaimana nilai koefisien korelasi dari metode korelasi pearson dan metode korelasi spearman untuk satu populasi pada berbagai ukuran sampel jika datanya terdapat pencilan. Bangkitan data pencilan dengan cara berikut:
Bangkitkan data berdasarkan distribusi multivariate Normal dengan vektor rataan (10,12) dan korelasi 0.85 serta vektor simpangan baku (2,3)
Bangkitan outlier pada kedua peubah dengan berdasarkan distribusi N(17,2) sebanyak 20% dari ukuran sampel pada langkah a.
Ganti amatan dari hasil bangkitan pada langkah a dengan hasil bangkitan pada langkah b secara acak.
Bangkitkan data dengan ulangan 10000 dan ukuran sampel yang diteliti adalah 10,20,30,40,50,60,70,80,90,100. Berikan kesimpulan dari hasil simulasi!
Jawaban
Langkah 1 abstraksi bangkitan data
langkah a
set.seed(1212)
# mendefinisikan vektor simpangan baku
sd0 <- c(2,3)
# mendefinisikan matrix korelasi
cor_mat <- matrix(c(1,0.85,0.85,1),nrow=2)
# mengubah matrix korelasi ke kovarians
cov_mat <- MBESS::cor2cov(cor_mat,sd0)
# mendefinisikan vektor rataan
mu_mvn <- c(10,12)
# membangkitkan data berdistribus multivariate normal
x_mvn <- rmvnorm(n=40, mean=c(10,12), sigma=cov_mat)langkah b
set.seed(1212)
#membangkitkan data berdasarkan distribusi normal
x1_out <- rnorm(8,17,2)
y1_out <- rnorm(8,17,2)
x1y1_out <-cbind(x1_out,y1_out)langkah c
# memilih secara acak amatan yang akan diganti dengan pencilan
index_out <- sample.int(nrow(x_mvn),8)
# mengganti amatan dari langkah a dengan amatan dari langkah b
x_mvn[index_out,] <- x1y1_outqplot(x = V1,y=V2,data = as.data.frame(x_mvn))
Langkah 2 ada mendefinisikan parameter-parameter yang digunakan untuk simulasi
# mendefinisikan ukuran sampel yang dikaji
n <- seq(10,100,10)
# mendefinisikan jumlah ulangan
num_rep2 <- 10000
# mendefinisikan vektor simpangan baku
sd0 <- c(2,3)
# mendefinisikan matrix korelasi
cor_mat <- matrix(c(1,0.85,0.85,1),nrow=2)
# mengubah matrix korelasi ke kovarians
cov_mat <- MBESS::cor2cov(cor_mat,sd0)
# mendefinisikan vektor rataan
mu_mvn <- c(10,12)
# mendefinisikan mu dan sigma untuk data pencilan
mu1 <- 17
sigma1 <- 2Langkah 3 mendefinisikan fungsi-fungsi yang dibutuhkan
# generate multivariate normal dengan pencilan
gen_mvn_out <- function(n0,mu0,sigma0,n1,mu1,sigma1){
x_mvn <- rmvnorm(n=n0, mean=mu0, sigma=sigma0)
x1_out <- rnorm(n1,mu1,sigma1)
y1_out <- rnorm(n1,mu1,sigma1)
x1y1_out <-cbind(x1_out,y1_out)
index_out <- sample.int(nrow(x_mvn),n1)
x_mvn[index_out,] <- x1y1_out
return(x_mvn)
}
# extract koefisien korelasi
cor_res <- function(x_mat,method){
cor(x_mat[,1],x_mat[,2],method = method)
}Langkah 4 Menjalankan simulasi
set.seed(1212)
# membuat objek data.frame dummy untuk menyimpan hasil simulasi
result_sim2 <- data.frame(n=0,correlation=0,type="")
for(i in seq_along(n)){
# membangkitan data berdistribusi multivariate normal dengan ulangan 10rb
sim_mc_mvn <- replicate(num_rep2, gen_mvn_out(n0 = n[i],mu_mvn,cov_mat,n1=0.2*n[i],mu1,sigma1),
simplify = F)
# menghitung keofisien korelasi pearson dari 10rb data hasil bangkitan
pearson_res <- sapply(sim_mc_mvn,cor_res,method="pearson")
# menghitung keofisien spearman pearson dari 10rb data hasil bangkitan
spearman_res <- sapply(sim_mc_mvn,cor_res,method="spearman")
# menghitung rata-rata koefisien
pearson_sim <- mean(pearson_res)
spearman_sim <- mean(spearman_res)
# menyimpan hasil simulasi
result_sim2 <- rbind(result_sim2,list(n[i],pearson_sim,"pearson"),
list(n[i],spearman_sim,"spearman")
)
}
result_sim2 <- result_sim2[-1,]Langkah 5 Membuat Grafik berdsarkan hasil simulasi
ggplot(data = result_sim2,aes(n,correlation))+
geom_point(aes(color=type))+geom_hline(yintercept = 0.85)+scale_y_continuous(limits = c(0.5,1))
Latihan 3
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
- Bangkitkan residual berdasarkan distribusi normal ϵ∼N(0,σ)
#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 1)2.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.59937298Pembangkitan peubah penjelas Xi dapat dilakukan juga dengan fungsi mvnorm jika ingin melakukan simulasi multikolinearitas.
- 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-16Berdasarkan 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 4
Berdasarkan Latihan 3 lakukan pengulangan sebanyak 100 ulangan dan kemudian hitung rata-rata dan simpangan baku dari masing-masing koefisien regresi.
Jawaban
#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.02379381Demikian, Terima Kasih
Satria June Adwendi, IPB University, sjadwendi@apps.ipb.ac.id↩︎