Simulación de Cadenas de Markov 1 de 10
David J. Santana
9 de octubre de 2016
Simulación de una trayectoria de una cadena de Markov.
Caminata aleatoria simple.
Damos una definición breve de una cadena de Markov con espacio de estados y espacio temporal discretos. Además simulamos el fácil ejemplo de la caminata aleatoria simple.
Definición 1.
Supongamos una sucesión de v.a.’s \(X_0, X_1, X_2, \ldots\) cuyo soporte es un subconjunto \(\mathcal{S}\) de los números naturales. Si para cualesquiera números enteros \(n\ge 1\) y \(m\ge 0\) se cumple que
- \[P(X_n = \alpha_n\mid X_{n-1}=\alpha_{n-1},\ldots, X_0=\alpha_0) = P(X_n = \alpha_n\mid X_{n-1}=\alpha_{n-1}),\]
- \[P(X_{n+m} =y\mid X_m = x) = P(X_n = y \mid X_0 = x),\]
donde \(\alpha_0, \alpha_1 , \ldots, \alpha_n, x, y \in \mathcal{S}\), entonces diremos que \(\left\{X_i\right\}_{i=0}^{\infty}\) es una cadena de Markov homogénea en el tiempo (CMHT). La primer condición es conocida como propiedad de Markov, y la segunda es la condición de homogeneidad en el tiempo.
Ejemplo 1.
Supongamos que \(\xi_1, \xi_2, \ldots\) es una sucesión de v.a.i.i.d. tales que \(P(\xi_k= 1)= p\) y \(P(\xi_k= -1)= 1 - p\) para \(k\ge 1\). Sea \(X_n = X_{n-1}+\xi_n\) para \(n\ge 1\) y definimos \(X_0=0\). Se puede demostrar que \(\left\{X_i\right\}_{i=0}^{\infty}\) es una CMHT. Esta cadena es conocida como caminata aleatoria simple.
## Simulación del valor de la cadena al tiempo n con probabilidad p de
## incrementar en uno el valor de la cadena después de un paso.
set.seed(123)
n <- 10
p <- 0.54
X0 <- 0
tornillos <- 2*rbinom(n,1,p)-1
X <- numeric(n)
X[1] <- X0 + tornillos[1]
for(i in 2:n){
X[i] <- X[i-1] + tornillos[i]
}
paste('Valor de la cadena al tiempo', 0:n , '=', c(X0,X))## [1] "Valor de la cadena al tiempo 0 = 0"
## [2] "Valor de la cadena al tiempo 1 = 1"
## [3] "Valor de la cadena al tiempo 2 = 0"
## [4] "Valor de la cadena al tiempo 3 = 1"
## [5] "Valor de la cadena al tiempo 4 = 0"
## [6] "Valor de la cadena al tiempo 5 = -1"
## [7] "Valor de la cadena al tiempo 6 = 0"
## [8] "Valor de la cadena al tiempo 7 = 1"
## [9] "Valor de la cadena al tiempo 8 = 0"
## [10] "Valor de la cadena al tiempo 9 = -1"
## [11] "Valor de la cadena al tiempo 10 = 0"## Otra forma más rápida.
set.seed(123)
n <- 10
p <- 0.54
X0 <- 0
tornillos <- 2*rbinom(n,1,p)-1
Xn <- c(X0, (X0 + cumsum(tornillos)))
paste('Valor de la cadena al tiempo', 0:n , '=', Xn)## [1] "Valor de la cadena al tiempo 0 = 0"
## [2] "Valor de la cadena al tiempo 1 = 1"
## [3] "Valor de la cadena al tiempo 2 = 0"
## [4] "Valor de la cadena al tiempo 3 = 1"
## [5] "Valor de la cadena al tiempo 4 = 0"
## [6] "Valor de la cadena al tiempo 5 = -1"
## [7] "Valor de la cadena al tiempo 6 = 0"
## [8] "Valor de la cadena al tiempo 7 = 1"
## [9] "Valor de la cadena al tiempo 8 = 0"
## [10] "Valor de la cadena al tiempo 9 = -1"
## [11] "Valor de la cadena al tiempo 10 = 0"Ahora un gráfico simple de la trayectoria.
## Simulación del valor de la cadena al tiempo n con probabilidad p de
## incrementar en uno el valor de la cadena después de un paso.
set.seed(123)
n <- 10
p <- 0.54
X0 <- 0
tornillos <- 2*rbinom(n,1,p)-1
X <- numeric(n)
X[1] <- X0 + tornillos[1]
for(i in 2:n){
X[i] <- X[i-1] + tornillos[i]
}
plot( 0:n , c(X0,X), type = 'b',col = 'red', lwd = 2, las = 2, ylab = '')
title(main = 'Caminata aleatoria simple')
abline(h=0,v=0)
## Otra forma más rápida.
set.seed(123)
n <- 10
p <- 0.54
X0 <- 0
tornillos <- 2*rbinom(n,1,p)-1
Xn <- c(X0, (X0 + cumsum(tornillos)))
plot( 0:n , Xn , type = 'o' ,col = 'blue', lwd = 2, las = 2, ylab = '')
title(main = 'Caminata aleatoria simple')
abline(h=0,v=0)
Ahora, vemos al mismo tiempo varias trayectorias simuladas. Primero en varias gráficas y después en una sola gráfica.
set.seed(123)
n <- 10
p <- 0.54
X0 <- 0
Xn <- matrix(0,(n+1),9)
for(i in 1:9){
tornillos <- 2*rbinom(n,1,p)-1
Xn[,i] <- c(X0, (X0 + cumsum(tornillos)))
}
par(mfrow = c(3,3),mar = rep(0, 4),cex.axis=0.8,
oma = rep(4.1, 4),bg = 'lightyellow')
miny <- min(Xn)
maxy <- max(Xn)
for(i in 1:9){
if(i<7){
plot( 0:n , Xn[,i] , type = 'o',
lwd = 2, axes = F,ylim = c(miny,maxy),
col=rainbow(9)[i],ylab = '', xlab = '')
abline(h=0,v=0)
box()
}else{
plot( 0:n , Xn[,i] , type = 'o',
lwd = 2, axes = F,ylim = c(miny,maxy),
col=rainbow(9)[i],ylab = '', xlab = '')
axis(1,0:n)
abline(h=0,v=0)
box()
}
}
mtext('Caminata aleatoria simple', side = 3, line = 3,font = 2, outer = TRUE)
par(mfrow=c(1,1), bg='white')
#############
#############
matplot(0:n,Xn, type = 'o',pch=16, lty=1,col = rainbow(9))
abline(h=miny:maxy,v=0:n,col = 'lightgray')
abline(h=0,v=0)
axis(2,c(miny,maxy))
mtext('Caminatas aleatorias simples', side = 3, line = 1,font = 3, outer = TRUE)
Por último, ahora simulamos las caminatas durante 50 pasos únicamente cambiando \(n\) a 50.
# set.seed(123) Sin semilla
n <- 50
p <- 0.54
X0 <- 0
Xn <- matrix(0,(n+1),9)
for(i in 1:9){
tornillos <- 2*rbinom(n,1,p)-1
Xn[,i] <- c(X0, (X0 + cumsum(tornillos)))
}
par(mfrow = c(3,3),mar = rep(0, 4),cex.axis=0.8,
oma = rep(4.1, 4),bg = 'lightyellow')
miny <- min(Xn)
maxy <- max(Xn)
for(i in 1:9){
if(i<7){
plot( 0:n , Xn[,i] , type = 'o',
lwd = 2, axes = F,ylim = c(miny,maxy),
col=rainbow(9)[i],ylab = '', xlab = '')
abline(h=0,v=0)
box()
}else{
plot( 0:n , Xn[,i] , type = 'o',
lwd = 2, axes = F,ylim = c(miny,maxy),
col=rainbow(9)[i],ylab = '', xlab = '')
axis(1,0:n)
abline(h=0,v=0)
box()
}
}
mtext('Caminata aleatoria simple', side = 3, line = 3,font = 2, outer = TRUE)
par(mfrow=c(1,1), bg='white')
#############
#############
matplot(0:n,Xn, type = 'o',pch=16, lty=1,col = rainbow(9))
abline(h=miny:maxy,v=0:n,col = 'lightgray')
abline(h=0,v=0)
axis(2,c(miny,maxy))
mtext('Caminatas aleatorias simples', side = 3, line = 1,font = 3, outer = TRUE)
No hay comentarios:
Publicar un comentario