simulate Boltzmann distribution

OsawaRyuTedukuri.R

#参考文献：永井佑紀『1週間で学べる！Julia数値計算プログラミング』（講談社）
library(animation)
library(reshape2)

make_initial <- function(n_people,n_chip){
  basket = integer(n_people)
  for(i in 1:n_chip){
    target = sample.int(n_people, 1)
    basket[target] = basket[target] + 1
  }
  return(basket)  
}

giveandtake <- function(basket){
  n_people = length(basket) 
  give = sample.int(n_people, 1)
  take = sample.int(n_people, 1)
  while (give==take) {
    give = sample.int(n_people, 1)
    take = sample.int(n_people, 1)
  }
  if(basket[give] > 0){
    basket[give] = basket[give]-1
    basket[take] = basket[take]+1
  }
  return(basket)
}

run_game <- function(n_time, n_people, n_chip, n_group){
  box_hist <- array(0L, dim = c(n_time+1, n_people, n_group))
  box = replicate(n_group, make_initial(n_people, n_chip))
  box_hist[1,,] = box
  for(i in 1:n_time){
    box = apply(box, 2, giveandtake)    
    box_hist[i+1,,] = box
  }
  return(box_hist)
} 
set.seed(1234);res = run_game(1000, 6, 100, 100)
brks = seq(0, 100, by=5)
y_ran = c(0, 400)

animation::saveGIF(
  for(i in 1:1001){
    hist(res[i,,], breaks = brks, ylim = y_ran, main=paste("time:", i-1), xlab="chip")
  }, interval=0.01)