The random walk that converges to exponential distribution

MCMC_exp.R

library(animation)

make_initial <- function(np, S){
  X = rep(S/np, np)
  return(X)
}

randomwalk <- function(X){
  np = length(X) 
  i = sample.int(np, 1)
  j = sample.int(np, 1)
  dx = runif(1, 0, 1)
  if(X[i] - dx > 0){
    X[i] = X[i] - dx
    X[j] = X[j] + dx
  }
  return(X)
}

run_rw <- function(n_time, np, S){
  x_hist <- array(0L, dim = c(n_time+1, np))
  x = make_initial(np, S)
  x_hist[1,] = x
  for(i in 1:n_time){
    x = randomwalk(x)    
    x_hist[i+1,] = x
  }
  return(x_hist)
}

set.seed(123456);res = run_rw(1000, 100, 100)
brk = seq(0, 10, 0.5)
yran = c(0, 2)
animation::saveGIF({
  for(i in 1:1001){
    hist(res[i,], breaks = brk, freq = FALSE, main=paste("step:", i-1), xlab="energy", ylim = yran)
    curve(dexp(x),　add=TRUE, col="royalblue", lwd=2, lty=2)    
  }}, interval=0.01)