ej4 tproba v2

library("markovchain")

matrix.power <- function(A, n) {
  if(n == 0)
    return(diag(dim(P)[1]))
  e <- eigen(A)
  M <- e$vectors
  d <- e$values

  if(length(M) == dim(P)[1])   # si es diagonalizable
    return(M %*% diag(d^n) %*% solve(M))
  else{
    newA = diag(dim(P)[1])
    for(i in 1:n)
      newA = newA %*% A
    return(newA)
  }
}

S = c("A", "B", "C", "D")

mu0 = c(1/8, 2/8, 1/8, 4/8)

P = matrix(
  c(0.1, 0.1, 0.8, 0,
    0.4, 0  , 0.5, 0.1,
    0.1, 0  , 0  , 0.9,
    0.4, 0.1, 0.5, 0),
  nrow = 4,
  ncol = 4,
  byrow = TRUE
)

if( sum(mu0) != 1 ) #el vector de distrib. inicial k debe sumar 1
  stop( "Las probabilidades del vector inicial no suman 1" )

for(i in 1:dim(P)[1]){
  if( sum (P[i,]) != 1) #las filas deben sumar 1
    stop( "Las filas individuales de P no suman 1" )

  if( P[i,i] == 1 ) #no debe haber estados absorbentes por enunciado
    stop( "Hay algún estado absorbente presente" )
}

Q = new("markovchain", transitionMatrix = P,
        states = S, name = "Buffer")

plot(Q, main = "Cadena de Markov", edge.arrow.size=0.5)

#simulate
tests=10000
n = 100

T = markovchainSequence(n, Q, t0 = sample(Q@states, prob = mu0, replace = TRUE, 1), include.t0 = TRUE)
plot(0:n, factor(T), yaxt = "n", main = "Simulación de una Trayectoria", ylab = "Estado", xlab = "Tiempo", type = "l", col = "blue")
axis(side=2, at = 1:4, labels = S)

R = c(0, 0, 0, 0)
names(R) = S
T
for(i in T)
  R[i] = R[i] + 1

print("Numero de visitas a los nodos (en n=100 pasos) simulada: ")
print(R)

#newP = P^n
newP = matrix.power(P, n)
muN = mu0 %*% newP  #probabilidad de estar en el i-esimo nodo luego de 'n' pasos