Tugas SML Regresi Nonparametrik

Tugas yang diberikan merujuk pada SPMbook halaman 126. Dalam regresi nonparametrik bentuk kurva regresi tidak diketahui, maka data diharapkan mencari sendiri bentuk estimasinya sehingga memiliki fleksibelitas yang tinggi. Kurva regresi hanya diasumsikan termuat dalam suatu ruang fungsi yang berdimensi tak hingga dan merupakan fungsi (smooth). Model umum regresi nonparametric adalah \(y_i=m(x_i)+ \varepsilon_i\) dengan \(i=1,2,3,\dotsc,n\). Estimasi fungsi \(m(x_i)\) dilakukan berdasarkan data pengamatan dengan menggunakan teknik smoothing tertentu. Fungsi regresi diestimasi menggunakan rumus \(\hat{m}(x_i)=\frac{1}{n}\sum_{i=1}^{n}W_{hi}(x)Y_i\). Estimator Nadaraya Watson : \(\frac{K_h{(x-X_i)}}{\hat{f}_h(x)}\)

Simulasi Data

Untuk mensiumulasikan kernel regresi nonparametrik, dilakukan pembangkitan data dengan kurva \(m(x)\)= \(\sin(2\pi x^3 )^3\), \(Y_i=m(X_i)+\varepsilon_i\), dan \(X_i\) berdistribusi \(U[0,1]\) serta \(\varepsilon_i\) berdistribusi \(N(0,0.1)\). Berikut merupakan bentuk plot random data

set.seed(12345)
n <- 500
eps <- rnorm(n, mean=0,sd = 0.1)
m <- function(x) (sin(2*pi*x^3))^3

X <- runif(n,min=0,max=1)
Y <- m(X) + eps
plot(X,Y,pch=20)

Dari plot menunjukkan bahwa pola nonlinier.Untuk mengestimasi data berbentuk nonlinier, digunakan kernel regression dengan teknik regresi nonparametrik. Dalam melakukan regresi kernel, terdapat dua konfigurasi yaitu fungsi kernel dan ukuran bandwidth. Bandwidth (h) adalah parameter (smoothing) yang berfungsi untuk mengontrol kemulusan dari kurva yang diestimasi. Bandwidth yang terlalu kecil akan menghasilkan kurva yang under-smoothing yaitu sangat kasar dan sangat fluktuatif, dan sebaliknya bandwidth yang terlalu lebar akan menghasilkan kurva yang over-smoothing yaitu sangat smooth, tetapi tidak sesuai dengan pola data (Hardle, 1994). Oleh karena itu perlu dipilih bandwidth yang optimal menggunakan Cross Validation.

Apabila dilakukan perbandingan beberapa bandwitdh dengan kernel yang sama, maka akan diperoleh plot sebagai berikut :

#perbandingan bandwidth
Kreg1 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth = 0.01)
Kreg2 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth = 0.03)
Kreg3 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth = 0.07)
plot(X,Y,pch=20,main="Nadaraya-Watson Estimator Perbandingan Bandwidth")
lines(Kreg1, lwd=3, col="orange")
lines(Kreg2, lwd=3, col="purple")
lines(Kreg3, lwd=3, col="limegreen")
legend("bottomleft", c("h=0.01","h=0.03","h=0.07"), lwd=3, col=c("orange","purple","limegreen"),cex=0.8)

Dari plot tersebut terdapat bias-variance tiap plot, ketika h kecil, variabilitasnya besar namun biasnya kecil, dan sebaliknya ketika h besar pada warna hijau, variabilitasnya kecil tapi biasnya besar. Dari perbandingan hasil bandwitdh secara visual, h=0.03 lebih optimum daripada 0.07.

Cross-Validation Estimator Nadaraya-Watson

Cross Validation digunakan untuk memperkirakan kesalahan prediksi dari estimasi regresi yang telah dibuat.Selain itu CV dapat digunakan untuk memilih nilai bandwidth(\(h\)) yang optimum.

#CV nadaraya-watson
n = length(X)
# n: sample size
h_seq = seq(from=0.01,to=1, by=0.01)
# smoothing bandwidths we are using
CV_err_h = rep(NA,length(h_seq))
for(j in 1:length(h_seq)){
  h_using = h_seq[j]
  CV_err = rep(NA, n)
  for(i in 1:n){
    X_val = X[i]
    Y_val = Y[i]
    # validation set
    X_tr = X[-i]
    Y_tr = Y[-i]
    # training set
    Y_val_predict = ksmooth(x=X_tr,y=Y_tr,kernel = "normal",bandwidth=h_using, x.points = X_val)
    CV_err[i] = (Y_val - Y_val_predict$y)^2
    # we measure the error in terms of difference square
  }
  CV_err_h[j] = mean(CV_err)
}
CV_err_h
##   [1] 0.01149638 0.01051951 0.01043847 0.01085760 0.01163032 0.01268788
##   [7] 0.01401589 0.01560004 0.01741004 0.01939943 0.02151497 0.02370473
##  [13] 0.02592676 0.02815073 0.03035872 0.03254349 0.03470433 0.03684583
##  [19] 0.03897666 0.04110508 0.04324026 0.04538905 0.04755743 0.04975032
##  [25] 0.05197002 0.05421688 0.05649080 0.05879093 0.06111434 0.06345762
##  [31] 0.06581686 0.06818814 0.07056724 0.07294942 0.07533046 0.07770622
##  [37] 0.08007237 0.08242522 0.08476102 0.08707638 0.08936814 0.09163335
##  [43] 0.09386957 0.09607440 0.09824576 0.10038185 0.10248117 0.10454227
##  [49] 0.10656393 0.10854517 0.11048524 0.11238340 0.11423915 0.11605214
##  [55] 0.11782215 0.11954903 0.12123275 0.12287347 0.12447134 0.12602666
##  [61] 0.12753972 0.12901100 0.13044092 0.13183004 0.13317892 0.13448817
##  [67] 0.13575844 0.13699041 0.13818478 0.13934230 0.14046371 0.14154979
##  [73] 0.14260132 0.14361910 0.14460392 0.14555660 0.14647794 0.14736876
##  [79] 0.14822987 0.14906208 0.14986619 0.15064299 0.15139328 0.15211784
##  [85] 0.15281744 0.15349285 0.15414480 0.15477404 0.15538130 0.15596728
##  [91] 0.15653268 0.15707818 0.15760445 0.15811214 0.15860187 0.15907428
##  [97] 0.15952995 0.15996949 0.16039345 0.16080240
plot(x=h_seq, y=CV_err_h, type="b", lwd=3, col="blue",
     xlab="Smoothing bandwidth", ylab="LOOCV prediction error")

h_opt=h_seq[which(CV_err_h == min(CV_err_h))]
h_opt
## [1] 0.03

Dari hasil CV yang telah dihitung, didapatkan nilai h optimum sebesar 0.03.

Nadaraya-Watson EStimator

Berikutnya, akan dicobakan h=0.03 pada estimator Nadaraya-Watson secara manual dari fungsi regresi \(\hat{m}(x_i)=\frac{1}{n}\sum_{i=1}^{n}W_{hi}(x)Y_i\).

mNW <- function(x, X, Y, h, K = dnorm) {
  # x: evaluation points
  # X: vector (size n) with the predictors
  # Y: vector (size n) with the response variable
  # h: bandwidth
  # K: kernel
  
  # Matrix of size n x length(x)
  Kx <- sapply(X, function(Xi) K((x - Xi) / h) / n*h)
    # Weights
  W <- Kx / rowSums(Kx) # Column recycling!
    # Means at x ("drop" to drop the matrix attributes)
  drop(W %*% Y)
  }
ps <- rnorm(n, mean=0,sd = 0.1)
m <- function(x) (sin(2*pi*x^3))^3

X <- runif(n,min=0,max=1)
Y <- m(X) + eps
xGrid <- seq(-5, 5, l = 700)
h2 <- 0.03
# Plot data
plot(X, Y,pch=20,col="#42b883",ylab="Y,m,mh",main="Plot Nadaraya Watson dan Fit Regression")
lines(xGrid, m(xGrid), col = 1,lwd=3)
lines(xGrid, mNW(x = xGrid, X = X, Y = Y, h = h2), col = 2,lwd=3)
legend("bottom", legend = c("Fit regression", "Nadaraya-Watson Estimator regression h=0.03"),
       lwd = 3, col = 1:2,cex=0.6)

Local Polynomial Regression Smoothing Method

Estimator Nadaraya-Watson dapat digunakan untuk kasus dengan estimator nonparametrik yang lebih luas, yang disebut estimator polinomial lokal.Simulasi menggunakan library KernSmooth.

library(KernSmooth)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
#localpolynomial regression
plot(X, Y,pch=20,main="Plot Local Polynomial Regression")
fit <- locpoly(X, Y, bandwidth = 0.02)
lines(fit,col="green",lwd=3)

Dari hasil perhitungan CV didapatkan bahwa bandwidth optimum sebesar 0.02. Bandwidth yang sesuai dapat meminimumkan bias.

Spline Regression

Smoothing regression yang lain adalah spline smoothing dengan mempertimbangkan jumlah residual kuadrat yang didefinisikan sebagai \(\sum_{i=1}^n\{Y_i-m(X_i)\}^2\), dengan fungsi \(m(X_i)=Y_i\) \(i=1,2,\dots,n\). \(m(X_i)\) diestimasi dengan \(\hat{m}_{\lambda}(x)=n^{-1}\sum_{i=1}^nW_{{\lambda_i}}(x)Y_i\). Pada smoothing spline dilakukan optimasi pada spar atau smoothing parameter. Ketika spar ditentukan, koefisien \(\lambda\) integral dari turunan kedua kuadrat dalam kriteria fit adalah fungsi monoton dari spar (R).

SS1 = smooth.spline(x=X,y=Y,spar=0.2)
SS2 = smooth.spline(x=X,y=Y,spar=0.7)
SS3 = smooth.spline(x=X,y=Y,spar=1.2)
plot(X,Y,pch=20)
lines(SS1, lwd=3, col="orange")
lines(SS2, lwd=3, col="purple")
lines(SS3, lwd=3, col="limegreen")
legend("bottomleft", c("spar=0.2","spar=0.7","spar=1.2"), lwd=6,
col=c("orange","purple","limegreen"),cex=0.8)

Cross Validation Spline Regression

Sama dengan regresi kernel, cross-validasi digunakan untuk menentukan nilai dari spar (smooth parameter).

n = length(X)
# n: sample size
sp_seq = seq(from=0.01,to=1.0, by=0.01)
# values of spar we are exploring
CV_err_sp = rep(NA,length(sp_seq))
for(j in 1:length(sp_seq)){
spar_using = sp_seq[j]
CV_err = rep(NA, n)
for(i in 1:n){
X_val = X[i]
Y_val = Y[i]
# validation set
X_tr = X[-i]
Y_tr = Y[-i]
# training set
SS_fit = smooth.spline(x=X_tr,y=Y_tr,spar=spar_using)
Y_val_predict = predict(SS_fit,x=X_val)
# we use the 'predict()' function to predict a new value
CV_err[i] = (Y_val - Y_val_predict$y)^2
}
CV_err_sp[j] = mean(CV_err)
}
CV_err_sp
##   [1] 0.01221694 0.01219839 0.01217773 0.01215486 0.01212972 0.01210226
##   [7] 0.01207248 0.01204040 0.01200606 0.01196956 0.01193101 0.01189055
##  [13] 0.01184837 0.01180465 0.01175960 0.01171347 0.01166646 0.01161882
##  [19] 0.01157074 0.01152241 0.01147396 0.01142550 0.01137708 0.01132871
##  [25] 0.01128038 0.01123207 0.01118375 0.01113540 0.01108703 0.01103870
##  [31] 0.01099048 0.01094251 0.01089496 0.01084801 0.01080185 0.01075670
##  [37] 0.01071275 0.01067019 0.01062918 0.01058990 0.01055249 0.01051707
##  [43] 0.01048375 0.01045264 0.01042379 0.01039726 0.01037309 0.01035128
##  [49] 0.01033185 0.01031481 0.01030019 0.01028806 0.01027852 0.01027177
##  [55] 0.01026810 0.01026794 0.01027184 0.01028056 0.01029503 0.01031640
##  [61] 0.01034603 0.01038548 0.01043655 0.01050126 0.01058186 0.01068090
##  [67] 0.01080124 0.01094613 0.01111926 0.01132477 0.01156725 0.01185161
##  [73] 0.01218288 0.01256597 0.01300525 0.01350418 0.01406499 0.01468839
##  [79] 0.01537351 0.01611806 0.01691869 0.01777161 0.01867328 0.01962132
##  [85] 0.02061521 0.02165694 0.02275146 0.02390666 0.02513315 0.02644359
##  [91] 0.02785167 0.02937089 0.03101317 0.03278746 0.03469850 0.03674604
##  [97] 0.03892441 0.04122286 0.04362654 0.04611809
plot(x=sp_seq, y=CV_err_sp, type="b", lwd=3, col="blue",
xlab="Value of 'spar'", ylab="LOOCV prediction error")

sp_seq[which(CV_err_sp == min(CV_err_sp))]
## [1] 0.56

Didapatkan spar optimum sebesar 0.55 dan dibuat plot seperti berikut.

#spline regression
plot(X,Y,pch=20)
xm <- data.frame(cbind(X,Y))
xm <- xm[order(X),]
mh <- smooth.spline(xm$X, xm$Y, spar=0.55)
lines(xm$X,fitted(mh),col="red",lwd=3)

Selanjutnya menunjukkan perbandingan estimasi dari kernel regresi, local polynomial regression dan spline regression.

plot(X,Y,pch=20)
lines(xGrid, mNW(x = xGrid, X = X, Y = Y, h = h2), col = "purple",lwd=3) #nadaraya
lines(xm$X,fitted(mh),col="red",lwd=3) #spline
lines(fit,col="blue",lwd=3) #localpoli
legend("bottom", legend = c("Nadaraya-Watson","Spline","Locpoly"),col=c("purple","red","blue"),
        lwd = 2,cex=0.8)