2020 基礎実習 R programming 2

N = 10
M = 10000
means = rep(0, M)
for (i_rep in 1:M){
  dat <- runif(N)
  means[i_rep] <- mean(dat)
}
hist(means, xlab = "Estmatedmeans", probability = T, 
  breaks = 30, xlim = c(0,1))
x.temp = seq(0, 1, length.out = 1000)
dens <- dnorm(x.temp, mean = 0.5, sd = (1/sqrt(12*N)))
lines(x.temp, dens, col = "red", lwd = 3)
legend("topright","theoretical desiity",lwd =3, col = "red")

X = 10
Y = 1
sumXY = X + Y

summation <- function(X,Y){
  sumXY = X + Y
  return(sumXY)
}
summation(X=1, Y=10)
summation(X=5, Y=2)

CLT_example <- function(N, M){
  means = rep(0, M)
  for (i_rep in 1:M){
    dat <- runif(N)
    means[i_rep] <- mean(dat)
  }
  hist(means, xlab = "Estmated means", probability = T, 
   breaks = 30, xlim =c(0,1))
  x.temp = seq(0, 1, length.out = 1000)
  dens <- dnorm(x.temp, mean = 0.5, sd = (1/sqrt(12*N)))
  lines(x.temp, dens, col = "red", lwd = 3)
  legend("topright","theoretical desiity",lwd =3, col = "red")
}

x0 = 1e5
y0 = 10
z0 = 0
dt = 0.01
T = 10000
alpha = 1e-5
beta = 0.1
x = y = z = rep(0,T)
x[1] = x0
y[1] = y0
z[1] = z0
for(t in 1:(T-1)){
  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
}
plot(x=1:T, x, type = "l", col = 2, lwd = 3, ylim = c(0,x0),
  xlab = "Time", ylab = "Number of people")
lines(x=1:T, y, col = 3, lwd = 3)
lines(x=1:T, z, col = 4, lwd = 3)
legend("right", c("not infected", "infected","no longer infected"), col=2:4,lwd=3)

2020 基礎実習 R programming1

# creating example vectors
students = paste("s", 1:20, sep = "")
suite = rep(c("s","h","c","d"),13)
ns = sort(rep(1:13,4))
cards  = paste(suite, ns, sep="")

# samples
rand_select = sample(1:10, 2)
rand_perm = sample(1:10)
with_replacemnt = sample(0:1, 10, replace = T)
temp = sample(1:5, 5, replace = T)

# concrete example

# random selection
pointed <- sample(students,2)

# random permutation
randID= sample(1:20)
h1 = students[randID[1:10]]
h2 = students[randID[11:20]]

h1order = sample(h1)
h2order = sample(h2)

shuffled = sample(cards)
cardPlay = data.frame(p1 = shuffled[1:5], p2 = shuffled[6:10], p3 = shuffled[11:15])

# with replacement

# for loop 
num = 1:5
chr = c("a","b","c","d","e")
for (i in num){
  print(i)
}
for (i in chr){
  print(i)
}
for (i in num){
  print(chr[i])
}

p1 = p2 = p3 = rep("joker",5)
Ncard = 5
Nplayer = 3
cardSeq = seq(1, Ncard*Nplayer, Nplayer)
cardN = 0
for (i_card in cardSeq){
  cardN = cardN + 1
  p1[cardN] = shuffled[i_card]
  p2[cardN] = shuffled[i_card+1]
  p3[cardN] = shuffled[i_card+2]
}

p1 = shuffled[seq(1,15,3)]
p2 = shuffled[seq(2,15,3)]
p3 = shuffled[seq(3,15,3)]

players = matrix("joker",nrow = Ncard, ncol = Nplayer)
colnames(players) <- c("p1","p2","p3")
cardN = 0
for (i_card in 1:Ncard){
  for (i_player in 1:Nplayer){
    cardN = cardN + 1
    players[i_card, i_player] = shuffled[cardN]
  }
}

players$p1 = shuffled[seq(1,15,3)]
players$p2 = shuffled[seq(2,15,3)]
players$p3 = shuffled[seq(3,15,3)]



p1WinC = p2WinC = tie = 0
Nplay = 10
for (i_play in 1:Nplay){
  p1 = sample(1:6,1)
  p2 = sample(1:6,1)
  if (p1 > p2){
    p1WinC = p1WinC + 1
  } else {
    if (p2 > p1){
      p2WinC = p2WinC + 1
    } else {
      tie = tie + 1 
    }
  }
}

result = sample(c(0,1,2),10, replace=T, prob = c(2/12,5/12,5/12))

win_thres = 5
p1WinC = p2WinC = tie = 0

while(max(c(p1WinC, p2WinC)) < win_thres){
  p1 = sample(1:6,1)
  p2 = sample(1:6,1)
  if (p1 > p2){
    p1WinC = p1WinC + 1
  } else {
    if (p2 > p1){
      p2WinC = p2WinC + 1
    }
  }
}
if (p1WinC>p2WinC){
  print("player 1 won")
} else {
    print("player 2 won")
}


win_thres = 5
p1WinC = p2WinC = tie = 0

repeat{
  p1 = sample(1:6,1)
  p2 = sample(1:6,1)
  if (p1 > p2){
    p1WinC = p1WinC + 1
  } else {
    if (p2 > p1){
      p2WinC = p2WinC + 1
    }
  }
  if(max(c(p1WinC, p2WinC)) >= win_thres){
    if (p1WinC>p2WinC){
      print("player 1 won")
    } else {
      print("player 2 won")
    }
    break
  }
}


vec = sample(1:15)
len_vec = length(vec)
#vec = 1:15
for (loop1 in 1:(len_vec-1)){
  for (loop2 in 2:(len_vec - loop1 + 1)){
    if (vec[loop2] < vec[(loop2 - 1)]){
      temp_num = vec[loop2]
      vec[loop2] = vec[(loop2-1)]
      vec[(loop2-1)] = temp_num
    }
  }
  print(paste("result after loop1 = ", loop1))
  print(vec)
}

vec = sample(1:15)
len_vec = length(vec)
#vec = 1:15
for (loop1 in 1:(len_vec-1)){
  if (all(1:15 == vec)) {
    print("sorted")
    print(vec)
    break
  }
  for (loop2 in 2:(len_vec - loop1 + 1)){
    if (vec[loop2] < vec[(loop2 - 1)]){
      temp_num = vec[loop2]
      vec[loop2] = vec[(loop2-1)]
      vec[(loop2-1)] = temp_num
    }
  }
  print(paste("result after loop1 = ", loop1))
  print(vec)
}

基礎実習 B02

bin2dec <- function(bin){
# convert binary to decimal
   b.rev = rev(bin)
   ones = which(b.rev == 1)
   ones.adj = ones - 1
   b.sq = 2^ones.adj
   dec = sum(b.sq)
   return(dec)
}


dec2bin <- function(dec, digit) {
# convert decimal to binary
bin = rep(0, digit)
r.counter  = 1 
   while(dec != 0){
     r = dec %% 2
     dec = dec %/% 2
     bin[r.counter] = r
     r.counter = r.counter + 1
  }
  return(rev(bin))
}

# cellur automata
# RULE 150
time = 100
n.cells =  101
cells = matrix(0, nrow = time, ncol = n.cells)
state = rep(0,n.cells)
state[51] = 1
cells[1,] = state
for (i.time in 2:time) {
   left = c(state[n.cells], state[1:(n.cells - 1)])
   right = c(state[2:n.cells], state[1])
   s = left + state + right
   state = s %% 2 
   cells[i.time, ] = state
}

time = 100
n.cells =  101
cells = matrix(0, nrow = time, ncol = n.cells)
cells[1,51] = 1
for (i.time in 2:time) {
   left = c(cells[(i.time-1), n.cells], cells[(i.time-1), 1:(n.cells - 1)])
   right = c(cells[(i.time-1), 2:n.cells], cells[(i.time-1), 1])
   s = left + cells[(i.time-1),] + right
   cells[i.time, ] = s %% 2 
}

transFUN<-function(st,ruleID){
  output=dec2bin(ruleID,8);
  a=matrix(c(1,1,1,
             1,1,0,
             1,0,1,
             1,0,0,
             0,1,1,
             0,1,0,
             0,0,1,
             0,0,0),nrow=8,byrow=T)
  newSt=output[which(apply(a,1,function(x) {all(x==st)}))]
  return(newSt)
}

ECA<-function(nCell, nGen,ruleID){
  res=matrix(0,nrow=nGen,ncol=nCell)
  res[1,ceiling(nCell/2)]=1;
  for (i_gen in 2:nGen) {
    for (i_cell in 2:(nCell-1)) {
      res[i_gen,i_cell]=transFUN(res[i_gen-1,(i_cell-1):(i_cell+1)],ruleID)
     }
   res[i_gen,1]=transFUN(c(res[i_gen-1,nCell],
                           res[i_gen-1,1],
                           res[i_gen-1,2]),ruleID)
   res[i_gen,nCell]=transFUN(c(res[i_gen-1,(nCell-1)],
                               res[i_gen-1,nCell],
                               res[i_gen-1,1]),ruleID)
  }
  return(res)
}
res<-ECA(200,100,99)
image(res) 


基礎実習 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, inits, max_time){
  dt=0.01
  max_timeR = max_time - 1
  x=c(inits[1], rep(0, max_timeR))
  y=c(inits[2], rep(0, max_timeR))
  z=c(inits[3], rep(0, max_timeR))
  # main loop
  for (t in 1:max_timeR) {
    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, c(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)