##package "leaps" library(leaps) head(swiss) #leaps : to select regression subsets leaps.model1=leaps(x=swiss[,-1] ,y=swiss[,1], nbest=3) leaps.model1 #nbest :Number of subsets of each size to report par(mfrow=c(1,1)) plot(leaps.model1$size,leaps.model1$Cp,ylim=c(0,40),xlab="model size",ylab="Cp") abline(0,1) #regsubsets sub.model1=regsubsets(Fertility ~ .,data=swiss) par(mfrow=c(1,2)) plot(sub.model1) plot(sub.model1,scale="Cp") #nbest: number of subsets of each size to record (default =1) #nvmax: maximum size of subsets to examine #force.in: index to columns of design matrix that should be in all models #force.out: index to columns of design matrix that should be in no models sub.model2=regsubsets(Fertility ~ .,data=swiss, nbest=2 , force.in=2) plot(sub.model2) plot(sub.model2,scale="Cp") summary(sub.model2) sub.model2=regsubsets(Fertility ~ .,data=swiss, force.out=c(2,3,7,8)) plot(sub.model2) ######################################################################################### #GenSa library(GenSA) ?GenSA # Try Rastrgin function (The objective function value for global minimum is 0 with all components of par are 0.) Rastrigin <- function(x) { sum(x^2 - 10 * cos(2 * pi * x)) + 10 * length(x) } #plot for 2-dimentional case z=matrix(0,ncol=100,nrow=100) x1=x2=seq(-5.12,5.12,length=100) for (i in 1:100){ for (j in 1:100){ x=c(x1[i],x2[j]) z[i,j]=Rastrigin(x) } } nx <- 100;ny <- 100;z=ras hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1,-1]) hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt)) persp(x1, x2, z, theta = 60, phi =30, expand = 0.5, col = cm.colors(10)[floor(9 * hgt + 1)]) persp(x1, x2, z, theta = 30, phi =10, expand = 0.5, col = cm.colors(10)[floor(9 * hgt + 1)]) #plot for 1-dimentional case a=seq(-5.12,5.12,by=0.05) res=vector() for (i in 1:length(a)){ res[i]= Rastrigin(a[i]) } res plot(a,res,ylab='',xlab='x',type='l') #use optim() with method SANN - perfoms poorly optim(par=c(5,5), fn=Rastrigin, gr = NULL, method = c("SANN"), lower = -Inf, upper = Inf, control = list(), hessian = T) optim(par=c(0.001,0.001), fn=Rastrigin, gr = NULL, method = c("SANN"), lower = -Inf, upper = Inf, control = list(), hessian = T) optim(par=c(0,0), fn=Rastrigin, gr = NULL, method = c("SANN"), lower = -Inf, upper = Inf, control = list(), hessian = T) # GenSA will stop after finding the targeted function value 0 with absolute tolerance 1e-13 set.seed(2014) dimension <- 2 global.min <- 0 tol <- 1e-13 lower <- rep(-5.12, dimension) upper <- rep(5.12, dimension) out <- GenSA(lower = lower, upper = upper, fn = Rastrigin, control=list(threshold.stop=global.min+tol,verbose=TRUE)) out[c("value","par","counts")]