Part 5 - Branching Processes and MCMC

Many real-world phenomena can be described as developing with a tree-like structure, for example,

• Growth of cells.

• Spread of viruses or other pathogens in a population.

• Nuclear chain reactions.

• Spread of funny cat videos on the internet.

• Spread of a surname over generations. 11This was the original motivation for studying branching processes by Francis Galton and Henry William Watson in the 19th century.

The process with which one node gives rise to “children” can be described as random. We will assume the probabilistic properties of this process is the same for all nodes.

Further, we will assume all nodes are organized into generatiomns, we are only concerned with the size of each generation.

Lastly, the important questions we are interested in are; How large are the generations? How much does the size vary? Will the process die out eventually (i.e., reach a generation with size 0)?

A (Galton-Watson) branching process is a discrete Markov chain $Z_0, Z_1, \ldots, Z_n, \ldots$ where,

• The state space is the non-negative integers.

• $Z_0 = 1$ (i.e., the process starts with one individual).

• $Z_n$ is the sum $X_1 + X_2 + \ldots + X_{Z_{n - 1}}$, where $X_i$ are independent random non-negative integers all with the same offspring distribution,

Connecting each of the $Z_n$ individuals in generation $n$ with their offspring in generation $n + 1$ we get a tree illustrating the branching process.

The offspring distribution is described by the probability vector $a = (a_0, a_1, \ldots)$ where $a_j = P(X_i = j)$.

We will mainly focus on the interesting cases where we assume $a_0 > 0$ (i.e., there is a positive probability of having no offspring) and $a_0 + a_1 < 1$ (i.e., there is a positive probability of having more than one offspring).

Then, we may compute that,

Thus, we directly get that,

Based on the value of $\mu$, we classify branching processes into three categories:

1. Subcritical branching processes: $\mu < 1$. The expected generation size decreases over time, $\lim_{n \to \infty} \mathbb{E}[Z_n] = 0$.

2. Critical branching processes: $\mu = 1$. The expected generation size remains constant over time, $\lim_{n \to \infty} \mathbb{E}[Z_n] = 1$.

3. Supercritical branching processes: $\mu > 1$. The expected generation size increases over time, $\lim_{n \to \infty} \mathbb{E}[Z_n] = \infty$.

Let $\mu$ be defined as above, and let $\sigma^2 = \mathrm{Var}(X_i)$ be the variance of the number of children.

Using the law of total variance we get,

From this, we can prove by induction that, for $n \geq 1$,

Let $X$ denote an offspring variable and $e$ the probability of extinction.

Using first-step analysis,

where we define the probability generating function (PGF) $G_X(s)$ as,

Let $e_n$ be the probability of extinction by generation $n$. Using first-step analysis again, we get,

Let $G_X$ be the probability generating function for the offspring distribution for a branching process. The probability of eventual extinction is the smallest positive root of the equation $G_X(s) = s$.

For any discrete random variable $X$ taking values in $\{0, 1, 2, \ldots\}$, define the probability function $G_X(s)$ as,

The series converges absolutely for $\Vert s \Vert \leq 1$. We assume $s$ is a real number in $[0, 1]$.

We get a 1-1 correspondence between probability vectors on $\{0, 1, 2, \ldots\}$ and functions represented by a series where the coefficients sum to 1.

Specifically, if $G_X(s) = G_Y(s)$ for all $s$ for random variables $X$ and $Y$, then $X$ and $Y$ have the same distribution.

Let $G_X$ be the probability generating function for the offspring distribution for a branching process. The probability of eventual extinction is the smallest positive root of the equation $G_X(s) = s$. Also, if the process is critical ($\mu = 1$) then the extinction probability is 1.

So far we have used Markov chains in the following way

1. Define a Markov chain

2. Find its limiting distribution, if it exists

Now, we want to,

1. Start with a desired distribution

Construct a Markov chain having this limiting distribution

In many contexts we want to compute the expecation $\mathbb{E}[r(X)]$ for some random variable $X$ and some function $r$.

One way to approximate this is to generate samples $x_1, \ldots, x_N$ from $X$ and use,

Thus, the idea is simulating from a Markov chain having $X$ as its limiting distribution generates an approximate sample. Average over it as above.

We know this works in classical statistics.

Strong law of large numbers of samples: If $Y_1, Y_2, \ldots, Y_m$ and $Y$ are i.i.d random variables from a distribution with finite mean, and if $\mathbb{E}[r(Y)]$ exists, then, with probability 1,

There is a corresponding result for Markov chains.

Strong law of large numbers for Markov chains: If $X_0, X_1, \ldots$ is an ergodic Markov chain with stationary distribution $\pi$, and if $\mathbb{E}_{X \sim \pi}[r(X)]$ exists, then, with probability 1,

where $X$ has distribution $\pi$.

Consider the Markov chain $X_0, X_1, \ldots$ with states $\{0, 1, 2\}$ and with,

One can find that the limiting distribution is $v = (\frac{20}{23}, \frac{2}{23}, \frac{1}{23})$.

Consider the function $r(x) = x^5$. If $X$ is a random variable with the limiting distribution,

If $Y_1, \ldots, Y_n$ are all i.i.d variables with the limiting distribution, we can check numerically that,

We also get for $X_0, X_1, \ldots$ that,

but in this case the limit is approached more slowly.

The Metropolis-Hastings algorithm is a method to construct a Markov chain having a desired limiting distribution. Let $\theta$ be a random variable with probability mass function $\pi( \theta)$. We also assume given a proposal distribution $T(\theta_{\text{new}} \mid \theta)$, which, for every given $\theta$, provides a PMF for a new $\theta_{\text{new}}$.

Finally, we define, for $\theta$ and $\theta_{\text{new}}$, the acceptance probability,

The Metropolis-Hastings algorithm is: Starting with some initial value $\theta_0$, generate $\theta_1, \theta_2, \ldots$ by, at each step, proposing a new $\theta$ based on the old using the proposal function and accepting it with probability $a$. If it is not accepted, the old value is used again.

If this defines an ergodic Markov chain, its unique stationary distribution is $\pi(\theta)$.

For any probability model over a vector $\theta = (\theta_1, \theta_2, \ldots, \theta_n)$, consider a Metropolis-Hastings proposal function changing only one coordinate, with the value of this coordinate simulated from the conditional distribution given the remaining coordinates. The acceptance probability is always 1, so the proposal is always accepted.

Putting together an algorithm updating different coordinates in different steps may create an ergodic Markov chain.

This is called Gibbs sampling. Sometimes the conditional distributions are easy to derive. Then this is an easy-to-use version of Metropolis-Hastings.

Uses a grid of vertices; We will assume an $n \times n$ grid. Two vertices $v$ and $w$ are neighbours, denoted $v \sim w$, if they are next to each other in the grid.

Each vertex $v$ can have value $+1$ or $-1$ (called its “spin”); We denote this by $\sigma_{v} = 1$ or $\sigma_{v} = -1$.

A configuration $\sigma$ consists of a choice of $+1$ or $-1$ for each vertex. Thus, the set $\Omega$ of possible configurations has $2^(n^2)$ elements. We define the energy of a configuration as $E(\sigma) = - \sum_{v \sim w} \sigma_v \sigma_w$.

The Gibbs distribution is the probability mass function on $\Omega$ defined by,

where $\beta$ is a parameter of the model; $\frac{1}{\beta}$ is called the temperature.

It turns out that when the temperature is high, samples from the model will show a chaotic pattern of spins, but when the temperature sinks below the phase transition value, samples will show chunks of neighbouring vertices with the same spin, i.e., the system will be “magnetized”.

For a vertex configuration $\sigma$ and a vertex $v$, let $\sigma_{-v}$ denote the part of $\sigma$ that does not involve $v$. Prpose a new configuration $\sigma^{\star}$ given an old configuration $\sigma$ by first choosing a vertex $v$, then, let $\sigma^{\star}$ be identical to $\sigma$ except possibly at $v$. Decide the spin at $v$ using the conditional distribution given $\sigma_{-v}$,

This works, but we will dicuss an even better approach “perfect sampling” later.

Given an ergodic Markov chain with finite state space of size $k$ and limiting distribution $\pi$. The idea is that, given $n$ prove that $X_n$ actually has reached the limit distribution $\pi$.

If we can rpove that the distribution at $X_n$ is independent of the starting state $X_0$, then we might be able to conclude that $X_n$ has the limiting distribution $\pi$.

Construct $k$ Markov chains that are dependent (“coupled”) but which are marginally Markov chains as above. If they start at the $k$ possible value at $X_0$ but have identical values at $X_n$, then we are done.

Consider the chains $X^{(j)}_{-n}, \ldots, X^{(j)}_0$ for $j = 1, 2, \ldots, k$. Instead of simulating $X^{(j)}_{i + 1}$ based on $X^{(j)}_i$ independently for each j, we define a function $g$ so that $X^{(j)}_{i + 1} = g(X^{(j)}_i, U_i)$ for all $j$ where $U_i \sim \mathrm{Uniform}(0, 1)$.

Thus, if two chains have identical values in $X_i$, they will also be identical in $X_{i + 1}$.

Thus, for a particular $n$, if all chains have not converged at $X_0$, we simulate $k$ chains from $X_{-2n}$ to $X_{-n}$. They might only hit a subset of the $k$ states at $X_{-n}$ and thus might coalesce to one state at $X_0$ using the old simulations. If not, double $n$ again.

One smart question is, do we need to keep track of all $k$ chains? We define a partial ordering on as et as a relation $x \leq y$ between some pairs $x$ and $y$ in the set, such that,

• If $x \leq y$ and $y \leq x$ then $x = y$.

• If $x \leq y$ and $y \leq z$ then $x \leq z$.

We will need that our partial ordering has a minimal element (an $m$ such that $m \leq x$ for all $x$) and a maximal element (an $M$ such that $x \leq M$ for all $x$).

If we have a partial ordering on the state space of the Markov chain and if $x \leq y$ implies $g(x, U) \leq g(y, U)$, then $g$ is monotone.

We can then prove that we only need to keep track of the chain starting at $m$ and the chain starting at $M$.

We saw above how to use Gibbs sampling to simulate from the Ising model. We then simulated from the conditional distributions, using,

Now, we extend this method into a perfect sampling.

Define the partial ordering on configurations $\sigma, \tau$ be defining,

We have minimal and maximal elements $m = (-1, \ldots, -1)$ and $M = (1, \ldots, 1)$.

Defining $g(\sigma, U)_v = 2 I (U \leq \pi(\sigma_v^{\star} = 1 \mid \sigma_{-v}) ) - 1$, makes $g$ monotone.

Given an Ising model with $\beta > 0$, start with choosing an integer $n > 0$ and setting configurations,