Tugas Praktikum Kelima STA561 Pemrograman Statistika
Klik disini untuk ke halaman rpubs.
Komputasi untuk Regresi (Praktikum Pertemuan Ke-11)
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
Diketahui data berikut
Lakukan prosedur Sweep
Jawab
#install.packages("fastmatrix")
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.3333333
Dekomposisi 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.083333Expectation-Maximization (EM) Algorithm (Praktikum Pertemuan Ke-12)
Algoritma EM
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 untuk missing value dalam rancop
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, toleransi = 1e-4, maksimal iterasi=100 (hasil= 26.99986)
emrancob(ftn, x,0,2,3, tol = 1e-4, max.iter = 100) #30x iterasi 26.99986## 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.99986emrancob(ftn, x,17.4,2,3, tol = 1e-4, max.iter = 100) #27x iterasi 26.99983## At iteration 1 value of x is: 20.6
## At iteration 2 value of x is: 22.73333
## At iteration 3 value of x is: 24.15556
## At iteration 4 value of x is: 25.1037
## At iteration 5 value of x is: 25.7358
## At iteration 6 value of x is: 26.1572
## At iteration 7 value of x is: 26.43813
## At iteration 8 value of x is: 26.62542
## At iteration 9 value of x is: 26.75028
## At iteration 10 value of x is: 26.83352
## At iteration 11 value of x is: 26.88901
## At iteration 12 value of x is: 26.92601
## At iteration 13 value of x is: 26.95067
## At iteration 14 value of x is: 26.96712
## At iteration 15 value of x is: 26.97808
## At iteration 16 value of x is: 26.98538
## At iteration 17 value of x is: 26.99026
## At iteration 18 value of x is: 26.9935
## At iteration 19 value of x is: 26.99567
## At iteration 20 value of x is: 26.99711
## At iteration 21 value of x is: 26.99808
## At iteration 22 value of x is: 26.99872
## At iteration 23 value of x is: 26.99914
## At iteration 24 value of x is: 26.99943
## At iteration 25 value of x is: 26.99962
## At iteration 26 value of x is: 26.99975
## At iteration 27 value of x is: 26.99983
## Algorithm converged## [1] 26.99983Contoh 2
Pendugaan parameter pada data genetika (Rao 1973), diasumsikan data phenotype y = (y1, y2, y3, y4) = (125, 18, 20, 34) memiliki distribusi multinomial dengan peluang
### Contoh_2
ftn2 <- function(x,yij,i,j){
x <- matrix(c(12,yij,18,22,23,24),ncol=3,byrow=TRUE)
mu <- mean(x)
alfai <- mean(x[i,])- mu
betaj <- mean(x[,j])- mu
yij <- mu + alfai + betaj
return(yij)
}
### Cari y12, dengan nilai awal 16.5, i=2, j=3
emrancob(ftn2, x,16.5,1,2, tol = 1e-4, max.iter = 100)## At iteration 1 value of x is: 16
## At iteration 2 value of x is: 15.66667
## At iteration 3 value of x is: 15.44444
## At iteration 4 value of x is: 15.2963
## At iteration 5 value of x is: 15.19753
## At iteration 6 value of x is: 15.13169
## At iteration 7 value of x is: 15.08779
## At iteration 8 value of x is: 15.05853
## At iteration 9 value of x is: 15.03902
## At iteration 10 value of x is: 15.02601
## At iteration 11 value of x is: 15.01734
## At iteration 12 value of x is: 15.01156
## At iteration 13 value of x is: 15.00771
## At iteration 14 value of x is: 15.00514
## At iteration 15 value of x is: 15.00343
## At iteration 16 value of x is: 15.00228
## At iteration 17 value of x is: 15.00152
## At iteration 18 value of x is: 15.00101
## At iteration 19 value of x is: 15.00068
## At iteration 20 value of x is: 15.00045
## At iteration 21 value of x is: 15.0003
## At iteration 22 value of x is: 15.0002
## At iteration 23 value of x is: 15.00013
## Algorithm converged## [1] 15.00013### Pendugaan data tersembunyi (sebaran multinomial)
### Contoh_1
x3<- 18
x4<- 20
x5<- 34
x2d1<- 0
x2d2<- 0
theta0<- 2
err<- 10
iterasi<- 1
while(err>10^-4){
theta1<- (x2d1+x5)/(x2d1+x3+x4+x5)
x1d2<- 125*(0.5/(0.5+0.25*theta1))
x2d2<- 125*((0.25*theta1)/(0.5+0.25*theta1))
err<- abs(theta1-theta0)
theta0<- theta1
x2d1<-x2d2
iterasi<- iterasi + 1
}
x1d2 # x1## [1] 95.17232x2d2 # x2## [1] 29.82768theta1## [1] 0.6268142iterasi## [1] 7Teladan 1
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.82795Resampling (Bootstrap dan Jackknife)(Praktikum Pertemuan Ke-13)
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.
Bootstrap
Bootstrap, diperkenalkan oleh Efron (1979), dapat melakukan pencarian informasi tentang parameter populasi berdasarkan hanya suatu data contoh dari populasi tersebut.Istilah ‘bootstrap’ dari “To pull oneself up by one’s bootstraps”.Bootstrap akan membentuk fungsi sebaran empirik sebagai gambaran fungsi sebaran sebenarnya.Tujuan dari bootstrap adalah mengumpulkan informasi secara cepat dan lebih murah karena tidak langsung dari populasi tetapi berdasarkan suatu data contoh saja yang dianggap populasi.
Bootstrap terutama digunakan untuk mempelajari suatu penduga parameter populasi, yang disebut statistik, berdasarkan data contoh.Bootstrap sering digunakan bila sebaran populasinya tidak diketahui dan hanya data contoh yang tersedia sebagai populasi semu (pseudo) dengan karakteristik yang serupa populasi sebenarnya.
Misalkan 𝜃duga adalah statistik contoh sebagai penduga parameter populasi yang diperoleh berdasarkan data contoh tersebut. Bila n besar (n ≥30), sebaran penarikan contohnya berbentuk simetris terhadap 𝜃; bila kasus n kecil, sebaran sulit diperoleh apalagi untuk statistik seperti median, korelasi, atau koefisien keragaman. Pendugaan suatu sebaran dengan bootstrap analog dengan pendugaan kepekatan. Bootstrap membangkitkan contoh-contoh acak dari suatu sebaran empirik. Suatu histogram dibentuk berdasarkan contoh-contoh tersebut untuk memperoleh suatu bentuk penduga fungsi kepekatan dan dipandang sebagai penduga kepekatan yang layak bagi sebaran empirik.
Data berukuran kecil dianggap sebagai populasi terhingga dan proses penarikan contoh dengan pengembalian (replacement) dilakukan secara berulang (resampling). Penarikan contoh dilakukan secara acak sebagaimana penarikan contoh acak dari populasi besar. Penduga parameter diperoleh dari setiap contoh acak sehingga pengulangan sebanyak B kali akan menghasilkan B penduga, yaitu 𝜃duga1, 𝜃duga2, … , 𝜃duga𝐵.Rata-rata dari B buah penduga, yaitu 𝜃dugabar.
Teladan 1
Menduga rata-rata, standard error, bias, dan selang kepercayaan bagi rata-rata
x <-c(8.26, 6.33, 10.4, 5.27, 5.35, 5.61, 6.12, 6.19, 5.2, 7.01, 8.74, 7.78, 7.02, 6, 6.5, 5.8, 5.12, 7.41, 6.52, 6.21, 12.28, 5.6, 5.38, 6.6, 8.74)
x <- as.matrix(x)
x## [,1]
## [1,] 8.26
## [2,] 6.33
## [3,] 10.40
## [4,] 5.27
## [5,] 5.35
## [6,] 5.61
## [7,] 6.12
## [8,] 6.19
## [9,] 5.20
## [10,] 7.01
## [11,] 8.74
## [12,] 7.78
## [13,] 7.02
## [14,] 6.00
## [15,] 6.50
## [16,] 5.80
## [17,] 5.12
## [18,] 7.41
## [19,] 6.52
## [20,] 6.21
## [21,] 12.28
## [22,] 5.60
## [23,] 5.38
## [24,] 6.60
## [25,] 8.74mean(x) # penduga rata-rata## [1] 6.8576b <- 5000
boot <- numeric(b)
n <- nrow(x)
for (i in 1:b) boot[i] <- mean(sample(x,n,replace=T))
mean(boot)## [1] 6.852311sqrt(var(boot)) #penduga salah baku bagi rata-rata## [1] 0.3339427bias<- mean(boot)-mean(x) #penduga bias
bias## [1] -0.00528936hist(boot)
# selang kepercayaan bagi rata-rata
quantile(boot,0.025)## 2.5%
## 6.26558quantile(boot,0.975)## 97.5%
## 7.57241Teladan 2
Menduga rata-rata, standard error, dan selang kepercayaan bagi median
### Teladan: Penduga Median dan Salah Baku
x <-c(8.26, 6.33, 10.4, 5.27, 5.35, 5.61, 6.12, 6.19, 5.2, 7.01, 8.74, 7.78, 7.02, 6, 6.5, 5.8, 5.12, 7.41, 6.52, 6.21, 12.28, 5.6, 5.38, 6.6, 8.74)
x <- as.matrix(x)
median(x) # penduga median## [1] 6.33n <- nrow(x)
b <- 5000
boot2 <-numeric(b)
for (i in 1:b) boot2[i] <- median(sample(x,n,replace=T))
mean(boot2)## [1] 6.370448sqrt(var(boot2)) #penduga salah baku bagi median## [1] 0.2841701hist(boot2)
# selang kepercayaan bagi median
quantile(boot2,0.025)## 2.5%
## 5.8quantile(boot2,0.975)## 97.5%
## 7.02Latihan 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.
# 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
')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, hitunglah galat baku (standard error) dari θduga, yang merupakan akar dari V(θduga), jika diketahui θduga=x_bar/u_bar!
Jawab
Cara 1 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(θduga) Ragam V(θduga) bisa didekati menggunakan ragam bootstrap V(θ^(B)).
# ragam
ragamTheta <- var(thetaBoost)
ragamTheta## [1] 0.001208174# galat baku
sqrt(ragamTheta)## [1] 0.03475879Cara 2 dengan package boot()
library(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.036673Jackknife
Pada awalnya metode jackknife atau leave-one-out bertujuan untuk menduga bias dari suatu penduga yang diperkenalkan oleh Maurice Quenouille pada tahun 1949 dan dikembangkan oleh John Tukey pada tahun 1957, yang dikenal dengan Quenouille-Tukey Jackknife.Metode jackknife serupa dengan metode bootstrap untuk salah baku. Metode ini seperti tipe validasi silang leave-one-out. Bila 𝒙=(𝑥1, 𝑥2, …, 𝑥𝑛) adalah contoh acak dan 𝒙((𝒊)) adalah subset dari 𝒙 tanpa pengamatan 𝑥𝑖, yaitu 𝒙((𝒊))=(𝑥1, …,𝑥((𝑖−1) ),𝑥((𝑖+1) ), …, 𝑥𝑛). Jika 𝜃duga=𝑓𝑛 (𝒙), maka definisikan bahwa 𝜃duga=𝑓(𝑛−1) ( 𝒙((𝒊) ) ), 𝑖=1, 2, …, 𝑛.
Teladan 1
x <-c(8.26, 6.33, 10.4, 5.27, 5.35, 5.61, 6.12, 6.19, 5.2, 7.01, 8.74, 7.78, 7.02, 6, 6.5, 5.8, 5.12, 7.41, 6.52, 6.21, 12.28, 5.6, 5.38, 6.6, 8.74)
x <- as.matrix(x)
n <- nrow(x)
theta.hat <- median(x)
print(theta.hat)## [1] 6.33# pendugaan median dengan jackknife (leave-one-out)
theta.jack <- numeric(n)
for (i in 1:n) theta.jack[i] <- median(x[-i])
medjack <- mean(theta.jack)
print(medjack) # penduga median dengan Jackknife## [1] 6.343# penduga bias
bias <- (n-1) * (mean(theta.jack) - theta.hat)
print(bias)## [1] 0.312
Teladan 2
x <-c(8.26, 6.33, 10.4, 5.27, 5.35, 5.61, 6.12, 6.19, 5.2, 7.01, 8.74, 7.78, 7.02, 6, 6.5, 5.8, 5.12, 7.41, 6.52, 6.21, 12.28, 5.6, 5.38, 6.6, 8.74)
x <- as.matrix(x)
n <- nrow(x)
theta.hat <- median(x)
print(theta.hat)## [1] 6.33# pendugaan median dengan jackknife (leave-one-out)
theta.jack <- numeric(n)
for (i in 1:n) theta.jack[i] <- median(x[-i])
medjack <- mean(theta.jack)
print(medjack) # penduga median dengan Jackknife## [1] 6.343# penduga salah baku
se <- sqrt((n-1)*mean(((theta.jack)-mean(theta.jack))^2))
print(se)## [1] 0.3482068
x <-c(8.26, 6.33, 10.4, 5.27, 5.35, 5.61, 6.12, 6.19, 5.2, 7.01, 8.74, 7.78, 7.02, 6, 6.5, 5.8, 5.12, 7.41, 6.52, 6.21, 12.28, 5.6, 5.38, 6.6, 8.74)
x <- as.matrix(x)
n <- length(x)
jack <- numeric(n-1)
pseudo <- numeric(n)
for (i in 1:n)
{for (j in 1:n) {if(j<i) jack[j]<-x[j] else if(j>i) jack[j-1]<- x[j]}
pseudo[i]<-n*median(x)-(n-1)*median(jack)}
mean(pseudo) #penduga rata2## [1] 6.018var(pseudo) #penduga ragam## [1] 3.0312var(pseudo)/n #ragam dari penduga rata2## [1] 0.121248# Penduga selang kepercayaan 95% bagi median #
mean(pseudo) + qt(0.975,n-1)*sqrt(var(pseudo)/n)## [1] 6.736664mean(pseudo) - qt(0.975,n-1)*sqrt(var(pseudo)/n)## [1] 5.299336Latihan 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.
# 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
')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 JACKNIFE, hitunglah galat baku (standard error) dari θduga, yang merupakan akar dari V(θduga), jika diketahui θduga=x_bar/u_bar!
Jawab
- 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(θduga). Ragam V(θduga) bisa didekati menggunakan ragam jacknife V(θ^(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.03453431Validasi Silang (Cross Validation-CV)
Validasi silang adalah metode partisi data yang dapat digunakan untuk mengetahui antara lain kestabilan penduga parameter, kesesuaian penduga model. Jackknife dapat dipandang sebagai kasus khusus dari validasi silang yang pada awalnya digunakan untuk pendugaan bias dan salah baku.
Ilustrasi 1
Pendugaan galat prediksi model dengan validasi silang
Andaikan ada empat model untuk memprediksi Y dari X. Linier : 𝑌=𝛽0+𝛽1 𝑋+𝜀 Kuadratik : 𝑌=𝛽0+𝛽1 𝑋+𝛽2 𝑋^2+𝜀 Eksponensial : log(𝑌)=log(𝛽_0)+𝛽_1 𝑋+𝜀 Log-Log : log(𝑌)=𝛽_0+𝛽_1 log(𝑋)+𝜀
## Data:
par(mfrow=c(2,2))
library(DAAG); attach(ironslag)## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:boot':
##
## melanomaa <- seq(10, 40, .1) #sequence for plotting fits
# Model linier
L1 <- lm(magnetic ~ chemical)
plot(chemical, magnetic, main="Linear", pch=16)
yhat1 <- L1$coef[1] + L1$coef[2] * a
lines(a, yhat1, lwd=2)
# Model Kuadratik
L2 <- lm(magnetic ~ chemical + I(chemical^2))
plot(chemical, magnetic, main="Quadratic", pch=16)
yhat2 <- L2$coef[1] + L2$coef[2] * a + L2$coef[3] * a^2
lines(a, yhat2, lwd=2)
# Model Exponensial
L3 <- lm(log(magnetic) ~ chemical)
plot(chemical, magnetic, main="Exponential", pch=16)
logyhat3 <- L3$coef[1] + L3$coef[2] * a
yhat3 <- exp(logyhat3)
lines(a, yhat3, lwd=2)
# Model LOG-LOG
L4 <- lm(log(magnetic) ~ log(chemical))
plot(log(chemical), log(magnetic), main="Log-Log", pch=16)
logyhat4 <- L4$coef[1] + L4$coef[2] * log(a)
lines(log(a), logyhat4, lwd=2) 
n <- length(magnetic) #in DAAG ironslag
e1 <- e2 <- e3 <- e4 <- numeric(n)
for (k in 1:n) {
y <- magnetic[-k]
x <- chemical[-k]
J1 <- lm(y ~ x) # model linier
yhat1 <- J1$coef[1] + J1$coef[2] * chemical[k]
e1[k] <- magnetic[k] - yhat1
J2 <- lm(y ~ x + I(x^2)) # model kuadratik
yhat2 <- J2$coef[1] + J2$coef[2] * chemical[k] + J2$coef[3] * chemical[k]^2
e2[k] <- magnetic[k] - yhat2
J3 <- lm(log(y) ~ x) # model eksponensial
logyhat3 <- J3$coef[1] + J3$coef[2] * chemical[k]
yhat3 <- exp(logyhat3)
e3[k] <- magnetic[k] - yhat3
J4 <- lm(log(y) ~ log(x)) # model log-log
logyhat4 <- J4$coef[1] + J4$coef[2] * log(chemical[k])
yhat4 <- exp(logyhat4)
e4[k] <- magnetic[k] - yhat4
}
c(mean(e1^2), mean(e2^2) , mean(e3^2), mean(e4^2)) # model kuadratik dengan galat prediksi paling kecil ## [1] 19.55644 17.85248 18.44188 20.45424Latihan 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!
Jawab:
# 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)*akar(Σ(y^−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.02602562Metode Monte Carlo (Praktikum Pertemuan Ke-14)
Mencakup sejumlah alat komputasi dalam statistika terapan, antara lain integrasi Monte Carlo. Mengacu pada metode penarikan kesimpulan (statisical inference) atau analisis numerik yang menggunakan simulasi. Ada ketidak pastian (uncertanty) pendugaan pada penarikan kesimpulan. Penarikan contoh secara berulang dari suatu sebaran peluang, seperti bootstrap parametrik, untuk mengetahui ketidak pastian. Penarikan contoh secara berulang dari data pengamatan, seperti bootstrap (nonparametrik). Dapat digunakan untuk pendugaan parameter, kuadrat tengah galat (MSE), peluang selang kepercaayaan dan salah jenis I.
Penduga bagi beda dua nilai tengah dan salah baku
m <- 1000
g <- numeric(m)
for (j in 1:m) {
x <- rnorm(2)
g[j] <- abs(x[1] - x[2])
}
est <- mean(g)
est## [1] 1.165277
m <- 1000
g <- numeric(m)
for (i in 1:m) {
x <- rnorm(2)
g[i] <- abs(x[1] - x[2])
}
estse <- sqrt(sum((g - mean(g))^2)) / m
estse## [1] 0.02594272Penduga bagi Kuadrat Tengah Galat (MSE)
Teladan
# trimmed mean untuk nilai tengah (rata-rata)
n <- 20
m <- 1000
tmean <- numeric(m) # trimmed mean
for (j in 1:m) {
x <- sort(rnorm(n)) # sebaran normal
tmean[j] <- sum(x[2:(n-1)]) / (n-2)
}
mse <- mean(tmean^2)
mse## [1] 0.04781138sqrt(sum((tmean - mean(tmean))^2)) / m #se## [1] 0.006905955# trimmed mean untuk median
# (median is actually a trimmed mean)
n <- 20
m <- 1000
tmean <- numeric(m) # trimmed mean
for (i in 1:m) {
x <- sort(rnorm(n))
tmean[i] <- median(x) # median (median berada di tengah)
}
mse <- mean(tmean^2)
mse## [1] 0.07614342sqrt(sum((tmean - mean(tmean))^2)) / m #se## [1] 0.008726004Penduga Bagi Tingkat Kepercayaan (Confidence Level)
# batas atas selang kepercayaan 95% bagi ragam
# contoh acak n=20 dari sebaran N ~ (0,4)
n <- 20
alpha <- 0.05
x <- rnorm(n, mean=0, sd=2)
UCL <- (n-1) * var(x) / qchisq(alpha, df=n-1) # upper confidence level
UCL## [1] 10.70307
# batas atas selang kepercayaan 95% bagi ragam
# contoh acak n=20 dari sebaran N ~ (0,4)
n <- 20
alpha <- .05
ucl<- numeric(10)
for (i in 1:10){x <- rnorm(n, mean=0, sd=2)
ucl[i]<-(n-1) * var(x) / qchisq(alpha, df=n-1)
}
ucl## [1] 10.217686 5.318907 5.391812 8.905902 5.570816 3.590928 10.480207
## [8] 8.909878 10.096666 8.047340
# µ = 0, σ = 2, n = 20, m = 1000 ulangan, dan α = 0.05
# Proporsi selang yang mencakup σ2 = 4 adalah penduga Monte Carlo
# bagi tingkat kepercayaan sebenarnya
# Versi 1
n <- 20
alpha <- .05
UCL <- replicate(1000, expr = {
x <- rnorm(n, mean = 0, sd = 2)
(n-1) * var(x) / qchisq(alpha, df = n-1)
} )
#count the number of intervals that contain sigma^2=4
sum(UCL > 4)## [1] 947#or compute the mean to get the confidence level
mean(UCL > 4)## [1] 0.947# µ = 0, σ = 2, n = 20, m = 1000 ulangan, dan α = 0.05
# Proporsi selang yang mencakup σ2 = 4 adalah penduga Monte Carlo
# bagi tingkat kepercayaan sebenarnya
# Versi 2
calcCI <- function(n, alpha) {
y <- rnorm(n, mean = 0, sd = 2)
return((n-1) * var(y) / qchisq(alpha, df = n-1))
}
UCL <- replicate(1000, expr = calcCI(n = 20, alpha = .05))
#count the number of intervals that contain sigma^2=4
sum(UCL > 4)## [1] 952#or compute the mean to get the confidence level
mean(UCL > 4)## [1] 0.952Metode Monte Carlo untuk Uji Hipotesis
Andaikan kita akan melakukan uji hipotesis berikut: H0: θ = θ0 vs H1: θ = θ1. Dalam pengujian hipotesis ada dua jenis kesalahan, yaitu salah jenis I yang terjadi jika H0 ditolak padahal H0 benar, dan salah jenis II yang terjadi jika H0 tidak ditolak padahal H0 salah.
n <- 20 # banyaknya data 𝑥_1, …, 𝑥_20
alpha <- 0.05
mu0 <- 500 # H0: µ = 500
sigma <- 100 # σ
m <- 10000 #number of replicates
p <- numeric(m) #storage for p-values
for (j in 1:m) {
x <- rnorm(n, mu0, sigma)
ttest <- t.test(x, alternative = "greater", mu = mu0)
p[j] <- ttest$p.value
}
p.hat <- mean(p < alpha) # proporsi
se.hat <- sqrt(p.hat * (1 - p.hat) / m) # salah baku
print(c(p.hat, se.hat))## [1] 0.046400000 0.002103498#install.packages("BSDA")
library(BSDA)##
## Attaching package: 'BSDA'## The following object is masked from 'package:datasets':
##
## Orangen <- 20 # banyaknya data 𝑥_1, …, 𝑥_20
alpha <- 0.05
mu0 <- 500 # H0: µ = 500
sigma <- 100 # σ
m <- 10000 #number of replicates
p <- numeric(m) #storage for p-values
for (j in 1:m) {
x <- rnorm(n, mu0, sigma)
ztest <- z.test(x, sigma.x=sigma, alternative = "greater", mu = mu0) # paket BSDA
p[j] <- ztest$p.value
}
p.hat <- mean(p < alpha) # proporsi
se.hat <- sqrt(p.hat * (1 - p.hat) / m) # salah baku
print(c(p.hat, se.hat))## [1] 0.047600000 0.002129184Metode Markov Chain Monte Carlo (MCMC)
Diperkenalkan oleh Metropolis dan Hastings untuk integrasi Monte Carlo. Integrasi untuk membangkitkan contoh-contoh (samples) dari fungsi kepekatan f(.).Pendekatan MCMC untuk penarikan contoh dari f(.) digunakan untuk mengkontruksi suatu Markov Chain dengan sebaran stasioner f(.), dan terus menerus untuk waktu yang lama sampai rantai konvergen (approximately) terhadap sebaran stasioner.Metode untuk mengkonstruksi rantai Markov antara lain Metropolis-Hastings dan Gibbs sampler.
#Sebaran Proposal Khi Kuadrat
# Fungsi untuk evaluasi kepekatan Rayleigh [versi a]
f <- function(x, sigma) {
if (any(x < 0)) return (0)
stopifnot(sigma > 0)
return((x / sigma^2) * exp(-x^2 / (2*sigma^2))) # sebaran target
}
# Pada simulasi ini suatu contoh dari sebaran Rayleigh dibangkitkan berdasarkan
# sebaran proposal khi kuadrat; pada setiap transisi, titik 𝑌 dibangkitkan dari 𝜒^2 (𝑑𝑓=𝑋_(𝑖−1) )
xt <- x[i-1]
y <- rchisq(1, df = xt)
m <- 10000
sigma <- 4
x <- numeric(m)
x[1] <- rchisq(1, df=1)
u <- runif(m)
for (i in 2:m) {
xt <- x[i-1]
y <- rchisq(1, df = xt)
num <- f(y, sigma) * dchisq(xt, df = y) # pembilang pada r(Xt,Y)
den <- f(xt, sigma) * dchisq(y, df = xt) # penyebut pada r(Xt,Y)
if (u[i] <= num/den) x[i] <- y
else x[i] <- xt
}
index <- 1:m
y1 <- x[index]
plot(index, y1, type="l", main="", ylab="x") # # plot sekuen Xt 
par(mfrow = c(1, 2))
y <- x[1:m] # xt yang diterima
a <- ppoints(100)
QR <- sigma * sqrt(-2 * log(1 - a)) # kuantil sebaran Rayleigh
Q <- quantile(x, a)
qqplot(QR, Q, main="", xlab="Rayleigh Quantiles", ylab="Sample Quantiles")
hist(y, breaks="scott", main="", xlab="", freq=FALSE)
lines(QR, f(QR, 4))
par(mfrow = c(1, 1))
mean(x) # penduga nilai tengah (E[X] ) ## [1] 4.963099var(x) # penduga ragam (Var) ̂[X]## [1] 6.992452sigma*sqrt(pi/2) # E[X]=σ√(π/2)## [1] 5.013257sigma^2*((4-pi)/2) # Var[X]=σ^2 (4-π)/2## [1] 6.867259#Sebaran Proposal Normal
m <- 50000
sigma <- 4
x <- numeric(m)
x[1] <- rnorm(1, 1, 1) # sebaran normal dengan mu=1 dan ragam1=1
u <- runif(m)
for (i in 2:m) {
xt <- x[i-1]
y <- rnorm(1, xt, 1)
num <- f(y, sigma) * dnorm(xt, y, 1) # pembilang pada r(X_t,Y)
den <- f(xt, sigma) * dnorm(y, xt, 1) # penyebut pada r(X_t,Y)
if (u[i] <= num/den) {x[i] <- y}
else {x[i] <- xt}
}
index <- 1:m
y1 <- x[index]
plot(index, y1, type="l", main="", ylab="x") # # plot sekuen Xt 
par(mfrow = c(1, 2))
y <- x[1:m]
a <- ppoints(100)
QR <- sigma * sqrt(-2 * log(1 - a)) # kuantil sebaran Rayleigh
Q <- quantile(x, a)
qqplot(QR, Q, main="", xlab="Rayleigh Quantiles", ylab="Sample Quantiles")
hist(y, breaks="scott", main="", xlab="", freq=FALSE)
lines(QR, f(QR, 4))
par(mfrow = c(1, 1))#Sebaran Proposal Khi Kuadrat
# khi kuadrat
m <- 50000
x <- numeric(m)
sigma <- 4
target <- function(x, sigma) {
if (any(x < 0)) return (0)
stopifnot(sigma > 0)
return((x / sigma^2) * exp(-x^2 / (2*sigma^2)))
}
sta561MCMC <- function(m, sigma, startval, proposalsd) {
x <- rep(0, m)
x[1] <- startval
for (i in 2:m) {
currentx <- x[i - 1]
proposedx <- rchisq(1, df = currentx)
num <- target(proposedx,sigma)*dchisq(currentx, df = proposedx)
den <- target(currentx,sigma)* dchisq(proposedx, df = currentx)
rasio <- num/den
if (runif(1) < rasio) {
x[i] <- proposedx
} else {
x[i] <- currentx
}
}
return(x)
}
z <- sta561MCMC(m,sigma,1,1)
plot(z, type = "l")
par(mfrow = c(1, 2))
y <- z[1:m]
a <- ppoints(100)
QR <- sigma * sqrt(-2 * log(1 - a)) # kuantil sebaran Rayleigh
Q <- quantile(z, a)
qqplot(QR, Q, main="", xlab="Rayleigh Quantiles", ylab="Sample Quantiles")
hist(y, breaks="scott", main="", xlab="", freq=FALSE)
lines(QR, f(QR, 4))
par(mfrow = c(1, 1))Latihan 1
library(ggplot2)
library(mvtnorm)
#install.packages("MBESS")
library(MBESS)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!
Jawab:
Langkah 1 ada mendefinisikan parameter-parameter yang digunakan untuk simulasi
# mendefinisikan ukuran sampel dari 5 sampai 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_out
qplot(x = V1,y=V2,data = as.data.frame(x_mvn))
Langkah 2 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
Jawab
- Bangkitkan residual berdasarkan distribusi normal ϵ∼N(0,σ)
#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 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.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.
Jawab
#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.02379381TERIMA KASIH
Referensi
Dito, Gerry Alfa dan Raharjo, Mulianto. 2021. EM Algorithm. Retrieved From https://newlms.ipb.ac.id
Dito, Gerry Alfa dan Raharjo, Mulianto. 2021. Metode Monte Carlo. Retrieved From https://newlms.ipb.ac.id
Dito, Gerry Alfa dan Raharjo, Mulianto. 2021. Resampling. Retrieved From https://newlms.ipb.ac.id
Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Expectation-Maximization (EM). Retrieved from https://newlms.ipb.ac.id/
Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Komputasi untuk Regresi. Retrieved from https://newlms.ipb.ac.id/
Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Metode Monte Carlo. Retrieved from https://newlms.ipb.ac.id/
Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Resampling (Bootstrap dan Jacknife). Retrieved from https://newlms.ipb.ac.id/
Wikibooks. 4 November 2018. Data Mining Algorithms In R/Clustering/Expectation Maximization (EM). Retrieved from https://en.wikibooks.org/wiki/Data_Mining_Algorithms_In_R/Clustering/Expectation_Maximization_(EM)
Tugas Praktikum 5 STA561 Pemrograman Statistika, Reni Amelia, Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, reniamelia@apps.ipb.ac.id↩︎