11. Stochastic processes#

“Stochastic, or random, processes is the dynamic side of probability. What differential equations is to calculus, stochastic processes is to probability.”

—Robert P. Dobrow, Introduction to Stochastic Processes with R

11.1. Notation#

English

General notation

Acceptance probability

\(A\)

Limiting distribution

Markov chain

Proposal distribution

\(S\)

Random walk

Reversibility

Stationary distribution (or equilibrium distribution)

Stochastic process

Transition density

\(T\)

11.2. Introduction#

A stochastic (random) process is simply one in which outcomes are uncertain. By contrast, a deterministic process always produces the same result for a given input. While functions and differential equations are often used to describe deterministic processes, we can use random variables and probability distribution functions to describe stochastic ones.

Stochastic processes are very useful for the modeling of systems that have inherent randomness. This is certainly true in physics, and not the least for quantum systems. We will not delve deep into this topic in these notes. You might encounter stochastic modeling in the future when studying, for example, statistical physics.

Instead, our specific goal is to become familiar with a special kind of stochastic processes known as Markov chains. Having defined and studied general properties, we will find that such chains be designed to produce samples from a stationary probability distribution of our choice.

Discuss

Why is it so useful in Bayesian inference to be able to produce samples from a user-specified probability distrbution?

Example 11.1 (An exponential growth model)

Consider a simple exponential growth model to describe a population of bacteria in a culture with a food source. In a deterministic model we assert that the population grows at a fixed rate \(\lambda\). The growth rate is then proportional to the population size as described by the differential equation

(11.1)#\[\begin{equation} \frac{dy}{dt} = \lambda y(t). \end{equation}\]

With the initial population \(y(0) = y_0\), this equation is solved by the exponential function

(11.2)#\[\begin{equation} y(t) = y_0 \exp(\lambda t), \end{equation}\]

from which the population is precisely determined. For example, with \(y_0 = 10\) and \(\lambda = 0.25\) per minute we find that the population has almost doubled four times after 11 minutes (since \(y(11) = 15.6 \approx 16\)).

This model, however, does not take into account the randomness of the reproduction of individual organisms. This can be considered in a stochastic growth model where we assign the probability \(\lambda \Delta t\) for a single bacteria to split once during a (very short) time interval \(\Delta t\). This random process can be simulated using a short python script (click below to show the code). Ten different realizations of the stochastic model are visualized in Fig. 11.1 and compared with the deterministic one. Questions that might be of scientific interest include:

  • What is the expected number of bacteria at time \(t\)?

  • What is the probability that the number of bacteria will exceed some threshold value after \(t\) minutes?

  • What is the distribution of the time it takes for the population to double in size?

Hide code cell content 
import numpy as np
import matplotlib.pyplot as plt
from myst_nb import glue

# Set the relevant time interval and the model parameters
t = np.linspace(0,16,num=20000)
lam=0.25
y0=10

# Deterministic model
def y_deterministic(t,y0=1,lam=0.2):
    return y0*np.exp(lam*t)
    
# Stochastic model
def y_stochastic(t,y0=1,lam=0.2):
    y = [y0]
    assert max(t[1:]-t[:-1])/lam < 0.01, "The time interval must be much smaller than lambda."
    t_old = t[0]
    for t_i in t[1:]:
        dt = t_i-t_old
        t_old=t_i
        y_i = y[-1]
        # Increase the number of bacteria by the number of splits
        y_i += np.sum(np.random.random(y_i)<lam*dt)
        y.append(y_i)
    return np.array(y)
    
# Create a visualization
fig, ax = plt.subplots(nrows=1, ncols=1, sharey=True, **{"figsize":(6,4)})
ax.plot(t,y_deterministic(t,lam=lam,y0=y0),label='deterministic model',lw=2,color='red')
y_max = 1.1 * y_deterministic(t[-1],lam=lam,y0=y0)
np.random.seed(seed=3)
num_runs=20
y_runs = np.zeros((num_runs,len(t)))
for irun in range(num_runs):
    y_run = y_stochastic(t,lam=lam,y0=y0)
    if irun==0: label='stochastic growth processes'
    else: label=''
    if irun%2==0: ax.plot(t,y_run,lw=1,color='k',label=label)
    y_runs[irun,:] = y_run
y_mean = y_runs.mean(axis=0)
ax.plot(t,y_mean,lw=1,ls='--',color='k',label=f'mean of {num_runs} stochastic runs')
ax.set_ylim((0,y_max))
ax.legend(loc='best')
ax.set_xlabel('Time (min)')
ax.set_ylabel('Population')
ax.set_title(rf'Bacterial growth ($\lambda={lam}$, $y_0={y0}$)')

glue("bacterial-growth_fig", fig, display=False)

# How large is the probability that the number of bacteria exceeds 100 after 10 mins?
i10 = np.argmin(np.abs(t-10))
fraction_100_after_10min = 100 * np.sum(y_runs[:,i10]>100) / num_runs
glue("growth_question", fraction_100_after_10min)

65.0
../../_images/293d21597f2947366eefeebf3ef3d74ad422b173f83ee8c84e6f7d8ef7d257c5.png
../../_images/293d21597f2947366eefeebf3ef3d74ad422b173f83ee8c84e6f7d8ef7d257c5.png

Fig. 11.1 Growth of a bacteria population. The deterministic exponential growth curve is plotted against a number of realizations of a stochastic growth process. We find for example that 65% of the processes give 100 bacteria after ten minutes.#

11.3. Definition of a stochastic process#

Using probability theory, a stochastic process \(X\) is a family

(11.3)#\[\begin{equation} \{X_t \; : \; t \in T\}, \end{equation}\]

of random variables \(X_t\) indexed by some set \(T\). All the random variables take values in a common state space, \(S\). The state space can be discrete or continuous. Although we might use discrete state spaces in some of the examples, one should realize that most examples in physics belong to the continuous type.

Furthermore, the family is known as discrete process if the indices take discrete values \(t \in \{0,1,2,\ldots\}\). This will in fact be the situation in the practical application of Markov chain Monte Carlo applications that we are aiming for. On the other hand, many physical stochastic processes are continuous and will have \(T \in \mathbb{R}\) or \(t \in [0,\infty)\). In either case we think of a stochastic process as a family of variables that evolve as time passes, although in general \(T\) might be any kind of ordering index and is not restricted to physical time.

Due to the randomness, different runs of the stochastic process will generate different sequences of outputs. An observer recording many such sequences generated from the same stochastic process would need to measure all the following conditional probabilities to fully describe the process:

(11.4)#\[\begin{equation} \begin{gathered} p_{X_0}(x_0), \\ p_{X_1 \vert X_0}(x_1 | x_0), \\ p_{X_2 \vert X_0, X_1}(x_2 | x_0, x_1), \\ \vdots\\ p_{X_n \vert X_0, X_1, \ldots, X_{n-1}} \left( x_n | x_0, x_1, \ldots, x_{n-1} \right). \end{gathered} \end{equation}\]

Example 11.2 (Conditional probabilities of a stochastic process)

Below is an example of a stochastic process (implemented using an abstract python class StochasticProcess).

Note that the process indeed is defined by:

  1. the distribution of the first random variable \(X_0\), which is included in the definition of class method start, and

  2. the conditional distributions for subsequent variables (as defined in the class method update).

Two relevant visualizations are also included in the example below

  • Plots of the first variables from a few runs.

  • A corner plot of bivariate and univariate marginal distributions.

import sys
import os
# Adding ../Utils/ to the python module search path
sys.path.insert(0, os.path.abspath('../Utils/'))

from StochasticProcess.StochasticProcess import StochasticProcess as SP

class StochasticProcessExample(SP):
    def start(self,random_state):
        return random_state.uniform()

    def update(self, random_state, history):
        return np.sum(history) + random_state.uniform()

# Create an instance of the class. It is possible to provide a seed
# for the random state such that results are reproducible.
example = StochasticProcessExample(seed=1)

# The following method call produces 4 runs, each of length 10.
example.create_multiple_processes(10,4)

# and here plotted
fig_runs, ax_runs = example.plot_processes()
../../_images/bfb20b1861937e1c178ba46116ac2083aceca4d8121695c761e7f05cfe82393b.png 
# Pretty corner plots generated with `prettyplease`
# https://github.com/svisak/prettyplease 
import prettyplease.prettyplease as pp

# Here we instead create 10,000 runs (to collect statistics), each of length 4.
example.create_multiple_processes(4,10000)
fig_corner = pp.corner(example.sequence.T,labels=[fr'$X_{{{i}}}$' for i in range(13)])
../../_images/146ee5c6bf705b905b3af54357153d35df95012a52d80d5b2f81689dd3104ebd.png

Exercise 11.1

Can you understand the initialization and update rules of Example 11.2?

Exercise 11.2

For each panel of the pairplot above:

  • What conditional distribution of random variables does each panel show?

  • What random variables are marginalized out, if any?

  • Explain how the pairplot answers the question “Are \(X_3\) and \(X_0\) independent?”?

Hints:

  • The top-left panel shows \(\p{X_0}\), with no marginalization since \(X_0\) starts the random sequence.

  • The independence of \(X_3\) and \(X_0\) would imply \(\pdf{X_3}{X_0}=\p{X_3}\).

Exercise 11.3

Define your own stochastic process by creating your own initial and update methods within the StochasticProcess class. Explore the conditional probabilities by creating many chains and making plots.

11.4. Example: Random walk#

One of the simplest random processes is the so called random walk. This process can arise in many ways. A standard example is the coin-flip casino game in which a gambler G starts with a fortune \(a\). The croupier tosses a (possibly biased) coin. Each time heads appears, the gambler gets one euro (let us assume that the game takes place in Monte Carlo). If tails appears, the player loses one euro.

The player’s fortune \(S_n\) after \(n\) tosses is

(11.5)#\[\begin{equation} S_n = a + \sum_{i=1}^n X_i, \end{equation}\]

where \(X_i\) is a random variable with possible outcomes \(\{-1,1\}\). The probability mass function is \(\{ p(X_i = -1) = p, p(X_i = 1) = 1-p \}\), where an unbiased coin would have \(p=0.5\). The fortune \(S_n\) can be considered a random walk stochastic process. A number of different scenarios can be considered by imposing external constraints on the gambler’s fortune (such as the termination of the game when \(S_n = 0\)).

Example 11.3 (Simple random walk)

In the following code example (again using the StochasticProcess class) we create a large number of such processes starting at \(X_0 = 0\). Here we investigate the questions:

  • What is the mean position of the random walk after \(n\) steps?

  • What is the average distance traveled after \(n\) steps?

Can you formulate these questions as some moments of a probability distribution?

class randomwalk(SP):
    def start(self,random_state):
        return 0.

    def update(self, random_state, history):
        step = np.sign(random_state.uniform()-0.5)
        return history[-1]+step

test=randomwalk(seed=1)
test.create_multiple_processes(20,5)
test.plot_processes()

test.create_multiple_processes(101,5000)

for step in [9,16,25,36,49,64,81,100]:
    mean = test.sequence[step,:].mean()
    std = test.sequence[step,:].std()
    print(f'After {step:>3} steps: mean = {mean:5.3f}, std = {std:6.2f}')

After   9 steps: mean = -0.021, std =   3.00
After  16 steps: mean = -0.004, std =   4.02
After  25 steps: mean = 0.017, std =   4.99
After  36 steps: mean = -0.043, std =   6.01
After  49 steps: mean = -0.015, std =   7.00
After  64 steps: mean = 0.047, std =   8.01
After  81 steps: mean = 0.112, std =   9.19
After 100 steps: mean = 0.013, std =  10.27
../../_images/8f5a3e59af2df730318df1fbc070fd0d162953834465608135650eef5fce2e7d.png

These walks can be made more general by allowing the steps \(X_n\) to have some continuous distribution on the reals, and it can be placed in a multi-dimensional space with both positions and steps in \(\mathbb{R}^p\). A particularly well-known random walk process involves the motion of a particle in a fluid where it is subjected to multiple collisions. This gives a so called Brownian motion.

11.5. Special case: Gaussian process#

The joint density function for a multivariate normal distribution was presented in Eq. (6.33) for \(\boldsymbol{x} = (x_1, x_2, \ldots, x_k)\) corresponding to random variables \(X_1, X_2, \ldots X_k\). This normal distribution is completely determined by its mean vector, \(\boldsymbol{\mu}\), and covariance matrix, \(\boldsymbol{\Sigma}\), where elements \(\mu_i = \expect{X_i}\) and \(\Sigma_{ij} = \cov{X_i}{X_j}\).

The multivariate normal distribution has the remarkable property that all marginal and conditional distributions are normal, and specified by the corresponding subsets of the mean vector and covariance matrix.

A Gaussian process extends the multivariate normal distribution to a stochastic process with a continuous time index \(T\).

Definition 11.1 (Gaussian process)

A Gaussian process \((X(t))_{t \geq 0}\) is a stochastic process with a continuous-time index \(t \in [0,\infty)\) if each finite-dimensional vector of random variables

\[ X(t_1), X(t_2), \ldots, X(t_k) \]

has a multivariate distribution for \(0 \leq t_1 < \ldots < t_k\).

A Gaussian process is completely determined by its mean function \(\mu(X(t))\) and covariance function \(C(X(s), X(t))\) for \(s, t \geq 0\).

A Gaussian process is stationary if the mean function \(\mu(X(t))\) is constant for all \(t\) and if the covariance function fulfils \(C(X(s), X(t)) = C(X(s+h), X(t+h))\) with \(h \geq 0\). Note that stationarity is not a requirement for a Gaussian process.

11.6. Solutions#

Here are answers and solutions to selected exercises.