#### Demmler & Reinsch baza #### #### smoothing spline primjer iz ESL 5.5.2 ### #### preuzeto od Feng Liang, UIUC #### library(splines) ### podaci i stvarna regresijska funkcija set.seed(1234) n = 50 err = 1.5 x = 1:50/50 y = sin(12*(x+0.2))/(x+0.2) + rnorm(n, 0, err); plot(x, y, col="red"); fx = 1:100/100; fy = sin(12*(fx+0.2))/(fx+0.2) lines(fx, fy, col=8, lwd=2); #### ss s razlicitim df par(mfrow=c(2,2)) plot(x,y, xlab='x', ylab='y'); lines(fx, fy, col=8, lwd=1.5); lines(predict(smooth.spline(x, y, df=5),fx), lty=2, col='blue', lwd=1.5); title('df=5'); plot(x,y,xlab='x', ylab='y'); lines(fx, fy, col=8, lwd=1.5); lines(predict(smooth.spline(x, y, df=9),fx), lty=2, col='blue', lwd=1.5); title('df=9'); plot(x,y,xlab='x', ylab='y'); lines(fx, fy, col=8, lwd=1.5); lines(predict(smooth.spline(x, y, df=15),fx), lty=2, col='blue', lwd=1.5); title('df=15'); plot(x,y,xlab='x', ylab='y'); lines(fx, fy, col=8, lwd=1.5); lines(predict(smooth.spline(x, y, df=20),fx), lty=2, col='blue', lwd=1.5); title('df=20') #### Demmler & Reinsch baza ## kod ss-a imamo y_hat=S_lambda y ## sljedeca funkcija vraca S_lambda ako su cvorovi x ## te stupnjevi slobode = df (u 1-1 vezi s lambda) smooth.matrix = function(x, df){ # pretpstavka je da su svi elementi iz x razliciti n = length(x); A = matrix(0, n, n); for(i in 1:n){ ### ako je y=(0,...,0,1,0,...,0) pri cemu je 1 na i-tom mjestu ### S_lambda y ce tocno dati i-ti redak matrice S_lambda y = rep(0, n); y[i]=1; yi = smooth.spline(x, y, df=df)$y; A[,i]= yi; } return(A) } S4 = smooth.matrix(x, df=4); S9 = smooth.matrix(x, df=9); eigen.S4 = eigen(S4); eigen.S9 = eigen(S9); u4 = eigen.S4$ve; ## svojstveni vektori -- vektori DR baze u9= eigen.S9$ve; ## u9=u4 (do na predznak) rho_4 = eigen.S4$va ## svojstvene vrijednosti rho_9 = eigen.S9$va ### vektore DR baze u_j shvacamo kao funkcije od x ### te gledamo S_y u_j = rho_j(lambda)*u_j = ### ---> smanjenje je vece za u_j za koje je rho_j(lambda) manje ## df=4 par(mfrow=c(2,3)); for(i in 1:12) { ## crtamo samo prvih 12 elemenata DR baze plot(c(min(x), max(x)),c(min(u4, u9), max(u4, u9)), xlab="x", ylab="", type="n"); lines(x, u4[,i], col=2,lty=1, lwd=2.5); lines(x, u4[,i]*rho_4[i], lty=2, lwd=2) } ## veci j odgovaraju zakrivljenijim funkcijama u_j (dakle, za ## za njih je smanjenje vece) rho_4[1:2] ## [1] 1 1 , tj. d_1=d_2=0 ## dakle, nema smanjenja u_1 i u_2 koji razapinju prostor linearnih funkcija ## (jer lin. funkcije imaju drugu derivaciju == 0) ##isto za df=9 par(mfrow=c(2,3)); for(i in 1:12) { plot(c(min(x), max(x)),c(min(u4, u9), max(u4, u9)), xlab="x", ylab="", type="n"); lines(x, u9[,i], col=2,lty=1, lwd=2.5); lines(x, u9[,i]*rho_9[i], lty=2, lwd=2) } ## ---> manje smanjenje zakrivljenijih funkcija #### CV metoda za odabir df df = seq(2, n-1, by=0.5) ### iz [2,n] m = length(df) mycv = rep(0,m) mygcv= rep(0,m) for (i in 1:m){ fit = smooth.spline(x, y, df=df[i], cv=T); mycv[i] = fit$cv fit = smooth.spline(x, y, df=df[i], cv=F); mygcv[i]=fit$cv } par(mfrow=c(1,1)) plot(df, mycv, col="red", type="l", xlab='df', ylab='CV scores'); points(df, mygcv, col="darkblue", type="l") ### slicno # optimalni df (koristeci LOOCV) optdf = df[mycv==min(mycv)] optdf ## [1] 8.5 # df[mygcv==min(mygcv)] ## [1] 8.5 # ss za optimalni df fitopt = smooth.spline(x, y, df=optdf); plot(x, y, xlab="x", ylab="y") lines(predict(fitopt, (1:100)/100),col="red", lwd=2) lines(fx, fy, col="gray", lwd=2)