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)
Category Archives: 基礎実習&演習
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 }
数理社会学: なぜ宣伝しなくても流行がおこるのか
社会を<モデル>でみる:数理社会学への招待
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)
読みやすい?例
## 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)