Simulación

Para poder simular la cadena de markov vamos a crear una función que simula una variable aleatoria discreta.

DiscretaN <- function(probas,valores) { # inicio funcion

  if (sum(probas) != 1) { #este if es un control para que no entreguen
    #valores de probas mayores a uno, o no representen una distribucion
    return(NA)
  }

  if (length(probas) != length(valores)) { #este if es un control que
    #indica que si mete diferente numero de probas que de valores nos
    #devuelva un na
    return(NA)
  }

  min = 0 #el inicio del primer intervalo
  max = probas[1] #fin del primer intervalo
  X = "no funciona" #control del variable final
  uniforme = runif(n = 1, min = 0,max = 1) #numero aleatorio enre 0 y 1

  for (n in 1:length(probas)) { #inicio del for

    if (uniforme >= min & uniforme <= max) { #si nuestro numero aleatorio
      #esta en el intervalo que nos devuelva el valor de ese intervalo
      X = valores[n]
      break
    }
    #si llegamos a esta parte del for es por que no encontramos al numero
    #aleatorio en los otros intervalos
    min = min + probas[n] #asi que el nuevo intervalo va ser el inicio del
    #intervalo anterior mas la proba del intervalo que acaba de pasar
    max = max + probas[n+1] #mientras que el max va ser el max anterior mas
    #la proba siguiente
  } #fin del for
  return(X)
} # fin funcion

Esta función nos va a servir, para poder simular el paso que dara dado el presente, y la información se dara gracias a la matriz de probabilidad.

5 trayectorias de 15 pasos

Pero antes de que podamos simular los pasos que vamos a dar, primero vamos a escoger una distribución inicial, la cual va a ser la siguiente:

\[ \pi = (1/8 \hspace{0.25 cm} 1/8 \hspace{0.25 cm} 1/8 \hspace{0.25 cm} 1/8 \hspace{0.25 cm}1/8 \hspace{0.25 cm} 1/8 \hspace{0.25 cm}1/8 \hspace{0.25 cm}1/8) \] y la matriz dada.

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
P = matrix(data =  c(0.1, 0.1, 0.2, 0.2, 0.1, 0.1, 0.1, 0.1,
                    0, 0.1, 0.1, 0.1, 0, 0.3, 0.2, 0.2,
                    0.6, 0, 0, 0.1, 0.1, 0.1, 0.1, 0,
                    0, 0, 0, 0.3, 0.7, 0, 0, 0, 
                    0, 0, 0, 0.7, 0.3, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0.3, 0.4, 0.3, 
                    0, 0, 0, 0, 0, 0.1, 0, 0.9, 
                    0, 0, 0, 0, 0, 0.8, 0.2, 0),
           nrow = 8, 
           ncol = 8, 
           byrow = TRUE) #matriz estocastica
E = c(0,1,2,3,4,5,6,7) # espacio de estados
Pi = c(1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8) #distribucion inicial

colnames(P) = E #nombre de los estados
rownames(P) = E #nombre de los estados

P
##     0   1   2   3   4   5   6   7
## 0 0.1 0.1 0.2 0.2 0.1 0.1 0.1 0.1
## 1 0.0 0.1 0.1 0.1 0.0 0.3 0.2 0.2
## 2 0.6 0.0 0.0 0.1 0.1 0.1 0.1 0.0
## 3 0.0 0.0 0.0 0.3 0.7 0.0 0.0 0.0
## 4 0.0 0.0 0.0 0.7 0.3 0.0 0.0 0.0
## 5 0.0 0.0 0.0 0.0 0.0 0.3 0.4 0.3
## 6 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.9
## 7 0.0 0.0 0.0 0.0 0.0 0.8 0.2 0.0
tra = -1 #inicio para guardar las trayectoria

for (j in 1:5) {
  S = 0 
  S[1] = DiscretaN(probas = Pi, valores = E) #uso de la distribucion incial
  
  for (i in 2:15) { #siguientes 14 pasos
    presente = S[i-1] #tomamos el paso anterior como presente
    indice = presente + 1 #indice dentro del lenguaje (R empieza en 1)
    paso = P[indice,] # la probabilidad del siguiente paso
    futuro = DiscretaN(probas = paso, valores = E) #simulamos el futuro
    S[i] = futuro #ya simulado se la asignamos el paso siguiente
  }
 tra = c(tra,S)#guardamos la simulacion
}
 tra = tra[-1]

 trayectorias = data.frame(paso = rep(1:15, 5),
                           numeroSim = factor(rep(1:5, each = 15)),
                           trayectoria = tra)#lo hacemos la tabla.
 
 #lo graficamos para observalo

 trayectorias %>% ggplot(aes(x = paso, 
                             y = trayectoria, 
                             col = numeroSim, 
                             type = numeroSim)) +
   facet_wrap(~numeroSim)+
   geom_line() +  
   scale_y_continuous(breaks = E) + 
   scale_x_continuous(breaks = 1:15) + 
   theme_bw() + 
   labs(title = "Simulación de 5 trayectorias",
        col = "Numero de simulación",
        y = "Trayectoria",
        x = "Paso") + theme(axis.text.x = element_text(angle = 90))

Aqui podemos ver la simulación de cada una de las trayectorias que toma la cadena de markov, hay que notar que en la cadena tiene clases cerradas por eso se identifican algunos patrones de trayectoria.

1000 trayectorias de 15 pasos

tra = -1 #inicio para guardar las trayectoria

for (j in 1:1000) {
  S = 0 
  S[1] = DiscretaN(probas = Pi, valores = E) #uso de la distribucion incial
  
  for (i in 2:15) { #siguientes 14 pasos
    presente = S[i-1] #tomamos el paso anterior como presente
    indice = presente + 1 #indice dentro del lenguaje (R empieza en 1)
    paso = P[indice,] # la probabilidad del siguiente paso
    futuro = DiscretaN(probas = paso, valores = E) #simulamos el futuro
    S[i] = futuro #ya simulado se la asignamos el paso siguiente
  }
 tra = c(tra,S)#guardamos la simulacion
}
 tra = tra[-1]

trayectorias = data.frame(paso = rep(1:15, 1000),
                           numeroSim = factor(rep(1:1000, each = 15)),
                           trayectoria = tra)#lo hacemos la tabla.

Tabla = trayectorias %>% 
  group_by(Estado = trayectoria) %>% 
  count() %>% 
  mutate(Proporcion = n/15000) %>%  select(Estado, Proporcion)

knitr::kable(Tabla)
Estado Proporcion
0 0.0196000
1 0.0106000
2 0.0144000
3 0.1994667
4 0.1954000
5 0.2365333
6 0.1357333
7 0.1882667

esta tabla es una recopilacion 1000 cadenas de quince pasos, por eso podemos ver que hay proporcion en cada uno de los estados, ya que nuestra distribucion inicial tiene la misma probabilidad empezar por cada uno de los estados.

5 trayectoria de mil pasos.

tablas = list()
for (j in 1:5) {
  

  S = 0 
  S[1] = DiscretaN(probas = Pi, valores = E) #uso de la distribucion incial
  
  for (i in 2:1000) { #siguientes 14 pasos
    presente = S[i-1] #tomamos el paso anterior como presente
    indice = presente + 1 #indice dentro del lenguaje (R empieza en 1)
    paso = P[indice,] # la probabilidad del siguiente paso
    futuro = DiscretaN(probas = paso, valores = E) #simulamos el futuro
    S[i] = futuro #ya simulado se la asignamos el paso siguiente
  }
  
  Trayectoria = data.frame(S)
  
  tabla = Trayectoria %>% 
    group_by(Estado = S) %>% 
    count() %>% 
    mutate(proporcion = n/1000)
  
  tablas[[j]] = tabla
  
}

Trayectoria 1

knitr::kable(tablas[[1]])
Estado n proporcion
5 422 0.422
6 237 0.237
7 341 0.341

Trayectoria 2

knitr::kable(tablas[[2]])
Estado n proporcion
5 408 0.408
6 242 0.242
7 350 0.350

Trayectoria 3

knitr::kable(tablas[[3]])
Estado n proporcion
0 1 0.001
2 1 0.001
3 498 0.498
4 500 0.500

Trayectoria 4

knitr::kable(tablas[[4]])
Estado n proporcion
5 427 0.427
6 231 0.231
7 342 0.342

Trayectoria 5

knitr::kable(tablas[[5]])
Estado n proporcion
5 422 0.422
6 239 0.239
7 339 0.339

Conclusión.

Ya con esta información se puede observar un poco mejor el comportamiento de las clases cerradas, por esa razon algunos estados no aparecen en algunas trayectorias.