7 April 2017

Random processes

Random processes are sequences of random variables:

\(X = \{ X_t: t \in T \}\) where \(T\) is some index set.

\(X\) is a discrete time process when \(T\) is discrete (e.g \(T = \{0,1,2,3, \ldots\}\))

\(X\) is a continuous time process when \(T\) is continuous (e.g. \(T = [0, \infty ]\))

Random Walk

Random Walk

  • Also known as the drunkard's walk
  • \(X_i\) records the position of a drunk after \(i\) steps.
  • Each step is either to the left or the right with equal proability.
  • Let \(X_0 = 0\) and \(X_{i+1} = X_i + S_{i}\) where \[S_{i} = \begin{cases} +1 & \mbox{with probability } 1/2 \\ -1 & \mbox{with probability } 1/2. \end{cases} \]
  • Equivalently, \(X_i = X_0 + \sum_{j = 0}^{i-1} S_j\)
  • Simple enough that it can be analysed analytically, by also very easy to simulate.

Example: random walk

nStep = 50
step = sample(c(-1,1), size = nStep, replace = T) # sample step direction
walk = cumsum(c(0,step)) # sum steps
plotStep(nStep,walk) # plot

While it is typically easier to view these things graphically, we can write down the random process as \(X =\) {0, 1, 2, 3, 4, 3, 2, 1, 0, -1, 0, 1, 0, 1, 0, 1, 0, -1, 0, 1, 0, 1, 0, -1, 0, -1, -2, -1, 0, 1, 0, 1, 0, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 0, 1, 0, -1, 0, 1, 2,…}

Again

step = sample(c(-1,1), size = nStep, replace = T)
walk = cumsum(c(0,step))
plotStep(nStep,walk)

Longer walk

nStep = 500
step = sample(c(-1,1), size = nStep, replace = T)
walk = cumsum(c(0,step))
plotStep(nStep,walk)

Still longer walk

nStep = 10000
step = sample(c(-1,1), size = nStep, replace = T)
walk = cumsum(c(0,step))
plotStep(nStep,walk)

Add in a bias: Prob of +1 step now 0.51

nStep = 10000
step = sample(c(-1,1), prob = c(0.49,0.51), size = nStep, replace = T)
walk = cumsum(c(0,step))
plotStep(nStep,walk,"Biased random walk with 10000 steps")

The Poisson process

The Poisson process

Models events that occur independently of each other and randomly in time (or space).

E.g. Genetic mutations, arrival times of customers, radioactive decay, occurrences of earthquakes, industrial accidents, errors in manufacturing processes (e.g. silicon wafers or cloth)

  • We'll consider events that occur in time
  • Events occur at rate \(\lambda\) (or with intensity \(\lambda\))
  • A Poisson process is a continuous time counting process \(\{N(t), \, t\geq 0 \}\) where \(N(t)\) is the number of events in the interval \((0,t)\).

A Poisson process is defined by the following properties

  1. \(N(0) = 0\) and \(N(s) < N(t)\) when \(s<t\).
  2. The number of events in any non-overlapping intervals are independent.
  3. In any sufficiently small time interval the probability of 2 events occurring is 0.

Example of Poisson process

Properties of the Poisson process

Let process have rate \(\lambda\).

  1. Number of events in an interval of length \(t\) is Poisson distributed with parameter \(t\lambda\).

  2. Given the number of events in a Poisson process in some interval \([s,t]\), the times of the events are uniform in that interval

  3. The times between events — interarrival-times — in a Poisson process are exponentially distributed with parameter \(\lambda\).

Properties of the Poisson process

Can see distribution of inter-arrival times is exponential when we plot them:

On left, inter-arrival times are time between red lines.

On right, histogram of these times along with exponential density with rate 3.

Simulating the Poisson Process

Properties suggest two ways to simulate a Poisson process with rate \(\lambda\) on interval \([0,T]\):

Either simulate exponential variables with rate \(\lambda\) until cumulative time exceeds \(T\):

t = randExp(rate = lambda)
eventCount = 0
eventTime = []
while (t < T)
  eventCount = eventCount + 1
  eventTime = [eventTime, t]
  t = t + randExp(1, rate = lambda)
end
return(count = 1:eventCount, time = eventTime)

Or simulate a Poisson variable with rate \(\lambda*T\) to get number of events and then distribute them uniformly across interval \([0,T]\).

eventCount = randPoisson(rate = lambda * T)
eventTime = randUniform(min = 0, max = T, number = eventCount)
eventTime = sort(eventTime)
return(count = 1:eventCount, time = eventTime)

But: easiest way to simulate a Poisson random variable is simulate a Poisson process with rate \(\lambda\) for time 1 and count the number of events.

Splitting Poisson processes

Suppose each event is one of \(k\) possible types.

An event is type \(i\) with probability \(p_i\), for \(i \in \{1,\ldots,k\}\)

Now observe only events of type \(j\).

They form a Poisson process with rate \(\lambda p_j\).

Example: Consider a Poisson process with rate \(\lambda = 2\). Suppose events are either blue (with probability \(0.75\)) or red (with probability 0.25).

Looking only at the blue events gives a Poisson process with rate \(\lambda p_{blue} = 2 \times 0.75 = 1.5\).

While including only the red events is a Poisson process with rate \(2 \times 0.25 = 0.5\)

Merging Poisson processes

If \(X\) is a Poisson process with rate \(\lambda\) and

\(Y\) is a Poisson process with rate \(\mu\),

Then \(Z = \{Z(t) = X (t) + Y (t), t \geq 0\}\) is a Poisson process with rate \(\lambda + \mu\).

Markov chains

Markov property

The Markov property is about memorylessness: future states depend only on the current state

To propagate a process with the Markov property forward, we need only be told the current state to generate the next state.

Formally, the sequence of random variables \(X_1,X_2, X_3,\ldots\) is a Markov chain if it has the Markov property
\[P(X_{n+1} = x | X_n = x_n, X_{n-1} = x_{n-1},\ldots,X_2 = x_2,X_1=x_1) = P(X_{n+1} = x |X_n = x_n).\]

Named after Andrei Andreyevich Markov a Russian mathematician.

Markov chain example: a simple Weather model

Suppose the weather tomorrow depends only on the weather today.

Weather can either be Sunny or Rainy, S or R for short. State space is \(\{S, R\}\)

If it is sunny today, tomorrow it is sunny with probability 0.7 and rainy otherwise.

If it is rainy today, tomorrow it is rainy with probability 0.6 and sunny otherwise.

Markov chain transition probabilites

Often write \(\Pr(X_{n+1} = j|X_n = i)\) as \(P_{ij}\) for short.

Define the transition matrix \[P = [P_{ij}]\] where \(i\)th row and \(j\)th column is probability of moving from state \(j\) given that currently in state \(i\).

Each entry is a probability (so between 0 and 1) and all rows sum to 1: \(\sum_j P_{ij} = 1\)

Example: The transition matrix for the weather example given above is \[\left[ \begin{array}{cc} 0.7 & 0.3 \\ 0.4 & 0.6 \end{array} \right] \] where rows and columns 1 and 2 are indexed by \(S\) and \(R\)

Markov chain transition probabilites

The \(m\)-step transition probability is the probability of going from state \(i\) to state \(j\) in exactly \(m\) steps, \[P_{ij}(m) = \Pr(X_{n+m} = j|X_n = m).\]

The \(m-\)step transition matrix is \(P_m = [P_{ij}(m)]\).

Result (Chapman-Kolmogorov equations): \(P_{m+n} = P_mP_n\) (where the right-hand side is just standard matrix multiplication).

In particular, \(P_n = P^n\). That is, the n-step transition matrix is just the \(n\)th power of the transition matrix.

Example: The probability that it is raining the day after tomorrow given that it is sunny today is the \((S, R)\)-entry of the two-step transition matrix \(P_2 = P^2\):

Since \[P = \left[ \begin{array}{cc} 0.7 & 0.3 \\ 0.4 & 0.6 \end{array} \right] \]

\[P^2 = \left[ \begin{array}{cc} 0.61 & 0.39 \\ 0.52 & 0.48 \end{array} \right] \]

so the probability of rain two days after a sunny day is 0.39.

An example that is NOT Markov: Weather

  • We have rainy and sunny days as before in our simple model but hurricanes are also possible:
  • A hurricane lasts exactly 7 days and it rains hard every day during a hurricane. Call this state \(H\).
  • This process is not Markov since can only leave state \(H\) if been there for seven days.
  • So probabilty of next state depends on not only on current day but also on previous six days.

Markov chain example: music composition

The transition matrix between pitches is constructed from empirical counts of the observed transitions.

From which we can generate random realisations:

A,G,F,G,F,G,A,B,G,
F,G,F,G,A,B,D,E,
B,C,A,F,G,
B\(\flat\),A,F,G,A,G,A,B,G,A.