rep_monte.R

library(animation)

#N <- 10
tmax <- 1e4
freq_ex <- tmax/100

theta <- 0.5

plotrange <- c(-1,1)
cols <- heat.colors(N)

Mu <- .3
mu1 <- c(Mu,Mu)
mu2 <- c(-Mu,-Mu)

vc <- .1
v.run <- 0.01

gen_gaussian_mix<- function(mu1,mu2,theta,v){
  function(x){
    theta*dnorm(dist(rbind(x,mu1)),0,v)+(1-theta)*dnorm(dist(rbind(x,mu2)),0,v)
  }
}
gen_gaussian_mix_v<- function(v){
  gen_gaussian_mix(mu1,mu2,theta,v)
}

gen_gaussian_mix_a<- function(mu){
  gen_gaussian_mix(mu,-mu,theta,vc)
}

metropolis<-function(x,pdf,v){
  xn <- x+rnorm(length(x),0,v)  
  if(runif(1,0,1)<pdf(xn)/pdf(x))
    x <- xn
  x
}

rep_monte <- function(tmax=tmax,N=1,show=F){
#  vars <- vc*(1:N)
#  pdf <- lapply(vars,gen_gaussian_mix_v)
  ms<-seq(Mu,0,length.out=N)
  mm<-cbind(ms,-ms)
  pdf <- apply(mm,1,gen_gaussian_mix_a)
  
  x <- matrix(runif(N*2,-0.4,-0.2),N,2)
  xx<- mean(x[,1])

  for(t in (1:tmax)){
    for(n in (1:N))
      x[n,]<-metropolis(x[n,],pdf[[n]],v.run)
  
#replica exhange
  if(t %% freq_ex==0 ){
    if(N>1){
    ind <- sample(seq(1,N-1),1)
    if(runif(1,0,1)<
       pdf[[ind+1]](x[ind,])*pdf[[ind]](x[ind+1,])/
       (pdf[[ind]](x[ind,])*pdf[[ind+1]](x[ind+1,]))){
      tmp <- x[ind+1,]
      x[ind+1,] <- x[ind,]
      x[ind,] <- tmp
    }
    }
    if(show){
      for(n in (1:N)){
        plot(x[n,1],x[n,2],xlim=plotrange,ylim=plotrange,col=cols[n])  
        par(new=T)
      }
      par(new=F)
    }
    xx<-rbind(xx,x[,1])
  }
}
  xx
}

xx<-rep_monte(tmax,N)
#plot(xx[,1])
#ts.plot(xx)
saveMovie({  rep_monte(tmax,N,T)}, movie.name = "test.gif", interval = 0.1, nmax = 30, ani.width = 600, ani.height = 600)