# ==============================================================================
# This code allows one to analysis some properties of a Markov chain of 3-states.
#
# Modelled & Coded by Ebert Brea
# ==============================================================================
require(markovchain)
## Loading required package: markovchain
## Package: markovchain
## Version: 0.8.6
## Date: 2021-05-17
## BugReport: https://github.com/spedygiorgio/markovchain/issues
library(latex2exp)
library(ggplot2)
library(shape)
library(diagram)
#===================================================
#
#===================================================
knowledge=TRUE # TRUE if the distribution function of first order PIn at n=0 is known
# FALSE if the distribution function of first order PIn at n=0 is unknown
#===================================================
P = matrix(c(0.4,0.1,0.5,
0.5,0 ,0.5,
0.5,0.5,0),nrow = 3,byrow = TRUE)
mc = new("markovchain",transitionMatrix=P,states=c("X1","X2","X3"),name="Markov Chain")
# For reporting respectively: Markov chain data frame; summary;
# and the recurrent, transient and absorbing states of the Markov chain
str(mc)
## Formal class 'markovchain' [package "markovchain"] with 4 slots
## ..@ states : chr [1:3] "X1" "X2" "X3"
## ..@ byrow : logi TRUE
## ..@ transitionMatrix: num [1:3, 1:3] 0.4 0.5 0.5 0.1 0 0.5 0.5 0.5 0
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:3] "X1" "X2" "X3"
## .. .. ..$ : chr [1:3] "X1" "X2" "X3"
## ..@ name : chr "Markov Chain"
summary(mc)
## Markov Chain Markov chain that is composed by:
## Closed classes:
## X1 X2 X3
## Recurrent classes:
## {X1,X2,X3}
## Transient classes:
## NONE
## The Markov chain is irreducible
## The absorbing states are: NONE
recurrentClasses(mc)
## [[1]]
## [1] "X1" "X2" "X3"
transientStates(mc)
## character(0)
absorbingStates(mc)
## character(0)
steadyStates(mc)
## X1 X2 X3
## [1,] 0.4545455 0.2121212 0.3333333
# For drawing the Markov chain graph
plot(mc,package="DiagrammeR")
#########################################################
R=30 # Number of steps %
St=3 # Number of states
# Defining matrix according to: the number of steps & the number of states
x <-matrix( nrow=R+1, ncol = St)
t <-matrix( nrow=R+1, ncol = 1)
if(knowledge){
# Distribution function for each state of Xn at n=0
x[1,]=c(1,0,0)
}else{
# unknowledged Distribution function for each state of Xn at n=0
p= runif(1,0,1)
x[1,]=c((1-p)^2,2*p*(1-p),p^2)
}
print(x[1,])
## [1] 1 0 0
for(n in 1:R){
x[n+1,]=x[n,]*mc
t[n]=n-1
}
# For reporting the distribution function for each state of x
print(x)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.0000000 0.0000000
## [2,] 0.4000000 0.1000000 0.5000000
## [3,] 0.4600000 0.2900000 0.2500000
## [4,] 0.4540000 0.1710000 0.3750000
## [5,] 0.4546000 0.2329000 0.3125000
## [6,] 0.4545400 0.2017100 0.3437500
## [7,] 0.4545460 0.2173290 0.3281250
## [8,] 0.4545454 0.2095171 0.3359375
## [9,] 0.4545455 0.2134233 0.3320313
## [10,] 0.4545455 0.2114702 0.3339844
## [11,] 0.4545455 0.2124467 0.3330078
## [12,] 0.4545455 0.2119585 0.3334961
## [13,] 0.4545455 0.2122026 0.3332520
## [14,] 0.4545455 0.2120805 0.3333740
## [15,] 0.4545455 0.2121416 0.3333130
## [16,] 0.4545455 0.2121110 0.3333435
## [17,] 0.4545455 0.2121263 0.3333282
## [18,] 0.4545455 0.2121187 0.3333359
## [19,] 0.4545455 0.2121225 0.3333321
## [20,] 0.4545455 0.2121206 0.3333340
## [21,] 0.4545455 0.2121215 0.3333330
## [22,] 0.4545455 0.2121211 0.3333335
## [23,] 0.4545455 0.2121213 0.3333333
## [24,] 0.4545455 0.2121212 0.3333334
## [25,] 0.4545455 0.2121212 0.3333333
## [26,] 0.4545455 0.2121212 0.3333333
## [27,] 0.4545455 0.2121212 0.3333333
## [28,] 0.4545455 0.2121212 0.3333333
## [29,] 0.4545455 0.2121212 0.3333333
## [30,] 0.4545455 0.2121212 0.3333333
## [31,] 0.4545455 0.2121212 0.3333333
# Making data frames
df <- as.data.frame(x)
xend=R-1
# Drawing distribution function
plot(ylim=c(0, 1),t,df$V1,main=TeX('$\tDistribution \tfunction \tfor \teach \tstate \tin \\{X_{n}^{(i)}$\\}_{i=1}^{3}'),xlab=TeX('$n$'),ylab=TeX('$p_{x}(n)$'),type="o",lwd = 2,col="blue")
axis(1,at=0:xend)
lines(t,df$V2,col="red",type="o",lwd = 2)
lines(t,df$V3,col="darkgreen",type="o",lwd = 2)
legend("topright",c(TeX('$\tState X^{(1)}$'),TeX('$\tState X^{(2)}$'),TeX('$\tState X^{(3)}$')), fill=c("blue","red","darkgreen"))
summary.data.frame(df)
## V1 V2 V3
## Min. :0.4000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.4545 1st Qu.:0.2121 1st Qu.:0.3333
## Median :0.4545 Median :0.2121 Median :0.3333
## Mean :0.4705 Mean :0.2033 Mean :0.3262
## 3rd Qu.:0.4545 3rd Qu.:0.2121 3rd Qu.:0.3333
## Max. :1.0000 Max. :0.2900 Max. :0.5000
print.data.frame(df)
## V1 V2 V3
## 1 1.0000000 0.0000000 0.0000000
## 2 0.4000000 0.1000000 0.5000000
## 3 0.4600000 0.2900000 0.2500000
## 4 0.4540000 0.1710000 0.3750000
## 5 0.4546000 0.2329000 0.3125000
## 6 0.4545400 0.2017100 0.3437500
## 7 0.4545460 0.2173290 0.3281250
## 8 0.4545454 0.2095171 0.3359375
## 9 0.4545455 0.2134233 0.3320313
## 10 0.4545455 0.2114702 0.3339844
## 11 0.4545455 0.2124467 0.3330078
## 12 0.4545455 0.2119585 0.3334961
## 13 0.4545455 0.2122026 0.3332520
## 14 0.4545455 0.2120805 0.3333740
## 15 0.4545455 0.2121416 0.3333130
## 16 0.4545455 0.2121110 0.3333435
## 17 0.4545455 0.2121263 0.3333282
## 18 0.4545455 0.2121187 0.3333359
## 19 0.4545455 0.2121225 0.3333321
## 20 0.4545455 0.2121206 0.3333340
## 21 0.4545455 0.2121215 0.3333330
## 22 0.4545455 0.2121211 0.3333335
## 23 0.4545455 0.2121213 0.3333333
## 24 0.4545455 0.2121212 0.3333334
## 25 0.4545455 0.2121212 0.3333333
## 26 0.4545455 0.2121212 0.3333333
## 27 0.4545455 0.2121212 0.3333333
## 28 0.4545455 0.2121212 0.3333333
## 29 0.4545455 0.2121212 0.3333333
## 30 0.4545455 0.2121212 0.3333333
## 31 0.4545455 0.2121212 0.3333333