23-Cadenas de Markov 5

Cadenas de Markov 5 de 10

David J. Santana

21 de octubre de 2016

Cadenas de Markov con espacio de estados infinito.

Ahora vamos a simular una cadena con espacio de estados infinito, con el fin de ver las probabilidades de regreso después de \(n\) pasos, a partir de un estado estado inicial determinado.

Cadena de rachas de éxitos y su simulación

La cadena de rachas de éxitos tiene espacio de estados \[\mathcal{S}=\left\{0,1,2,\ldots \right\},\] y probabilidades de transición dadas por \begin{eqnarray} p_{i,i+1}&=& p,\\ p_{i,0}&=& 1-p, \end{eqnarray}

para \(i = 0,1,2,\ldots\). Con el siguiente código podemos ver los estados que la cadena visita durante \(n\) pasos iniciando en el estado \(X_0\) y con parámetro \(p\).

trayectoria_rachas <- function(n,X0,p = 1/2){
        sigue_o_regresa <- rbinom(n,1,p) #1-sigue la racha, 0-vuelve
        trayectoria <- X0
        for(i in 1:n){
                trayectoria[i+1] <- (trayectoria[i]+1)*(sigue_o_regresa[i] == 1)
        }
        trayectoria
}

## Ejemplo

trayectoria_rachas(50,0)

##  [1] 0 0 0 1 2 3 4 5 0 0 1 2 0 1 0 1 2 3 4 5 6 7 0 1 0 0 0 1 2 3 0 1 2 3 0
## [36] 0 1 2 0 1 0 1 2 0 1 2 3 4 0 1 2

Si además queremos ver una gráfica de la trayectoria, entonces podemos usar el siguiente código.

## Ejemplo 2
set.seed(10232016)
path <- trayectoria_rachas(50,0,0.65)
plot(0:50,path,type = 'o',main = 'Trayectoria de la cadena de Rachas de éxitos',
     col = 2, pch = 16, las = 2)
abline(h = 0, v = 0)
abline(h = seq_along(path), col = 'lightgray', lty = 2)

## Ejemplo 3
path2 <- trayectoria_rachas(500,0,0.75)
plot(0:500,path2,type = 'l',main = 'Trayectoria de la cadena de Rachas de éxitos',
     col = 2, pch = 16, las = 2)
abline(h = 0, v = 0)
abline(h = seq_along(path2), col = 'lightgray', lty = 2)

Ahora nos interesa hacer simulaciones para aproximar la probabilidad de ir del estado \(i\) al estado \(i\) en \(n\) pasos, \(p_{i,i}(n)\), y la probabilidad de ir por primera vez al estado \(i\) justo en el \(n\)-ésimo paso, \(f_{i,i}(n)\), comenzando en \(i\). Teóricamente \[p_{0,0}(n)=1-p\] y \[f_{0,0}(n)=(1-p)p^{n-1}\] para el estado \(i = 0\) pero es dífícil calcular estas probabilidades para otros estados. Con simulación podemos calcular fácilmente aproximaciones para ambas cantidades.

## Función para calcular p_{ii}(n)

piin <- function(i,n,p, sim = 20000){
        pii_n <- numeric(sim)
        for(k in 1:sim){
                sigue_o_regresa <- rbinom(n,1,p)
                aux <- i
                for(j in 1:n){
                        aux <- (aux+1)*(sigue_o_regresa[j] == 1)
                }
        pii_n[k] <- 1*(aux == i)
        }
        mean(pii_n)
}

## Ejemplo 3
set.seed(123)
piin(0,7,0.6) # P(X_7 = 0 / X_0 = 0) = p_{0,0}(7) = 1 - p = 0.4

## [1] 0.39695

piin(0,20,0.6) # P(X_20 = 0 / X_0 = 0) = p_{0,0}(20) = 1 - p = 0.4

## [1] 0.39815

piin(0,50,0.6) # P(X_50 = 0 / X_0 = 0) = p_{0,0}(50) = 1 - p = 0.4

## [1] 0.40265

piin(1,7,0.6) # P(X_7 = 1 / X_0 = 1) = p_{1,1}(7)

## [1] 0.23865

piin(2,7,0.6) # P(X_7 = 2 / X_0 = 2) = p_{2,2}(7)

## [1] 0.1422

## Función para calcular f_{ii}(n)

fiin <- function(i,n,p, sim = 20000){
        fii_n <- numeric(sim)
        for(k in 1:sim){
                sigue_o_regresa <- rbinom(n,1,p)
                aux <- (i+1)*(sigue_o_regresa[1] == 1)
                j <- 1
                while(aux!=i && j<n){
                        aux <- (aux+1)*(sigue_o_regresa[j+1] == 1)
                        j <- j+1
                }
                if(j < n){exito <- 0}else{exito <- 1*(aux == i)}
                fii_n[k] <- exito
        }
        mean(fii_n)
}

## Ejemplo 4
set.seed(2016)
fiin(0,7,0.6) # P(X_7 = 0, X_j != 0; j=1,...,6 / X_0 = 0) = f_{0,0}(7) = 0.0187

## [1] 0.0185

fiin(1,7,0.6) # P(X_7 = 1, X_j != 1; j=1,...,6 / X_0 = 1) = f_{1,1}(7)

## [1] 0.0536

fiin(2,7,0.6) # P(X_7 = 2, X_j != 2; j=1,...,6 / X_0 = 2) = f_{2,2}(7)

## [1] 0.1013