Le catene di Markov sono dei processi casuali discreti che compaiono spesso in fisica, ad esempio in fisica statistica e nelle simulazioni Monte Carlo. In questo articolo vedremo il più semplice esempio di catena di Markov, in cui il sistema ha solo due stati possibili.
Come abbiamo detto, il nostro sistema ha solo due stati possibili, ${A,B}$. Possiamo pensarli come “acceso” e “spento”, “0” e “1”, o come volete, però qui li chiamiamo $A$ e $B$. Lo stato al tempo $t$ della catena lo chiamiamo $\psi_t$. Poiché questo è un processo di Markov, $\psi_t$ dipende solo dallo stato precedente $\psi_{t-1}$. Questa è la proprietà che definisce un processo di Markov: non ha memoria di ciò che è avvenuto in passato. Se sono al tempo $t$ e devo andare al tempo $t+1$, l’unica cosa che conta è lo stato in cui mi trovo adesso.
La matrice di transizione
Pertanto, ad ogni tempo $t$, se $\psi_t = A$ ci sarà una certa probabilità $p$ che $A \to B$, cioè $\psi_{t+1} = B$. Poiché non c’è nessun’altra possibilità, la probabilità di $A \to A$ sarà $1-p$. Alla stessa maniera, se al tempo $t$ siamo in $B$, avremo una certa probabilità $q$ di passare in $A$, e una certa probabilità $1-q$ di rimanere in $B$. In linea di principio $p \neq q$, dato che non c’è nessuna ragione per cui debbano essere uguali.
Possiamo quindi salvare queste informazioni in una matrice di transizione $T$. Se chiamiamo $A=1$ e $B=2$ allora $T_{ij}$ è la probabilità di $j \to i$ ad ogni passo temporale. Questo non è un refuso, ma ricordiamo che la matrice di transizione è davvero definita “al contrario”. Pertanto
$$T = \begin{pmatrix} 1-p & q \\ p & 1-q \end{pmatrix}$$
La matrice di transizione è definita in questa maniera perché ci permette di fare quanto segue. Se volessimo sapere qual è la probabilità di $j \to i$ in $k$ passi, cioè partendo da $\psi_t = j$, vogliamo sapere la probabilità che $\psi_{t+k} = i$, ci basta moltiplicare la matrice di transizione con se stessa ed estrarre l’elemento corrispondente, $\pqty{T^k}_{ij}$.
La distribuzione stazionaria
Immaginiamo di partire da uno stato $\mathbf{x}_0$, cioè un vettore che ci dice la probabilità di trovarci inizialmente in $A$ e in $B$. Allora dopo $k$ passi della catena di Markov le probabilità di trovarci in $A$ e in $B$ rispettivamente saranno date da $T^k \mathbf{x}_0$. Supponiamo che ad un certo punto il processo diventi stazionario, cioè le probabilità diventino invariate ad ogni iterazione. Tipicamente ciò avviene nel limite $k\to \infty$. La distribuzione stazionaria $\boldsymbol{\pi}$ a questo punto soddisfa
$$T \boldsymbol{\pi} = \boldsymbol{\pi}$$
che è esattamente l’affermazione che $\boldsymbol{\pi}$ rimane invariata ad ogni iterazione. Ovvero $\boldsymbol{\pi}$ è l’autovettore con autovalore $1$ della matrice di transizione. Un calcolo diretto ci mostra che $T$ ha autovalori $1$ e $\lambda = 1-p-q$ e che quindi
$$\boldsymbol{\pi} = \frac{1}{p+q}\begin{pmatrix} q \\ p \end{pmatrix}$$
dove abbiamo normalizzato $\boldsymbol{\pi}$ in modo tale che la somma dei suoi componenti (le probabilità stazionarie) faccia $1$ come dovrebbe.
La distribuzione stazionaria può anche essere ottenuta come limite del processo di Markov, ed è così che viene spesso calcolata in pratica:
$$\boldsymbol{\pi} = \lim_{k\to \infty} T^k \mathbf{x}_0$$
dove $\mathbf{x}_0$ è una distribuzione iniziale qualsiasi. Infatti, nelle situazioni più interessanti di utilizzo lo spazio degli stati è enorme e quindi diagonalizzare $T$ può essere impossibile in pratica. Quest’ultima formula è valida poiché abbiamo
$$T^k \mathbf{x}_0 = T T^{k-1} \mathbf{x}_0$$
e quindi prendendo il limite $k \to \infty$ sia $T^k \mathbf{x}_0 \to \boldsymbol{\pi}$ sia $T^{k-1} \mathbf{x}_0 \to \boldsymbol{\pi}$, così riproduciamo la definizione di $\boldsymbol{\pi}$.
In questo caso semplificato possiamo però trovare esplicitamente $T^k$ e capire come ci avviciniamo al limite $k \to \infty$. Diagonalizzando $T= S D S^{-1}$ anche tramite Wolfram Alpha, otteniamo $T^k=S D^k S^{-1}$ e quindi
$$T^k = \frac{1}{p+q} \begin{pmatrix} q & q \\ p & p \end{pmatrix} + \frac{\lambda^k}{p+q} \begin{pmatrix} p & -q \\ -p & q \end{pmatrix}$$
dove $\lambda = 1-p-q$ è l’altro autovalore di $T$. Poiché $0 < p, q <1$, $\abs{\lambda} < 1$ e quindi nel limite $k \to \infty$ rimane solo il primo termine. Tuttavia, in ogni caso pratico non possiamo davvero prendere il limite $k\to \infty$ ma solo $k$ grande; ciò significa in particolare che il valore di $\lambda$ è molto importante. Se $\lambda$ è vicino ad $1$ ci vorrà un sacco di tempo prima di raggiungere lo stato stazionario; se invece $\lambda$ è piccola, la convergenza sarà rapida.
Vediamo cosa succede applicando $T^k$ ad un qualsiasi vettore $\mathbf{x}_0$ delle probabilità iniziali. Queste devono sommare ad $1$, per cui scriviamo $\mathbf{x}_0 = \begin{pmatrix} a \\ 1-a \end{pmatrix}$. Pertanto otteniamo
$$T^k \mathbf{x}_0 = \underbrace{\frac{1}{p+q} \begin{pmatrix} q \\ p \end{pmatrix}}_{\boldsymbol{\pi}} + \lambda^k\pqty{a-\frac{q}{p+q}} \begin{pmatrix} 1 \\ -1 \end{pmatrix}$$
Il primo termine è semplicemente la distribuzione stazionaria, per cui come avevamo detto $T^k \mathbf{x}_0 \to \boldsymbol{\pi}$ per $k \to \infty$. Notiamo che questo primo termine è indipendente dallo stato iniziale $\mathbf{x}_0$. Il secondo termine, invece, dipende da $a$, cioè dalla distribuzione iniziale di probabilità, e misura di quanto ci discostiamo dalla distribuzione stazionaria per $k$ finito.