6.4 Distribuzioni invarianti

Gli esempi delle sezioni precedenti, sia nel caso di catene di Markov che per processi a salti, mostrano che, partendo da una certa densità marginale al tempo iniziale, il processo raggiunge (anche piuttosto rapidamente) un “equilibrio” in cui le densità marginali sono costanti nel tempo. Tale fenomeno è utile in svariate applicazioni, e lo studio delle possibili densità limite è quindi particolarmente rilevante. Per definire tali densità, basta considerare rispettivamente l’equazione ricorsiva (6.4) o l’equazione differenziale (6.6) e imporre che la densità marginale non cambi nel tempo e sia pertanto invariante. Diamo quindi due definizioni, nel caso a tempi discreti (catene) e a tempi continui (processi a salti).

Definizione 6.4 Sia $Q$ una matrice di transizione. Si dice che un vettore riga $\pi \in \R^E$ corrispondente ad una densità discreta sull’insieme degli stati, quindi tale che $\pi_x \in [0,1]$ per ogni $x \in E$, e $\sum_{x \in E} \pi_x = 1$ è una distribuzione invariante per $Q$ se vale \[\tag{6.8}Q. (#eq:inv-Q)\]

Definizione 6.5 Sia $L$ una matrice di intensità di salto. Si dice che un vettore riga $\pi \in \R^E$ corrispondente ad una densità discreta sull’insieme degli stati, quindi tale che $\pi_x \in [0,1]$ per ogni $x \in E$, e $\sum_{x \in E} \pi_x = 1$ è una distribuzione invariante per $L$ se vale \tag{6.9}L. (#eq:inv-L)\]

Remark. In entrambi i casi, si può quindi affermare che $\pi$ è invariante se e solo se, qualora si consideri un processo $X$ (catena o a salti) tale che la legge marginale al tempo iniziale sia $\pi$, allora tutte le leggi marginali coincidono con $\pi$.

Remark. La condizione di invarianza si può riscrivere anche come (nel caso delle catene) segue: per ogni $x \in E$, \[ \sum_{y \neq x} \pi_x Q_{x \to y} = \sum_{y \neq x} \pi_y Q_{ y \to x}.\] In questa formulazione il membro a sinistra si interpreta come flusso (di probabilità) uscente dallo stato $x$, mentre il membro a destra è un flusso entrante. L’equazione esprime quindi un bilancio di flusso. Nel caso di processi a salti, l’equazione diventa \[ \sum_{y \neq x} \pi_x L_{x y} = \sum_{y \neq x} \pi_y L_{ yx}.\]

La seguente proposizione collega il concetto di stazionarietà con il fatto che la densità delle marginali sia invariante. È ovvio che se la catena è stazionaria, allora la densità delle marginali deve essere invariante. Il viceversa richiede qualche osservazione in più, che qui non riportiamo per brevità.

Proposizione 6.1 Sia $X$ una catena di Markov o un processo di Markov a salti. Allora $X$ è stazionario se e solo se la marginale al tempo iniziale $X_0$ ha come densità una distribuzione invariante.

La domande teoriche che ci poniamo ora sono: data $Q$ (oppure $L$) le distribuzioni invarianti esistono sempre? se sì, quante sono?

Dal lato pratico è invece importante disporre di algoritmi efficienti per poter calcolare, almeno in modo approssimato, le distribuzioni invarianti.

Tali problemi si possono affrontare anche con tecniche di algebra lineare, perché le equazioni (6.8) o (6.9) sono dei sistemi di equazioni lineari omogenei nelle incognite date dalle componenti del vettore $\pi$. In questa sezione ci limitiamo ad esporre i risultati principali, accennando alle dimostrazioni.

6.4.1 Esistenza

Il primo risultato riguarda l’esistenza di (almeno) una distribuzione invariante. La risposta è affermativa se l’insieme degli stati $E$ è finito.

Teorema 6.1 Se l’insieme degli stati $E$ è finito, e quindi la matrice di transizione $Q$ oppure delle intensità di salto $L$ sono matrici con un numero finito di righe e colonne, allora esiste sempre almeno una distribuzione invariante $\pi$.

Proof. Vediamo prima la dimostrazione nel caso delle catene di Markov. Si consideri una qualsiasi densità discreta $\pi_0$ sull’insieme degli stati – intepretata come vettore riga. Sappiamo che se $\pi_0$ è la densità marginale di una catena di Markov $(X_n)_n$ al tempo $n=0$, la densità marginale al tempo $k=0,1,2,\ldots$, è $\pi_0 Q^k$. Se esiste il limite per $k \to \infty$ delle $\pi_0 Q^k$, esso è una distribuzione invariante, ma l’esistenza in generale non è garantita (anzi in generale è falsa). Tuttavia, se si considerano le medie aritmetiche \[ \bar{\pi}_n = \frac 1 n \sum_{i=1}^n \pi_0 Q^k,\] allora ciascuna $\bar{\pi}_n$ è una densità discreta di probabilità (quindi un vettore a componenti in $[0,1]$ e a somma $1$), e per l’estensione al caso vettoriale del teorema di Bolzano-Weierstrass, esiste una sottosuccessione $\bar{\pi}_{n_k}$ con $n_k \to \infty$ che converge ad un limite \[ \bar{\pi}_\infty = \lim_{k \to \infty} \bar{\pi}_{n_k},\] ossia tale che ogni componente del vettore $\bar{\pi}_{n_k}$ converge alla corrispondente componente di $\bar{\pi}_\infty$. Anche il limite è una densità discreta di probabilità sugli stati $E$, perché ciascuna componente del vettore è in $[0,1]$, essendo limiti di valori compresi tra $0$ e $1$, e la somma dei limiti delle componenti coincide con il limite della somma, che vale $1$ (qui si usa che $E$ è finito, quindi la somma non è una serie). Quindi, basta dimostrare che vale \[ \bar{\pi}_\infty = \bar{\pi}_\infty Q.\] Per ogni $n$, si ha l’identità \[ \begin{split} \bar{\pi}_n Q & = \bra{ \frac 1 n \sum_{k=1}^n \pi_0 Q^k} Q \\ & = \frac 1 n \sum_{k=1}^n \pi_0 Q^{k+1} \\ & = \frac{1}{n}\sum_{k=1}^n \pi_0 Q^{k} + \frac{1}{n}\bra{ \pi_0 Q^{n+1} - \pi_0 Q}\\ & = \bar{\pi}_n + \frac{1}{n}\bra{ \pi_0 Q^{n+1} - \pi_0 Q}.\end{split}\] Al tendere di $n \to \infty$, il termine \[ \frac{1}{n}\bra{ \pi_0 Q^{n+1} - \pi_0 Q} \to 0\] è infinitesimo al tendere di $n \to \infty$, perché le componenti del vettore $\pi_0 Q^{n+1}$ sono comprese tra $[0,1]$, e si divide per $n$. Ponendo $n = n_k \to \infty$, concludiamo quindi che $\bar{\pi}_\infty$ è una distribuzione invariante.

Nel caso di processi di Markov a salti, l’idea è analoga, ma si sostituisce alla media aritmetica una media integrale sui tempi, ponendo \[ \bar{\pi}_t = \frac 1 t \int_0^t \pi_0 \exp\bra{ s L} d s. \] Si può estrarre anche in questo caso una successione $t_k\to \infty$ in modo che $\bar{\pi}_{t_k}$ converga a un vettore $\bar{\pi}_\infty$, che è una densità discreta su $E$. Per mostrare che è invariante, si calcola \[ \begin{split} \bar{\pi}_t L & = \frac 1 t \int_0^t \pi_0 \exp\bra{ s L} L d s\\ & = \frac 1 t \int_0^t \pi_0 \frac{d}{ds} \exp\bra{ s L} d s \\ & = \frac 1 t \pi_0 \int_0^t \frac{d}{ds} \exp\bra{ s L} d s\\ & = \frac 1 t \pi_0 \bra{ \exp\bra{t L} - Id},\end{split}\] che al tendere di $t \to \infty$ tende a zero (si noti l’analogia con il caso a tempi discreto).

Questo risultato teorico garantisce quindi che i sistemi lineari (6.8) o (6.9) ammettano almeno una soluzione, e pertanto è sufficiente usare un qualsiasi metodo risolutivo per determinarla (ad esempio la riduzione a gradini di Gauss). Tuttavia, la dimostrazione stessa suggerisce un metodo approssimato basato sul calcolo delle medie $\bar{\pi}_n$ nel caso di catene o $\bar{\pi}_{t}$ nel caso di processi di Markov a salti.

Esempio 6.10 Si consideri la matrice $Q$ definita in (6.2). Per risolvere il sistema omogeneo, conviene passare alla trasposta e determinare tutti i vettori colonna $x \in \R^2$ tali che \[ (Id - Q^T) x = 0.\] Nel caso di interesse, si trova che sono tutti del tipo \[ x = u (7,1 )\] dove $u \in \R$ è un parametro da determinare imponendo che la somma delle componenti di $x$ sia $1$. Concludendo si trova quindi $u=8$ e quindi \[ \pi = (\frac 7 8, \frac 1 8 ).\] Vediamo come procedere tramite R. Diamo prima una soluzione “esatta” mediante il comando eigen(), che determina autovalori e autovettori: in questo caso ci interessano infatti gli autovettori di $Q^T$ con autovalore $1$, o equivalentemente gli autovettori di $Id - Q^T$ con autovalore $0$ (il nucleo di $Id - Q^T$).

Q <- matrix(c(0.9, 0.1, 0.7, 0.3), nrow = 2,
  byrow = TRUE)

# calcoliamo autovalori e autovettori
# della matrice Q trasposta

sol <- eigen(t(Q))

sol$values

## [1] 1.0 0.2

sol$vectors

##           [,1]       [,2]
## [1,] 0.9899495 -0.7071068
## [2,] 0.1414214  0.7071068

# siamo interessati all'autovettore di
# autovalore 1 che è il primo

pi_inv <- sol$vectors[, 1]
pi_inv <- pi_inv/sum(pi_inv)

# la distribuzione invariante è quindi

pi_inv

## [1] 0.875 0.125

# rappresentiamola con un barplot

alternative <- c("OFF", "ON")
colori <- miei_colori[1]

barplot(pi_inv, col = colori, names.arg = alternative,
  ylab = "probabilità", xlab = "stato")

Esempio 6.11 Determiniamo le distribuzioni invarianti nel caso della matrice di intensità di salto (6.5).

L <- matrix(c(-15, 5, 10, 1, -4, 3, 0, 4,
  -4), nrow = 3, byrow = TRUE)


sol <- eigen(t(L))

sol$values

## [1] -1.514005e+01 -7.859945e+00  8.881784e-16

sol$vectors

##            [,1]        [,2]        [,3]
## [1,]  0.7539667 -0.09196364 -0.04908437
## [2,] -0.1055968 -0.65662548 -0.73626560
## [3,] -0.6483699  0.74858912 -0.67491013

# siamo interessati all'autovettore di
# autovalore 0, che è il terzo

pi_inv <- sol$vectors[, 3]
pi_inv <- pi_inv/sum(pi_inv)

# la distribuzione invariante è quindi

pi_inv

## [1] 0.03361345 0.50420168 0.46218487

# rappresentiamola con un barplot

alternative <- c("Off", "Standby", "On")
colori <- miei_colori[1]

barplot(pi_inv, col = colori, names.arg = alternative,
  ylab = "probabilità", xlab = "stato")

Remark. È fondamentale non confondere l’equazione $\pi (Id - Q) = 0$, oppure $\pi L = 0$ con le “trasposte” $(Id - Q) v = 0$ e $L v = 0$. Infatti queste hanno sempre come soluzione la densità uniforme $v_x = 1/n$, dove $n$ è il numero degli elementi di $E$. Questo perché la somma su ciascuna riga delle matrici di transizione vale $1$ (mentre vale $0$ nel caso delle matrici di intensità di salto).

Se non si specifica l’operazione di trasposizione si trova sempre tale soluzione, che però non è quella cercata.

Vi è tuttavia un caso in cui la densità uniforme è sicuramente una distribuzione invariante: si tratta del caso in cui la matrice $Q$ sia bistocastica, ossia anche la somma sulle colonne dia $1$. Un caso particolare è quando $Q$ sia simmetrica.

Esempio 6.12 L’potesi che l’insieme degli stati $E$ sia finito è necessaria. Si consideri una catena sugli stati $E = \mathbb{N}$ e $Q_{x \to x+1} = 1$ per ogni $x \in \mathbb{N}$ e $0$ altrimenti. Allora ogni distribuzione invariante deve essere uniforme, ma essendo gli stati infiniti una tale distribuzione non esiste.

6.4.2 Unicità

Affrontiamo il problema dell’unicità, mostrando in quali condizioni la distribuzione invariante è unica e che è possibile classificare tutte le distribuzioni invarianti associate ad una matrice (di transizione $Q$ oppure di intensità di salto $L$). Cominciamo da un esempio:

Esempio 6.13 Un gioco tra Alice e Bruno consiste nel lanciare un dado a sei facce ripetutamente fintanto che non esca il numero $1$ (e in tal caso vince Alice) oppure il numero $6$ (e e in tal caso vince Bruno). Possiamo rappresentare una partita tramite una catena di Markov sugli stati $E = \cur{A, G, B}$, dove $A$ indica che Alice ha vinto, $B$ Bob ha vinto, e $G$ il gioco continua (si deve lanciare nuovamente il dado).

Figura 6.6: Rappresentazione grafica del gioco tra Alice e Bruno

È chiaro che vi sono almeno due distribuzioni stazionarie, corrispondenti rispettivamente al caso in cui Alice vinca, $\pi_A = (1,0,0)$ oppure Bruno vinca, $\pi_B = (0,0,1)$. Ma in realtà sono infinite, perché ogni combinazione \[ \alpha \pi_A + (1-\alpha) \pi_B = (\alpha, 0, 1-\alpha),\quad \text{con $\alpha \in [0,1]$}\] è pure una distribuzione invariante (corrispondente al fatto che Alice vinca con probabilità $\alpha$ e Bob con probabilità $1-\alpha$). In effetti, è intuitivamente chiaro che, se la catena di Markov al tempo iniziale si trova in $G$, la densità limite è quella corrispondente ad $\alpha = 1/2$, perché le regole del gioco non favoriscono né Alice né Bruno.

L’esempio sopra evidenzia un fatto generale: se vi sono almeno due distribuzioni invarianti (diverse), allora le distribuzioni invarianti sono infinite perché tutte le combinazioni come sopra, al variare del parametro $\alpha\in [0,1]$ sono pure invarianti. Se l’insieme degli stati $E$ è finito, è possibile determinare un numero finito di distribuzioni invarianti “di base” in modo da poter rappresentare tutte le altre distribuzioni invarianti come combinazioni di esse. La procedura è “geometrica” e si basa sulla decomposizione dell’insieme degli stati in sottoinsiemi più piccoli, dette classi irriducibili. Diamo alcune definizioni.

Definizione 6.6 Data una matrice di transizione $Q$, si dice che lo stato $y \in E$ è accessibile (o raggiungibile) da $x \in E$ se esiste un cammino $\gamma = (x_0 \to x_1 \to \ldots \to x_n)$ con $x_0 = x$ e $x_n = y$ e peso strettamente positivo: \[ Q_\gamma = \prod_{k=1}^n Q_{x_{k-1} x_k} >0. \] Nel caso di una matrice di intensità di salto $L$, il peso $Q_\gamma$ va sostituito con \[ L_\gamma = \prod_{k=1}^n L_{x_{k-1} x_k}.\]

Ricordando che $Q^n_{xy}$ è la somma dei pesi di tutti i cammini lunghi $n$ che collegano $x$ a $y$, si può equivalentemente dire che $y$ è raggiungibile da $x$ se esiste $n \ge 1$ tale che $Q^n_{xy}>0$.

Osserviamo che $y \in E$ è raggiungibile da $x\in E$, e $z \in E$ è raggiungibile da $y$, allora $z$ è raggiungibile da $x$. Tuttavia non è detto che $x$ sia raggiungibile da $y$ (o da $z$): se questo appunto non accade, lo stato $x$ è detto transitorio.

Definizione 6.7 Data una matrice di transizione $Q$ (o una matrice di intensità di salto $L$) uno stato $x \in E$ è detto transitorio se esiste $y \in E$ tale che $y$ è raggiungibile da $x$, ma $x$ non lo è da $y$. Se $x$ non è transitorio, è detto ricorrente.

Equivalentemente, $x$ è ricorrente se, per ogni $y\in E$ tale che $y$ è raggiungibile da $x$, anche $x$ lo è da $y$. Se l’insieme degli stati $E$ è finito, non possono essere tutti transitori e deve esserci almeno uno stato ricorrente.

Esempio 6.14 Consideriamo una catena di Markov di cui rappresentiamo graficamente la matrice di transizione come segue (completare per esercizio le probabilità mancanti, in modo che la somma sulle righe, ossia gli archi uscenti da ciascuno stato, sia $1$): Allora stato $1$ è raggiungibile dallo stato $2$, mentre lo stato $2$ non è raggiungibile dallo stato $1$. Lo stato $2$ è quindi transitorio. Anche lo stato $5$ lo è. Tutti i rimanenti stati $\cur{1,3,4,6}$ sono ricorrenti.

Introduciamo ora il concetto di classe chiusa irriducibile.

Definizione 6.8 Data una matrice di transizione $Q$ (o una matrice di intensità di salto $L$), un sottoinsieme $C \subseteq E$ di stati è detto classe chiusa se, per ogni $x \in C$ e $y \in E$ raggiungibile da $x$, anche $y \in C$. Una classe chiusa $C$ è detta irriducibile se non contiene altre classi chiuse $C' \subseteq C$ (diverse dai casi banali $C' = \emptyset$ oppure $C' = C$ stessa). La matrice $Q$ (oppure $L$) è detta irriducibile se tutto l’insieme degli stati $E$ è una classe chiusa irriducibile.

Equivalentemente, da uno stato in una classe chiusa non è possibile raggiungere stati al di fuori di essa (mentre è possibile entrarvi), e una classe chiusa è irriducibile quando da ogni stato in essa si può raggiungere qualsiasi altro stato in essa. Data una classe chiusa è ben definita la restrizione della matrice $Q$ su $C \times C$, perché $Q_{x \to y}=0$ per $x \in C$ e $y \notin C$, e quindi $(Q_{x \to y})_{y \in C}$ sono densità discrete di probabilità (la somma delle righe vale ancora $1$). Un ragionamento analogo vale nel caso di matrici di intensità di salto $L$.

Esempio 6.15 Riprendiamo l’Esempio 6.14. Gli stati $\cur{1,4}$ sono una classe chiusa (irriducibile), mentre ad esempio gli stati $\cur{1,2}$ non lo sono, perché si può “uscirne” visitando ad esempio lo stato $4$.

L’importanza della condizione di irreducibilità è dovuta del seguente risultato.

Teorema 6.2 Sia $Q$ una matrice di transizione (oppure $L$ di intensità di salto) irriducibile su un insieme di stati $E$ finito. Allora esiste una e una sola distribuzione invariante.

Proof. Abbiamo già visto che almeno una distribuzione invariante $\pi$ esiste, quindi basta mostrare l’unicità. Nel caso di matrice $Q$ di transizione, introduciamo la matrice \[ R = \frac 1 2 \sum_{n=0}^\infty 2^{-n} Q^n,\] che è pure una matrice di transizione, in particolare la somma sulle righe vale $1$ perché \[ \frac 1 2 \sum_{n=0}^\infty 2^{-n} = 1.\] Grazie all’ipotesi di irriducibilità è tale che $R_{xy}>0$ per ogni $x$, $y \in E$. Nel caso di $L$ matrice di intensità di salto, introduciamo invece $R = \exp(L)$ che pure è una matrice di transizione e similmente ha tutte entrate strettamente positive. Se $\pi$ è una distribuzione invariante per $Q$ (o per $L$), vale l’identità \[ \pi R = \frac 1 2 \sum_{n=0}^\infty 2^{-n} \pi Q^ n = \pi \frac 1 2 \sum_{n=0}^\infty 2^{-n} = \pi,\], ossia $\pi$ è distribuzione invariante anche per la matrice di transizione $R$ (nel caso di Markov a salti il ragionamento è anche più semplice). Consideriamo quindi una seconda distribuzione invariante $\tilde{\pi}$ e per mostrare che $\pi = \tilde \pi$ scriviamo la seguente diseguaglianza: \[\begin{split} \sum_{x \in E} | \pi_x - \tilde \pi_x| & = \sum_{x \in E} | \sum_{y \in E} \pi_y R_{yx} - \sum_{y \in E} \tilde \pi_y R_{yx}| \\ & = \sum_{x \in E} | \sum_{y \in E} (\pi_y - \tilde \pi_y) R_{yx} | \\ & \le \sum_{x \in E} \sum_{y \in E} | \pi_y-\tilde \pi_y| R_{yx} \\ & = \sum_{y \in E} | \pi_y-\tilde \pi_y| \sum_{x \in E} R_{yx} \\ & = \sum_{y \in E} | \pi_y-\tilde \pi_y|.\end{split}\] Poiché la prima e l’ultima espressione coincidono (è solamente cambiato l’indice di somma), devono essere tutte uguaglianze, in particolare nel passaggio in cui si è stimato, per ogni $x \in E$, \[ | \sum_{y \in E} (\pi_y - \tilde \pi_y) R_{yx} | \le \sum_{y \in E} | \pi_y-\tilde \pi_y| R_{yx}\] deve in realtà valere l’uguaglianza. È noto tuttavia che la diseguaglianza triangolare tra numeri reali \[ | \sum_i z_i | \le \sum_i |z_i|\] è una uguaglianza se e solo se hanno tutti lo stesso segno, ossia $z_i \ge 0$ per ogni $i$ oppure $z_i \le 0$ per ogni $i$. Se quindi vale per ogni $y \in E$ \[ (\pi_y - \tilde \pi_y )R_{yx} \ge 0,\] essendo $R_{yx}>0$ si ottiene (dividendo) che $\pi_y \ge \tilde \pi_y$ per ogni $y \in E$. Ma poiché sono entrambe densità di probabilità, \[ \sum_{y \in E} \pi _y = 1 = \sum_{y \in E}\tilde{\pi}_y,\] ne segue che deve valere $\pi_y = \tilde \pi_y$ (se la diseguaglianza fosse stretta per qualche $y$ allora non potrebbe valere l’uguaglianza nella somma).

Remark. La diseguaglianza principale nella dimostrazione sopra si può usare per mostrare anche che, data una qualsiasi densità discreta $\pi_0$ su $E$ la quantità detta variazione totale tra $\pi_0 Q^n$ e la distribuzione invariante $\pi$, \[ \sum_{x\in E} | \pi - \pi_0 Q^n |\] decresce al crescere di $n$. Questa è una quantità utile per stimare la convergenza delle densità marginali di una catena verso la distribuzione invariante.

Nel caso in cui $Q$ non sia irriducibile, bisogna preliminarmente decomporre $E$ in classi chiuse irriducibili. Precisamente, ogni classe chiusa irriducibile non contiene stati transitori (intuitvamente basterebbe infatti rimuoverli per ottenere una classe chiusa più piccola), e due classi chiuse irriducibili diverse sono necessariamente disgiunte (altrimenti l’intersezione sarebbe più piccola di entrambe). Ne segue che, se $E$ è finito, si può partizionare in insiemi a due a due disgiunti \[ E = \cur{\text{stati transitori}} \cup C^1 \cup C^2 \cup \ldots \cup C^k,\] dove le $C^i$ sono classi chiuse irriducibili ($k \ge 1$ perché non tutti gli stati sono transitori).

Se consideriamo la restrizione di $Q$ su ciascuna classe chiusa irriducibile $C^i$, allora esiste una e una sola distribuzione invariante associata $\pi^i$, che può essere vista come una distribuzione invariante su tutto l’insieme degli stati $E$ (semplicemente la probabilità assegnata al di fuori di $C^i$ è nulla). Si può dimostrare (non lo faremo) che tutte le distribuzioni invarianti sono ottenibili come combinazioni \[ \pi = \alpha_1 \pi^1 + \alpha_2 \pi^2 + \ldots, + \alpha_k \pi^k,\] dove $\alpha_1$, $\alpha_2$, …, $\alpha_k \in [0,1]$ sono tali che \[\alpha_1+\alpha_2 +\ldots + \alpha_k = 1,\] ossia $(\alpha_i)_{i=1}^k$ sono una (qualsiasi) densità discreta di probabilità. In particolare, ogni distribuzione invariante è nulla sugli stati transitori (che quindi si possono subito trascurare se in un problema è richiesto di calcolare tutte le distribuzioni invarianti).

Esempio 6.16 Riprendiamo l’Esempio 6.14. Per determinare tutte le distribuzioni invarianti, basta calcolare l’unica distribuzione per la classe chiusa irriducibile $C_1 = \cur{1,4}$, che si trova impostando il bilancio di flusso, ad esempio in $1$: \[ \pi_1 \frac 2 3 = \pi_4 \frac 1 3 \quad \text{da cui} \quad (\pi_1, \pi_4) = (\frac 1 3 , \frac 2 3 ).\] Similmente, per la classe chiusa irriducibile $C_2 = \cur{3,6}$, si trova \[ \pi_3 \frac 2 3 = \pi_6 \quad \text{da cui} \quad (\pi_3, \pi_6) = (\frac 3 5, \frac 2 5).\] Di conseguenza tutte le distribuzioni invariante della catena originaria sono parametrizzate nel seguente modo: \[ (\alpha \frac 1 3, 0, (1-\alpha) \frac 3 5, \alpha \frac 2 3, 0, (1-\alpha)\frac 2 5 ),\] dove $\alpha \in [0,1]$ è un parametro (invece di usare $2$ parametri $\alpha_1$, $\alpha_2$ che sommano a $1$ ne indichiamo solo uno).

6.4.3 Sul limite delle potenze della matrice di transizione

In questa sezione, principalmente utile per svolgere alcuni esercizi, analizziamo più nel dettaglio sotto quali condizioni, data una catena di Markov¹⁵ con matrice di tranzione $Q$ su un insieme di stati finito, il limite \[ \lim_{n \to \infty} (Q^n)_{ij}\] esiste per ogni coppia di stati $i$, $j \in E$. Inoltre presentiamo un metodo per calcolare tale limite mediante sistemi di equazioni lineari.

Osserviamo che l’intepretazione di $(Q^n)_{ij}$ in termini di probabilità è, per via del risultato generale sulle densità marginali, ponendo $\pi_0$ la densità discreta che vale $1$ nello stato $i$ e $0$ altrimenti, \[ (Q^n)_{ij} = P(X_n = j | X_0 =i).\] D’altra parte le potenze $Q^n$ intervengono sia nel teorema di esistenza delle distribuzioni invarianti, sia nel teorema di unicità (un po’ implicitamente, nelle ipotesi per via della condizione di irriducibilità, e nella dimostrazione per la costruzione della matrice $R$, con tutte le entrate positive).

L’osservazione di base è che, se il limite esiste, \[ Q^\infty_{ij} = \lim_{n \to \infty} (Q^n)_{ij},\] allora è una matrice le cui righe, ossia i vettori $(Q^\infty_{ij})_{j \in E}$, sono tutte distribuzioni invarianti. Infatti vale \[ \begin{split} Q^\infty & = \lim_{n \to \infty }Q^n = \lim_{n \to \infty } Q^{n+1} = \lim_{n \to \infty }Q^n \cdot Q\\ & = Q^\infty\cdot Q.\end{split}\] Che scritta componente per componente diventa il sistema definente le distribuzioni invarianti: \[ Q^\infty_{ij} = \sum_{k} Q^{\infty}_{ik} Q_{k j}.\] Inoltre poiché le righe di $Q^n$ sono densità discrete di probabilità, anche il limite lo è.

Si potrebbe pensare che, se la catena è irriducibile, quindi esiste una sola distribuzione invariante $\pi$, allora necessariamente il limite $Q^\infty$ esiste ed è la matrice in cui tutte le righe sono identiche e uguali a $\pi$. Il ragionamento è corretto, purché il limite esista, cosa che non accade sempre, come il prossimo esempio mostra.

Esempio 6.17 Sia $E = \cur{\text{Off}, \text{On}}$ e sia (ordinando gli stati nell’ordine in cui sono scritti), \[ Q = \bra{ \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}}.\] Il grafo associato è il seguente:

La “dinamica” della catena è molto semplice: se inizia nello stato Off, allora passerà a tempi alterni da Off in On. Questo si vede anche calcolando le potenze di $Q$, \[ Q^2 = \bra{ \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}} = Id \quad Q^3 = \bra{ \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}}= Q, \quad Q^4 = Id\ldots\] quindi il limite $Q^\infty$ non esiste perché le componenti si alternano tra i valori $0$ e $1$. In questo esempio tuttavia è facile vedere che la catena è irriducibile con l’unica distribuzione invariante (uniforme).

La condizione di irriducibilità non è quindi sufficiente a garantire che il limite $Q^\infty$ esista. La seguente nozione è invece il concetto corretto.

Definizione 6.9 Sia $Q$ una matrice di transizione di una catena di Markov su un insieme $E$ di stati finito. Diciamo che $Q$ (o la catena) è regolare se esiste $n\in \mathbb{N}$ tale che la potenza $Q^n$ ha tutte entrate strettamente positive, ossia \[ (Q^n)_{ij} >0 \quad \text{per ogni $i$, $j \in E$.}\]

La definizione somiglia molto a quella di catena irriducibile, ma a ben vedere è diversa: nel caso di catena irriducibile si richiede che per ogni $i$, $j \in E$ esista un cammino di una qualsiasi lunghezza $n$ che li collega, ossia $(Q^n)_{ij}>0$. Nel caso di catena regolare, si richiede che la lunghezza $n$ sia la stessa per tutti gli $i$, $j \in E$ (anche quando $i=j$). Questo ragionamento ci dice quindi che se la catena è regolare, allora è anche irriducibile. Ma non necessariamente il viceversa, come mostra l’esempio di prima.

L’importanza del concetto è data dal seguente teorema, che non dimostriamo, ma è utile negli esercizi.

Teorema 6.3 Sia $Q$ una matrice di transizione irriducibile su un insieme di stati $E$ finito. Allora $Q$ è regolare se e solo se esiste il limite \[ Q^\infty_{ij} = \lim_{n \to \infty}(Q^{n})_{ij} \quad \text{per ogni $i$, $j \in E$.}\] Più in generale, anche se $Q$ non è irriducibile, il limite esiste (per ogni $i$, $j \in E$) se e solo se la catena ristretta a ciascuna classe chiusa irriducibile è regolare.

Appoggiandosi a questo risultato possiamo garantire l’esistenza di $Q^\infty$ in molti problemi concreti: basta prima classificare gli stati e le classi chiuse irridicubili e controllare che la matrice ristretta a ciascuna classe sia regolare.

Ci sono però due aspetti pratici da non sottovalutare: come verificare che $Q$ sia regolare? Una strategia è di moltiplicare $Q$ per se stessa, finché le componenti non sono tutte positive. Un metodo più semplice è di considerare solo le potenze di $2$, ossia calcolare \[Q^2 = Q \cdot Q, \quad Q^4 = Q^2 \cdot Q^2, \quad Q^8 = Q^4\cdot Q^4, \quad \text{ecc.,}\] In questo modo l’esponente $n$ cresce più rapidamente e l’algoritmo termina prima (se termina). Senza un calcolatore, tuttavia tale algoritmo non è molto pratico. Ci sono però diversi criteri per la regolarità, e suggeriamo il seguente (senza dimostrazione) che può essere utile negli esercizi.

Proposizione 6.2 Sia $Q$ una matrice di transizione irriducibile su un insieme di stati $E$ finito. Se esiste (almeno) uno stato $i\in E$ tale che $Q_{ii}>0$, allora $Q$ è anche regolare.

Basta quindi osservare, dopo aver controllato che $Q$ sia irriducibile, che la diagonale della matrice $Q$ non sia identicamente nulla. Attenzione tuttavia, la condizione è solo sufficiente, esistono catene $Q$ irriducibili con diagonale tutta nulla ma comunque regolari. Questo criterio comunque risolve la maggior parte dei casi concreti (in particolare negli esercizi).

Appurato quindi che $Q^\infty$ esista (oppure ignorando il problema temporaneamente), la domanda successiva è come calcolarlo. L’osservazione chiave è di ripartire dall’argomento che mostrava tutte le righe invarianti, e trovare una seconda equazione (o sistema di equazioni):

\[ \begin{split} Q^\infty & = \lim_{n \to \infty }Q^n = \lim_{n \to \infty } Q^{1+n} = \lim_{n \to \infty }Q\cdot Q^n\\ & = Q \cdot Q^\infty.\end{split}\]

Scrivendo la relazione per ciascun elemento della matrice, diventano le equazioni \[ Q^{\infty}_{ij} = \sum_{k\in E} Q_{ik} Q^\infty_{kj}, \quad \text{per $i$, $j\in E$.}\]

Queste equazioni, aggiunte all’informazione che le righe di $Q^\infty$ sono distribuzioni invarianti, permettono di ottenere un sistema che, una volta risolto, determina completamente la matrice $Q^\infty$.

Esempio 6.18 Riprendiamo l’Esempio 6.13 del gioco tra Alice e Bruno. Sappiamo che tutte le distribuzioni invarianti sono della forma $(\alpha, 0, 1-\alpha)$, con $\alpha \in [0,1]$. È ovvio che se lo stato iniziale della catena è $A$ (vince Alice), che corrisponde ad $\alpha = 1$, allora la catena rimane in quello stato (è assorbente), e quindi la riga corrispondente nella matrice $Q^\infty$ è \[(Q^\infty_{A A}, Q^\infty_{AG}, Q^\infty_{AB}) = (1,0,0).\] Similmente, se lo stato iniziale è $B$, \[(Q^\infty_{B A}, Q^\infty_{BG}, Q^\infty_{BB}) = (0,0,1).\] Resta da determinare la riga corrispondente allo stato iniziale transitorio $G$. Abbiamo già osservato che per simmetria della situazione deve essere $\alpha =1/2$, però possiamo anche procedere scrivendo le equazioni sopra (non serve scriverle tutte, basta trovarne una che permetta di determinare $\alpha$). Ad esempio con $i= G$, $j=A$, troviamo \[ \begin{split} \alpha & = Q^\infty_{GA} = Q_{GA}Q^{\infty}_{AA} + Q_{GG}Q^{\infty}_{GA} + Q_{GB}Q^{\infty}_{BA}\\ &= \frac 1 6 \cdot 1 + \frac 4 6 \cdot \alpha + \frac 1 6 \cdot 0\end{split}\] dove abbiamo usato le righe di $Q^\infty$ precedentemente calcolate (quelle relative agli stati $A$ e $B$). Concludiamo che $\alpha = \frac 1 2$ come avevamo già osservato.

Remark. Esiste in realtà una formula generale per determinare $Q^\infty$. Tuttavia per ottenerla conviene effettuare un passaggio intermedio per semplificare la struttura della catena e riportarsi in un certo senso all’esempio sopra, in cui gli stati sono o transitori o assorbenti.

A meno di scegliere opportunamente un ordinamento degli stati, possiamo supporre che che gli stati transitori siano $\cur{1, \dots, \ell}$ e vi siano $k$ classi chiuse irriducibili $C^1$, $C^2$, …, $C^k$. Dal teorema di classificazione, le distribuzioni invarianti sono determinate da una scelta dei parametri $(\alpha_{c^j})_{j=1}^k$, dove $\alpha_{c^j}$ indica la probabilità di “entrare” nella classe $C^j$. Il nostro obiettivo è quindi determinare solamente tali probabilità, a partire dagli stati transitori $\cur{1, \ldots, \ell}$ (perché se partiamo da uno stato ricorrente, esso apparterrà ad una classe chiusa irriducibile $C^j$ e quindi sarà solo $\alpha_{c^j} =1$ e gli altri nulli). Per determinare $Q^\infty$ la vera incognita sono quindi le probabilità \[ \alpha := (\alpha_{ic^j})_{i=1,\ldots, \ell}^{j=1, \ldots, k} \in [0,1]^{\ell \times k}\] che abbiamo organizzato in una matrice. L’idea è ora di “collassare” le classi chiuse in singoli stati assorbenti, definendo una matrice di transizione più semplice da trattare. Precisamente, introduciamo l’insieme degli stati \[ \bar E = \cur{1, \ldots,\ell} \cup\cur{c^1, c^2, \ldots, c^{k} }\] in cui ogni classe chiusa $C^j$ corrisponde ora ad un singolo stato. Definiamo una nuova matrice di transizione $Q'$ nel seguente modo: lasciamo invariate le probabilità di transizione tra coppie di stati transitori, ossia poniamo per $i,j=1, \ldots, \ell$, \[\bar Q _{i j} = Q_{ij},\] mentre definiamo tutti gli stati come $c^j$ assorbenti, ossia \[ \bar Q_{c^jc^j} = 1 \quad \text{e} \quad \bar Q_{i c^j} = 0 \text{se $i\neq c^j$,} \] e infine la probabilità di transizione da uno stato transitorio $i =1, \ldots, \ell$ ad uno stato $c^j$, $j=1, \ldots, k$, è data dalla somma di tutte le probabilità (secondo $Q$) di passare da $i$ ad uno qualsiasi degli stati in $C^j$, ossia \[ \bar Q_{i c^j} = \sum_{x \in C^j} Q_{ix},\] che è semplicemente la probabilità di entrare dallo stato $i$ nella classe $C^j$ (in un singolo passo della catena).

Tale catena è in pratica molto semplice da definire, e ha naturalmente una struttura “a blocchi”, per via delle definizioni date: \[ Q' = \bra{ \begin{array}{cc} Q_{TT} & \bar{Q}_{TC} \\ 0 & Id\end{array}},\] dove $Q_{TT} = \bar{Q}_{TT}$ indica il blocco $\ell \times \ell$ delle probabilità di transizione tra stati transitori, $\bar{Q}_{TC}$ indica il blocco $\ell \times k$ corrispondente alle probabilità di transzione da transitori agli stati assorbenti corrispondenti alle classi chiuse irriducibili, $0$ e $Id$ sono rispettivamente una matrice $\ell \times k$ di tutti zeri e la matrice identità di dimensione $k\times k$. Tenendo conto di questa decomposizione in blocchi, si può mostrare che la matrice $\alpha$ (che contiene i parametri da determinare nel problema originario) è data da \[ \alpha = (Id-Q_{TT})^{-1} \bar{Q}_{TC}.\]

Questo perché si mostra che la matrice $\bar{Q}^\infty = \lim_{n \to \infty} \bar{Q}^n$ esiste ed è data da \[ \bar Q^\infty = \bra{ \begin{array}{cc} 0 & \alpha \\ 0 & Id\end{array}}\] e quindi il sistema da risolvere diventa, in forma matriciale \[ \begin{split} \bra{ \begin{array}{cc} 0 & \alpha \\ 0 & Id\end{array}} & = \bra{ \begin{array}{cc} Q_{TT} & \bar{Q}_{TC} \\ 0 & Id\end{array}} \cdot \bra{ \begin{array}{cc} 0 & \alpha \\ 0 & Id\end{array}} \\ & = \bra{ \begin{array}{cc} 0 & Q_{TT} \alpha +\bar{Q}_{TC}\\ 0 & Id\end{array}}\end{split}\] e quindi troviamo che \[ \alpha = Q_{TT} \alpha + \bar{Q}_{TC} \quad \text{ossia} \quad (Id-Q_{TT}) \alpha = \bar{Q}_{TC}\] e invertendo $Id-Q_{TT}$ (si mostra che è possibile farlo) si trova la formula.

Esempio 6.19 Riprendiamo l’Esempio 6.14. Ponendo $C_1 =\cur{1,4}$, $C_2= \cur{3,6}$ troviamo che la catena di Markov su $\bar{E}$ è data graficamente da In questo caso, le matrici $\bar{Q}_{TT} = Q_{TT}$ e $\bar{Q}_{TC}$ sono rispettivamente \[ Q_{TT} = \bra{ \begin{array}{cc} 0 & 2/7 \\ 5/7 & 0 \end{array}},\] mentre \[ \bar{Q}_{TC} = \bra{ \begin{array}{cc} 3/7 & 2/7 \\ 1/7 & 1/7 \end{array}}. \] Osserviamo che $\det(Id - Q_{TT}) = 1-10/49 = 39/49>0$, da cui \[ (Id - Q_{TT})^{-1} = \frac {49}{39} \bra{ \begin{array}{cc} 1 & 2/7 \\ 5/7 & 1 \end{array}}\] Si conclude quindi che $\alpha$ vale \[ \begin{split} \alpha & = (Id - Q_{TT})^{-1}\bar{Q}_{TC} = \frac {49}{39} \bra{ \begin{array}{cc} 1 & 2/7 \\ 5/7 & 1 \end{array}} \bra{ \begin{array}{cc} 3/7 & 2/7 \\ 1/7 & 1/7 \end{array}} \\ & = \frac {1}{39} \bra{ \begin{array}{cc} 23 & 16 \\ 12 & 17 \end{array}}.\end{split} \]

6.4.4 Esercizi

Esercizio 6.3 Rita e Bruno effettuano il seguente “gioco”: da un’urna contenente $R$ palline rosse e $B= N-R$ palline blu, si effettuano estrazioni con rimpiazzo fintanto che non si osservano o due palline rosse estratte consecutivamente (e in tal caso vince Rita) o due palline blu estratte consecutivamente (e in tal caso vince Bruno). Si può modellizzare tale gioco tramite una catena di Markov sull’insieme degli stati $E = \cur{RR, RB, BR, BB}$, in cui si tiene conto delle ultime due estrazioni effettuate. In particolare lo stato $RR$ rappresenta la vittoria di Rita, lo stato $BB$ quella di Bruno. Calcolare tutte le distribuzioni invarianti. Visualizzare con un grafico la variazione totale tra la densità marginale al tempo $t$ e il tempo successivo $t+1$, per $t=0,1,2,\ldots, 10$.