Bootstrap
Misalkan \(x=(x_{1},x_{2},\ldots,x_{n})\) adalah sampel/contoh yang digunakan untuk menduga suatu parameter \(\theta\) dengan menggunakan penduga \(\hat{\theta}\). Prosedur bootstraping berdasarkan kondisi ini adalah sebagai berikut :
- Bangkitan sampel baru \(x^{*}\) yang berukuran \(n\) dengan menggunakan metode penarikan sampel/contoh dengan pengembalian (with replacement) sebanyak B kali.
- Hitung penduga \(\hat{\theta}^{*}\) untuk setiap \(x^{*}\), sehingga diperoleh \(\tilde{\theta}^{*}=(\hat{\theta}_{1}^{*},\hat{\theta}_{2}^{*},\ldots, \hat{\theta}_{B}^{*})\)
- Gunakan \(\tilde{\theta}^{*}\) untuk menghitung mean \(\bar{\tilde{\theta}}^{*} = \frac{1}{B} \Sigma_{b=1}^{B}\hat{\theta}_{b}^{*}\) dan ragam \(V(\tilde{\theta}^{*}) = \frac{1}{B-1} \Sigma_{b=1}^{B}(\theta_{b}-\bar{\tilde{\theta}}^{*})^{2}\) dari hasil resampling. Kemudian, penduga bootstrap \[\hat{\theta}^{(B)}= \bar{\tilde{\theta}}\] dan ragam dari penduga bootstrap \[V(\hat{\theta}^{(B))}=V(\tilde{\theta}^{*})\].
Jackknife
Misalkan \(x=(x_{1},x_{2},\ldots,x_{n})\) adalah sampel/contoh yang digunakan untuk menduga suatu parameter \(\theta\) dengan menggunakan penduga \(\hat{\theta}\). Prosedur Jacknife berdasarkan kondisi ini adalah sebagai berikut :
- Bangkitan sampel baru \(x^{*}\) yang berukuran \(n-1\) dengan menghilangkan satu amatan pada setiap ulangan. Sehingga banyaknya ulangan yang dapat dilakukan adalah sebanyak \(n\).
- Hitung penduga \(\hat{\theta}^{*}\) untuk setiap \(x^{*}\), sehingga diperoleh \(\tilde{\theta}^{*}=(\hat{\theta}_{1}^{*},\hat{\theta}_{2}^{*},\ldots, \hat{\theta}_{B}^{*})\)
- Gunakan \(\tilde{\theta}^{*}\) untuk menghitung mean \(\bar{\tilde{\theta}}^{*} = \frac{1}{n} \Sigma_{j=1}^{n}\hat{\theta}_{n}^{*}\) dan ragam \(V(\tilde{\theta}^{*}) = \frac{1}{n-1} \Sigma_{n=1}^{n}(\theta_{b}-\bar{\tilde{\theta}}^{*})^{2}\) dari hasil resampling. Kemudian, penduga jacknife \[\hat{\theta}^{(Jack)}= n \hat{\theta}-(n-1)\bar{\tilde{\theta}}\] dan ragam dari penduga jacknife \[V(\hat{\theta}^{(Jack)})=\frac{(n-1)^2}{n}V(\tilde{\theta}^{*})\]
Validasi silang (CV)
- Bagi amatan data menjadi k bagian. Misal k=10 dan n=250 maka tiap bagian berisi 25.
- Bagian ke-k digunakan untuk menduga respons, sedangkan bagian lainnya digabungkan untuk menduga model.
- Hitung galat prediksi: \(\text{cv}(-k) = \frac{1}{nk} \sqrt{\Sigma(\hat{y}-y)^{2}}\)
- Ulangi sampai setiap bagian pernah digunakan untuk menduga galat prediksi
- Hitung RMSEP
Sebagai ilustrasi terdapat 10 bagian yang tiap bagian berisi 25 amatan. Kemudian, misalkan saja bagian ke-6 digunakan untuk menduga respons dan galat prediksi, sehingga bagian 1-5 dan 7-10 digabungkan untuk selanjutnya digunakan menduga model. Begitu seterusnya sampai setiap bagian pernah digunakan untuk menduga respon dan galat prediksi.
Latihan
Latihan 1
Tabel di bawah ini melaporkan jumlah populasi sebanyak \(n=49\) amatan dari sepasang kota di US pada tahun 1920 dan 1930, yang dinotasikan dengan \(u\) dan \(x\).
| 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 |
Berdasarkan data ini, peniliti tertarik untuk menghitung rasio dari mean (rata-rata/rataan) yang nantinya dapat digunakan untuk menduga populasi total USA dalam tahun 1930 dari 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 \(\theta = \frac{E(X)}{E(U)}\), yang merupakan parameter rasio yang akan diduga. Dengan menggunakan bootstrap hitunglah galat baku (standard error) dari \(\hat{\theta}\), yang merupakan \(\sqrt{V(\hat{\theta})}\), jika diketahui \(\hat{\theta}=\frac{\bar{x}}{\bar{u}}\) !
Jawaban
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 membuat sintaks 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], ]))
}
length(thetaBoost)## [1] 1000head(thetaBoost)## [1] 1.332714 1.189381 1.248475 1.251276 1.323737 1.228723- Gunakan hasil pada b untuk menghitung ragam \(V(\hat{\theta})\)
Ragam \(V(\hat{\theta})\) bisa didekati menggunakan ragam bootstrap \(V(\theta^{(B)})\).
# ragam
ragamTheta <- var(thetaBoost)
ragamTheta## [1] 0.001208174- 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.036673Latihan 2
Gunakanlah metode jackknife pada soal Latihan 1
Jawaban
- 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], ]))
}
length(thetaJack)## [1] 49head(thetaJack)## [1] 1.244711 1.241282 1.240336 1.231179 1.235917 1.235599- Gunakan hasil pada b untuk menghitung ragam \(V(\hat{\theta})\)
Ragam \(V(\hat{\theta})\) bisa didekati menggunakan ragam jacknife \(V(\hat{\theta}^{(Jack)})\).
# ragam
ragamTheta2 <- (((n-1)^2)/n) * var(thetaJack)
ragamTheta2## [1] 0.001192619# standard error
sqrt(ragamTheta2)## [1] 0.03453431Latihan 3
Gunakanlah data swiss yang ada di R, 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
sisanya merupakan peubah prediktor!.
Jawaban
Importing data
data(swiss)- 5-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 47class(indexfolds)## [1] "list"- Bagian ke-k digunakan untuk menduga respons, sedangkan bagian lainnya digabungkan untuk menduga model.
- Hitung galat prediksi: \(\text{cv}(-k) = \frac{1}{nk} \sqrt{\Sigma(\hat{y}-y)^{2}}\)
- Ulangi sampai setiap bagian pernah digunakan untuk menduga galat prediksi
cv <- 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
cv <- c(cv, 1/(n * k) * sqrt(c(cv, sum(galat ^ 2))))
}- Hitung RMSEP
RMSEP <- mean(cv)
RMSEP## [1] 0.01718263Untuk LOO-CV tetapkan nilai k=n
Statistika dan Sains Data, IPB University, alfanugraha@apps.ipb.ac.id↩︎