基礎実習 MA02

select.act <- function(state){
  poss.act = which(state=="N")
  if (length(poss.act)==1){
    act = poss.act
  } else {
    act = sample(poss.act, 1)
  }
  return(act)
}

ck.term <- function(state) {
  term = F
  result = "undecided"
  term.idx = matrix(c(1,2,3,4,5,6,7,8,9,1,4,7,
                      2,5,8,3,6,9,1,5,9,3,5,7),ncol=3,byrow=T)
  if (sum(state == "N")==0) {
    term=T
    result = "tie"
  }                      
  for (i.ck in 1:8) {
    O.won = all(state[term.idx[i.ck,]]=="O")
    X.won = all(state[term.idx[i.ck,]]=="X")
    if (O.won ==1 ) {
      term= T
      result = "O won"
      break
    }
    if (X.won ==1){
      term=T
      result = "X won"
      break
    }
  }      
  return(list(term = term, result=result))
}

ck.term2 <- function(state) {
  term = F
  result = "undecided"
  term.idx = matrix(c(1,2,3,4,5,6,7,8,9,1,4,7,
                      2,5,8,3,6,9,1,5,9,3,5,7),ncol=3,byrow=T)
  if (sum(state == "N")==0) {
    term=T
    result = "tie"
  }                      
  O.won = apply(term.idx,1,function(x) all(state[x]=="O"))
  X.won = apply(term.idx,1,function(x) all(state[x]=="X"))
  if (any(O.won) == T ) {
    term= T
    result = "O won"
  } else if (any(X.won) ==1){
    term=T
    result = "X won"
  }
  return(list(term = term, result=result))
}

tictactoe <-function(){ 
  coord = matrix(c(1,3,1,2,1,1,2,3,2,2,2,1,3,3,
                   3,2,3,1), nrow = 9, byrow=T)
  plot(0,0, type="n", xlim= c(0.5,3.5), ylim=c(0.5,3.5))
  abline(h = 1.5); abline(h = 2.5)
  abline(v = 1.5); abline(v = 2.5)
  state = rep("N",9); marker = c("O","X")
  repeat { 
    for (i.player in 1:2) {
      act = select.act(state)
      state[act] = marker[i.player]
      text(coord[act,1],coord[act,2], marker[i.player], cex=4)
      res = ck.term2(state)
      if (res$term == T) { break }
      Sys.sleep(0.5)
    }
    if (res$term == T) { 
      print(paste("result:", res$result))
      break 
    }
  }
}

# 実行例
> tictactoe()
[1] "result: X won"

# 対戦版 - 無作為に動くので非常に弱いです。
# 上のselect.actを変更することでマシになると思います。
play.tictactoe.R <-function(){ 
  coord = matrix(c(1,3,1,2,1,1,2,3,2,2,2,1,3,3,
                   3,2,3,1), nrow = 9, byrow=T)
  plot(coord, pch=paste(1:9), col='gray', cex=3, xlim= c(0.5,3.5), ylim=c(0.5,3.5))
  abline(h = 1.5); abline(h = 2.5)
  abline(v = 1.5); abline(v = 2.5)
  state = rep("N",9); 
  repeat { 
    for (i.player in 1:2){
      if (i.player == 1){
        act = as.numeric(readline(prompt="Enter box ID: "))
      } else {act = select.act(state)}
      state[act] = marker[i.player]
      text(coord[act,1],coord[act,2], marker[i.player], cex=4)
      res = ck.term2(state)
      if (res$term == T) { break }
      Sys.sleep(0.5)
    }
    if (res$term == T) { 
      print(paste("result:", res$result))
      break 
    }
  }
}

# 実行例:
> play.tictactoe.R()
Enter box ID: 5
Enter box ID: 1
Enter box ID: 9
[1] "result: O won"

2019 基礎実習MA01

# random number generators
x=rnorm(n=1,mean=100,sd=15)
y=runif(n=3,min=1,max=10)

N = 10000
# N = 1000
random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
  probability = T, main = "Histogram of Random Data")

# density of generated data
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)

# theoretical density
x = seq(-4,4,0.1)  
true.norm = dnorm(x, mean = 0, sd = 1)
lines(x,true.norm, col = "green", lty = 3, lwd = 4)
legend("topleft",c("empirical", "theoretical"), lty = c(1,3),
  col = c('orange','green'),lwd=4)

# random sampling
sample(1:10,3)   
sample(c(“gu”,“choki”,“pa”),1)
sample(1:10)
sample(0:1, 10, replace=T)
sample(c("Head","Tail"), 10, replace=T)
sample(c("Head","Tail"), 10, replace=T, prob=c(0.9,0.1)) 

# flow control
for (i_loop in 1:5){print(i_loop)}
counter <- 1

while(counter<=10){
  print(counter)
  counter<-counter+1
}

counter <- 1
while(counter^2 <= 10){
  print(c(counter, counter^2))
  counter<-counter+1
}

affil<-"cogsci"
if (affil=="cogsci") {
    print("you are wonderful")
}

affil<-"phil"
if (affil=="cogsci") {
   print("you are wonderful")
 } else {
   print("still, you are wonderful")
}

v1=1633;v2=355;
repeat {
  r=v1%%v2
  print(paste('v1 =',v1,v2 = ',v2, remainder = ',r))
  v1=v2;v2=r
  if (r==0){ break}
}

counter=6
repeat{ 
  print(counter)
  counter = counter + 1
  if(counter>5){break}
}

counter=6
repeat{
  if(counter>5){break}
  print(counter)
  counter+counter+1
}

数理社会学: なぜ宣伝しなくても流行がおこるのか

vogue2015plot

社会を<モデル>でみる:数理社会学への招待
18章:なぜ宣伝しなくても流行がおこるのか。
モデルの説明:
x:流行に影響されていない人
y:流行に影響されている人
z:以前流行に影響されてが、もう影響されていない
t:時間
\delta x = -\alpha * x[t] * y[t] * \delta t
\delta y = (\alpha * x[t] * y[t] – \beta * y[t]) * \delta t
\delta z = \beta * y[t] * \delta t

# initialization
vogue<-function(alpha,beta,x.init,y.init,z.init,max_time){
  dt=0.01
  x=c(x.init,rep(0,max_time))
  y=c(y.init,rep(0,max_time))
  z=c(z.init,rep(0,max_time))
 # main loop
  for (t in 1:max_time) {
    x[(t+1)]=x[t]-alpha*x[t]*y[t]*dt
    y[(t+1)]=y[t]+(alpha*x[t]*y[t]-beta*y[t])*dt
    z[(t+1)]=z[t]+beta*y[t]*dt
  }
 return(data.frame(x,y,z))
}
vogue.plot<-function(dat) {
 plot(dat$x,type="l",lwd=4,col="blue", main="Typical Time Series of \"Vogue\" ",
  xlab="time", ylab="Number of People",cex.lab=1.3,ylim=c(0,dat$x[1]))
 lines(dat$y,type="l",lwd=4,col="red")
 lines(dat$z,type="l",lwd=4,col="green")
 legend("right", c("affected","not affected", "no longer affected"), 
  col=c("red","blue","green"),lty=rep(1,3),lwd=4)
}

# running functions
res<-vogue(0.0001,0.1,10000,100,0,5000)
vogue.plot(res)

CLT, LLN, GCD

# CLT
nSample=10;nRep=10^5;
x=matrix(runif(nSample*nRep),nrow=nSample);
x.means<-colMeans(x)
hist(x.means,50,main='Distribution of Means of Uniform Distribution', 
 xlab='Means', probability=T)
x.means.dens<-density(x.means)
lines(x.means.dens,lwd=3,lty=1,col='red')
x2=seq(0,1,0.001);CLT=dnorm(x2,mean=0.5,sd=(sqrt(1/12))/(sqrt(nSample)))
lines(x2,CLT,lwd=3,lty=3,col='cyan')
legend("topright",c("Density of Actual Means","Normal Distribution"), 
 lty=c(1,3), col=c('red','cyan'),lwd=3)

# LAW of LARGE NUMBER
dat<-sample(1:6,1000,replace=T)
prob6<-cumsum(dat==6)/(1:1000)
par(mfrow=c(2,1))
plot(dat,type='p',col="red",main="results of 1000 rolls of a die",xlab="rolls",
 ylab="outcome", pch=20)
plot(prob6,type='l',main="Cumulative probability that rolls of a die come out SIX",
 xlab="rolls", ylab="probability",lwd=3,ylim=c(0.0,0.5))
abline(h=1/6,col='red',lwd=2,lty=2)


# WHILE
r=-99;v1=1633;v2=355
while (r!=0){
  r=v1%%v2
  print(paste('v1 =',v1,', v2 = ',v2,',remainder = ',r))
  v1=v2
  v2=r
}

# REPEAT
v1=1633;v2=355;
repeat {
  r=v1%%v2
  print(paste("v1 =",v1,",v2 = ",v2,",remainder = ",r))
  v1=v2;v2=r
  if (r==0){ break}
 }

タカハトゲーム 

タカ・ハトゲーム
タカ戦略とハト戦略のみが存在する世界において進化的に安定な状態を探るゲーム。
### 設定
# pay-offは以下の通り:
タカ vs. ハト、餌(food)を総取り
タカ vs. タカ、餌からコストを引いたものを2分割
ハト vs. タカ、無し(0)
ハト vs. ハト、餌を2分割
# 引数:
個体数の初期値([タカ、ハト])、世代数、コスト、餌
# 個体数の推移
p.hawk[i_gen]=(p.hawk[i_gen-1]*(p.hawk[i_gen-1]*eHawkHawk+p.dove[i_gen-1]*eHawkDove))
 /mean.W;

HDgame<-function(N,generation,cost,food){
  N.hawk=N[1]; N.dove=N[2]; 
  p.hawk=rep(0,generation);p.hawk[1]=N.hawk/(N.dove+N.hawk)
  p.dove=rep(0,generation);p.dove[1]=1-p.hawk[1]
  eHawkHawk=0.5*(food-cost);eHawkDove=food
  eDoveHawk=0;eDoveDove=0.5*food
  for (i_gen in 2:generation){
    mean.W=p.hawk[i_gen-1]^2*eHawkHawk+p.hawk[i_gen-1]*p.dove[i_gen-1]*
     (eHawkDove+eDoveHawk)+p.dove[i_gen-1]^2*eDoveDove
    p.hawk[i_gen]=(p.hawk[i_gen-1]*(p.hawk[i_gen-1]*eHawkHawk
      +p.dove[i_gen-1]*eHawkDove))/mean.W;
    p.dove[i_gen]=1-p.hawk[i_gen]
  }
  plot(1:generation,p.hawk,type='o',col='red',ylim=c(0,1),pch=19,
    main="Result of Hawk-Dove Game",ylab="Proportion",xlab="Generation")
  lines(1:generation, p.dove,type='o',col='blue',pch=17)
  legend("topright",c("Hawk","Dove"),pch=c(19,17),col=c('red','blue'))
}
# 実行例
HDgame(c(20,80),50,10,6)

# 実装例、その2
HDgameDE<-function(N,generation,cost,food){
  N.hawk=N[1]; N.dove=N[2]; 
  tStep=0.01;ts=seq(0,generation,tStep);Nts=length(ts);
  p.hawk=rep(0,Nts);p.hawk[1]=N.hawk/(N.dove+N.hawk)
  p.dove=rep(0,Nts);p.dove[1]=1-p.hawk[1]
  eHawkHawk=0.5*(food-cost);eHawkDove=food
  eDoveHawk=0;eDoveDove=0.5*food
  for (i_gen in 2:Nts){
    mean.W=p.hawk[i_gen-1]^2*eHawkHawk+p.hawk[i_gen-1]*p.dove[i_gen-1]*
     (eHawkDove+eDoveHawk)+p.dove[i_gen-1]^2*eDoveDove
    p.hawk[i_gen]=p.hawk[i_gen-1]+(p.hawk[i_gen-1]*(p.hawk[i_gen-1]*eHawkHawk
     +p.dove[i_gen-1]*eHawkDove-mean.W)/mean.W)*tStep;
    p.dove[i_gen]=1-p.hawk[i_gen]
  }
  plot(1:Nts,p.hawk,type='l',col='red',ylim=c(0,1),lwd=5,
   main="Result of Hawk-Dove Game",ylab="Proportion",xlab="Generation")
  lines(1:Nts, p.dove,type='l',col='blue',lwd=5)
  legend("topright",c("Hawk","Dove"),lwd=5,col=c('red','blue'))
}
# 実行例
HDgameDE(c(20,80),50,10,6)

D_H01

読みやすい?例

## a simple version ##
HDgame<-function(N,food,cost,n_gen){
  p.hawk=N[1]/sum(N);p.dove=N[2]/sum(N)
  e.HH=1/2*(food-cost);e.HD=food;e.DH=0;e.DD=1/2*(food)
  p.vec=c(p.hawk,p.dove);p.mat<-outer(p.vec,p.vec)
  payMat<-matrix(c(e.HH,e.DH,e.HD,e.DD),ncol=2)
  hist=matrix(0,n_gen,2);hist[1,]=p.vec
  for (i_gen in 2:n_gen){
    w.H=sum(p.vec*payMat[1,])
    w.D=sum(p.vec*payMat[2,])
    w.mean=sum(p.mat*payMat)
    p.vec=c(p.vec[1]*w.H/w.mean,p.vec[2]*w.D/w.mean)
    p.mat<-outer(p.vec,p.vec)
    hist[i_gen,]=p.vec
  }
  plot(1:n_gen,hist[,1],pch=20,type='o',lwd=2,cex=2,col='red',
   ylim=c(0,1),ylim=c(0,1.25),xlab="time",ylab="proportion",cex.lab=2)
  lines(1:n_gen,hist[,2],pch=20,type='o',lwd=2,cex=2,col='blue')
  legend("topright",c("Hawk","Dove"),col=c('red','blue'),lwd=2,pch=20)
}

## DE version ##
HDgameDE<-function(N,food,cost,n_gen){
  p.hawk=N[1]/sum(N);p.dove=N[2]/sum(N)
  e.HH=1/2*(food-cost);e.HD=food;e.DH=0;e.DD=1/2*(food)
  p.vec=c(p.hawk,p.dove);p.mat<-outer(p.vec,p.vec)
  payMat<-matrix(c(e.HH,e.DH,e.HD,e.DD),ncol=2)
  ts=0.01;Nts=length(seq(1,n_gen,ts));
  hist=matrix(0,Nts,2);hist[1,]=p.vec
  for (i_gen in 2:Nts){
    w.H=sum(p.vec*payMat[1,])
    w.D=sum(p.vec*payMat[2,])
    w.mean=sum(p.mat*payMat)
    p.vec[1]=p.vec[1]+ts*p.vec[1]*(w.H-w.mean)/w.mean
    p.vec[2]=p.vec[2]+ts*p.vec[2]*(w.D-w.mean)/w.mean
    p.mat<-outer(p.vec,p.vec)
    hist[i_gen,]=p.vec
  }
  plot(1:Nts,hist[,1],type='l',lwd=4,cex=2,col='red',ylim=c(0,1.25),
    xlab="time",ylab="proportion",cex.lab=2)
  lines(1:Nts,hist[,2],type='l',lwd=4,cex=2,col='blue')
  legend("topright",c("Hawk","Dove"),col=c('red','blue'),lwd=4)
  return(hist)
}

最適化法の比較

認知情報科学基礎実習 課題03
提出期限:2013.01.31 23:59

勾配法、単純確率的最適化法、焼き鈍し法を比較してみましょう。
目的関数:x^4-16*x^2-5*x+y^4-16*y^2-5*y
初期値:一様分布(-0.5~0.5)に従う乱数
繰り返し数:1000回
各最適化法の結果の平均を可視化して(下の図のようなもの)傾向を考察してみてください。
各最適化は下のものを参考にしてください。
各最適化法のパラメターを変えてみてましょう。

3optims

 #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;
x<-seq(xmin, xmax,length.out=n_len);y<-x;z<-outer(x,y,funZ)
 #gradient descent 
minGD<-function(initXY=c(0,0), maxItr=1000, stepSize=0.01){
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    x=x-stepSize*(4*x^3-32*x-5);
    y=y-stepSize*(4*y^3-32*x-5);
    z=funZ(x,y);
    opHist[counter,]=c(x,y,z);
  }
  return(opHist[1:counter,])
}

res<-minGD(c(0.5,0.5),1000,0.01)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(res[,1:2],col='cyan',type='o',pch=19)
#Simple Stochastic Optimization
minSTOCH<-function(initXY=c(0,0),maxItr=1000,width=0.1) {
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    xTemp=x+rnorm(1,mean=0,sd=width);
    yTemp=y+rnorm(1,mean=0,sd=width);
    zTemp=funZ(xTemp,yTemp);
    if (z > zTemp) {
  	x=xTemp;y=yTemp;z=zTemp;
    }
    opHist[i_loop,]=c(x,y,z);
  }
  return(opHist)
}

opHist<-minSTOCH(c(0,0),1000,0.1)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(opHist[,1:2],type='o',lwd=2,col='green',pch=20)
#Simulated Annealing
minSA<-function(initXY=c(0,0),maxItr=1000,C=1,eta=0.99,width=10) {
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    xTemp=x+rnorm(1,mean=0,sd=width)
    yTemp=y+rnorm(1,mean=0,sd=width)
    zTemp=funZ(xTemp, yTemp);
    delta=zTemp-z;
    prob=1/(1+exp(delta/(C*width)))
    if (runif(1) < prob) {
      x=xTemp;y=yTemp;z=zTemp;
     } 
    opHist[i_loop,]=c(x,y,z);
    width=width*eta
  }
  return(opHist)
}

res<-minSA(c(0,0),1000,1,0.99,10)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(res[,1:2],type='o',lwd=2,col='green',pch=20)