simple example of using a gillespie algorithm to simulate a ctmc with 2 events

simpleGillespie.R

events = c(
  birth = 1/1.5,
  death = 1/2.5
)

state = c(
  animals = 5
)

tmax = 1e3
t = 1
i = 2

while(t[length(t)] <= tmax){
  
  birth_rate = state[i-1] * events[["birth"]]
  death_rate = state[i-1] * events[["death"]]
  
  overall_rate = birth_rate + death_rate
  
  # sample when the thing happens
  t[i] = t[i-1] + rexp(n = 1,rate = overall_rate)
  
  # sample what happens
  what = sample(x = c(1,2),size = 1,prob = c(birth_rate,death_rate))
  
  if(what==1){
    state[i] = state[i-1] + 1
  }
  if(what==2){
    state[i] = state[i-1] - 1
  }
  
  if(state[i]==0){
    break()
  }
  
}

# plot(state,type="l")
plot(x = t,y = state,type="l")