COMPLETE Code BY RSTUDIO of Markov Model

MarkovModelCodeinR.txt

library(markovchain)
weatherStates <- c("sunny", "cloudy", "rain")
byRow <- TRUE
weatherMatrix <- matrix(data = c(0.50, 0.25, 0.25,
0.5, 0.1, 0.4,
0.1, 0.7, 0.2), byrow = byRow, nrow = 3,
dimnames = list(weatherStates, weatherStates))
#transition matrix and probabilities
mcWeather <- new("markovchain", states = weatherStates, byrow =
byRow,transitionMatrix = weatherMatrix, name = "Weather")
defaultMc <- new("markovchain") #Default markov created
mcList <- new("markovchainList", markovchains = list(mcWeather, defaultMc),
name = "A list of Markov chains")
initialState <- c(0, 1, 0)
after2Days <- (t(mcWeather) * t(mcWeather)) * initialState
after7Days <- (t(mcWeather) ^ 7) * initialState
after2Days

fvals<-function(mchain,initialstate,n) {
out<-data.frame()
names(initialstate)<-names(mchain)
for (i in 0:n)
{
iteration<-initialstate*mchain^(i)
out<-rbind(out,iteration)
}
out<-cbind(out, i=seq(0,n))
out<-out[,c(4,1:3)]
return(out)
}
fvals(mchain=mcWeather,initialstate=c(80,10,10),n=4)

states(mcWeather)
names(mcWeather)
dim(mcWeather)


name(mcWeather)
name(mcWeather) <- "New Name"
name(mcWeather)

markovchain:::sort(mcWeather)

transitionProbability(mcWeather, "cloudy", "rain")

mcWeather[2,3]

print(mcWeather)

show(mcWeather)

plot(mcWeather, package="diagram", box.size = 0.05)

mcDf <- as(mcWeather, "data.frame")
mcNew <- as(mcDf, "markovchain")
mcDf

mcIgraph <- as(mcWeather, "igraph")
require(msm)
## Loading required package: msm
Q <- rbind ( c(0, 0.25, 0, 0.25),
+ c(0.19, 0, 0.19, 0.19),
+ c(0, 0.2, 0, 0.2),
+ c(0, 0, 0, 0) )
cavmsm <- msm(state ~ years, subject = PTNUM, data = cav, qmatrix = Q, death
= 4)
msmMc <- as(cavmsm, "markovchain")
msmMc

library(etm)
data(sir.cont)
sir.cont <- sir.cont[order(sir.cont$id, sir.cont$time), ]
for (i in 2:nrow(sir.cont)) {
if (sir.cont$id[i]==sir.cont$id[i-1]) {
if (sir.cont$time[i]==sir.cont$time[i-1]) {
sir.cont$time[i-1] <- sir.cont$time[i-1] - 0.5
}
}
}
tra <- matrix(ncol=3,nrow=3,FALSE)
tra[1, 2:3] <- TRUE
tra[2, c(1, 3)] <- TRUE
tr.prob <- etm(sir.cont, c("0", "1", "2"), tra, "cens", 1)
tr.prob

etm2mc<-as(tr.prob, "markovchain")
etm2mc

myMatr<-matrix(c(.1,.7,.2,.2,.7,.1,.4,.4,.2), byrow=TRUE, ncol=3)
myMc<-as(myMatr, "markovchain")
myMc

plot(myMc)

stateNames = c("H", "I", "D")
Q0 <- new("markovchain", states = stateNames,
transitionMatrix =matrix(c(0.6, 0.3, 0.1,0.3, 0.4, 0.3,0, 0, 1),
byrow = TRUE, nrow = 3), name = "state t0")
Q1 <- new("markovchain", states = stateNames,
transitionMatrix = matrix(c(0.1, 0.6, 0.3,0, 0.4, 0.6,0, 0, 1),
byrow = TRUE, nrow = 3), name = "state t1")
Q2 <- new("markovchain", states = stateNames,
transitionMatrix = matrix(c(0.2, 0.3, 0.5,0.2, 0, 0.8,0, 0, 1),
byrow = TRUE,nrow = 3), name = "state t2")
Q3 <- new("markovchain", states = stateNames,
transitionMatrix = matrix(c(0, 0, 1, 0, 0, 1, 0, 0, 1),
byrow = TRUE, nrow = 3), name = "state t3")
mcCCRC <- new("markovchainList",markovchains = list(Q0,Q1,Q2,Q3),
name = "Continuous Care Health Community")
print(mcCCRC)

mcCCRC[[1]]

dim(mcCCRC)

conditionalDistribution(mcWeather, "sunny")

steadyStates(mcWeather)

gamblerRuinMarkovChain <- function(moneyMax, prob = 0.5) {
require(matlab)
matr <- zeros(moneyMax + 1)
states <- as.character(seq(from = 0, to = moneyMax, by = 1))
rownames(matr) = states; colnames(matr) = states
matr[1,1] = 1; matr[moneyMax + 1,moneyMax + 1] = 1
for(i in 2:moneyMax)
{ matr[i,i-1] = 1 - prob; matr[i, i + 1] = prob }
out <- new("markovchain",
transitionMatrix = matr,
name = paste("Gambler ruin", moneyMax, "dim", sep = " ")
)
return(out)
}
mcGR4 <- gamblerRuinMarkovChain(moneyMax = 4, prob = 0.5)

steadyStates(mcGR4)
absorbingStates(mcGR4)
absorbingStates(mcWeather)
.commclassesKernel <- function(P){
m <- ncol(P)
stateNames <- rownames(P)
T <- zeros(m)
i <- 1
while (i <= m) {
a <- i
b <- zeros(1,m)
b[1,i] <- 1

old <- 1
new <- 0
while (old != new) {
old <- sum(find(b > 0))
n <- size(a)[2]
matr <- matrix(as.numeric(P[a,]), ncol = m,
nrow = n)
c <- colSums(matr)
d <- find(c)
n <- size(d)[2]
b[1,d] <- ones(1,n)
new <- sum(find(b>0))
a <- d}
T[i,] <- b
i <- i+1 }
F <- t(T)
C <- (T > 0)&(F > 0)
v <- (apply(t(C) == t(T), 2, sum) == m)
colnames(C) <- stateNames
rownames(C) <- stateNames
names(v) <- stateNames
out <- list(C = C, v = v)
return(out)
}
P <- matlab::zeros(10)
P[1, c(1, 3)] <- 1/2;
P[2, 2] <- 1/3; P[2,8] <- 2/3;
P[3, 2] <- 1;
P[4, 3] <- 1;
P[5, c(4, 5, 8)] <- 1/3;
P[6, 6] <- 1;
P[7, 7] <- 1/4; P[7,9] <- 3/4;
P[8, c(3, 4, 8, 9)] <- 1/4;
P[9, 3] <- 1;
P[10, c(2, 5, 9)] <- 1/3;
rownames(P) <- letters[1:10]
colnames(P) <- letters[1:10]
probMc <- new("markovchain", transitionMatrix = P,
name = "Probability MC")
.commclassesKernel(P)


summary(probMc)

transientStates(probMc)

probMcCanonic <- canonicForm(probMc)
probMc

is.accessible(object = probMc, from = "a", to = "c")
is.accessible(object = probMc, from = "g", to = "c")
E <- matrix(0, nrow = 4, ncol = 4)
E[1, 2] <- 1
E[2, 1] <- 1/3; E[2, 3] <- 2/3
E[3,2] <- 1/4; E[3, 4] <- 3/4
E[4, 3] <- 1
mcE <- new("markovchain", states = c("a", "b", "c", "d"),
transitionMatrix = E,
name = "E")
is.irreducible(mcE)

period(mcE)
require(matlab)
mathematicaMatr <- zeros(5)
mathematicaMatr[1,] <- c(0, 1/3, 0, 2/3, 0)
mathematicaMatr[2,] <- c(1/2, 0, 0, 0, 1/2)
mathematicaMatr[3,] <- c(0, 0, 1/2, 1/2, 0)
mathematicaMatr[4,] <- c(0, 0, 1/2, 1/2, 0)
mathematicaMatr[5,] <- c(0, 0, 0, 0, 1)
statesNames <- letters[1:5]


mathematicaMc <- new("markovchain", transitionMatrix = mathematicaMatr,
name = "Mathematica MC", states = statesNames)
.firstpassageKernel <- function(P, i, n){
G <- P
H <- P[i,]
E <- 1 - diag(size(P)[2])
for (m in 2:n) {
G <- P %*% (G * E)
H <- rbind(H, G[i,])
}
return(H)
}
firstPassagePdF <- firstPassage(object = mcWeather, state = "sunny",
n = 10)
firstPassagePdF[3, 3]
plot(myMc)