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)
