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)

# 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)
}

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) 

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"

# 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
}

# 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. タカ、無し（０）
ハト 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)

# 実装例、その２
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)

読みやすい？例

## 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回
各最適化法の結果の平均を可視化して（下の図のようなもの）傾向を考察してみてください。
各最適化は下のものを参考にしてください。
各最適化法のパラメターを変えてみてましょう。

 #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)