A <- 5 # player A's starting amount
B <- 10
jua <- function(){
n <- 0
while( A != 0 & B != 0){
X = rbinom(1, 1, 0.5)
if(X == 1) {
A = A+1
B = B-1
}
else{
A = A-1
B = B+1
}
n= n+1
}
if(B == 0) return(1)
else return(c(0, n))
}Markov Chains
An introduction
Warning: This is not a primer, only a computational exploration. These are not the most efficient implementations; the markovchain package is there for that. This is where I try to implement things I learned in the course MTH212 from scratch using matrix and element wise multiplication, and some logical computation.
A Markov chain is a sequence of Random Variables \(X_0, X_1 … X_n\) where the probability of \(X_k\) taking any value depends only on the value of \(X_{k-1}\).
We deal with finite state spaces i.e. the number of values our random variables of interest can take is finite.
Each element \(P_{i,j}\) of a transition probability matrix \(P\) represents the probability of going from state i to j in one step (i.e. \(P(X_{k+1} = j | X_k = i)\) ).
The Monte Carlo Way
We can estimate the probability of a Markov Chain ending in a particular state by simulating it and then calculating the proportion of times it ends up at that state.
For simulating one simulation of a zero sum game:
Now, to get probability of A winning, and expected length of game:
mtx <- jua()
for(i in 2:1e3){
mtx <- cbind(mtx, jua())
}
print(mean(mtx[1,])) # average of outcomes (1 means A wins)[1] 0.32print(mean(mtx[2,])) # average of game length[1] 28.354Stochastic matrix
A stochastic matrix is by one that can be a transition probability matrix of a markov chain. Since one row represents the probabilities of going to all the n states in the next step from a fixed current one, they should add up to 1. Also, a good stochastic matrix should be square and a 2 D array.
mock <- matrix(c(1, 0, 0, 0,
.5, 0, .5, 0,
0, .5, 0, .5,
0, 0, 0, 1 ), 4,4,
byrow = T)
print(mock) [,1] [,2] [,3] [,4]
[1,] 1.0 0.0 0.0 0.0
[2,] 0.5 0.0 0.5 0.0
[3,] 0.0 0.5 0.0 0.5
[4,] 0.0 0.0 0.0 1.0is_stochastic <- function(P){
dp <- dim(P)
if(length(dp) ==2) if(dp[1] == dp[2]) if( setequal(rowSums(P),numeric(length = dp[1])+1)) return(TRUE)
else return(FALSE)
}
print(is_stochastic(mock))[1] TRUESimulating a Markov chain given the TPM and initial state
simulate_MC <- function(P, init){
n <- dim(P)[1]
count <- 0
current <- init
while(P[current, current] != 1) { # stop if we're in an absorbing state
count <- count +1
current <- sample(1:n, 1, prob = P[current,])
}
return(c(current, count))
}
simulate_MC(mock, 2)[1] 1 1Absorbing states
The easiest ones to catch. If a state can only lead to itself, it’s absorbing:
absorbing_states <- function(P, mode = "binarr"){
n <- dim(P)[1]
ans <- ((P * diag(1, n, n)) %*% matrix(1, n, 1) == 1)
if(mode == "binarr" ) return(ans)
else if (mode == "indices") return(which(ans))
}
print(absorbing_states(mock, "indices"))[1] 1 4Reachable states
Given a state i, if there is a finite probability of starting from it and going to some state j, j is called reachable from i. By definition, a state is reachable from itself.
The probabilities \(P(X_{k+2} = j | X_k = i)\) can be represented by the elements of the matrix \(P^2\) i.e. the matrix product of P with itself. Similarly, \(P(X_{k+m} = j | X_k = i)\) are arranged in the matrix \(P^m\).
We can prove (and intuitively reason) that in a finite state space Markov chain with n elements, if an element is reachable at all, it must be reachable in at most n-1 steps.
We can therefore reason that if the (i,j)th element of \(\sum_{i=1}^n P^i\) is non zero, j must be accessible or reachable from i.
access <- function(P){
if(is_stochastic(P))
{
n <- dim(P)[1]
access_mtx <- matrix(0,n,n)
temp_mtx <- (P + diag(1,n,n))
# a state that is accessible can be reached in at most n - 1 steps
# check for states accessible in 1st, 2nd...(n-1)th steps
for(i in 1:n-1){
access_mtx <- (access_mtx + temp_mtx)
temp_mtx <- temp_mtx %*% P
}
return((access_mtx>0) +0)
}
else return("Invalid stochastic matrix")
}
print(access(mock)) [,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 1 1 1 1
[3,] 1 1 1 1
[4,] 0 0 0 1The access function returns a binary matrix \(Q\) where \(Q_{i,j}\) is the Boolean truth value of whether j is accessible from i.
Communication Classes
A Markov chain’s state space is inevitably split into multiple equivalence classes such that all constituent states communicate (i.e. are mutually accessible).
To find out which ones are, we take the Boolean accessibility matrix \(Q\) (where \(Q_{i,j}=1\) means \(i \rightarrow j\)) and multiply it elementwise with its transpose (\(Q^T_{i,j}=1\) implies \(Q_{j,i}=1\) and thus \(j \rightarrow i\)). We get the Boolean matrix for which classes communicates with others.
The rows for all members of an equivalent class will be identical (they all communicate with the same states, ie each other). Thus, taking the unique rows, we get each distinct class as a row, with 1s denoting the member states.
comm_classes <- function(P, mode = "bin_arr"){
n <- dim(P)[1]
step_mtx <- ((P + diag(1,n,n))>0) + 0 # classes accessible in first step
access_mtx <- access(P)
commun_mtx <- access_mtx * t(access_mtx) # states which are accessible from each other communicate
eq_class <- unique(commun_mtx) # states belonging to same equivalent class communicate with same set of states
# the rest is dressing up the output
foo <- character(0)
for(i in 1:dim(eq_class)[1]) foo <- c(foo, paste("class", i, sep= ""))
rownames(eq_class) <- foo
if(mode == "bin_arr") return(eq_class) # jth row means jth equivalent class
df <- as.data.frame(matrix(0, dim(eq_class)[1],1))
chars <- character(dim(eq_class)[1])
for(i in 1:dim(eq_class)[1]){
nom <- character(0)
for(j in 1:n){
if(eq_class[i,j] == 1) nom <- paste(nom, j)
}
chars[i] <- nom
}
rownames(df) <- foo
df[,1] <- chars
colnames(df) <- "elements"
if(mode == "visual") return(df)
}
comm_classes(mock) [,1] [,2] [,3] [,4]
class1 1 0 0 0
class2 0 1 1 0
class3 0 0 0 1comm_classes(mock, "visual") elements
class1 1
class2 2 3
class3 4Classifying States as Recurrent or Transient
For a given state \(i\), if all the states that are accessible from \(i\) also communicate with it (i.e. i is also accessible from all those states), that implies we will eventually get back to \(i\) with probability 1 if we start at \(i\). If it’s not recurrent, it’s called transient.
That means that if we subtract the Boolean communication matrix from the accessibility matrix and get all zeros in a row, that state is recurrent. (which is equivalent to saying that the rowsum is zero since it will always be non-negative- no communication without accessibility is possible)
classify_states <- function(P, mode = "bin_arr"){
n <- dim(P)[1]
access_mtx <- access(P)
commun_mtx <- access_mtx * t(access_mtx)
good_mtx <- access_mtx - commun_mtx
rt_state <- matrix((rowSums(good_mtx) ==0)+0,n,1)
if(mode == "bin_arr") return(rt_state)
else if(mode == "visual"){
vis <- as.data.frame(matrix(as.character(factor(rt_state, levels = c(0,1), labels = c("transient", "recurrent"))), n,1))
colnames(vis) <- "type"
return(vis)
}
}
classify_states(mock, "visual") type
1 recurrent
2 transient
3 transient
4 recurrentClassifying Classes as Recurrent/Transient
Relatively easy. All states in a communication class are either recurrent or transient. A class is recurrent if all its states are recurrent and vice versa.
classify_classes <- function(P, mode = "bin_arr"){
n <- dim(P)[1]
eq_class <- comm_classes(P)
rt_state <- classify_states(P)
ans <- (eq_class %*% rt_state > 0)+0
y<- factor(ans, levels = c(0,1), labels = c("transient", "recurrent"))
y <- as.character(y)
df <- as.data.frame(matrix(0, dim(eq_class)[1],2))
chars <- character(dim(eq_class)[1])
for(i in 1:dim(eq_class)[1]){
nom <- character(0)
for(j in 1:n){
if(eq_class[i,j] == 1) nom <- paste(nom, j)
}
chars[i] <- nom
}
df[,1] <- chars
df[,2] <- y
colnames(df) <- c("elements", "type")
foo <- character(0)
for(i in 1:dim(eq_class)[1]) foo <- c(foo, paste("class", i, sep= ""))
rownames(df) <- foo
if(mode == "bin_arr") return(ans)
else if(mode == "visual") return(df)
}
classify_classes(mock, "visual") elements type
class1 1 recurrent
class2 2 3 transient
class3 4 recurrentTo do
Stable configurations i.e. probability of ending up in state i after infinite time
Infinite state space ones, maybe in a separate