# 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 programming１

# 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)  # 基礎実習 MA0２ 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 }  # 数理社会学: なぜ宣伝しなくても流行がおこるのか 社会を＜モデル＞でみる：数理社会学への招待 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. タカ、無し（０）
ハト 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)
}


# 最適化法の比較

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