禁忌演算法
禁忌演算法是啟發式演算法對個體的應用的一種。由於在運用最速下降或者最速上升區域性搜尋最值的時候可能會因為到了區域性最小值後停止搜尋。這裡禁忌演算法是一種可以look back的演算法，但是需要設定一些禁忌目錄來保證搜尋不是無限的。
例如：在Baseball salary案例中，我們希望尋找幾個與salary最connected的features去做multiple regression所以，需要在p個features中尋找s個features。我們在這裡通過線性迴歸的衡量引數最大擬合模型的AIC來做一個需要優化的function。。我們需要尋找的是哪s個feature使得AIC最大。所以這是一個NP問題。
下面進行程式碼的說明:

cun.current <- rbinom(m,1,0.5)#在m個feature裡邊隨便挑選一半的feature進行初始的train
cun.var=baseball.sub[,cun.current==1]
g <-lm(salary.log~.,cun.var)#用選中的feature進行aic值的訓練。
cun.aic.current <- extractAIC(g)[2]


for (i in 1:m)
{
cun.next = cun.current#將我們選中的feature的序號賦值給cun.next
cun.next[i]  <- !cun.current[i]#挨個換我們之前選得feature比如之前在123456789中選擇了13579現在由於將第一個取反變成了3579.第二次迴圈將變成選擇123579
 

cun.var=baseball.sub[,cun.next==1]
g <- lm(salary.log~.,cun.var)
cun.aic <- extractAIC(g)[2]
if(cun.aic<cun.aic.current)
{
  cun.current <- cun.next#我們把取得最大aic值的順序記錄下來放在下一次訓練中。其他的小於的序列捨去
  cun.aic.current <- cun.aic
}
}

將本段程式迴圈15次（我們最初設定的最大迭代次數。）應該是找到aic不再動的時候停止。
這樣我們找到了最優的aic值。這時，選中的feature序列數就是我們需要找到最優的多元線性迴歸需要用到的最優features.
主要程式就是這樣子了~

baseball.dat = read.table("D:\\sunyao\\baseball.dat",header=TRUE)
baseball.dat$freeagent = factor(baseball.dat$freeagent)
baseball.dat$arbitration = factor(baseball.dat$arbitration)
baseball.sub = baseball.dat[,-1]
salary.log = log(baseball.dat$salary)
n = length(salary.log)
m = length(baseball.sub[1,])#去掉第一列的salary
num.starts = 5
runs = matrix(0,num.starts,m)
itr = 15
runs.aic = matrix(0,num.starts,itr)

# INITIALIZES STARTING RUNS
set.seed(19676) 
for(i in 1:num.starts){runs[i,] = rbinom(m,1,.5)}

## MAIN
for(k in 1:num.starts){
    run.current = runs[k,]#隨機選取要被訓練的樣本。

    # ITERATES EACH RANDOM START
    for(j in 1:itr){
        run.vars = baseball.sub[,run.current==1]#選取一行中的的一半的feature隨機 進行train
        g = lm(salary.log~.,run.vars)
        run.aic = extractAIC(g)[2]#選取本很函式中的第二個value一共兩個：第一個自由度。第二個是AIC值
        run.next = run.current

        # TESTS ALL MODELS IN THE 1-NEIGHBORHOOD AND SELECTS THE
        # MODEL WITH THE LOWEST AIC
        for(i in 1:m){
            run.step = run.current
            run.step[i] = !run.current[i]#給第一個feature取反
            run.vars = baseball.sub[,run.step==1]
            g = lm(salary.log~.,run.vars)
            run.step.aic = extractAIC(g)[2]
            if(run.step.aic < run.aic){
                run.next = run.step
                run.aic = run.step.aic
            }
        }
        run.current = run.next
        runs.aic[k,j]=run.aic
    }
    runs[k,] = run.current
}

## OUTPUT
runs        # LISTS OF PREDICTORS
runs.aic    # AIC VALUES

我們可以用影象表示我們迭代過程中aic值的變化情況。

plot(1:(itr*num.starts),-c(t(runs.aic)),xlab="Cumulative Iterations",
  ylab="Negative AIC",ylim=c(360,420),type="n")
for(i in 1:num.starts) {
  lines((i-1)*itr+(1:itr),-runs.aic[i,]) }

我發現R包一個tabuSearch可以成功的找到我們需要的函式的最小值。
下面對tabuSearch函式進行簡介。
tabuSearch用來搜尋一個二進位制陣列的最佳形式。其中我們需要自己寫一個evaluate函式來確定演算法的評判標準。
下面是引數說明：
tabusearch()
一種優化二進位制字串的禁忌搜尋演算法。 它需要使用者定義的目標函式
並報告在整個搜尋過程中找到的最佳二進位制配置，即最高的二進位制配置
目標函式值。 該演算法可用於變數選擇。
size:二進位制配置的長度。
iters:迭代次數
objFunc：引數為01的目標函式用來評判選擇變數
config：初始值
neigh;設定每次迭代時候鄰居的位置：初始值是all size
listSize：禁忌列表的長度
nRestarts：演算法強化的時候最多可以重新啟動的次數
repeatAll：可以重複搜尋的次數
verbose：T\F T：列印所在階段的名稱：初步階段。集約化階段。多元化階段
對於上一道題。我不明白為啥算出來的aic值成了最小的了····

evaluate <- function(th){
  if(sum(th)==0)return(0)
  run.vars <- as.data.frame(cbind(salary.log~.,baseball.sub[,th==1]))
  g <- lm(salary.log~.,run.vars)
  cun.aic <- extractAIC(g)[2]
  return(cun.aic)
}
res <- tabuSearch(size = m, iters = 50, objFunc = evaluate, config = matrix(1,1,m),
                  listSize = 5, nRestarts = 4)