En este capítulo estudiaremos sistemas estocásticos que serán modelados haciendo uso de Cadenas de Markov

Proceso de Markov. Proceso aleatorio donde el valor de la variable aleatoria en el instante \(n\) depende solamente de su valor pasado inmediato en el instante \(n-1\)

En un proceso de Markov, la variable aleatoria representa el estado del sistema en un instante dado \(n\)

Los procesos de Markov son ejemplos de procesos que salen en situaciones de la vida real:

• Protocolos de telecomunicación y sistemas de maquinaria
• Llegadas y salidas de clientes en bancos
• Pasar por caja en los supermercados
• Mutación de un virus o molécula de ADN
• Paseo aleatorio como por ejemplo el movimiento Browniano
• Llegada de coches a una intersección
• El estado del tiempo diario

Si el espacio de estados (espacio muestral) de un proceso de Markov es discreto, el proceso de Markov se denomina Cadena de Markov.

En este caso, los estados se etiquetan con los números \(0,1,2,\dots\)

En este tema, estudiaremos cadenas de Markov de tiempo discreto, ya que son las que salen en la mayoría de sistemas de comunicación.

Tiempo de espera. En una cadena de Markov, es el tiempo que estamos en un estado determinado

Dependiendo de como midamos este tiempo de espera, tenemos dos tipos de cadenas de Markov:

• Cadenas de Markov de tiempo discreto. Los tiempos de espera toman valores enteros. Por tanto, los cambios de estados pasan en tiempos discretos \(t = T_0,T_1,T_2,\dots\). Si estos tiempos están equiespaciados, tenemos \(T_n = nT\), con \(n=0,1,2,\dots\)
• Cadenas de Markov de tiempo continuo. En este caso, los cambios de estados pueden tener lugar en cualquier instante de tiempo

En una cadena de Markov de tiempo discreto el valor de la variable aleatoria \(S(n)\) representa el estado en el tiempo \(n\)

La variable aleatoria \(S(n)\) solo depende de la variable aleatoria \(S(n-1)\). Esta propiedad se denomina falta de memoria de las cadenas de Markov.

La probabilidad de que \(S(n)\) tome el valor \(s_i\) solo depende de los valores \(s_j\) que toma \(S(n-1)\). Matemáticamente

\[p\{S(n) = s_i\} = f(s_j)\] para todos \(s_i,s_j\in S\) donde \(S\) es el conjunto de estados del sistema (el espacio muestral)

Definimos \(p_{ij}(n)\) la probabilidad de que en el tiempo actual \(n\), nuestro sistema se encuentre en el estado \(i\) suponiendo que en el tiempo anterior estábamos en el estado \(j\) (tiempo \(n-1\)).

Matemáticamente

\[p_{ij}(n) = p\{S(n) = i\ |\ S(n-1) = j\}\]

En el caso de que la probabilidad anterior \(p_{ij}(n)\) sea independiente del tiempo \(n\), diremos que la cadena de Markov es homogénea.

En este caso, podemos escribir

\[p_{ij} = p\{S(n) = i\ |\ S(n-1) = j\}\]

Definimos \(s_i(n)\) la probabilidad de que en el tiempo \(n\) estemos en el estado \(i\)

\[s_i(n) = p\{S(n) = i\}\]

Haciendo uso del Teorema de las probabilidades totales y suponiendo que la cadena de Markov es homogénea, podemos escribir

\[s_i(n) = \sum_{j = 1}^mp_{ij}s_j(n-1)\] suponiendo que los estados de nuestra cadena de Markov son \(\{1,2,\dots,m\}\)

Matricialmente, podemos escribir la expresión anterior como \[s(n) = P\cdot s(n-1)\]

donde \(P\) es la llamada matriz de transición de la cadena de Markov homogénea:

\[P = \begin{pmatrix} p_{11} & p_{12} & \cdots & p_{1m}\\ p_{21} & p_{22} & \cdots & p_{2m}\\ \vdots & \vdots & \ddots & \vdots\\ p_{m1} & p_{m2} & \cdots & p_{mm} \end{pmatrix}\]

Definiremos \(s(n)\) el vector de probabilidades de la cadena de markov en el tiempo \(n\):

\[s(n) = (s_n(1), s_n(2),\dots,s_n(m))^t\]

Como que \(s(n)\) es un vector de probabilidades, se debe verificar

\[\sum_{i = 1}^ms_n(i) = 1\]

El estudio de la matriz de transición \(P\) y más concretamente sus valores y vectores propios, nos determinará el comportamiento de nuestra cadena de Markov

Cada columna \(j-\)ésima de la matriz de transición \(P\) representa las probabilidades de ir hacia cada uno de los estados \(i\) desde el estado \(j\), para \(i = 1,\dots, m\). Por tanto, \[\sum_{i = 1}^mp_{ij} = 1\]

Las matrices que cumplen que la suma de sus columnas valen 1, se llaman matrices estocásticas

Matriz estocástica. Diremos que una matriz \(P\) es estocástica o una matriz de Markov si verifica las siguientes propiedades:

• La matriz \(P\) es cuadrada
• Los elementos \(p_{ij}\) de \(P\) están entre 0 y 1 (\(0\le p_{ij}\le 1\))
• La suma de cada columna de \(P\) vale 1: si \(P\) es una matriz de orden \(m\), entonces \[\sum_{i = 1}^mp_{ij} = 1\qquad \forall j = 1,\dots,m\]

A continuación veremos que estas matrices verifican

• Tienen el valor propio 1
• Todos los valores propios \(\lambda_i\) de \(P\) verifican \(|\lambda_i|\le 1\)

Teorema. Sea \(P\) una matriz estocástica. Entonces tiene el valor 1 como valor propio

Teorema. Sea \(A\) una matriz de dimensiones \(m\times n\). Sea \(P\) una matriz estocástica de orden \(m\). Entonces, la suma de cada columna de la matriz \(A\) y de la matriz \(B = PA\) es la misma. Es decir, si \(\sigma_j(M)\) simboliza la suma de la columna \(j-\)ésima de una matriz \(M\), tenemos que \(\sigma_j(A) = \sigma_j(B)\quad \forall j=1,\dots,n\)

En pocas palabras, las sumas de las columnas de una matriz cualquiera no cambian cuando la multiplicamos por la izquierda por una matriz estocástica

Vayamos a interpretar las diagonales de la matriz \(P\) estudiando un caso particular de cadenas de Markov: las colas

En una cola, los clientes llegan y se van y los estados de la cadena de Markov son el número de clientes que hay en la cola en cada momento

Una cola se caracteriza por el número medio de clientes que llegan por unidad de tiempo, el número de servidores, el tamaño de la cola y el número medio de clientes que abandonan la cola por unidad de tiempo.

Las diagonales de \(P\) reflejan las características de la cola

Teorema. Sea \(P\) una matriz estocástica de orden \(m\). Entonces,

• Cualquier valor propio \(\lambda\) de \(P\) verifica \(|\lambda|\le 1\)
• Si las componentes de la matriz \(P\) son estrictamente positivas (\(p_{ij}>0\)) y \(\lambda\ne 1\), entonces \(|\lambda|<1\)

Teorema. Sea \(P\) una matriz estocástica y regular. Si \(P\) es diagonalizable y todas sus componentes son estrictamente positivas, \(p_{ij}>0\), entonces el valor propio 1 tiene multiplicidad 1.

Teorema. Sea \(P\) una matriz estocástica. Sea \(x\) un vector propio del valor propio \(\lambda\ne 1\). Entonces la suma de las componentes de \(x\) es 0: \(\sigma(x) =0\)

Teorema. Sea \(P\) una matriz estocástica y diagonalizable de orden \(m\). Sea \(x\) un vector propio de valor propio 1. Entonces, la suma de las componentes de \(x\) no es 0: \(\sigma(x) \ne0\)

Corolario. Sea \(P\) una matriz estocástica y diagonalizable. Siempre podemos escoger un vector propio \(x\) de valor propio 1 tal que la suma de sus componentes sea \(1\): \(\sigma(x) = 1\)

Calculamos los valores y vectores propios

P = matrix(c(0.3, 0.1, 0, 0.2,
             0.1, 0.6, 0.3, 0.1,
             0.4, 0.2, 0.4, 0.5,
             0.2, 0.1, 0.3, 0.2),
           nrow = 4, byrow = T)
round(eigen(P)$values, 4) #Vaps

[1]  1.0000  0.3655  0.2211 -0.0866

round(eigen(P)$vectors, 4) #Veps

        [,1]    [,2]    [,3]    [,4]
[1,] -0.1965  0.1002 -0.5887 -0.4286
[2,] -0.6309  0.7720 -0.3983  0.1644
[3,] -0.6516 -0.5190  0.5555 -0.4821
[4,] -0.3724 -0.3532  0.4315  0.7463

round(apply(eigen(P)$vectors, MARGIN = 2, sum), 4)

[1] -1.8515  0.0000  0.0000  0.0000

Para construir la matriz de transición \(P\), se deben seguir los siguientes pasos:

• Verificar que nuestro sistema verifica las propiedades de Markov
• Todos los estados del sistema deben de estar identificados y etiquetados
• Todas las transiciones entre los estados deben de estar indicadas (diagrama de transición)
• Se deben calcular las probabilidades entre todas las transiciones
• Al final, tenemos que construir la matriz de transición
• Podemos reetiquetar los estados de forma que la matriz de transición sea más simple

Recordemos que \(s_n\) representaba el vector de probabilidades de la cadena de Markov en el instante \(n\)

Por tanto, el vector \(s_1\) será en función del vector incial de probabilidades \(s_0\): \(s_1 = Ps_0\)

Del mismo modo, \(s_2 = Ps_1 = P^2s_0\) y, en general, tendremos que el vector de probabilidades en el instante \(n\) será

\[s_n = P^ns_0\qquad n = 1,2,\dots\]

\(P^n\) es una matriz estocástica para cualquier valor de \(n\) ya que hemos visto que multiplicar por una matriz estocástica no cambia la suma de columnas de una matriz

Un elemento de \(P\) nunca puede ser negativo en \(P^n\)

Recordemos que el vector de probabilidades \(s_n\) en el instante \(n\) vale en función del vector de probabilidades inicial \(s_0\)

\[s_n = P^ns_0\]

Para calcular el vector de probabilidades, necesitamos técnicas para calcular \(P^n\). Estas pueden ser

• Multiplicaciones sucesivas de \(P\). Muy costoso computacionalmente
• Expandir el vector de probabilidades inicial \(s_0\) en una base adecuada
• Diagonalizar (si es posible) la matriz \(P\)
• Encontrar la forma normal de Jordan de la matriz \(P\)
• Hacer uso de la transformada-z

Sea \(P\) la matriz de transición de orden \(m\). Supongamos que \(P\) es diagonalizable. Esto quiere decir que existe una base de vectores propios de \(P\), \(x_1,\dots,x_m\) asociados a valores propios \(\lambda_1,\dots,\lambda_m\)

En este caso, consideremos la expansión del vector de probabilidades inicial \(s_0\) en la base de vectores propios \(x_i,\quad i=1,\dots,m\)

\[s_0 = c_1x_1+c_2x_2+\cdots+c_mx_m\]

Definimos la matriz \(X\) de columnas los vectores propios \(x_i\quad i = 1,\dots,m\)

\[X = [x_1,x_2,\dots,x_m]\]

Definimos \(c\) el vector de coeficientes de la expansión de \(s_0\)

\[c = (c_1,c_2,\dots,c_m)^t\]

De cara al cálculo de \(s_n\), tenemos que tener en cuenta que

\[Px_i = \lambda_ix_i\]

Así tendremos que

\[s_1 = Ps_0\] \[=P(c_1x_1+c_2x_2+\cdots+c_mx_m)\] \[=c_1Px_1+c_2Px_2+\cdots+c_mPx_m\] \[=c_1\lambda_1x_1+c_2\lambda_2x_2+\cdots+c_m\lambda_mx_m\]

y, en general

\[s_n = P^ns_0\] \[=P^n(c_1x_1+c_2x_2+\cdots+c_mx_m)\] \[=c_1P^nx_1+c_2P^nx_2+\cdots+c_mP^nx_m\] \[=c_1\lambda_1^nx_1+c_2\lambda_2^nx_2+\cdots+c_m\lambda_m^nx_m\]

Si llamamos a la matriz \(D\) la matriz diagonal formada por los vectores propios

\[D = \begin{pmatrix} \lambda_1 & 0 & \cdots & 0\\ 0 & \lambda_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \lambda_m \end{pmatrix}\]

podemos escribir las expansiones matricialmente como

\[s_1 = XDc\] \[s_n = XD^nc\]

donde \(D^n\) es la matriz

\[D^n = \begin{pmatrix} \lambda_1^n & 0 & \cdots & 0\\ 0 & \lambda_2^n & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \lambda_m^n \end{pmatrix}\]

P = matrix(c(0.3, rep(0.2, 3), 0.3, 0.4, rep(0.2, 8), rep(0.4, 2)),
           nrow = 4, byrow = T)
round(eigen(P)$values, 3) #Vaps

[1] 1.0 0.2 0.1 0.0

round(eigen(P)$vectors, 3) #Veps

       [,1]   [,2]   [,3]   [,4]
[1,] -0.439  0.000 -0.707  0.000
[2,] -0.548 -0.707  0.707  0.000
[3,] -0.395  0.000  0.000  0.707
[4,] -0.592  0.707  0.000 -0.707

X = round(eigen(P)$vectors, 3)
X

       [,1]   [,2]   [,3]   [,4]
[1,] -0.439  0.000 -0.707  0.000
[2,] -0.548 -0.707  0.707  0.000
[3,] -0.395  0.000  0.000  0.707
[4,] -0.592  0.707  0.000 -0.707

s0 = c(1,0,0,0)
c = round(solve(X)%*%s0, 3)
c

       [,1]
[1,] -0.507
[2,] -0.707
[3,] -1.100
[4,] -0.283

D = diag(round(eigen(P)$values, 3))
s1 = round(X%*%D%*%c, 3)
s1

     [,1]
[1,]  0.3
[2,]  0.3
[3,]  0.2
[4,]  0.2

s2 = round(X%*%D^2%*%c, 3)
s2

     [,1]
[1,] 0.23
[2,] 0.29
[3,] 0.20
[4,] 0.28

s5 = round(X%*%D^5%*%c, 3)
s5

      [,1]
[1,] 0.223
[2,] 0.278
[3,] 0.200
[4,] 0.300

s20 = round(X%*%D^(20)%*%c, 3)
s20

      [,1]
[1,] 0.223
[2,] 0.278
[3,] 0.200
[4,] 0.300

s_n = function(X, D, s0, n){
  c = solve(X)%*%s0
  sn = X%*%D^n%*%c
  return(round(sn, 3))
}
s20 = s_n(X,D,s0,20) #Comprobamos que funciona
s20

      [,1]
[1,] 0.222
[2,] 0.278
[3,] 0.200
[4,] 0.300

P = matrix(c(0.1, 0.4, 0.2,
             0.1, 0.4, 0.6,
             0.8, 0.2, 0.2),
           nrow = 3, byrow = T)
round(eigen(P)$values, 3)

[1]  1.00+0.000i -0.15+0.312i -0.15-0.312i

round(eigen(P)$vectors, 3)

          [,1]          [,2]          [,3]
[1,] -0.432+0i -0.170+0.354i -0.170-0.354i
[2,] -0.673+0i -0.510-0.354i -0.510+0.354i
[3,] -0.601+0i  0.679+0.000i  0.679+0.000i

X = round(eigen(P)$vectors, 3)
X

          [,1]          [,2]          [,3]
[1,] -0.432+0i -0.170+0.354i -0.170-0.354i
[2,] -0.673+0i -0.510-0.354i -0.510+0.354i
[3,] -0.601+0i  0.679+0.000i  0.679+0.000i

s0 = c(1,0,0)
c = round(solve(X)%*%s0, 3)
c

             [,1]
[1,] -0.586+0.00i
[2,] -0.259-0.93i
[3,] -0.259+0.93i

D = diag(round(eigen(P)$values, 3))
s1 = s_n(X, D, s0, 1)
s1

         [,1]
[1,] 0.100+0i
[2,] 0.100+0i
[3,] 0.799+0i

s5 = s_n(X, D, s0, 5)
s5

         [,1]
[1,] 0.251+0i
[2,] 0.399+0i
[3,] 0.350+0i

Si la matriz de transición \(P\) es diagonalizable, significa que existen matrices \(X\) de vectores propios y una matriz diagonal \(D\) de valores propios tales que

\[P = XDX^{-1}\]

Para calcular el vector de probabilidades \(s_n\) en el instante \(n\), debemos hacer

\[s_n = P^ns_0 = XD^nX^{-1}s_0\]

ya que \(P^n = XD^nX^{-1}\)

Nos damos cuenta de que este método es equivalente al método anterior de expandir el vector de probabilidades inicial \(s_0\), donde

\[c = X^{-1}s_0\]

Estudiemos el comportamiento de la matriz \(X^{-1}\)

Lema. Supongamos que la primera columna de la matriz de vectores propios \(X\) corresponde al vector propio normalizado (suma 1) de valor propio 1. Además, supongamos que este valor propio tiene multiplicidad 1. Entonces, la matriz \(X^{-1}\) tiene unos en su primera fila. O sea,

\[X^{-1} = \begin{pmatrix} 1 & 1 & \cdots & 1\\ \vdots & \vdots & \ddots & \vdots\end{pmatrix}\]

Recordemos que si \(P\) es diagonalizable, existen matrices \(X\) de vectores propios y una matriz diagonal \(D\) de valores propios tales que

\[P = XDX^{-1}\]

Definimos la matriz \(A_j\) de la forma siguiente: Primero, consideramos la matriz \(1_{jj}\), que vale 0 para todo elemento excepto en \((j,j)\), que vale 1. Después definimos \(A_j\) como el producto de la matriz \(X\) por la matriz anterior por la inversa de \(X\)

\[A_j = X\begin{pmatrix} 0 & \cdots & 0 & \cdots & 0\\ \vdots & \ddots & \vdots & \ddots & \vdots\\ 0 & \cdots & 1 & \cdots & 0\\ \vdots & \ddots & \vdots & \ddots & \vdots\\ 0 & \cdots & 0 & \cdots & 0 \end{pmatrix}X^{-1}:= X\cdot1_{jj}\cdot X^{-1}\]

Supongamos que los valores propios de \(P\) son \[\lambda_1 = 1,\lambda_2,\dots,\lambda_m\]

Escribimos la matriz diagonal \(D\) como

\[D = \lambda_1 1_{11}+\lambda_2 1_{22}+\cdots +\lambda_{m}1_{mm}\]

Sustituyendo la expresión anterior en \(P=XDX^{-1}\), nos queda

\[P=\lambda_1X\cdot 1_{11}\cdot X^{-1}+\cdots \lambda_m X\cdot 1_{mm}\cdot X^{-1}\] \[=\lambda_1A_1+\lambda_2A_2 +\cdots +\lambda_mA_m\] \[=A_1+\lambda_2A_2+\cdots+\lambda_mA_m\]

Del mismo modo, como \(P^n = XD^nX^{-1}\) nos queda

\[P^n=\lambda_1^nX\cdot 1_{11}\cdot X^{-1}+\cdots \lambda_m^n X\cdot 1_{mm}\cdot X^{-1}\] \[=\lambda_1^nA_1+\lambda_2^nA_2 +\cdots +\lambda_m^nA_m\]

\[=A_1+\lambda_2^nA_2+\cdots+\lambda_m^nA_m\]

Estudiemos la matriz \(A_1\)

Teorema. La matriz \(A_1\) es estocástica (la suma de sus columnas vale 1) y todas sus columnas son idénticas

En el caso de que existan, definimos el vector de probabilidades en el estado del infinito \(s_{\infty}\) como

\[s_{\infty} = \lim_{n\to\infty}s_n\]

El siguiente teorema nos dice como calcular \(s_{\infty}\)

Teorema. Consideremos una cadena de Markov donde la matriz de transición \(P\) tiene el valor propio 1 con multiplicidad 1. Entonces, existe el vector de probabilidades en el estado de equilibrio, es independiente del vector de probabilidades \(s_0\) y es igual a cualquier columna de la matriz \(A_1\)

Teorema. Las matrices \(A_i\) para \(i\ne 1\) son diferenciales. O sea, la suma de sus columnas valen 0

Teorema. Dada una matriz \(A\) y una matriz diferencial \(V\) (la suma de sus columnas vale 0), la matriz \(B = VA\) será una matriz diferencial

P = rbind(c(0.1, 0.3, 1),
          c(0.2, 0.3, 0),
          c(0.7, 0.4, 0))
round(eigen(P)$values, 3) #Vaps

[1]  1.000 -0.769  0.169

round(eigen(P)$vectors, 3) #Veps

       [,1]   [,2]   [,3]
[1,] -0.757  0.768  0.526
[2,] -0.216 -0.144 -0.804
[3,] -0.616 -0.624  0.277

X = round(eigen(P)$vectors, 3)

A1 = round(X%*%diag(c(1,0,0))%*%solve(X), 3)
A1

      [,1]  [,2]  [,3]
[1,] 0.476 0.476 0.477
[2,] 0.136 0.136 0.136
[3,] 0.388 0.387 0.388

A2 = round(X%*%diag(c(0,1,0))%*%solve(X), 3)
A2

       [,1]   [,2]   [,3]
[1,]  0.495  0.102 -0.645
[2,] -0.093 -0.019  0.121
[3,] -0.403 -0.083  0.524

A3 = round(X%*%diag(c(0,0,1))%*%solve(X), 3)
A3

       [,1]   [,2]   [,3]
[1,]  0.028 -0.578  0.168
[2,] -0.043  0.883 -0.257
[3,]  0.015 -0.304  0.088

P = rbind(c(0.2, 0.4, 0.2),
          c(0.1, 0.4, 0.6),
          c(0.7, 0.2, 0.2))
round(eigen(P)$values, 3) #Vaps

[1]  1.0+0.0i -0.1+0.3i -0.1-0.3i

round(eigen(P)$vectors, 3) #Veps

         [,1]          [,2]          [,3]
[1,] 0.475+0i -0.129+0.387i -0.129-0.387i
[2,] 0.660+0i -0.516-0.387i -0.516+0.387i
[3,] 0.581+0i  0.645+0.000i  0.645+0.000i

X = round(eigen(P)$vectors, 3)

A1 = round(X%*%diag(c(1,0,0))%*%solve(X), 3)
A1

         [,1]     [,2]     [,3]
[1,] 0.277+0i 0.277+0i 0.277+0i
[2,] 0.385+0i 0.385+0i 0.385+0i
[3,] 0.339+0i 0.339+0i 0.339+0i

A2 = round(X%*%diag(c(0,1,0))%*%solve(X), 3)
A2

              [,1]          [,2]          [,3]
[1,]  0.362+0.008i -0.138-0.159i -0.138+0.174i
[2,] -0.192+0.539i  0.308-0.128i -0.192-0.295i
[3,] -0.169-0.546i -0.169+0.287i  0.331+0.120i

A3 = round(X%*%diag(c(0,0,1))%*%solve(X), 3)
A3

              [,1]          [,2]          [,3]
[1,]  0.362-0.008i -0.138+0.159i -0.138-0.174i
[2,] -0.192-0.539i  0.308+0.128i -0.192+0.295i
[3,] -0.169+0.546i -0.169-0.287i  0.331-0.120i

Algunas veces, si reetiquetamos los estados, la matriz de transición se simplifica en el sentido de que tiene una estructura más sencilla. Veamos un ejemplo

En este último capítulo estudiaremos el comportamiento del vector de probabilidades \(s_n\) en el infinito. O sea,

\[s_{\infty} = \lim_{n\to \infty}s_n\]

Este vector de probabilidades nos permitirá estudiar el rendimiento de la cadena de Markov.

Anteriormente vimos que si nuestra cadena de Markov es homogénea, el vector \(s_{\infty}\) es un vector propio de valor propio 1 de la matriz de transición \[Ps_{\infty} = s_{\infty}\]

Si el valor propio 1 tiene multiplicidad 1, este vector es único

Vector de probabilidades en el estado de equilibrio. Es el vector \(s_{\infty}\) de la cadena de Markov

Cuando estamos en el estado de equilibrio, quiere decir que la cadena de Markov ya no cambia en el tiempo. Se mantiene estable

La matriz de transición \(P^n\), cuando \(n\to\infty\) también llega a alcanzar el estado de equilibrio

Llamaremos \(s\) al vector de probabilidades en el estado de equilibrio

Estudiaremos un conjunto de técnicas para calcular \(s\). Estas técnicas dependerán de la estructura de la matriz de transición \(P\)

P = matrix(c(0.2, 0.4, 0.5,
             0.8, 0, 0.5,
             0, 0.6, 0), 
           ncol = 3, byrow = T)
round(P%^%2, 3)

     [,1] [,2] [,3]
[1,] 0.36 0.38  0.3
[2,] 0.16 0.62  0.4
[3,] 0.48 0.00  0.3

round(P%^%10, 3)

      [,1]  [,2]  [,3]
[1,] 0.353 0.354 0.354
[2,] 0.404 0.404 0.404
[3,] 0.242 0.243 0.242

round(P%^%20, 3)

      [,1]  [,2]  [,3]
[1,] 0.354 0.354 0.354
[2,] 0.404 0.404 0.404
[3,] 0.242 0.242 0.242

1. Multiplicaciones repetidas de \(P\) para calcular \(P^n\) para valores grandes de \(n\)
2. Cálculo del vector propio de valor propio 1 de la matriz \(P\)
3. Ecuaciones en diferencias
4. Transformada \(Z\) (funciones generatrices)
5. Métodos directos para resolver un sistema lineal de ecuaciones
6. Métodos iterativos para resolver un sistema lineal de ecuaciones
7. Técnicas iterativas para poner los estados de \(P\) en función de los otros estados

¿Qué técnica es la más adecuada? Depende de la estructura de la matriz \(P\)

• La técnica de multiplicaciones sucesivas de \(P\) es muy sensible a errores numéricos
• La técnica 2 es adecuada cuando el tamaño de la matriz \(P\) no es demasiado grande y \(P\) está expresada numéricamente
• La técnica de ecuaciones en diferencia es adecuada cuando \(P\) es una matriz banda con pocas subdiagonales
• La técnica 4 es adecuada cuando \(P\) es triangular inferior o Hessengberg inferior donde cada subdiagonal tiene los mismos elementos
• Las técnicas 5 y 6 son adecuadas cuando \(P\) se expresa numéricamente y no tiene ninguna estructura en particular

Para calcular \(s\), tenemos que resolver el sistema lineal indeterminado siguiente

\[Ps = s\]

Una vez tengamos todas las soluciones del sistema anterior, imponemos que \(\sum_{i = 1}^ms_i = 1\) y encontramos el único \(s\) que verifica el sistema anterior y que la suma de sus componentes es 1

P = rbind(c(0.8, 0.7, 0.5),
          c(0.15, 0.2, 0.3),
          c(0.05, 0.1, 0.2))
round(eigen(P)$values, 3)

[1] 1.0 0.2 0.0

round(eigen(P)$vectors, 3)

      [,1]   [,2]   [,3]
[1,] 0.973  0.816  0.535
[2,] 0.215 -0.408 -0.802
[3,] 0.088 -0.408  0.267

X = round(eigen(P)$vectors, 3)
sum(X[,1])

[1] 1.276

Esta técnica es útil cuando la matriz de transición es una matriz banda

Esta técnica es útil si la matriz de transición es Hessenberg superior y los elementos en una diagonal son idénticos

Este tipo de cadenas de Markov tienen lugar cuando modelamos colas con múltiples llegadas y una salida por unidad de tiempo

Más concretamente, supongamos que en nuestra cadena de Markov podemos ir desde un estado \(j\) a cualquier otro estado \(i\) con \(i\ge j-1\)

Esto correspondería a un buffer o cola, donde su tamaño solamente podría disminuir en un paquete pero podría aumentar en uno o más paquetes en cada unidad de tiempo

Supongamos que la probabilidad de que lleguen \(k\) paquetes vale \(a_k\)

\[p\{k\text{ llegadas}\} = a_k,\qquad k=0,1,2,\dots\]

Supongamos que la probabilidad de que un paquete abandone la cola vale \(c\)

La condición para que la cola sea estable es

\[\sum_{k=1}^{\infty}ka_k<c\]

lo que indica que el número medio de llegadas es inferior al número medio de salidas por unidad de tiempo

La matriz de transición será

\[P = \begin{pmatrix} a_0 & b_0 & 0 & 0 & \cdots\\ a_1 & b_1 & b_0 & 0 & \cdots\\ a_2 & b_2 & b_1 & b_0 & \cdots\\ a_3 & b_3 & b_2 & b_1 & \cdots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix}\]

donde \(b_i = a_ic+a_{i-1}d\) y \(d=1-c\) es la probabilidad de que un paquete no abandone la cola por unidad de tiempo. Supongamos que \(a_i = 0\) para \(i<0\)

El vector de probabilidades en el equilibrio cumple \(Ps = s\). En nuestro caso, esto es equivalente a resolver el siguiente sistema de ecuaciones

\[s_i = a_is_0+\sum_{j=1}^{i+1}b_{i-j+1}s_j\quad i=0,1,2,\dots\]

De cara a resolver el sistema lineal anterior, definimos las funciones siguientes

\[A(z) = \sum_{i=0}^{\infty}a_iz^{-i}\qquad B(z)=\sum_{i=0}^{\infty}b_iz^{-i}\]

Definimos la transformada \(Z\) del vector de probabilidades en el equilibrio \(s\) como

\[S(z) = \sum_{i = 0}^{\infty}s_iz^{-i}\]

donde \(s_i\) son las componentes de \(s\)

\[s = (s_0,s_1,s_2,\dots)\]

El objetivo es calcular \(S(z)\) y a partir de esta función, calcular \(s\)

Haciendo uso del sistema de ecuaciones lineal que verificaban las componentes \(s_i\), tenemos

\[S(z) = s_0\sum_{i=0}^{\infty}a_iz^{-i}+\sum_{i=0}^{\infty}\sum_{j=1}^{i+1}b_{i-j+1}s_jz^{-i}\] \[=s_{0}A(z)+\sum_{j= 1}^{\infty}s_j\sum_{i=j-1}^{\infty}b_{i-j+1}z^{-i}\]

donde hemos cambiado el orden de sumación en la última suma

Si continuamos haciendo cálculos

\[S(z)=s_0A(z)+\sum_{j\ge1}s_jz^{-j}\sum_{i\ge j-1}b_{i-j+1}z^{-(i-j)} = s_0A(z)+\sum_{j\ge 1}s_jz^{-j}\sum_{i\ge 0}b_iz^{-(i-1)}\] \[=s_0A(z)+(S(z)-s_0)zB(z)\]

de donde despejando \(S(z)\) tenemos

\[S(z) = \frac{s_0(z^{-1}A(z)-B(z))}{z^{-1}-B(z)}\]

Los valores \(s_i\) salen de desarrollar la expresión anterior como potencias de \(z^{-1}\)

\[S(z) = \frac{s_0(z^{-1}A(z)-B(z))}{-b_0+z^{-1}\left(1-\sum_{i\ge0}b_{i+1}z^{-i}\right)} = \frac{s_0(z^{-1}A(z)-B(z))}{-b_0\left(1-\left(\frac{z^{-1}}{b_0}\right)\left(1-\sum_{i\ge 0}b_{i+1}z^{-i}\right)\right)}\]

Sean ahora las funciones

\[C(z) = \left(\frac{s_0}{-b_0}\right)(z^{-1}A(z)-B(z))\] \[D(z) = \left(\frac{z^{-1}}{b_0}\right)\left(1-\sum_{i\ge 0}b_{i+1}z^{-i}\right)\]

La función \(S(z)\) vale

\[S(z) = \frac{C(z)}{1-D(z)}\]

Desarrollando dicha expresión, tendremos

\[S(z) = C(z)(1+D(z)+D(z)^2+\cdots)\]

De esta forma iremos encontrando los coeficientes \(s_i\) del vector de probabilidades \(s\)

Calculemos el valor \(s_0\). Sustituyendo en la expresión inicial por \(S(1)\), obtenemos

\[S(1) = 1= \frac{s_0(A(1)-B(1))}{1-B(1)}\]

Ahora bien, como \(A(1)=B(1)=1\), tenemos una expresión indeterminada de la forma \(\frac{0}{0}\)

Para resolverla, aplicamos la regla del Hôpital

\[S(1) = 1 = \lim_{z\to 1}S(z) = \lim_{z\to 1}\frac{s_0(z^{-1}A(z)-B(z))}{z^{-1}-B(z)}\] \[=\lim_{z\to 1}\frac{s_0(-z^{-2}A(z)+z^{-1}A'(z)-B'(z))}{-z^{-2}-B'(z)} = \frac{s_0(-A(1)+A'(1)-B'(1))}{-1-B'(1)}\] \[= \frac{s_0(-1+A'(1)-B'(1))}{-1-B'(1)}\]

de donde deducimos que

\[s_0 = \frac{1+B'(1)}{1+B'(1)-A'(1)}\]

Observación. El término \(-A'(1)\) representa el número medio de paquetes que llegan a la cola por unidad de tiempo cuando la cola está vacía y \(-B'(1)\) representa el número medio de paquetes que llegan a la cola cuando ésta no está vacía ya que

\[-A'(1) = \sum_{i\ge 1}ia_i<0\qquad -B'(1)=\sum_{i\ge1}ib_i<0\]

Esta técnica se utiliza cuando la matriz de transición \(P\) es Heissemberg inferior y los elementos en cada diagonal no son iguales. Por ejemplo, para \(m =6\), \(P\) vale

\[P = \begin{pmatrix} h_{11} & h_{12} & 0 & 0 & 0 & 0\\ h_{21} & h_{22} & h_{23} & 0 & 0 & 0\\ h_{31} & h_{32} & h_{33} & h_{34} & 0 & 0\\ h_{41} & h_{42} & h_{43} & h_{44} & h_{45} & 0\\ h_{51} & h_{52} & h_{53} & h_{54} & h_{55} & h_{56}\\ h_{61} & h_{62} & h_{63} & h_{64} & h_{65} & h_{66} \end{pmatrix}\]

Esta matriz de transición correspondería a una cola donde pueden llegar como máximo \(m\) paquetes por unidad de tiempo

Para calcular el vector de probabilidades en el equilibrio, debemos resolver \(Ps = s\)

\[\begin{pmatrix} h_{11} & h_{12} & 0 & 0 & \cdots\\ h_{21} & h_{22} & h_{23} & 0 & \cdots\\ h_{31} & h_{32} & h_{33} & h_{34} & \cdots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix}\begin{pmatrix} s_1\\ s_2\\ s_3\\ \vdots\end{pmatrix} = \begin{pmatrix} s_1\\ s_2\\ s_3\\ \vdots \end{pmatrix}\]

donde \(s = (s_1,s_2,\dots)^t\)

La primera fila da

\[s_1 = h_{11}s_1+h_{12}s_2\]

Supongamos que \(s_1=1\). Al final, cuando normalicemos, ya calcularemos el verdadero valor de \(s_1\).

El valor \(s_2\) será, por tanto \[s_2 = \frac{1-h_{11}}{h_{12}}\]

La segunda fila da

\[s_2 = h_{21}s_1+h_{22}s_2+h_{23}s_3\]

de donde

\[s_3 = \frac{(1-h_{11})(1-h_{22})}{h_{12}h_{23}}-\frac{h_{21}}{h_{23}}\]

Continuando con las otras filas, iríamos encontrando \(s_i\), \(i>3\).

Al final tenemos que normalizar. Sea \(b\) la suma \[b = \sum_{i=1}^ms_i\] Entonces, el vector de probabilidades en el equilibrio será

\[s = \frac{1}{b}(s_1,s_2,\dots,s_m)\]

Las técnicas directas se utilizan cuando la matriz de transición no tiene ninguna estructura en particular, pero su tamaño no es suficientemente grande para que este tipo de técnicas afecte a los errores de redondeo de cálculo.

Para calcular \(s\), tenemos que resolver el sistema de ecuaciones homogéneas indeterminado de \(m\) ecuaciones y \(m\) incógnitas

\[(P-I_m)s=0\]

De cara a resolverlo, escogemos una submatriz de orden \(m-1\) de la matriz \(P-I_m\) con rango \(m-1\)

Sean \(f_1,\dots,f_{m-1}\) las filas de la submatriz anterior.

Consideremos el sistema de ecuaciones lineales formado por las filas anteriores con una ecuación extra, la condición de normalización

\[\sum_{i=1}^ms_i = 1\]

De esta forma, tendremos un sistema de ecuaciones compatible y determinado de \(m\) ecuaciones y \(m\) incógnitas. Por ejemplo, supongamos que la submatriz anterior de rango \(m-1\) está formada por las \(m-1\) primeras filas y columnas. El sistema a resolver será

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1m}\\ a_{21} & a_{22} & \cdots & a_{2m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{m-1,1} & a_{m-1,2} & \cdots & a_{m-1,m}\\ 1 & 1 & \cdots & 1 \end{pmatrix}\begin{pmatrix} s_1\\ s_2\\ \vdots\\ s_{m-1}\\s_m\end{pmatrix} = \begin{pmatrix}0\\0\\\vdots\\0\\1\end{pmatrix}\]

Las técnias iterativas son útiles cuando la matriz de transición no tiene ninguna estructura especial pero su tamaño hace que las técnicas iterativas no sean adecuadas

Estas técnicas son las siguientes

• Sobrerelajación
• Método de Jacobi
• Método de Gauss-Seidel

El vector de probabilidades \(s = (s_1,s_2,\dots)\) en el estado de equilibrio verifica

\[s_i = \sum_{j}p_{ij}s_j\qquad i = 1,2,\dots\]

Las ecuaciones anteriores se denominan ecuaciones de balance. Veamos el por qué

Sabemos que las probabilidaes \(p_{ij}\) de la matriz de transición verifican

\[\sum_jp_{ji} = 1,\quad i = 1,2,\dots\]

Multiplicando por \(s_i\) la ecuación anterior, obtenemos

\[\sum_{j}p_{ji}s_i = s_i = \sum_{j}p_{ij}s_j\]

donde hemos utilizado las ecuaciones de balance.

La parte de la izquierda de la expresión anterior representa todas las probabilidades de abandonar el estado \(i\) y la parte de la derecha representa todas las probabilidades de llegar al estado \(i\). En pocas palabras, tenemos que el flujo con respecto al estado \(i\) está balanceado