We will choose p_0 = 0.1, p_1 = 0.6, and p_2 = 0.3. This implies mu = 1.2.
p_0 <- 0.1
p_1 <- 0.6
p_2 <- 0.3
mu <- p_0 * 0 + p_1 * 1 + p_2 * 2
mu
## [1] 1.2
Since mu = 1.2, we should have 0 < a < 1. It turns out a = 1/3.
roots <- polyroot(c(p_0, p_1 - 1, p_2))
a <- Re(roots[1])
a
## [1] 0.3333
Running the simulation
p <- c(p_0, p_1, p_2) # probability of 0, 1, or 2 offspring
# graph sample chain -- run multiple times to observe extinction vs
# population explosion
nsteps <- 100
X <- matrix(0, nsteps, 1)
X[1] <- 1
for (n in 2:nsteps) {
ifelse(X[n - 1] > 0, X[n] <- sum(sample(0:(length(p) - 1), X[n - 1], replace = TRUE,
p)), X[n] <- 0)
}
plot(X, type = "l", xlab = "Generation", ylab = "Population", main = "Branching process")
We will choose the probabilities as
p_0 <- 0.4
p_1 <- 0.4
p_2 <- 0.2
p <- c(p_0, p_1, p_2)
sum((0:(length(p) - 1)) * p) # mu
## [1] 0.8
Calculating the extinction probability a
min(sapply(polyroot(c(p[1], p[2] - 1, p[3])), Re)) # a (extinction probability)
## [1] 1
Running the simulation
nsteps <- 100
X <- matrix(0, nsteps, 1)
X[1] <- 1
for (n in 2:nsteps) {
ifelse(X[n - 1] > 0, X[n] <- sum(sample(0:(length(p) - 1), X[n - 1], replace = TRUE,
p)), X[n] <- 0)
}
plot(X, type = "l", xlab = "Generation", ylab = "Population", main = "Branching process")
Calculating a numerically as the limit of the sequence a_n(1) as n -> Infty.
G <- function(s) p[1] + p[2] * s + p[3] * s^2
niteration <- 100
a <- matrix(0, niteration + 1, 1)
a[1] <- G(0)
for (n in 1:niteration) a[n + 1] <- G(a[n])
plot(1:(niteration + 1), a, xlab = "n", ylab = "a n(1)", ylim = c(0, 1))
Computing the partial sums of E(T), the expected time to extinction.
plot(1:(niteration + 1), 1 + cumsum((1 - a[1:(niteration + 1)])), xlab = "n",
ylab = "Partial sum")
E(T) appears to scale linearly with n and is therefore likely unbounded.
Filling in the probabilities
p_0 <- 0.15
p_1 <- 0.7
p_2 <- 0.15
p <- c(p_0, p_1, p_2)
sum((0:(length(p) - 1)) * p) # mu
## [1] 1
Computing extinction probability a
min(sapply(polyroot(c(p[1], p[2] - 1, p[3])), Re)) # a (extinction probability)
## [1] 1
Running simulation
nsteps <- 100
X <- matrix(0, nsteps, 1)
X[1] <- 1
for (n in 2:nsteps) {
ifelse(X[n - 1] > 0, X[n] <- sum(sample(0:(length(p) - 1), X[n - 1], replace = TRUE,
p)), X[n] <- 0)
}
plot(X, type = "l", xlab = "Generation", ylab = "Population", main = "Branching process")
It exhibits transient behavior, sometimes dying out initially like in the subcritical case while other times growing for a few generations and eventually dying out. The behavior is in between the subcritical and supercritical cases.
Computing partial sums of expected extinction time E(T)
G <- function(s) p[1] + p[2] * s + p[3] * s^2
niteration <- 1000
a <- matrix(0, niteration + 1, 1)
a[1] <- G(0)
for (n in 1:niteration) a[n + 1] <- G(a[n])
plot(1:(niteration + 1), 1 + cumsum((1 - a[1:(niteration + 1)])), xlab = "n",
ylab = "Partial sum")
The partial sums of the infinite series E(T) = sum_{n=1}^{\infty} (1-a_n(1)) also
are increasing in n but in this case take on a shape similar to logarithmic growth. It is appears that the partial sums still do not converge but this time grow at a much slower rate than previously.
Determining C ~ 1 - a_n(1)
plot(1:niteration, (1:niteration) * (1 - a[2:(niteration + 1)]), xlab = "n",
ylab = "n(1-a n(1))")
We know that the sequence n(1-a_n(1)) should approach C in the limit n -> Infty. From this plot, it appears that C = 6.583.