Posted on June 20, 2018 by matsukat

認知情報解析演習　W10

n.var = 10
n.child = 10
child = matrix(sample(0:1,n.var*n.child, replace =T),nrow =n.child)

GA_mutate <- function(child, p){
  mut.mat = matrix(sample(0:1,length(child),replace = T, 
                          prob = c(1-p, p)),nrow=nrow(child))
  return(abs(child-mut.mat))
}

GA_crossover<-function(parent){
  n.parent = nrow(parent)
  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[p1.idx,]
  c1.copy = c1
  c2 = parent[p2.idx,]
  c1[co.idx] = c2[co.idx]
  c2[co.idx] = c1.copy[co.idx]
  return(child = rbind(c1,c2))
}


GA_survive<-function(parent, child, fitP, fitC){
  nPop=nrow(parent);
  fitT=c(fitP,fitC);
  fitMax=sort(fitT,index.return=TRUE,decreasing = TRUE)
  tempX=rbind(parent,child);
  fittest=tempX[fitMax$ix[1:nPop],]
  return(fittest)
}

GA<-function(parent, nGen,p,X,Y){
  optHist=matrix(0,nrow=nGen,ncol=1)
  fittest.hist = matrix(0,nrow=nGen,ncol = ncol(parent))
  nPop = nrow(parent)
  nVar = ncol(parent)
  fitP = rep(0,nPop)
  fitC = fitP
  for (i.pop in 1:nPop){
    dat.temp = X[,which(parent[i.pop,]==1)]
    fitP[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
  }
  optHist[1]=max(c(fitP))
  fittest.hist[1,] = parent[which.max(fitP),]
  for (i_gen in 2:nGen) {
    child<-GA_crossover(parent)
    child<-GA_mutate(child,p)
    # assuming same numbers of set of genes
    for (i.pop in 1:nPop){
      dat.temp = X[,which(parent[i.pop,]==1)]
      fitP[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
      if (sum(child[i.pop,])==0){
        child[i.pop,sample(1:ncol(child),sample(1:nVar,1))]=1
      }
      dat.temp = X[,which(child[i.pop,]==1)]
      fitC[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
    }
    parent<-GA_survive(parent, child, fitP, fitC)
    optHist[i_gen]=max(c(fitP,fitC))
    fittest.hist[i_gen, ] = parent[1,]
  }
  return(list(fittest = parent, optHist = optHist, fittest.hist = fittest.hist))
}

set.seed(19)
nVar = 10
nData = 50
X=matrix(rnorm(nVar*nData,mean=10,sd=5),ncol=nVar);
y=X%*%c(-2,0,2,0,2,0,-3,0,-2,0)+rnorm(nData,mean=0,sd=8);
summary(lm(y~X[,seq(1,9,2)]))
summary(lm(y~X))
parent = matrix(sample(0:1,300,replace = T),nrow=30)
res<-GA(parent,100,0.1,X,y)
head(res$fittest)

より複雑なモデルへ

C = combn(3,2)
paste("X",C,sep="")
temp.lab =toString(paste("X",C[ , 1], sep = ""))
gsub(",", ":", temp.lab) 

n.comb = ncol(C)
var.labels=c()
for (i.comb in 1: n.comb){
  temp.label = toString(paste("X",C[ , i.comb], sep = ""))
  var.labels = c(var.labels, gsub(",", ":", temp.label) )
}


mk.labels <- function(C){
  var.labels = c()
  n.comb = ncol(C)
  for (i.comb in 1: n.comb){
    temp.label = toString(paste("X",C[ , i.comb], sep = ""))
    var.labels = c(var.labels, gsub(",", ":", temp.label) )
  }
  return(var.labels)
}
                                    
n.var = 3
max.interaction = 3
# single variable
var.labels = paste("X", 1:n.var, sep = "")


# interaction terms
for (i.interaction in 2:max.interaction ){
  combination = combn(n.var, i.interaction)
  var.labels = c(var.labels, mk.labels(combination))
}

model.def = paste("Y ~", gsub(",", "+",   toString(var.labels[c(1,1,1,1,1,0,0) == 1]))) 

set.seed(20)
nVar = 3
nData = 50
X=matrix(rnorm(nVar*nData,mean=10,sd=5),ncol=nVar);
y=X%*%c(1,1,-2)+rnorm(nData,mean=0,sd=3);
dat = data.frame(Y=y,X1=X[,1],X2=X[,2],X3=X[,3])

# example
model.def1 = paste("Y ~", gsub(",", "+",  toString(var.labels[c(1,1,1,0,0,0,0) == 1]))) 
model.lm = lm(model.def1, dat)

Evolutionary Strategy

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)
}

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))
}
# 評価の使用例
fitC = apply(child$theta,1,reg.error, X, y)

Related