---
title: "Branching Processes"
author: "Math 365 Tanya Leise"
date: "February 22, 2016"
output: html_document
---

Email the pdf knit from your R markdown to tleise@amherst.edu by Friday 2/26.

The objective of this lab is to examine some examples of Galton-Watson branching processes more closely to determine their long-term behavior and properties. In the 1870s Francis Galton posed a question about the probability of aristocratic surnames becoming extinct, and a solution was proposed by the Reverend Henry William Watson. Together they published a paper introducing a branching process to model patrilineal propagation of aristocratic family names. The same process was reinvented by Leo Szilard (who spurred the creation of the Manhattan Project) in the 1930s to model the proliferation of free neutrons in a nuclear fission reaction.

Assume $X_n$ is a branching process with $X_0=1$, where $p_k$ is the probability of an individual having $k$ offspring.

### Supercritical case $\mu>1$

Choose positive values for $p_0$, $p_1$, and $p_2$ satisfying $p_0+p_1+p_2=1$ (so $p_n=0$ for $n\ge3$) such that the expected number of offspring of an individual is greater than 1.

```{r}
p <- c(.5,.3,.2) # probability of 0, 1, or 2 offspring (change these values to your choice)
mu <- sum((0:(length(p)-1))*p) # expected number of offspring
```

(1) State your values of $p_0$, $p_1$, and $p_2$ and the resulting value of $\mu$.

(2) Calculate the value of $a$ (extinction probability).

(3) Run several simulations using the code below and your values of $p_0$, $p_1$, and $p_2$ entered above. You should observe ``boom or bust" behavior (rapid extinction or population explosion).

```{r}
# graph sample chain -- run multiple times to observe
# extinction vs population explosion
nsteps<-20
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")
```

### Subcritical case $\mu<1$

Choose positive values for $p_0$, $p_1$, and $p_2$ satisfying $p_0+p_1+p_2=1$ (so $p_n=0$ for $n\ge3$) such that the expected number of offspring of an individual is less than 1. 

```{r}
p <- c(.5,.3,.2) # probability of 0, 1, or 2 offspring (change these values to your choice)
mu <- sum((0:(length(p)-1))*p) # expected number of offspring
```

(1) State your values of $p_0$, $p_1$, and $p_2$ and the resulting value of $\mu$.

(2) Calculate the value of $a$ (extinction probability).

(3) Run several simulations using these values of $p_0$, $p_1$, and $p_2$. You should observe eventual extinction with probability 1.

```{r}
nsteps<-20
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")
```

(4) The expected time to extinction $T$ can be calculated via
$\mathbb{E}(T)=1+\sum_{n=1}^\infty{(1-a_n(1))}$, where we can recursively calculate the probability of extinction at time $n$ using $a_{n}(1)=G(a_{n-1}(1))$ with $a_1(1)=G(0)$. Here $G(s)=\sum_{k=0}^\infty{p_ks^k}$ is the generating function and $a_n(1)=G_n(0)$ is the probability of extinction by step $n$. We can calculate these values using a for-loop, then examine a plot to check convergence:

```{r}
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),type="l")
```

You should observe the $a_n(1)$ converging to $a$, which is a fixed point of the iteration. 

Next compute partial sums of $\mathbb{E}(T)=1+\sum_{n=1}^\infty{(1-a_n(1))}$ and plot them to check convergence of the infinite series:

```{r}
plot(1:(niteration+1),1+cumsum((1-a[1:(niteration+1)])),
     xlab="n", ylab="Partial sum",type="l")
```

What do you estimate $\mathbb{E}(T)$ to be, based on your graph?

### Critical case $\mu=1$

Choose positive values for $p_0$, $p_1$, and $p_2$ satisfying $p_0+p_1+p_2=1$ (so $p_n=0$ for $n\ge3$) such that the expected number of offspring of an individual equals 1. 

```{r}
p <- c(.5,.3,.2) # probability of 0, 1, or 2 offspring (change these values to your choice)
mu <- sum((0:(length(p)-1))*p) # expected number of offspring
```

(1) State your values of $p_0$, $p_1$, and $p_2$ and the resulting value of $\mu$.

(2) Calculate the value of $a$ (extinction probability).

(3) Run several simulations using these values of $p_0$, $p_1$, and $p_2$. What do you observe for this borderline case?

```{r}
nsteps<-20
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")
```

(4) Plot the partial sums of the infinite series $\mathbb{E}(T)=1+\sum_{n=1}^\infty{(1-a_n(1))}$ for this case. What do you observe?

(5) The expected extinction time is infinite for this critical case, even though the extinction probability $a$ is 1. The reason is that $1-a_n(1) \sim C/n$ for large $n$. Verify that this occurs for your example by estimating the value of $C$ from a plot of $n(1-a_n(1))$:

```{r}
plot(1:niteration,(1:niteration)*(1-a[2:(niteration+1)]), xlab="n", ylab="n(1-a_n(1))",type="l")
```

The theoretically derived value is $C=2/\sigma^2$, where $\sigma^2$ is the variance $\sum_{k=0}^\infty{p_k(k-\mu)^2}$.
