Analisis de tiempos promedio

# Creando el vector de estados --------------------------------------------

estados = c()

for (i in 0:3) {
  for (j in 0:2) {
    for (k in 0:2) {
      #la funcion paste une los indices i,j,k y los separa con una coma
      estados = c(estados, paste(i,j,k, sep = ","))
    }
  }
}
estados

 [1] "0,0,0" "0,0,1" "0,0,2" "0,1,0" "0,1,1" "0,1,2" "0,2,0" "0,2,1" "0,2,2" "1,0,0"
[11] "1,0,1" "1,0,2" "1,1,0" "1,1,1" "1,1,2" "1,2,0" "1,2,1" "1,2,2" "2,0,0" "2,0,1"
[21] "2,0,2" "2,1,0" "2,1,1" "2,1,2" "2,2,0" "2,2,1" "2,2,2" "3,0,0" "3,0,1" "3,0,2"
[31] "3,1,0" "3,1,1" "3,1,2" "3,2,0" "3,2,1" "3,2,2"

# ¿Cuantos estados tengo?
length(estados)

[1] 36

# Creando la matriz de tasas ----------------------------------------------

matrizQ = matrix(data = 0, nrow = length(estados), ncol = length(estados),
                 dimnames = list(estados,estados))

#Noten que debemos separar mi string para conocer el valor de cada varible y convertirlo a #
estados[1]

[1] "0,0,0"

as.numeric(strsplit(estados[1],split = ",")[[1]][3])

[1] 0

for (fila in estados){
  for (columna in estados){
    
    #Para conocer los valores de cada variable de estado, vamos a dividir cada estado en sus tres componentes, utilizando strsplit. Para mayor facilidad, convertimos los estados de cada una de las variables en un valor numérico.
    
    i = as.numeric(strsplit(fila, ",")[[1]][1])
    j = as.numeric(strsplit(fila, ",")[[1]][2])
    k = as.numeric(strsplit(fila, ",")[[1]][3])
    
    #Una vez más, realizamos la división de cada estado en sus tres componentes.
    
    l = as.numeric(strsplit(columna, ",")[[1]][1])
    m = as.numeric(strsplit(columna, ",")[[1]][2])
    n = as.numeric(strsplit(columna, ",")[[1]][3])
    
    #Definición de la matriz de tasas de transición, de acuerdo con la formulación general mostrada antes.
    
    #Llegadas a la primera zona (prensas presión)
    if(i<3 & l==i+1 & m==j & n==k){
      matrizQ[fila,columna] = 0.05
    }
    
    #Procesamiento máquina satinadora (zona 3)
    if(k>0 & l==i & m==j & n==k-1){
      matrizQ[fila,columna] = (11/12)*0.2
    }
    
    #Procesamiento zona recubrimiento (zona 2)
    if(j>0 & k<2 & l==i & m==j-1 & n==k+1){
      matrizQ[fila,columna] = 0.125*min(j,2)
    }
    
    #Procesamiento prensa presión (zona 1)
    if(i>0 & j<2 & l==i-1 & m==j+1 & n==k){
      matrizQ[fila,columna] = 0.104
    }
    
  }
}
#Llenar la diagonal 
for (i in 1:length(estados)) {
  matrizQ[i,i] = -sum(matrizQ[i,])
}

rowSums(matrizQ)

        0,0,0         0,0,1         0,0,2         0,1,0         0,1,1         0,1,2 
 0.000000e+00  1.387779e-17  1.387779e-17  1.387779e-17  1.387779e-17  1.387779e-17 
        0,2,0         0,2,1         0,2,2         1,0,0         1,0,1         1,0,2 
 1.387779e-17  1.387779e-17  1.387779e-17  0.000000e+00 -2.775558e-17 -2.775558e-17 
        1,1,0         1,1,1         1,1,2         1,2,0         1,2,1         1,2,2 
-2.775558e-17 -2.775558e-17 -2.775558e-17  1.387779e-17  1.387779e-17  1.387779e-17 
        2,0,0         2,0,1         2,0,2         2,1,0         2,1,1         2,1,2 
 0.000000e+00 -2.775558e-17 -2.775558e-17 -2.775558e-17 -2.775558e-17 -2.775558e-17 
        2,2,0         2,2,1         2,2,2         3,0,0         3,0,1         3,0,2 
 1.387779e-17  1.387779e-17  1.387779e-17  0.000000e+00  1.387779e-17  1.387779e-17 
        3,1,0         3,1,1         3,1,2         3,2,0         3,2,1         3,2,2 
 1.387779e-17  1.387779e-17  1.387779e-17  0.000000e+00  0.000000e+00  0.000000e+00 

library(markovchain)

cadenaContinua = new(Class = "ctmc", states = estados, byrow = T, 
                     generator = matrizQ)


#Para calcular el tiempo de primera pasada en cadenas continuas usamos la funcion
#expected tie, que recibe como parametros la cadena, el indice del estado inicial y el
#indice del estado final.

ExpectedTime(cadenaContinua, 14, 36)

[1] 44224.58

# a veces tenemos muchos estaods y puedeo ser laborioso encontrar los indices correpondientes
# para eso utilizamos la funcion which

estado_inicial = which(estados == "1,1,1")
estado_final = which(estados == "3,2,2")

ExpectedTime(cadenaContinua, estado_inicial, estado_final)

[1] 44224.58