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
Monthly Archives: July 2015
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)
}
# 古いのバージョン
# 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章：なぜ差別しなくても外国人居住区ができるのか
○と♯の２つのグループが存在し、以下の条件で他の場所へ移動する。
○: 近隣に２人以下○の場合、
♯: 近隣の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の内、１度のみ移動する場合
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")
}