欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

模拟退火算法的R语言实现

程序员文章站 2024-02-27 23:47:27
...

模拟退火算法:

  • 原理:
    固体在慢慢冷却下来的时候某内部的分子热能一般随机减少。但是在也有可能以boltzmann的概率随机增加。即在温度τ,能量增加幅度为ΔE 的概率密度为exp{ΔE/kτ} 其中kτ 是常数。但长时间后一定会以不变的能量存在。

  • 伪代码:
    1.initialize the x value
    2.随机寻找x_new
    3.calculate Δ=f(xnew)f(x)
    4.如果Δ<0 (求最小值),则令x=x_new
    5.如果不满足4,则以概率p=exp{ΔE/kτ} (其中k是变的)更新 x=x_new
    6.更新温度函数
    7.按事先指定的次数重复2-6
    说明:在这里,在不同的温度平台可以更新不同的次数。比如一共迭代15次,在220度更新60次,在120度更新120次,在60度更新240次。

  • R语言代码求解棒球队员薪资问题。这也是一个01变量求最优的故事。
    1.先进行数据处理,包括导入数据,把薪资这一列变成log形式做多元回归。把数据中的分类变量变成因子形式。

baseball.da <- read.csv("baseball.dat",head=T)
salary.log <- log(baseball.dat[,1])
baseball.dat$freeagent = factor(baseball.dat$freeagent)
baseball.dat$arbitration = factor(baseball.dat$arbitration)
baseball.sub <-as.data.frame(cbind(salary.log,baseball.da[,-1]))
m=length(baseball.sub[1,-1])
n=length(salary.log)

2.下面进行退火算法初始值的设置
设置更新初值,求目标函数的初始值。设置温度初值和温度变化。设置冷却平台的迭代次数。设置k

set.seed(19495)
run <- rbinom(m,1,0.5)
run.current <- run
run.vars <- baseball.sub[,c(1,run.current)]
g <- lm(salary.log~., run.vars)
run.aic <- extractAIC(g)[2]
best.aic <- run.aic
aic <- NULL
temp <- matrix(1000,1,m)
for(j in 1:15) {temp[j] <- temp[j-1]*0.95}
cooling <- c(rep(60,5),rep(120,5),rep(240,5))

3.开始写主要的退火过程。主要体现的是x变量的更新过程

for( j in 1:15){
for( i in 1:cooling[j]){
  pos <- sample(1,1:m)
  run.step <- run.current
  run.step[pos] <- !run.current[pos]
  run.vars <- baseball.sub[,c(1,run.step)]
  g <- lm(salary.log~.,run.vars)
  aic.step <- extractAIC(g)[2]
  p <-min(1,exp((run.aic-extractAIC(g)[2])/tau[j]))
  if(aic.step>run.aic)
  {
  run.aic <-aic.step
  run.current <- run.step
  }
 if(rbinom(1,1,p))
 {
  run.aic <-aic.step
  run.current <- run.step
 }
 if(aic.step>best.aic)
 {
 best.aic <-aic.step
 run <- run.step
 }
}
}

4.输出

## OUTPUT
run         # BEST LIST OF PREDICTORS FOUND
best.aic    # AIC VALUE
aics        # VECTOR OF AIC VALUES

## PLOT OF AIC VALUES
plot(aics,ylim=c(-420,-360),type="n",ylab="AIC", xlab="Iteration")
lines(aics)