# Thomas Stibor, 20.06.2008 library(e1071); ############################################## create.sin.cos.data <- function(m,noise=0.2) { ############################################## x1 <- c(1:2*m); x2 <- c(1:2*m); for (i in 1:m) { x1[i] <- (i/m) * pi; x2[i] <- sin(x1[i]) + rnorm(1,0,noise); } for (j in 1:m) { x1[m+j] <- (j/m + 1/2) * pi; x2[m+j] <- cos(x1[m+j]) + rnorm(1,0,noise); } typ <- c(rep(+1,m),rep(-1,m)) return(data.frame(x1,x2,typ)) } dataset <- create.sin.cos.data(200); # number of rows n <- nrow(dataset); # number of folds folds <- 10 samps <- sample(rep(1:folds, length=n), n, replace=FALSE) # test error vector test.error <- rep(0,folds); for (fold.iter in 1:folds) { train <- dataset[samps!=fold.iter,] # fit the model test <- dataset[samps==fold.iter,] # predict # features x1 and x2 from training set x_train <- train[,1:2]; # class labels y_train <- train[,3]; # determine the hyperplane model <- svm(x_train,y_train,type="C-classification", cost=10,kernel="radial", probability = TRUE, scale = FALSE); # features x1 and x2 from test set (no class labels) x_test <- test[,1:2]; # predict class labels for x_test pred <- predict(model,x_test); # prediction kernlab package ## pred.kernlab <- predict(model.kernlab, x_test); # get true class labels y_true <- test[,3]; # check accuracy: table(pred, y_true); # determine test error test.error[fold.iter] <- sum(pred != y_true)/length(y_true); } average.test.error <- sum(test.error)/length(test.error); print(average.test.error); # determine the hyperplane for complete set for generating a fancy image model <- svm(dataset[,1:2],dataset[,3],type="C-classification", cost=10,kernel="radial",scale = FALSE); # evaluate resol * resol many points, required for obtaining a detailed density # plot (see image command) x <- cbind(dataset$x1,dataset$x2); resol <- 200; rng <- apply(x,2,range); tx <- seq(rng[1,1],rng[2,1],length=resol); ty <- seq(rng[1,2],rng[2,2],length=resol); pnts <- matrix(nrow=length(tx)*length(ty),ncol=2); k <- 1 for(j in 1:length(ty)){ for(i in 1:length(tx)){ pnts[k,] <- c(tx[i],ty[j]); k <- k+1; } } # SVM predictition pred <- predict(model,pnts,decision.values = TRUE) z <- matrix(attr(pred,"decision.values"),nrow=length(tx),ncol=length(ty)) # uncomment if you want encapsulated postscript e.g. for LaTeX documents postscript("svmimage.eps", width = 10.0, height = 10.0, horizontal = FALSE, onefile = FALSE, paper = "special"); # png("svm_image.png", width = 600.0, height = 600.0, # horizontal = FALSE, onefile = FALSE, paper = "special"); image(tx,ty,z,xlab="",ylab="",axes=FALSE, xlim=c(rng[1,1],rng[2,1]),ylim=c(rng[1,2],rng[2,2]), # uncomment if you prefere a different color schema # col = heat.colors(100)); col = rainbow(200, start=0, end=.25)); # plot minus, zero and plus lines contour(tx,ty,z,add=TRUE, drawlabels=TRUE, level=0, lwd=3) contour(tx,ty,z,add=TRUE, drawlabels=TRUE, level=1, lty=1, lwd=1, col="grey") contour(tx,ty,z,add=TRUE, drawlabels=TRUE, level=-1, lty=1, lwd=1, col="grey") # uncomment if you want a title at the top # title(main = "SVM classification", font.main = 4) # plot points from -1 and +1 class points(dataset[dataset$typ==1,1:2],pch=21,col=1,cex=1) points(dataset[dataset$typ==-1,1:2],pch=19,col=4,cex=1) # plot highlighted SV points sv <- dataset[c(model$index),]; sv1 <- sv[sv$typ==1,]; sv2 <- sv[sv$typ==-1,]; points(sv1[,1:2],pch=13,col=1,cex=2) points(sv2[,1:2],pch=13,col=4,cex=2) dev.off();