Geometric and binomial distribution -> exponential and Poisson distribution

sim_exponential.R

library(animation)

sim_pois_1d <- function(t_n, g_n, p, seed){
  set.seed(seed)
  counts <- integer(t_n)
  x <- seq(0, 1, length.out=g_n)
  delta_list <- vector("list", t_n)
  hit <- matrix(0, t_n, g_n)
  pud = 10
  for(i in 1:t_n){
    y <- rbinom(g_n,1,p)
    hit[i,] <- y
    delta = diff(x[y==1L]) # geometric
    delta_list[[i]] <- delta
    counts[i] <- sum(y)
  }
  list(x=x,
       hit=hit,
       delta=delta_list,
       counts=counts)
}

res = sim_pois_1d(t_n=100, g_n=100, p=0.1, 1234)

tab = table(res$counts)
ran_x = range(res$counts)
prob_theo = dpois(ran_x[1]:ran_x[2],10)

ran_y = c(0, 0.3)
h2 = hist(unlist(res$delta), breaks = "scott", freq=FALSE)

saveGIF({
  layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
  for(i in 1:100){
    plot(res$x, res$hit[i,], type = "h", xlab = "x", ylab = "", main=paste("time:", i))
    
    tab1 = table(res$counts[1:i])
    plot(tab1/sum(tab1), xlim=ran_x, ylim=ran_y, 
         xlab="count", ylab="prob", type = "h", main = "poisson dist.")
    lines(ran_x[1]:ran_x[2], prob_theo, type = "b")
    
    hist(unlist(res$delta[1:i]), h2$breaks, freq=FALSE,
         main = "exponential dist.", xlab = "diff")
    curve(dexp(x,10), add=TRUE)  
  }
},movie.name = "sim_pois.gif", interval=0.2)