######################################################
# Simulation of a birth-death process
#
# To use, you must supply three items
# 1. Enter an expression for the total rate of
# births and death
# It could be a function of N, t, lambda and mu.
# Replace [enter expression] with the expression
# for the total rate.
# 2. Enter an expression for the probability
# that an event is a birth.
# It could be a function of lambda and mu.
# Replace [enter expression] with the expression
# for the probability of a birth
# 3. If you wish to plot the deterministic solution
# as well, then set plot_deterministic=TRUE
# and enter the deterministic solution
# inside the function N_deterministic.
# This solution should be a function of time t,
# initial condition N0, lambda, and mu.
#
# You can also set the parameters
# N0: the initial condition
# t0: the initial time
# tmax: the approximate ending time
# n_samples: the number of samples
# Once you have filled in these expressions,
# you can
# you can source this file to simulate the
# birth-death process and plot the results.
######################################################
t0=0
N0=10
tmax=100
lambda=0.01
mu=0.02
n_samples=20
# if an event occurs
# determine the probability that the event is a birth
# (This doesn't depend on population size
# so can compute once at the beginning)
prob_birth = [enter expression]
# lists to keep track of history of times and birth counts
times=list();
Ns = list();
# keep track of maximum number so can set axis scale
max_N = 0
# number of times that the population went extinct
n_extinct = 0
plot_deterministic = TRUE
# solution to deterministic approximation to mutations
N_deterministic = function(t) {
# enter an expression here that is a function of time t
}
for(samp in 1:n_samples) {
t=t0
N=N0
t_vector = t
N_vector = N
while(t max_N) {
max_N = N
}
# if N is zero, the population went extinct.
# If this extinction occured before the time tmax
# record that an extinction occured.
# Either way, stop simulation
if(N==0) {
if(t <= tmax) {
n_extinct = n_extinct+1
}
break
}
}
# finished simulation for sample
# record results in lists
times[[samp]] = t_vector
Ns[[samp]] = N_vector
}
cat("number of populations that went extinct:", n_extinct)
# plot all n_sample lines on the same graph
plot(times[[1]],Ns[[1]],'l',ylim=c(0,max_N), xlim=c(0,tmax),ylab="count", xlab="time", col=1)
if(n_samples>1) {
for(samp in 2:n_samples) {
lines(times[[samp]],Ns[[samp]],'l',col=samp)
}
}
if(plot_deterministic) {
curve(N_deterministic, from=0, to=tmax, add=TRUE, col= "black", lwd=3)
}