認知情報解析 GA/ES preview

ES_recomb<-function(G,child) {
  nParent=nrow(G$b);nChild=nrow(child$b);nVar=ncol(G$b)
  for (i_child in 1:nChild) {
    parentID=sample(1:nParent,2)
    coID=sample(c(0,1),nVar,replace=T)
    child$b[i_child,]=G$b[parentID[1],]
    child$b[i_child,which(coID==1)]=G$b[parentID[2],which(coID==1)]
    child$sigma[i_child,]=0.5*(G$sigma[parentID[1],]+G$sigma[parentID[2],])
  }
  return(child)
}
ES_mutate<-function(child,tau){
  nChild=nrow(child$b);nVar=ncol(child$b)
  child$sigma<-child$sigma*exp(matrix(rnorm(nChild*nVar)*tau,nrow=nChild))
  child$b=child$b+child$sigma*matrix(rnorm(nChild*nVar),nrow=nChild,ncol=nVar)
  return(child)
}
ES_er<-function(b,x,y){
  yhat<-x%*%b
  return(sum((y-yhat)^2))
}

x=matrix(rnorm(4*50,mean=10,sd=2),ncol=4);x=cbind(rep(1,50),x)
y=x%*%c(10,2,5,-3,-5)+rnorm(50,mean=0,sd=2);
    
G=list();child=list();
nGen=1000;nParent=30;nChild=60;tau=1;nVar=5
G$b=matrix(rnorm(nVar*nParent,mean=0,sd=1),ncol=nVar)
G$sigma=matrix(runif(nVar*nParent),ncol=nVar)+0.5
child$b=matrix(0,nrow=nChild,ncol=nVar)
child$sigma=matrix(0,nrow=nChild,ncol=nVar)
optHist=matrix(0,nrow=nGen,ncol=1)
for (i_gen in 1:nGen) {
  child<-ES_recomb(G,child)
  child<-ES_mutate(child,tau)
  fitG=apply(G$b,1,ES_er,x,y);fitC=apply(child$b,1,ES_er,x,y)
  fitT=c(fitG,fitC);fitMin=sort(fitT,index.return=T)
  tempB=rbind(G$b,child$b);tempS=rbind(G$sigma,child$sigma)
  G$b=tempB[fitMin$ix[1:nParent],]
  G$sigma=tempS[fitMin$ix[1:nParent],]
  optHist[i_gen]=fitMin$x[1]
}

>  head(G$b)
         [,1]     [,2]     [,3]      [,4]      [,5]
[1,] 8.597338 2.118279 5.081328 -2.982018 -5.062286
[2,] 8.597347 2.118279 5.081328 -2.982018 -5.062286
[3,] 8.597336 2.118280 5.081328 -2.982018 -5.062286
[4,] 8.597333 2.118280 5.081328 -2.982018 -5.062286
[5,] 8.597343 2.118280 5.081328 -2.982018 -5.062286
[6,] 8.597336 2.118280 5.081328 -2.982018 -5.062287
> summary(lm(y~x[,2:5]))

Call:
lm(formula = y ~ x[, 2:5])

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1866 -1.7043  0.3254  1.3890  4.1851 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.5973     3.9641   2.169   0.0354 *  
x[, 2:5]1     2.1183     0.2235   9.476 2.73e-12 ***
x[, 2:5]2     5.0813     0.2077  24.464  < 2e-16 ***
x[, 2:5]3    -2.9820     0.1606 -18.566  < 2e-16 ***
x[, 2:5]4    -5.0623     0.2120 -23.874  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.342 on 45 degrees of freedom
Multiple R-squared:  0.9737,	Adjusted R-squared:  0.9713 
F-statistic: 415.7 on 4 and 45 DF,  p-value: < 2.2e-16

Traveling Sales Man Problem

# route inversion
TSP_inv<-function(route,len) {
  len_route=length(route)
  invP=sample(1:len_route,1)
  route[invP:min(len_route,(invP+len-1))]
    =rev(route[invP:min(len_route,(invP+len-1))])
  return(route)
}
> TSP_inv(1:10,4)
 [1]  1  2  3  4  5  9  8  7  6 10

# route switching
TSP_switch<-function(route){
  len_route=length(route)
  switchP=sample(1:len_route,2)
  route[switchP]=route[rev(switchP)]
  return(route)
}
> TSP_switch(1:10)
 [1]  1  2  3  4  5  6  7 10  9  8

# route translation
TSP_trans<-function(route){
  len_route=length(route)
  transP=sample(1:len_route,2);
  tempR=route[-transP[1]]
  if (transP[2]==1){
    tempR=c(route[transP[1]],tempR)
  } else if (transP[2]==len_route) {
    tempR=c(tempR,route[transP[1]])
  } else {
    tempR=c(tempR[1:(transP[2]-1)],route[transP[1]],tempR[(transP[2]):  
     (len_route-1)])
  }
  return(tempR)
}
> TSP_trans(1:10)
 [1]  1  3  4  5  6  7  8  9  2 10

# initialize cities' locations
FUN_initTSP<-function(n_city=10) {
  return(matrix(runif(n_city*2,1,100),nrow=n_city,ncol=2))
}
> print(loc<-FUN_initTSP(10))
          [,1]     [,2]
 [1,] 36.78996 40.33464
 [2,] 91.67296 97.33156
 [3,] 87.15730 58.82665
 [4,] 56.19289 44.84425
 [5,] 46.06971 43.22932
 [6,] 50.91508 30.14634
 [7,] 84.18521 56.17124
 [8,] 77.68784 43.77393
 [9,] 42.48220 11.74252
[10,] 12.99097 21.11037

# calc. total distance
FUN_costTSP<-function(route,cities) {
  route=c(route,route[1]);n_cities=nrow(cities);totDist=0;
  for (i_cities in 1:n_cities) {
    totDist=totDist+dist(cities[route[i_cities:(i_cities+1)],])
  }
  return(totDist)
}
> FUN_costTSP(1:10,loc)
         1
2 341.4417

# TSP demo
TSP_demo<-function(n_city=20, maxItr=1000) {
  loc=FUN_initTSP(n_city);route=1:n_city
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=50; 
  ##
  optDist=FUN_costTSP(route,loc)
  optHist=matrix(0,nrow=maxItr,ncol=(length(route)+1))
  optHist[1,]=c(route,optDist)
  for (i_loop in 2:maxItr) {
    rand.op=sample(c('inv','sw','trans'),1,prob=c(0.75,0.125,0.125))
    if (rand.op=='inv') {
      new.route=TSP_inv(route,sample(2:n_city,1))
    } else if (rand.op=='sw') {
      new.route=TSP_switch(route)
    } else { new.route=TSP_trans(route)}
    new.dist=FUN_costTSP(new.route,loc)
    delta=new.dist-optDist
    prob=1/(1+exp(delta/(C*TEMP)))
    if (runif(1) < prob) {
      route=new.route;optDist=new.dist;
    } 
    optHist[i_loop,]=c(route,optDist);
    TEMP=TEMP*eta
  }
  par(mfrow=c(1,2)) 
  plot(rbind(loc,loc[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(loc[optHist[1000,c(1:n_city,1)],],type='o',pch=20, col='magenta',
    cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')
  return(optHist)
}

TSP_plot

# 古いのバージョン
# initializing cities
FUN_initTSP<-function(n_city=10) {
  return(matrix(runif(n_city*2,1,100),nrow=n_city,ncol=2))
}
# calculating total (Euclidean) distance
FUN_costTSP<-function(route,cities) {
  route=c(route,route[1]);n_cities=nrow(cities);totDist=0;
  for (i_cities in 1:n_cities) {
    totDist=totDist+dist(cities[route[i_cities:(i_cities+1)],])
  }
  return(totDist)
}

# route inversion
TSP_inv<-function(route,len) {
  len_route=length(route)
  invP=sample(1:len_route,1)
  route[invP:min(len_route,(invP+len-1))]=rev(route[invP:min(len_route,(invP+len-1))])
  return(route)
}

# route switching
TSP_switch<-function(route){
  len_route=length(route)
  switchP=sample(1:len_route,2)
  route[switchP]=route[rev(switchP)]
  return(route)
}

# route translation / insertion
TSP_trans<-function(route){
  len_route=length(route)
  transP=sample(1:len_route,2);
  tempR=route[-transP[1]]
  if (transP[2]==1){
    tempR=c(route[transP[1]],tempR)
  } else if (transP[2]==len_route) {
    tempR=c(tempR,route[transP[1]])
  } else {
    tempR=c(tempR[1:(transP[2]-1)],route[transP[1]],tempR[(transP[2]):(len_route-1)])
  }
  return(tempR)
}

## main TSP script ###
n_city=20;loc=FUN_initTSP(n_city);route=1:n_city
maxItr=5000;C=1;eta=0.99;TEMP=20;
optDist=FUN_costTSP(route,loc)
optHist=matrix(0,nrow=maxItr,ncol=(length(route)+1))
optHist[1,]=c(route,optDist)
for (i_loop in 2:maxItr) {
  rand.op=sample(c('inv','sw','trans'),1,prob=c(0.75,0.125,0.125))
  if (rand.op=='inv') {
    new.route=TSP_inv(route,sample(2:n_city,1))
  } else if (rand.op=='sw') {
    new.route=TSP_switch(route)
  } else { new.route=TSP_trans(route)}
  new.dist=FUN_costTSP(new.route,loc)
  delta=new.dist-optDist
  prob=1/(1+exp(delta/(C*TEMP)))
  if (runif(1) < prob) {
    route=new.route;optDist=new.dist;
  } 
  optHist[i_loop,]=c(route,optDist);
  TEMP=TEMP*eta
}
par(mfrow=c(1,2)) 
plot(rbind(loc,loc[1,]),type='b',pch=20,cex=2,
  xlab='location @X',ylab='location @Y',main='Initial route')
plot(loc[optHist[1000,c(1:n_city,1)],],type='b',pch=20,
  cex=2,xlab='location @X',ylab='location @Y',main='Optimized route')

## DEMO function ###
TSP_demo<-function(n_city=20, maxItr=1000) {
  loc=FUN_initTSP(n_city);route=1:n_city
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=20; 
  ##
  optDist=FUN_costTSP(route,loc)
  optHist=matrix(0,nrow=maxItr,ncol=(length(route)+1))
  optHist[1,]=c(route,optDist)
  for (i_loop in 2:maxItr) {
    rand.op=sample(c('inv','sw','trans'),1,prob=c(0.75,0.125,0.125))
    if (rand.op=='inv') {
      new.route=TSP_inv(route,sample(2:n_city,1))
    } else if (rand.op=='sw') {
      new.route=TSP_switch(route)
    } else { new.route=TSP_trans(route)}
    new.dist=FUN_costTSP(new.route,loc)
    delta=new.dist-optDist
    prob=1/(1+exp(delta/(C*TEMP)))
    if (runif(1) < prob) {
      route=new.route;optDist=new.dist;
    } 
    optHist[i_loop,]=c(route,optDist);
    TEMP=TEMP*eta
  }
  par(mfrow=c(1,2)) 
  plot(rbind(loc,loc[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(loc[optHist[1000,c(1:n_city,1)],],type='o',pch=20, col='magenta',
    cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')
  return(optHist)
}

最適化問題

最適化問題

# 勾配法
tol=0.0001;grad=100;
x=10;hist.x=x;lambda=0.1;
while(grad>tol) {
 grad=(2*x+2)
 x=x-lambda*grad
 hist.x=c(hist.x,x)
}
par(mfrow=c(1,2));
xs=seq(-10,10,0.1)
plot(xs,xs^2+2*xs+1)
plot(hist.x)

# naive stochastic optimization
sx=10;hist.sx=sx
for (i_loop in 1:100){
  sx.temp=sx+rnorm(1);
  if ((sx.temp^2+2*sx.temp+1)<(sx^2+2*sx+1)){
    sx=sx.temp
  }
  hist.sx=c(hist.sx,sx)
}

# simulated annealing
x=10;maxIt=1000;
c=1;s=1;temp=1;eta=0.99;
hist.x=matrix(0,nrow=maxIt);
for (i_loop in 1:maxIt) {
  x.new=x+rnorm(1,mean=0,sd=temp*s)
  E.old=x^2+2*x+1;
  E.new=x.new^2+2*x.new+1;
  Paccept=1/(1+exp((E.new-E.old)/(c*temp)))
  if (Paccept>runif(1)){x=x.new}
  temp=temp*eta;
  hist.x[i_loop]=x;
}

数理社会学:なぜ差別しなくても外国人居住区ができるのか

社会を<モデル>でみる:数理社会学への招待
29章:なぜ差別しなくても外国人居住区ができるのか
○と♯の2つのグループが存在し、以下の条件で他の場所へ移動する。
○: 近隣に2人以下○の場合、
♯: 近隣の1/3以上が♯でない場合、

# 1 epochの内、数回移動する可能性のある場合
FUNmigration<-function(field_size=6, Nsharp=10, Ncircle=10) {
  
Nempty=(field_size^2-Nsharp-Ncircle)
society<-matrix(sample(c(rep(0,Nempty),rep(1,Ncircle),rep(2,Nsharp))),ncol=field_size)
# plotting initial config.
par(mfrow=c(1,2))
idx1<-which(society==1,arr.ind=T);idx2<-which(society==2,arr.ind=T)
plot(idx1[,1],idx1[,2],type="n",xlim=c(0.5,field_size+0.5),
  ylim=c(0.5,field_size+0.5),main="Initial Arrangement",
  xlab="location @ x", ylab="location @ y")
text(idx1[,1],idx1[,2],"o",cex=4,col="red");text(idx2[,1],idx2[,2],"#",cex=4,col="blue")

# main 
moved=1;counter=0
while (moved>0&counter<1000) {
  counter=counter+1;moved=0
  for (i_row in 1:field_size) {
    for (i_col in 1:field_size) {
      # checking sharps
      if (society[i_row,i_col]==2) {
        neR_IDX=max((i_row-1),1):min((i_row+1),field_size);
        neC_IDX=max((i_col-1),1):min((i_col+1),field_size);
        n_sharp=sum(society[neR_IDX,neC_IDX]==2)-1;
        n_circle=sum(society[neR_IDX,neC_IDX]==1);
        if (n_sharp+n_circle==0 | n_sharp/(n_sharp+n_circle) < (1/3)) {
          moved=moved+1;
          loc_mov=sample(which(society==0),1)
          society[i_row,i_col]=0
          society[loc_mov]=2
        }
      }
      # checking circles
     if (society[i_row,i_col]==1) {
       neR_IDX=max((i_row-1),1):min((i_row+1),field_size);
       neC_IDX=max((i_col-1),1):min((i_col+1),field_size);
       n_circle=sum(society[neR_IDX,neC_IDX]==1)-1;
       if (n_circle < 3) {
         moved=moved+1;
         loc_mov=sample(which(society==0),1)
         society[i_row,i_col]=0
         society[loc_mov]=1
        }
      }
    }
  }
}
# plotting final config.
idx1<-which(society==1,arr.ind=T)
idx2<-which(society==2,arr.ind=T)
plot(idx1[,1],idx1[,2],type="n",xlim=c(0.5,field_size+0.5),ylim=c(0.5,field_size+0.5),
  main="Arragement After Migration",,xlab="location @ x", ylab="location @ y")
text(idx1[,1],idx1[,2],"o",cex=4,col="red")
text(idx2[,1],idx2[,2],"#",cex=4,col="blue")
}
# 1 epochの内、1度のみ移動する場合
FUNmigration <- function(soc.size = 6, n.circle = 10, n.sharp = 10) {
  r.sample = sample(soc.size ^ 2)
  society = matrix(0, soc.size, soc.size)
  society[r.sample[1:n.circle]] = 1
  society[r.sample[(n.circle + 1):(n.circle + n.sharp)]] = 2
  
  # plotting initial config.
  par(mfrow = c(1, 2))
  idx1 <- which(society == 1, arr.ind = T)
  idx2 <- which(society == 2, arr.ind = T)
  plot(idx1[, 1], idx1[, 2], type = "n", xlim = c(0.5, soc.size + 0.5),
    ylim = c(0.5, soc.size + 0.5), main = "Initial Arrangement",
    xlab = "location @ x", ylab = "location @ y")
  text(idx1[, 1], idx1[, 2], "o", cex = 3, col = "red")
  text(idx2[, 1], idx2[, 2], "#", cex = 3, col = "blue")
  move = 1
  while (move != 0 ) {
    # circles
    move = 0
    c.pos = which(society == 1, arr.ind = T)
    for (i.c in 1:n.circle) {
      r.idx = c(max(1, c.pos[i.c, 1] - 1), min(soc.size, c.pos[i.c, 1] + 1))
      c.idx = c(max(1, c.pos[i.c, 2] - 1), min(soc.size, c.pos[i.c, 2] + 1))
      neighber  = society[r.idx[1]:r.idx[2], c.idx[1]:c.idx[2]]
      neighber.c = sum(neighber == 1) - 1
      if (neighber.c < 3) {
        move = 1
        move.idx = which(society == 0)
        society[sample(move.idx, 1)] = 1
        society[c.pos[i.c, 1], c.pos[i.c, 2]] = 0
      }
    }
    # sharps
    s.pos = which(society == 2, arr.ind = T)
    for (i.s in 1:n.sharp) {
      r.idx = c(max(1, s.pos[i.s, 1] - 1), min(soc.size, s.pos[i.s, 1] + 1))
      c.idx = c(max(1, s.pos[i.s, 2] - 1), min(soc.size, s.pos[i.s, 2] + 1))
      neighbor  = society[r.idx[1]:r.idx[2], c.idx[1]:c.idx[2]]
      neighbor.s = sum(neighbor == 2) - 1
      neighbor.all = sum(neighbor != 0) - 1
      prop.s = max(0, neighbor.s / neighbor.all, na.rm = T)
      if (prop.s < 1 / 3) {
        move = 1
        move.idx = which(society == 0)
        society[sample(move.idx, 1)] = 2
        society[s.pos[i.s, 1], s.pos[i.s, 2]] = 0
      }
    }
  }
  idx1 <- which(society == 1, arr.ind = T)
  idx2 <- which(society == 2, arr.ind = T)
  plot(idx1[, 1], idx1[, 2], type = "n",
    xlim = c(0.5, soc.size + 0.5), ylim = c(0.5, soc.size + 0.5),
    main = "Arragement After Migration", xlab = "location @ x", ylab = "location @ y")
  text(idx1[, 1], idx1[, 2], "o", cex = 3, col = "red")
  text(idx2[, 1], idx2[, 2], "#", cex = 3, col = "blue")
}