Posted on June 27, 2018 by matsukat

認知情報解析演習A W12

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

Related