miércoles, 12 de octubre de 2016

21-Cadenas de Markov 3

Simulación de Cadenas de Markov 3 de 10

Simulación de m trayectorias de tamaño n de una cadena de Markov cuando conocemos su matriz de probabilidades de transición.

Caso cuando el espacio de estados es finito.

## Primero la función que simula una sola trayectoria.

simula_1_Cadena <- function(n,X0,P){
        dim <- length(P[1,])
        Xn <- numeric((n+1))
        Xn[1] <- X0
        for(i in 2:(n+1)){
                aux <- Xn[i-1]
                Xn[i] <- sample(1:dim,1,T,P[aux,])
        }
        Xn
}
## Ahora simulamos m trayectorias
## Los argumentos serán
## m <- número de cadenas a simular
## n <- número de pasos
## X0 <- estado inicial
## P <- matriz de probabilidades de transición

simula_m_Cadenas <- function(m,n,X0,P){
        MatResultados <- matrix(0,m,(n+1))
        for(i in 1:m){
        MatResultados[i,] <- simula_1_Cadena(n,X0,P)      
        }
        estados <- numeric(n+1)
        for(i in 0:n) estados[i+1] <- paste('X',i)
        colnames(MatResultados) <- estados
        MatResultados
}

Ejemplo 1.

\[P= \begin{bmatrix} 1/2 & 1/2 & 0 \\ 1/2 & 1/3 & 1/6\\ 0 & 1 & 0 \end{bmatrix}\]

p0 <- c(1/2,1/2,0)
p1 <- c(1/2,1/3,1/6)
p2 <- c(0,1,0)
P <- rbind(p0,p1,p2)
simula_m_Cadenas(16,6,2,P)
##       X 0 X 1 X 2 X 3 X 4 X 5 X 6
##  [1,]   2   2   2   1   2   2   1
##  [2,]   2   2   2   1   1   2   1
##  [3,]   2   1   1   1   2   2   1
##  [4,]   2   1   1   2   2   3   2
##  [5,]   2   1   1   2   3   2   2
##  [6,]   2   1   2   2   2   1   1
##  [7,]   2   1   2   1   2   1   2
##  [8,]   2   1   2   3   2   1   1
##  [9,]   2   2   1   2   1   2   2
## [10,]   2   1   2   1   2   1   1
## [11,]   2   2   1   1   1   2   2
## [12,]   2   1   1   1   1   1   2
## [13,]   2   1   1   1   1   2   1
## [14,]   2   1   2   2   3   2   2
## [15,]   2   1   1   1   2   2   3
## [16,]   2   2   2   3   2   1   2

Ejemplo Cadena de Erhenfest

#Ejemplo Erhenfest N=9
p0=c(0,1,0,0,0,0,0,0,0,0)
p1=c(1/9,0,8/9,0,0,0,0,0,0,0)
p2=c(0,2/9,0,7/9,0,0,0,0,0,0)
p3=c(0,0,3/9,0,6/9,0,0,0,0,0)
p4=c(0,0,0,4/9,0,5/9,0,0,0,0)
p5=c(0,0,0,0,5/9,0,4/9,0,0,0)
p6=c(0,0,0,0,0,6/9,0,3/9,0,0)
p7=c(0,0,0,0,0,0,7/9,0,2/9,0)
p8=c(0,0,0,0,0,0,0,8/9,0,1/9)
p9=c(0,0,0,0,0,0,0,0,1,0)
P=rbind(p0,p1,p2,p3,p4,p5,p6,p7,p8,p9)
simula_m_Cadenas(15,10,3,P)
##       X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 10
##  [1,]   3   4   5   4   5   4   3   4   3   4    5
##  [2,]   3   2   3   4   5   4   5   6   7   6    7
##  [3,]   3   4   5   6   5   6   5   4   5   6    5
##  [4,]   3   4   3   4   5   4   5   4   5   6    7
##  [5,]   3   4   5   6   5   6   7   6   5   4    3
##  [6,]   3   4   3   4   3   4   5   4   3   4    5
##  [7,]   3   4   5   4   5   6   5   6   7   8    7
##  [8,]   3   4   3   4   5   4   5   6   5   6    5
##  [9,]   3   4   3   4   3   4   3   4   5   6    5
## [10,]   3   4   5   4   5   4   5   4   5   6    7
## [11,]   3   4   3   2   3   4   5   6   5   6    5
## [12,]   3   2   3   4   3   2   3   2   1   2    3
## [13,]   3   2   3   4   5   4   3   4   5   4    5
## [14,]   3   2   3   4   3   2   3   4   5   6    5
## [15,]   3   4   3   4   3   4   3   4   5   6    5

Graficando las trayectorias

#Ejemplo Erhenfest N=9
trayectorias <- t(simula_m_Cadenas(15,1000,3,P))
matplot(0:1000,trayectorias,type = 'l', pch = 16)

No hay comentarios:

Publicar un comentario