Practice Markov Chain in Two Ways: Excel & Python

Simulate Markov Chain in Excel & Python, using probability theories. Which way is easier and why?

Yung Codes

5 min read

·

Jul 5, 2020

--

1

Listen

Share

In this flash-card on Markov Chain, I will summarize a basic understanding of Markov Chain and show how to implement this method in the context of a simple problem using Excel and Python.

Background: What is Markov Chain?

• A stochastic (random) model that simulates the sequence of possible events where the next event is determined by the state of the previous event.
• Note: “Memoryless” property of Markov Chains means a process for which predictions of the future outcome is solely based on its present state. So the past and future is independent of each other. However, the future is conditional upon the present state.
• Discrete (finite state or space): Examples of discrete Markov Chain are random walks based on an integer.
• Continuous (over time): Examples of continuous Markov Chain are Wiener process (Brownian motion) and Poisson process.

Motivation: Why should I care?

• An important and fundamental concept that has many applications.
• It sounds fancy and has historical significance. For example, Einstein's Nobel Prize can be inferred to relate to this probability method (via Brownian motion, Wiener process). You also use it every day whenever you google anything.

Application: When to use Markov Chain?

• Famous applications: Google uses Markov Chain to rank webpages. Heavily used in forecasting such as certain weather pattern and the stock market trading algorithm (Hidden Markov Model (HMM): assumes that there is another process Y whose behavior depends on X, an unobservable state), and more. Chances are your life’s savings that you threw into investment funds are utilizing automated trading using Markov Chain principles.

• Built Environment applications: Transportation (traffic) forecasting, solar irradiance variability
• Building Energy: Used in occupancy simulation and predict energy usage — usually associate occupancy state (1, 2, 3, 4).

Practice: How can I apply the Markov Chain method in practice?

Simple Problem: Suppose the weather is divided into two states: “Sunny” and “Rain”. We will denote “R” as “Rain” and “S” as “Sunny”. We are given the probability matrix of the transition states:

Suppose we initialize the weather state on Monday or X_0 as “Rain”.

Question:

• What is the probability that it will rain on Wednesday, given that it is raining on Monday?
• In other words, what is Pr(X_2 = ‘R’ | X_0 = ‘R’)?

Practice: Excel

Logic used:

P(R|Xi-1)== If(Previous Day R/S = "R", 0.8 or R-> R, 0.4 or S -> R)
Ui == RAND()
Ui < P0 == If (Ui < P(R|Xi-1), “Yes”, “No”)
R/S == If (Ui<P0 == “Yes”, “R”, “S”) 
Count == COUNTIF(Range of R/S, “R”)

Press enter or click to view image in full size

• What I learned: The trickiest part about implementing this in Excel was the sheer amount of columns for the X_n desired event (X_2 in our case).
• Pro: Excel provides a good drawing board for simple tests and helps visually map out the logic. Good for beginners.
• Con: If n is large for X_n, Excel computation time may be slower and cumbersome.

Practice: Python

First, import numpy and input the 2 x 2 events.

import numpy as np 
import random as rm#define the states
states = np.array(["R", "S"])#Position sequence of events
transition_name = np.array([["RR", "RS"],["SR", "SS"]])#Probability matrix of transition states
P = np.array([[0.3, 0.7], [0.1, 0.9]])try:
  sum(P[0]) + sum(P[1]) == 2
  print("Probability Matrix input looks good to go!")
except:
  print("Error: Check your Probability Matrix Input!")

The try/except gatekeeper is kind of excessive, but it’s good to flex debug muscles and also points out that each ROW of the transition matrix needs to add up to 1. On the other hand, it doesn’t matter if the columns add up to 1.

Next, I built a function to evaluate the next state based on the current state and the state in question (what we care about).

def get_next_state(current_state, question_state): 
#Return next state  #Generate random variable
  U_i = rm.uniform(0,1)  #Get P(R|X_i-1)
  next_state = current_state + question_state
  P_index = np.where(transition_name == next_state)
  Pr_next_state = P[P_index[0], P_index[1]]  #Condition U_i < P(R|X_i-1)
  if U_i < Pr_next_state:
     next_state = question_state
  else:
     next_state = states[states != question_state].item(0)  return next_state

Next, it’s time to simulate this many times.

def get_prob_next_state(days, n_sample, initial_state, question_state):  assert initial_state == "R" or initial_state == "S" , "Need to     assert R or S for the initial state"  assert question_state == "R" or question_state == "S" , "Need to assert R or S for the question state"  state_list = []
  i = 0
  j = 0  for j in range(n_sample):
    current_state = initial_state
    for i in range(days):  
      next_state = get_next_state(current_state,question_state)
      current_state = next_state  
    state_list.append(current_state)prob = state_list.count(question_state)/(len(state_list) + 1)print("The probability of " , question_state , "on " , days , "day(s) since the initial day with condition " , initial_state , ", for " , n_sample , " samples, is: " , prob)

Now… let’s test it out with this:

get_prob_next_state(2, 1000, “R”, “R”)

The result:

The probability of  R on  2 day(s) since the initial day with condition  R , for  1000  samples, is:  0.17882117882117882

This shows a similar result from the Excel version. Here is the Google Colab notebook.

• What I learned: Numpy array can be faster than Python’s list, especially easier to maneuver than “list of list”.
• Pro: Scaleable. Flexible to answer more than the original question.
• Con: Take more time to develop compared to the Excel tool initially. Easier to make mistakes, thus need to conduct tests to make sure the code performs as intended.

Verify with Probabilities

Alternatively, you can also do a probability calculation in this problem. There are two possibilities:

• Mon (R) -> Tues (R) -> Wed (R) : (0.3) * (0.3) = 0.09
• Mon (R) -> Tues (S) -> Wed (R) : (0.7) * (0.1) = 0.07

Therefore, the total possibilities that it can rain on Wednesday is the sum (0.09+0.07 = 0.16). That is very close to 0.17 we got from our simulations.

How to Level Up?

• Try to solve more complicated examples (1, 2, 3).
• Recreate interesting papers using Markov Chain (1, 2, 3).

In this series, I’m going to write a flash-card version of interesting topics I learned while being a student at Georgia Tech OMSA. I find that writing these summaries to do “actively recall” (aka studying). Last update: July 2020.