ES
n.parent = 10
n.theta = 10
parent = list(theta = matrix(rnorm(n.parent*n.theta),nrow=n.parent),
sigma = matrix(runif(n.parent*n.theta)+1,nrow = n.parent))
ES_crossover<-function(parent){
n.parent = nrow(parent$theta)
rand.order = sample(1:n.parent)
p1.idx= rand.order[1:(n.parent/2)]
p2.idx= rand.order[(n.parent/2+1):n.parent]
co.idx = sample(1:(length(parent)/2),length(parent)/4)
c1 = parent$theta[p1.idx,]
c1.copy = c1
c2 = parent$theta[p2.idx,]
c1[co.idx] = c2[co.idx]
c2[co.idx] = c1.copy[co.idx]
sigma = (parent$sigma[p1.idx,]+parent$sigma[p2.idx,])/2
return(child = list(theta = rbind(c1,c2), sigma = rbind(sigma,sigma)))
}
ES_mutate<-function(child, tau){
child$sigma = child$sigma*
exp(tau*matrix(rnorm(length(child$sigma)),nrow=nrow(child$sigma)))
child$theta = child$theta + matrix(rnorm(length(child$sigma),0,child$sigma),nrow=nrow(child$sigma))
return(child)
}
ES_survive<-function(parent, child, fitP, fitC){
nPop=nrow(parent$theta);
fitT=c(fitP,fitC);
fitMax=sort(fitT,index.return=TRUE,decreasing = FALSE)
tempTheta=rbind(parent$theta, child$theta)
tempSigma=rbind(parent$sigma, child$sigma)
fittest=list(theta = tempTheta[fitMax$ix[1:nPop],], sigma = tempSigma[fitMax$ix[1:nPop],])
return(fittest)
}
set.seed(20); nData = 100
X=matrix(rnorm(9*nData,mean=10,sd=2),ncol=9);X=cbind(rep(1,nData),X)
y=X%*%c(10,2,5,-3,-5,0,0,0,0,0)+rnorm(nData,mean=0,sd=2);
#summary(lm(y~X[,2:10]))
reg.error<-function(b,X,y){
yhat<-X%*%b
return(sum((y-yhat)^2))
}
apply(child$theta,1,reg.error, X, y)
ES<-function(parent, nGen,tau,X,Y){
optHist=matrix(0,nrow=nGen,ncol=1)
fittest.hist = matrix(0,nrow=nGen,ncol = ncol(parent$theta))
nPop = nrow(parent$theta)
nVar = ncol(parent$theta)
fitP = apply(parent$theta, 1, reg.error, X, Y)
#fitC = fitP
optHist[1]=max(c(fitP))
fittest.hist[1,] = parent$theta[which.max(fitP),]
for (i_gen in 2:nGen) {
child<-ES_crossover(parent)
child<-ES_mutate(child,tau)
fitP = apply(parent$theta,1, reg.error, X, Y)
fitC = apply(child$theta, 1, reg.error, X, Y)
parent<-ES_survive(parent, child, fitP, fitC)
optHist[i_gen]=max(c(fitP,fitC))
fittest.hist[i_gen, ] = parent$theta[1,]
}
return(list(fittest = parent, optHist = optHist, fittest.hist = fittest.hist))
}
> res<-ES(parent,50000,1,X,y)
> res$fittest$theta[1,]
[1] 9.508492513 2.074884770 4.929865865 -2.989781723 -5.129436719 0.050919574 0.024012318
[8] -0.003250697 0.081489175 -0.028049091
> coef(lm(y~X[,2:9]))
(Intercept) X[, 2:9]1 X[, 2:9]2 X[, 2:9]3 X[, 2:9]4 X[, 2:9]5 X[, 2:9]6 X[, 2:9]7
9.424585305 2.072861468 4.929182784 -2.993701821 -5.135528730 0.048416459 0.023772334 -0.005105199
X[, 2:9]8
0.078469949
ES & PSO
funZ<-function(x,y) {x^4-16*x^2-5*x+y^4-16*y^2-5*y}
xmin=-5;xmax=5;n_len=100;
x<-seq(xmin, xmax,length.out=n_len);
y<-x;
z<-outer(x,y,funZ)
contour(x,y,z,nlevels= 50, drawlabels = F,col='blue')
# useful memo
funZ<-function(x,y) {x^4-16*x^2-5*x+y^4-16*y^2-5*y}
funZ<-function(x) {
x[1]^4-16*x[1]^2-5*x[1]+x[2]^4-16*x[2]^2-5*x[2]
}
n.theta = 10;
n.iter = 100;
Wp = 1;
Wg = 1;
theta = matrix(rnorm(n.theta*2), nrow=n.theta)
v = matrix(rnorm(n.theta*2), nrow=n.theta)
theta.hist = array(0,c(n.theta,2,n.iter))
theta.hist[,,1]=theta
p.best.v <- apply(theta,1,funZ)
p.best = theta
g.b.v = min(p.best.v)
g.b.idx = which.min(p.best.v)
g.best <- theta[g.b.idx,]
v = v + Wp*runif(n.theta)*(p.best - theta)+ Wg*runif(n.theta)*t(g.best - t(theta))
theta = theta + v