Markov Chains - Portfolio

Markov Chains #

Custom markov chain functions along with computing absorbing states and graph regularity.

.rmd knitted with execution results available for download here

Set up

1
if (!require("markovchain")) {
2
    install.packages("markovchain")
3
    library(markovchain)
4
}

1
## Loading required package: markovchain

1
## Warning: package 'markovchain' was built under R version 4.3.2

1
## Package:  markovchain
2
## Date:     2023-09-24 09:20:02 UTC
3
## BugReport: https://github.com/spedygiorgio/markovchain/issues

1
if (!require("igraph")) {
2
    install.packages("igraph")
3
    library(igraph)
4
}

1
## Loading required package: igraph

1
##
2
## Attaching package: 'igraph'

1
## The following objects are masked from 'package:stats':
2
##
3
##     decompose, spectrum

1
## The following object is masked from 'package:base':
2
##
3
##     union

Q1#

1
# set up helper functions
2

3
# limit of 'p^k' k->inf (actually, 1e5), is p^k positive
4
reg = function(p){
5
  pk = p
6
  for(i in c(1, 1e5)){
7
    if(reg1(pk)){
8
      # compute stationary probability
9
      # normalize eigen vector associated w/ eigen value '1' (always e$vectors[,1])
10
      e = eigen(p)
11
      # normalize eigen vector
12
      norm_e = abs(e$vectors[,1]) / sum(abs(e$vectors[,1]))
13
      # return(c(norm_e, e$vectors[,j]))
14
      return(norm_e)
15

16
    }
17
    pk = p %*% p
18
  }
19
  # didn't reach limit in reasonable time period
20
  return(FALSE)
21
}
22

23
# if matrix positive, return true
24
reg1 = function(p_in){
25
  p = p_in %*% p_in
26
  # loop thru rows
27
  for(i in c(1, length(p[1,]))){
28
    # loop thru cols
29
    for(j in c(1, length(p[,1]))){
30
      # entry non-positve
31
      if(p[i,j] <= 0){
32
        # entire matrix non-positive -> non-regular markov-chain
33
        return(FALSE)
34
      }
35
    }
36
  }
37

38
  return(TRUE)
39
}
40

41
# a markov chain has absorbing states if:
42
# some value on the main diagonal is '1' -> stuck at that state
43
# returns false if 'p' is not absorbing,
44
# returns the absorbing states otherwise
45
absorb = function(p, l=TRUE){
46
  # assuming p is valid (cols sum up to 1)
47
  val = FALSE
48
  count = 1
49
  for(i in c(1:length(p[,1]))){
50
    if(p[i,i] == 1){
51
      # absorbing state: 'i'
52
      if(typeof(val) == "logical"){ # first absorbing state
53
        val = c()
54
      }
55
      if(l){
56
        val[count] = paste("s",toString(i),sep="")
57
      } else {
58
        val[count] = i
59
      }
60
      count = count + 1
61
    }
62
  }
63
  return(val)
64
}

1
# Define matrices
2
p1 = matrix(c(c(1,0),c(0,1)), nrow=2, ncol=2)
3
p2 = matrix(c(c(0.2, 0.5, 0.3), c(0, 0.5, 0.5), c(0.5, 0.5, 0)), nrow=3, ncol=3)
4
p3 = matrix(c(c(0.1, 0.9, 0, 0), c(0.8, 0.1, 0.1, 0), c(0, 0.1, 0.3, 0.6), c(0, 0, 0.5, 0.5)), nrow=4, ncol=4)
5
p1

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
p2

1
##      [,1] [,2] [,3]
2
## [1,]  0.2  0.0  0.5
3
## [2,]  0.5  0.5  0.5
4
## [3,]  0.3  0.5  0.0

1
p3

1
##      [,1] [,2] [,3] [,4]
2
## [1,]  0.1  0.8  0.0  0.0
3
## [2,]  0.9  0.1  0.1  0.0
4
## [3,]  0.0  0.1  0.3  0.5
5
## [4,]  0.0  0.0  0.6  0.5

Plot the graphs.

1
G1 = make_graph(c(c(1,1), c(2,2)))
2
E(G1)$weights = c(1,1)
3
par(mar=c(0,0,0,0)+.1)
4
plot(G1)

1
G2 = make_graph(c(
2
  c(1,1), c(1,2),
3
  c(1,3), c(2,2), c(2,3),
4
  c(3,1), c(3,2)))
5
E(G2)$weights = c(0.2, 0.5, 0.3, 0.5, 0.5, 0.5, 0.5)
6
par(mar=c(0,0,0,0)+.1)
7
plot(G2, edge.arrow.size=E(G2)$weights*10, layout=layout_in_circle(G2))

1
G3 = make_graph(c(
2
  c(1,1), c(2,2), c(3,3), c(4,4),
3
  c(1,2), c(2,1), c(2,3), c(3,2),
4
  c(3,4), c(4,3)))
5
E(G3)$weights = c(0.1, 0.1, 0.3, 0.5, 0.9, 0.8, 0.1, 0.1, 0.6, 0.5)
6
par(mar=c(0,0,0,0)+.1)
7
plot(G3, edge.arrow.size=E(G3)$weights*10, layout=layout_in_circle(G3))

1
# is p1 regular?
2
reg(p1)

1
## [1] FALSE

1
# is p2 regular?
2
reg(p2)

1
## [1] 0.1923077 0.5000000 0.3076923

1
# is p3 regular?
2
reg(p3)

1
## [1] 0.2173913 0.2445652 0.2445652 0.2934783

So p1 is not regular, but p2 and p3 are, and their stationary distribution vectors are the normalized eigen vectors associated with the eigen value ‘1’.

1
# is p1 absorbing?
2
absorb(p1)

1
## [1] "s1" "s2"

1
# is p2 absorbing?
2
absorb(p2)

1
## [1] FALSE

1
# is p3 absorbing?
2
absorb(p3)

1
## [1] FALSE

So only ‘p1’ is absorbing, with the absorbing states: {1, 2}

Verify Q1 with ‘markovchain’ Library#

1
m1 = new('markovchain',
2
         transitionMatrix=p1,
3
         states=c('s1','s2'), byrow=FALSE)
4
m1

1
## Unnamed Markov chain
2
##  A  2 - dimensional discrete Markov Chain defined by the following states:
3
##  s1, s2
4
##  The transition matrix  (by cols)  is defined as follows:
5
##    s1 s2
6
## s1  1  0
7
## s2  0  1

1
is.regular(m1)

1
## [1] FALSE

1
absorbingStates(m1)

1
## [1] "s1" "s2"

1
# so P1 is not regular, but is absorbing, with the absorbing states {1, 2}, so my Q1 is right
2

3
m2 = new('markovchain',
4
         transitionMatrix=p2,
5
         states=c('s1','s2', 's3'), byrow=FALSE)
6
m2

1
## Unnamed Markov chain
2
##  A  3 - dimensional discrete Markov Chain defined by the following states:
3
##  s1, s2, s3
4
##  The transition matrix  (by cols)  is defined as follows:
5
##     s1  s2  s3
6
## s1 0.2 0.0 0.5
7
## s2 0.5 0.5 0.5
8
## s3 0.3 0.5 0.0

1
is.regular(m2)

1
## [1] TRUE

1
steadyStates(m2)

1
##         [,1]
2
## s1 0.1923077
3
## s2 0.5000000
4
## s3 0.3076923

1
absorbingStates(m2)

1
## character(0)

1
# so P2 is regular w/ same stationary probability, and isn't absorbing, so my Q1 is right
2

3
m3 = new('markovchain',
4
         transitionMatrix=p3,
5
         states=c('s1','s2', 's3', 's4'), byrow=FALSE)
6
m3

1
## Unnamed Markov chain
2
##  A  4 - dimensional discrete Markov Chain defined by the following states:
3
##  s1, s2, s3, s4
4
##  The transition matrix  (by cols)  is defined as follows:
5
##     s1  s2  s3  s4
6
## s1 0.1 0.8 0.0 0.0
7
## s2 0.9 0.1 0.1 0.0
8
## s3 0.0 0.1 0.3 0.5
9
## s4 0.0 0.0 0.6 0.5

1
is.regular(m3)

1
## [1] TRUE

1
steadyStates(m3)

1
##         [,1]
2
## s1 0.2173913
3
## s2 0.2445652
4
## s3 0.2445652
5
## s4 0.2934783

1
absorbingStates(m3)

1
## character(0)

1
# so P3 is regular w/ same stationary probability, and isn't absorbing, so my Q1 is right

So all my Q1 answers are correct.

Q2#

Create standard form function

1
# converts markov chain matrix to standard form
2
sForm = function(p){
3

4
  states = absorb(p, FALSE) # list of absorbing states
5
  if(typeof(states) == "logical"){ # no absorbing states
6
    return(FALSE)
7
  }
8

9

10
  # col swap
11
  # move the 'states' columns into the front
12
  M1 = p[,states]
13

14
  # now add non-absorbing columns
15
  M1 = cbind(M1, p[,-states])
16

17
  # print(M1)
18

19
  # row swap
20
  # move the '1's in the absorbing state columns to the top to create an identity matrix
21
  M2 = M1[states,] # grab row from M1
22

23
  # now add non-absorbing rows
24
  M2 = rbind(M2, M1[-states,])
25

26
  return(M2)
27
}
28

29
val_comp = function(sf, a){
30
  sf_l = length(sf[1,])
31
  x = sf_l - length(a)
32
  # print(x)
33
  # print(sf_l)
34
  I = diag(nrow=x ,ncol=x)
35
  Q = sf[(x+1):sf_l, (x+1):sf_l]
36

37
  N = solve(I-Q) # (I-Q)^-1
38

39
  Ti = list()
40
  for(r in c(1:length(N[,1]))){ # for each row
41
    Ti[r] = sum(N[r,])
42
  }
43

44
  R = sf[1:length(a), (x+1):sf_l] # gets sf: rows [1,length(a)), cols [length(a), length(sf[,1]))
45

46
  B = R %*% N
47

48
  print(R)
49
  print(N)
50

51
  # print('---- N is ---- ')
52
  # print(N)
53
  # print('---- Ti is ---- ')
54
  # print(Ti)
55
  # print('---- B is ---- ')
56
  # print(B)
57

58
  return(c(as.data.frame(N), as.data.frame(Ti), as.data.frame(B)))
59
}

1
P = matrix(
2
  c(0.2, 0.3, 0.4, 0.1, 0, 1, 0, 0, 0.5, 0.2, 0.2, 0.1, 0, 0, 0, 1),
3
  nrow=4, ncol=4)
4
P

1
##      [,1] [,2] [,3] [,4]
2
## [1,]  0.2    0  0.5    0
3
## [2,]  0.3    1  0.2    0
4
## [3,]  0.4    0  0.2    0
5
## [4,]  0.1    0  0.1    1

1
a = absorb(P)
2
a

1
## [1] "s2" "s4"

1
sf = sForm(P)
2
sf

1
##      [,1] [,2] [,3] [,4]
2
## [1,]    1    0  0.3  0.2
3
## [2,]    0    1  0.1  0.1
4
## [3,]    0    0  0.2  0.5
5
## [4,]    0    0  0.4  0.2

1
vc = val_comp(sf, a)

1
##      [,1] [,2]
2
## [1,]  0.3  0.2
3
## [2,]  0.1  0.1
4
##           [,1]     [,2]
5
## [1,] 1.8181818 1.136364
6
## [2,] 0.9090909 1.818182

1
# N
2
as.data.frame(vc[1:2])

1
##          V1       V2
2
## 1 1.8181818 1.136364
3
## 2 0.9090909 1.818182

1
# Ti
2
as.data.frame(vc[3:4])

1
##   X2.95454545454545 X2.72727272727273
2
## 1          2.954545          2.727273

1
# B
2
as.data.frame(vc[5:6])

1
##          V1        V2
2
## 1 0.7272727 0.7045455
3
## 2 0.2727273 0.2954545

1
M = new("markovchain",
2
        transitionMatrix=P,
3
        states=c('s1','s2','s3','s4'), byrow=FALSE)
4
M

1
## Unnamed Markov chain
2
##  A  4 - dimensional discrete Markov Chain defined by the following states:
3
##  s1, s2, s3, s4
4
##  The transition matrix  (by cols)  is defined as follows:
5
##     s1 s2  s3 s4
6
## s1 0.2  0 0.5  0
7
## s2 0.3  1 0.2  0
8
## s3 0.4  0 0.2  0
9
## s4 0.1  0 0.1  1

1
canonicForm(M)

1
## Unnamed Markov chain
2
##  A  4 - dimensional discrete Markov Chain defined by the following states:
3
##  s2, s4, s1, s3
4
##  The transition matrix  (by cols)  is defined as follows:
5
##    s2 s4  s1  s3
6
## s2  1  0 0.3 0.2
7
## s4  0  1 0.1 0.1
8
## s1  0  0 0.2 0.5
9
## s3  0  0 0.4 0.2

Which is the same as my standard form matrix, so my Q2 computes the right standard form matrix.

Markov Chains#

Q1#

Verify Q1 with ‘markovchain’ Library#

Q2#

Markov Chains #