サイコロを３つなげ、その和が何であるかをかけるギャンブルがあるそうです。
和が、4~10ならsmall, 11~17ならlargeとおおよそ２択のギャンブルだそうです。（3をx、18をyと呼ぶことにします）
Kくんが体験したのは「過去50回の結果がlargeが70％と表示されていたので、次はsmallじゃないかと思いsmallに賭けて、実際にsmallが出た」ということでした。この出来事がどの程度よく起きるのかシミュレーションを用いて検証しまししょう。

# 非効率だけど、直感的に分かりやすいと思われる例

# 1試行の関数
die <- function() {
  die1<-sample(1:6,1)
  die2<-sample(1:6,1)
  die3<-sample(1:6,1)
  die.sum = die1 + die2 + die3
  if (die.sum==3) {
     result="X"
  } else { 
    if (die.sum==18) {
      result="Y"
   } else {
     if (die.sum<11){ 
       result="small"
       } else {
         result="large"
       }
     }
   }
  return(result)
}

# 51試行の結果を表す関数
SL51=function( ) {
  n.rep=51
  res=rep("N",51)
  for (i_rep in 1:n.rep) {
    res[i_rep]=die()
  }
  p.large=sum(res[1:50]=="large")/50
  return(list(p.large=p.large,last.result=res[51]))
}
  
# 51試行をN.REP回繰り返す関数
ck.kurimoto <- function(n.rep) {
  p.hist = rep(0,n.rep)
  lr.hist = rep("N",n.rep)
  for (i_rep in 1:n.rep) {
     result=SL51()
     p.hist[i_rep] = result$p.large
     lr.hist[i_rep] = result$last.result
  }
  return(list(p.hist=p.hist,lr.hist=lr.hist))
}

# 50試行の内、Lが出た割合が70%以上のものを選択
index70=which(m.result$p.hist>=0.7)
sum(m.result$lr.hist[index70]=="small")/length(index70)

# 上の３つの関数をまとめたもの
n.rep=1e6
die.mat = matrix(sample(c("X","Y","S","L"),
 n.rep*51,prob=c(1/16,1/16,7/16,7/16),replace=T),nrow=n.rep)
p.large=rowSums(die.mat[,1:50]=="L")/50
index70=which(p.large>=0.7)
ck51L70=die.mat[index70,51]
table(ck51L70)/length(index70)

# 直感的にはわかりやすいけど、非効率な例
bet.example<-function(n.trial,p.win,max.bet) {
  rew.history=rep(0,n.trial)
  bet.history=rep(0,n.trial)
  WL.history=rep("W",n.trial)
  loss.counter=0
  i_trial=0;
  WL="L"
  loss.counter.reset="F"
  while (i_trial < n.trial | WL=="L" ) {
    i_trial=i_trial+1
    bet=2^loss.counter
    if (bet > max.bet) {
      bet = max.bet
      loss.counter.reset = "T"
    }
    bet.history[i_trial]=bet
    WL=sample(c("W","L"),1,prob=c(p.win, 1-p.win))
    WL.history[i_trial]=WL
    if (WL=="L") {
      loss.counter=loss.counter+1
      reward=-bet
    } else {
      loss.counter=0
      reward=bet
    }
    if (loss.counter.reset=="T"){loss.counter=0}
    rew.history[i_trial]=reward
  }
  return(list(reward=sum(rew.history),money.required=max(bet.history)))
}

# より効率的な例 - max.betは未導入
bet.example02<-function(n.trial,p.win) {
  WL=sample(c("W","L"),n.trial,replace=T,prob=c(p.win, 1-p.win))
  if (WL[n.trial]!="W") {
    WL.extra=sample(c("W","L"),1,prob=c(p.win,1-p.win))
    while (tail(WL.extra,1)=="L") {
      WL.extra=c(WL.extra,sample(c("W","L"),1,prob=c(p.win,1-p.win)))
     }
    WL=c(WL,WL.extra)
  }
  reward=sum(WL=="W")
  result=rle(WL)
  loss.index=result$values=="L"
  money.required=2^max(result$lengths[loss.index])
  return(list(reward,money.required))
}

traveling salesman problem
来週までに、inversionとtranslationを行う関数を作っておいてください。

# initialisation, etc.
n.city=20
location=matrix(runif(n.city*2,min=0,max=100),nrow=n.city)
route=sample(n.city) 
calc.Dist<-function(location,route) {
  n.city=length(route)
  total.Dist=0
  route=c(route,route[1])
  for (i_city in 1:n.city){
    total.Dist=total.Dist+dist(location[c(route[i_city],route[i_city+1]),])
  }
  return(total.Dist)
}

# plotting results
plot(rbind(location[route,],location[route[1],]),type='o',pch=20,cex=2.5, col='red',
  xlab='location @X',ylab='location @Y',main='Initial route')
plot(rbind(location[best.route,],location[best.route[1],]),type='o',pch=20, 
  col='magenta',cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')

#switching 
TSP.switch<-function(route) {
  two.cities=sample(route,2,replace=F)
  route[two.cities]=route[rev(two.cities)]
  return(route)
}

# inversion
TSP.inv<-function(route,inv.length) { 
  inv.begin=sample(route,1)
  inv.end=min(inv.begin+inv.length-1,length(route))
  route[inv.begin:inv.end]=rev(route[inv.begin:inv.end])
  return(route)
}

# translation
TSP.trans<-function(route,tr.length) {
  trP1=sample(route,1)
  tr.vec=route[trP1:min(length(route),trP1+tr.length-1)]
  temp.vec=route[-(trP1:min(length(route),trP1+tr.length-1))]
  trP2=sample(1:length(temp.vec),1)
  if (trP2==length(temp.vec)) {
    new=c(temp.vec,tr.vec)
  } else {
    new=c(temp.vec[1:trP2],tr.vec,temp.vec[(trP2+1):length(temp.vec)])
  }
return(new)
}

# demo
TSP.demo<-function(n.city=20, maxItr=1000) {
  location=matrix(runif(n.city*2,min=0,max=100),nrow=n.city)
  route=sample(n.city) 
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=50; 
  ##
  optDist=calc.Dist(location,route)
  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,sample(2:(round(n.city/2)),1))}
    new.dist=calc.Dist(location,new.route)
    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(location,location[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(location[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)
}

最適化問題

 
#Objective Function 
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;
xs=ys<-seq(xmin, xmax,length.out=n_len)
z<-outer(xs,ys,funZ)
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))

勾配法の実装例

 
# initialization
x=x.hist=runif(1,min=-5,max=5) # initial value of x
y=y.hist=runif(1,min=-5,max=5) # initial value of y
lambda=0.01                    # step size
zeta=1e-8                      # tol (stopping criterion)
t=0                            # time
g=100                          # maximum gradient
# end initialization

### note ###
# obj fun: z=x^4-16*x^2-5*x+y^4-16*y^2-5*y
# part. deriv. x: 4*x^3-32*x-5
# part. deriv. y: 4*y^3-32*y-5
#############

# main 
while (g>zeta) {
  t=t+1
  g.x=4*x^3-32*x-5
  g.y=4*y^3-32*y-5
  x=x-lambda*g.x
  y=y-lambda*g.y
  g=max(abs(g.x),abs(g.y))
  x.hist=c(x.hist,x)
  y.hist=c(y.hist,y)
}

# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

単純確率的最適化法の実装例

 
# initialization
x=runif(1,min=-5,max=5)             # initial value of x
y=runif(1,min=-5,max=5)             # initial value of y
f.xy=x^4-16*x^2-5*x+y^4-16*y^2-5*y  # initial value of obj. fun.
w=0.1                               # search width
t=0                                 # time
max.t=1000                          # max iter (stopping criterion)
stop.crit=F                         # stop criterion
x.hist=rep(0,max.t)                 # history of x
y.hist=rep(0,max.t)                 # history of x
x.hist[1]=x
y.hist[1]=y
# end initialization

# main
while (stop.crit==F) {
  t=t+1
  x.temp=x+rnorm(1,mean=0,sd=w)
  y.temp=y+rnorm(1,mean=0,sd=w)
  f.temp=funZ(x.temp,y.temp)
  if (f.temp < f.xy) {
    x=x.temp
    y=y.temp
    f.xy=f.temp
  }
  x.hist[t]=x
  y.hist[t]=y
  if (t==max.t){stop.crit=T}
}

# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

焼きなまし法の実装例

# initialization
x=runif(1,min=-5,max=5)             # initial value of x
y=runif(1,min=-5,max=5)             # initial value of y
f.xy=x^4-16*x^2-5*x+y^4-16*y^2-5*y  # initial value of obj. fun.
t=0                                 # time
max.t=1000                          # max iter (stopping criterion)
stop.crit=F                         # stop criterion
x.hist=rep(0,max.t)                 # history of x
y.hist=rep(0,max.t)                 # history of x
x.hist[1]=x
y.hist[1]=y
tau=1
c=1
g=0.999
# end initialization
 
# main
while (stop.crit==F) {
  t=t+1
  x.temp=x+rnorm(1,mean=0,sd=tau)
  y.temp=y+rnorm(1,mean=0,sd=tau)
  f.temp=funZ(x.temp,y.temp)
  deltaE=f.temp-f.xy
  p=1/(1+exp(deltaE/(c*tau)))
  if (runif(1) < p) {
    x=x.temp
    y=y.temp
    f.xy=f.temp
  }
  x.hist[t]=x
  y.hist[t]=y
  tau=tau*g
  if (t==max.t){stop.crit=T}
}
 
# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

３つの手法の比較

 

GD<-function(x.init,y.init,n.iter,lambda){
x=x.init;y=y.init;z.hist=rep(0,n.iter)
# main 
for (i_iter in 2:n.iter) {
  g.x=4*x^3-32*x-5
  g.y=4*y^3-32*y-5
  x=x-lambda*g.x
  y=y-lambda*g.y
  g=max(abs(g.x),abs(g.y))
  x.hist=c(x.hist,x)
  y.hist=c(y.hist,y)
  z.hist[i_iter]=funZ(x,y)
}
return(z.hist)
}

stochOpt<-function(x.init,y.init,n.iter,w) {
x=x.init;y=y.init;f.xy=funZ(x,y)
z.hist=rep(0,n.iter);z.hist[1]=f.xy
# main
for (i_iter in 2:n.iter) {
  x.temp=x+rnorm(1,mean=0,sd=w)
  y.temp=y+rnorm(1,mean=0,sd=w)
  f.temp=funZ(x.temp,y.temp)
  if (f.temp < f.xy) {
    x=x.temp;y=y.temp;f.xy=f.temp
  }
  z.hist[i_iter]=f.xy
}
return(z.hist)
}
 
simAnn<-function(x.init,y.init,n.iter,tau,c,g) {
x=x.init;y=y.init;f.xy=funZ(x,y)
z.hist=rep(0,n.iter);z.hist[1]=f.xy
# main
for (i_iter in 2:n.iter) {
  x.temp=x+rnorm(1,mean=0,sd=tau)
  y.temp=y+rnorm(1,mean=0,sd=tau)
  f.temp=funZ(x.temp,y.temp)
  deltaE=f.temp-f.xy
  p=1/(1+exp(deltaE/(c*tau)))
  if (runif(1) < p) {
    x=x.temp;y=y.temp;f.xy=f.temp
  }
  z.hist[i_iter]=f.xy
  tau=tau*g
}
return(z.hist)
}

n.rep=500
n.iter=500
xs=runif(n.rep,-3,3);ys=runif(n.rep,-3,3)
GD.hist=SA.hist=SO.hist=matrix(0,n.iter,n.rep)
for (i.rep in 1:n.rep) {
  GD.hist[,i.rep]=GD(xs[i.rep],ys[i.rep],n.iter,0.025)
  SO.hist[,i.rep]=stochOpt(xs[i.rep],ys[i.rep],n.iter,2)
  SA.hist[,i.rep]=simAnn(xs[i.rep],ys[i.rep],n.iter,3,1,0.999)
}

plot(rowMeans(GD.hist),type='l',ylim=c(-160,0),lwd=3)
lines(rowMeans(SO.hist),col='red',lwd=3)
lines(rowMeans(SA.hist),col='blue',lwd=3)

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

例１
#
# 1 epochで複数回移動するエージェントがいてもいいバージョン
# 次回は移動は1epoch内で最大１回とするものを作りましょう。
#

# function to check circles
ck.circle <- function(world,i_row,i_col,max.row,max.col) {
  neighbor=world[max(i_row-1,1):min(i_row+1,max.row),
                 max(i_col-1,1):min(i_col+1,max.col)]
  circles=length(which(neighbor==2))
  act=ifelse(circles<3,"move","stay")
  return(act)
}

# function to check sharps
ck.sharp <- function(world,i_row,i_col,max.row,max.col) {
  neighbor=world[max(i_row-1,1):min(i_row+1,max.row),
                 max(i_col-1,1):min(i_col+1,max.col)]
  all.neighb=length(which(neighbor!=0))-1
  if (all.neighb==0){
    return("move")
  } else {
    sharps=length(which(neighbor==1))-1
    act=ifelse(sharps/all.neighb<(1/3),"move","stay")
    return(act)
  }
}

FUNmigration<-function(N.row=10, N.col=10, N.sharp=10, N.circle=10) {
   
# initialization
world=matrix(0,nrow=N.row,ncol=N.col)
temp.vec=c(rep(1,N.sharp),rep(2,N.circle))
world[sample(1:(N.row*N.col),N.sharp+N.circle,replace=F)]=temp.vec
moved=1;counter=0

# plotting initial config.
par(mfrow=c(1,2))
idx1<-which(world==2,arr.ind=T);idx2<-which(world==1,arr.ind=T)
plot(idx1[,1],idx1[,2],type="n",xlim=c(0.5,N.row+0.5),
  ylim=c(0.5,N.col+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 loop, max.iter=1000
while (moved>0 & counter<1000) {
  moved=0
  counter=counter+1
  for (i_row in 1:N.row) {
    for (i_col in 1:N.col) {
      if (world[i_row,i_col]==1) {
        act=ck.sharp(world,i_row,i_col,N.row,N.col)
        if (act=="move") {   
          dest=sample(which(world==0),1)
          world[dest]=1
          world[i_row,i_col]=0
          moved=1  
        }
      }
      if (world[i_row,i_col]==2) {
        act=ck.circle(world,i_row,i_col,N.row,N.col)
        if (act == "move") {
         dest=sample(which(world==0),1)
         world[dest]=2
         world[i_row,i_col]=0
         moved=1
        }
}}}}

# plotting result
idx1<-which(world==1,arr.ind=T)
idx2<-which(world==2,arr.ind=T)
plot(idx1[,1],idx1[,2],type="n",xlim=c(0.5,N.row+0.5),ylim=c(0.5,N.col+0.5),
  main="Arragement After Migration",,xlab="location @ x", ylab="location @ y")
text(idx1[,1],idx1[,2],"#",cex=4,col="blue")
text(idx2[,1],idx2[,2],"o",cex=4,col="red")
return(world)
}

例2
#
# 移動は1epoch内で最大１回とするバージョン
#

FUNmigration2<-function(N.row=10, N.col=10, N.sharp=10, N.circle=10) {
   
# initialization
world=matrix(0,nrow=N.row,ncol=N.col)
temp.vec=c(rep(1,N.sharp),rep(2,N.circle))
world[sample(1:(N.row*N.col),N.sharp+N.circle,replace=F)]=temp.vec
moved=1;counter=0

# plotting initial config.
par(mfrow=c(1,2))
idx1<-which(world==2,arr.ind=T);idx2<-which(world==1,arr.ind=T)
plot(idx1[,1],idx1[,2],type="n",xlim=c(0.5,N.row+0.5),
  ylim=c(0.5,N.col+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")

ck.circle <- function(world,i_row,i_col,max.row,max.col) {
  neighbor=world[max(i_row-1,1):min(i_row+1,max.row),
                 max(i_col-1,1):min(i_col+1,max.col)]
  circles=length(which(neighbor<0))
  act=ifelse(circles<3,"move","stay")
  return(act)
}
 
ck.sharp <- function(world,i_row,i_col,max.row,max.col) {
  neighbor=world[max(i_row-1,1):min(i_row+1,max.row),
                 max(i_col-1,1):min(i_col+1,max.col)]
  all.neighb=length(which(neighbor!=0))-1
  if (all.neighb==0){
    return("move")
  } else {
    sharps=length(which(neighbor>0))-1
    act=ifelse(sharps/all.neighb<(1/3),"move","stay")
    return(act)
  }
}

# ここから違う
moved=1;counter=0;N.row=10;N.col=10;
while (moved>0 & counter<1000) {
  moved=0
  counter=counter+1
  for (i_sharp in 1:N.sharp) {
    idx=which(world==i_sharp,arr.ind=T)
    act=ck.sharp(world,idx[1],idx[2],N.row,N.col);
    if (act=="move") {   
      dest=sample(which(world==0),1)
      world[dest]=i_sharp
      world[idx[1],idx[2]]=0
      moved=1  
    }
  }
  for (i_circle in -1:-N.circle) {
    idx=which(world==i_circle,arr.ind=T)
    act=ck.circle(world,idx[1],idx[2],N.row,N.col);
    if (act=="move") {   
      dest=sample(which(world==0),1)
      world[dest]=i_circle
      world[idx[1],idx[2]]=0
      moved=1  
    }  
  }
} 

# plotting result
idx1<-which(world==1,arr.ind=T)
idx2<-which(world==2,arr.ind=T)
plot(idx1[,1],idx1[,2],type="n",xlim=c(0.5,N.row+0.5),ylim=c(0.5,N.col+0.5),
  main="Arragement After Migration",,xlab="location @ x", ylab="location @ y")
text(idx1[,1],idx1[,2],"#",cex=4,col="blue")
text(idx2[,1],idx2[,2],"o",cex=4,col="red")
return(world)
}