In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import networkx as nx
Definim o functie care genereaza graful de tranzitie al unui lant Markov cu $m$ stari, codificate $0,1, \ldots, m-1$, avand matricea de tranzitie $Q$:
In [2]:
def GraphTr(Q):
G=nx.from_numpy_matrix(Q, create_using=nx.DiGraph())
nx.draw(G, node_color='b', alpha=0.3)
Functia EquilD calculeaza distributia de echilibru a unui lant ireductibil si aperiodic
ca vector propriu corespunzator valorii proprii 1 a matricii de tranzitie, transpusa:
In [3]:
def EquilD(Q):
Lam, V=np.linalg.eig(Q.T)
absLam=np.abs(Lam)
j=np.argmax(absLam)#indicele valorii proprii ce are valoarea absoluta maxima
v=V[:,j] # extragem vectorul propriu corespunzator lui 1;
# calculam distributia de echilibru impartind fiecare coordonata a lui v la
#suma coordonatelor:
return v/np.sum(v)
In [4]:
#Exemplu de definire a grafului de tranzitie a unui Lant Markov
#Constructia distributiei de echilibru ca vector propriu corespunzator valorii 1
# Generarea sirului pi_n, cu n suficient de mare
pia=np.array([0.3, 0.45, 0.25]) #distributia initiala de probabilitate
Q=[[0, 3./4, 1./4],
[1./3, 0, 2./3],
[0, 1.0, 0]]
Q=np.array(Q)# am introdus elementele lui Q
#intr-o lista de liste si apoi am convertit la array
#Distributia de echilibru calculata ca un vector propriu al matricii Q transpus:
GraphTr(Q)
pie=EquilD(Q)
#Distributia de echilibru ca o aproximatie a limitei sirului pi_n
n=200
for i in range (n):
pia=np.dot(pia,Q)
print 'Distributia de echilibru calculata ca vector proriu este:\n', pie
print 'Distributia aproximata de sirul pi_n este:\n',pia.round
print 'norma diferentei dintre cele doua:\n', np.linalg.norm(pie-pia)
Un lant Markov cu 4 stari are matricea de tranzitie $Q$.
In [5]:
Q=[[0.4, 0, 0.3, 0.3], [0, 0.5, 0.3, 0.2], [0, 0.3, 0.3, 0.4], [0.2, 0.2, 0.45, 0.15]]
Q=np.array(Q)
print Q
GraphTr(Q)
Precizam ca desi matricea de tranzitie are elemente nenule pe diagonala principala, networkx nu
afiseaza si buclele $(0,0)$, $(1,1)$, $(2,2)$, $(3,3)$ (nu este implementata afisarea grafica a buclelor)
In [6]:
pii=np.array([0.25, 0.25, 0.25, 0.25])#distributia initiala de probabilitate
Calculam si afisam distributiile de probabilitate $\pi_n$ cu $n=\overline{0,9}$. Apoi calculam fara a afisa $\pi_n$ cu $n=\overline{10,49}$. Incepand de la $\pi_{50}$ afisam din nou:
In [7]:
for i in range(10):
print pii.round(4)
pii=np.dot(pii, Q)
Observam ca pe masura ce $n$ creste, se schimba probabilitatile de vizita a nodurilor 0,1,2,3. Initial, fiecare nod avea aceeasi probabilitate de a fi vizitat. Dupa 9 pasi probabilitatea ca nodul 0 sa fie vizitat a scazut de la 0.25 cat era initial, la 0.0879, in timp ce pentru nodurile 1,2, 3 a crescut.
In [8]:
for i in range(10,50):
pii=np.dot(pii, Q)
In [9]:
for i in range(50, 60):
print pii.round(4)
pii=np.dot(pii, Q)
Se observa ca distributia de probabilitate a vizitelor nodurilor s-a stabilizat si deci distributia de echilibru experimentala calculata cu 4 zecimale este $\pi=[0.0878, 0.3091, 0.3395, 0.2635]$
Suma coordonatelor este:
In [10]:
print np.sum(pii)
si deci $\pi$ este vector probabilist. $\pi(k)$, $k=0,1,2,3$, este probabilitatea de vizita asimptotica a nodului $k$.
Sa calculam acum distributia de echilibru ca vector propriu corespunzator valorii proprii 1, a matricii $Q^T$:
In [11]:
pit= EquilD(Q)
print pit
In [12]:
print np.sum(pit)
Norma diferentei dintre cele doua distributii de echilibru, pii, dedusa in urma plimbarii aleatoare pe graf, conform matricii de tranzitie $Q$ si cea teoretica (ca vector propriu)
este:
In [13]:
print np.linalg.norm(pii-pit)
adica foarte mica.
Luand o alta distributie initiala de probabilitate si afisand sirurile asociate ca in cazul precedent avem:
In [14]:
pii=np.array([0.2, 0.3, 0.4, 0.1])
In [15]:
for i in range(10):
print pii.round(4)
pii=np.dot(pii, Q)
In [16]:
for i in range(10,50):
pii=np.dot(pii, Q)
In [17]:
for i in range(50, 60):
print pii.round(4)
pii=np.dot(pii, Q)
Remarcam ca si in acest caz pentru $n\geq 50$, $\pi_n$ s-a stabilizat la $\pi_n=[ 0.0878, 0.3091, 0.3395, 0.2635]^T$, adica la fel ca in cazul in care am avut distributia uniforma ca distributie initiala de probabilitate.
Acest exemplu ilustreaza ca daca un lant Markov contine traiectorii periodice, atunci desi matricea $Q^T$ admite valoarea proprie 1, vectorul propriu corespunzator nu e sigur ca are toate coordonatele nenule si de acelasi semn. In plus in acest caz fixand o distributie initiala de probabilitate, sirul corespunzator $\pi_n$, cu $\pi_n^T=\pi_{n-1}^TQ$, nu e convergent.
In [18]:
Q=[[0, 0, 1, 0, 0, 0, 0, 0], [0,0,0,1, 0,0,0,0], [0,1, 0,0, 0, 0, 0, 0],
[1,0,0,0, 0, 0, 0, 0], [0.6, 0, 0, 0, 0, 0.4, 0, 0], [0.2, 0, 0, 0, 0, 0, 0.8, 0],
[0.25, 0, 0, 0, 0.2, 0.4, 0, 0.15], [0, 0, 0, 0.65, 0, 0, 0.35, 0]]
Q=np.array(Q)
In [19]:
print Q
Graful asociat este:
In [23]:
GraphTr(Q)
Remarcam din graf ca exista un drum periodic de perioada 4, si anume: $0\to2\to1\to3\to 0$.
Fie $\pi=[0,0,0,0,0.25, 0.25, 0.25,0.25]$ distributia initiala de probabilitate. Probabilitatea ca miscarea aleatoare pe graf sa inceapa din nodurile 0,1,2,3 este 0 si respectiv in mod uniform poate incepe din nodurile 4,5,6,7.
Graful (matricea de tranzitie) nu este ireductibil(a) deoarece nu exista drum de arce intre nodul $3$ de exemplu si nodul $5$.
Valorile si vectorii proprii ale matricii $Q^T$ sunt:
In [24]:
Lam, V=np.linalg.eig(Q.T)
print Lam.round(2)
Mai precis valorile proprii sunt $-1, i, -i, 1, 0.68, -0.49, -0.19, 0$ iar valorile lor absolute:
In [25]:
absLam=np.abs(Lam)
print absLam.round(3)
Remarcam ca nu exista un unic maxim al valorilor absolute, ci patru.
Matricea $Q$ fiind stochastica pe linii are sigur valoarea proprie 1 si vectorul ei propriu este:
In [26]:
v=V[:,3]
print v.round(3)
Contrar cazului cand matricea este ireductibila si aperiodica, in acest caz vectorul propriu corespunzator valorii 1, nu are toate coordonatele nenule.
Sa investigam acum distributiile de probabilitate la momentele $n=1,2,\ldots$, cand distributia initiala este distributia uniforma:
In [27]:
pii=0.125*np.ones(8)
print pii
In [28]:
for i in range(35):
print pii.round(3)
pii=np.dot(pii, Q)
Remarcam ca in loc ca sirul $\pi_n$ sa tinda la o limita el are 4 subsiruri ce tind la limite diferite (ultimele 4 distributii se repeta dupa un rang $n$)
Luam acum o alta distributie initiala de probabilitate: $\pi=[0,0,0,0,1,0,0,0]^T$ si calculam sirul distributiilor $\pi_n$ asociat:
In [29]:
pii=np.array([0,0,0,0,1,0,0,0])
for i in range(35):
print pii.round(3)
pii=np.dot(pii,Q)
Remarcam din nou ca ultimele 4 distributii se repeta incepand de la un anumit rang. Prin urmare exista 4 subsiruri convergente ale sirului $(\pi_n)$ si limitele depind de distributia initiala.
Acest exemplu ilustreaza ce importante sunt ipotezele de matrice ireductibila si aperiodica in asigurarea unicitatii distributiei de echilibru si a coordonatelor nenule si de acelasi semn pentru vectorul propriu al lui $Q^T$, corespunzator valorii 1.
Analizam doua lanturi Markov absorbante ce modeleaza doua metode de forwardare a pachetelor de informatie in relele wireless (vezi descrierea in lista de probleme relativ la lanturi Markov).
In [30]:
from IPython.display import Image
Image(filename='Imags/MarkovAbsC.png')
Out[30]:
Matricea de tranzitie pentru o astfel de retea cu $n$ noduri si probabilitatea forwardarii cu succes dintr-un nod i spre nodul i+2, $i=\overline{0, n-3}$, egala cu p, este generata de functia:
In [31]:
def MatrixC(n, p):
Q=[]
for i in range(n):
Q.append([])
for i in range(n):
[Q[i].append(0) for j in range(n)]
for i in range(n-2):
Q[i][i:i+3]=[(1-p)**2, p*(1-p), p]
Q[n-2] [n-2:n]=[(1-p),p]
Q[n-1][n-1]=1
return np.array(Q)
Pentru cazul $n=5$ vizualizat in imaginea de mai sus si $p=0.85$, matricea de tranzitie este:
In [32]:
Q=MatrixC(5, 0.85)
print Q
Se observa si din graf si din matricea de tranzitie ca nodul $n=4$, destinatia pachetului, este absorbant. Aceasta este o forma standard a matricii de tranzitie oarecum opusa celei discutate la curs. In acest exemplu primele noduri sunt tranzitorii si ultimul este absorbant. In aceasta scriere matricea tranzitorie $T$ este formata din elementele de intersectie ale liniilor si coloanelor $0,1,2,3$:
In [33]:
T=Q[:4,:4]
print T
Matricea fundamentala $N$ este:
In [34]:
I=np.eye(4)
N=np.linalg.inv(I-T)
print N.round(4)
Suma elementelor de pe linia 0 a matricii N este timpul mediu petrecut de un pachet in retea, inainte de ajunge la destinatie (nodul 4).
In [35]:
timp=np.sum(N[0])
print timp
Protoculul de forwardare este in acest caz ilustrat pe o retea de 5 noduri, indexate $0,1,2,3,4$:
In [36]:
Image(filename='Imags/MarkovAbsNon.png')
Out[36]:
Matricea de tranzitie este in acest caz:
In [37]:
def MatrixNonC(n,p):
Q=[]
for i in range(n):
Q.append([])
for i in range(n):
[Q[i].append(0) for j in range(n)]
for i in range(n-2):
Q[i][i:i+3]=[(1-p), 0, p]
Q[n-2][n-2:n]=[(1-p),p]
Q[n-1][n-1]=1.0
return np.array(Q)
In [38]:
print MatrixNonC(5,0.85)
In continuare comparam performanta teoretica si experimentala a celor doua metode de forwardare a pachetelor intr-o retea de $n=100$ noduri si probabilitatea, $p$, de forwardare de la nodul $i$ la nodul $i+2$, $i=\overline{0,n-3}$, mai apropiat de destinatie, egala cu $p=0.85$.
O masura a performantei teoretice este timpul mediu petrecut in retea inainte de ajunge la destinatie. Deci dintre cele doua protocoale va fi mai performant cel pentru care suma elementelor de pe linia $0$ a matricii fundamentale, $N$, mai mica.
Performanta experimentala se evalueaza simuland mai multe forwardari (nr) si dupa fiecare calculand numarul de noduri disticte sau nu, vizitate de pachet inainte de ajunge la destinatie.
Dupa cele nr forwardari simulate, se calculeaza media aritmetica a numarului de noduri vizitate.
Fiecare nod este vizitat intr-o unitate de timp si deci aceasta medie este timpul mediu petrecut in retea.
Pentru a calcula matricea fundamentala a lantului, $N=(I-T)^{-1}$, apelam functia
inv din scipy.linalg, deoarece este mai performanta decat functia np.linalg.inv.
Pentru a simula forwardarea, adica lantul Markov, definim functia:
In [39]:
def simul(pr):
k=0
F=pr[0]
u=np.random.random()
while(u>F):
k+=1
F=F+pr[k]
return k
In [40]:
import scipy.linalg as spl
n=100
p=0.85
Q=MatrixC(n,p)#matricea de tranzitie conform primului protocol de forwardare
T=Q[:n-1,:n-1]
I=np.eye(n-1)
NC=spl.inv(I-T)
perfC=np.sum(NC[0])
S=MatrixNonC(n,p)
T=S[:n-1,:n-1]
NN=spl.inv(I-T)
perfNC=np.sum(NN[0])
print 'timpul mediu petrecut in retea de un pachet forwardat cf protocolului CARQ', perfC
print 'timpul mediu petrecut in retea de un pachet forwardat cf protocolului NCARQ', perfNC
Deci protocolul CARQ este mai performant, din punct de vedere teoretic. Sa verificam si experimental. mai intai generam o traiectorie a unui pachet conform fiecarui protocol si o afisam:
In [41]:
def forward(Q):
n=Q.shape[0]
nod=0
pr=Q[0]# pachetul porneste din nodul 0
traj=[nod]
while(nod!=n-1):
nod=simul(pr)
traj.append(nod)
pr=Q[nod]
return traj
In [42]:
traj1=forward(Q)
traj2=forward(S)
In [43]:
print 'traiectoria pachetului conform protocolului CARQ\n', traj1
print 'nr de noduri vizitate in fiecare unitate de timp:', len(traj1)-1
In [44]:
print 'traiectoria pachetului conform protocolului NCARQ\n', traj2
print 'nr de noduri vizitate in fiecare unitate de timp:', len(traj2)-1
Forwardam 50 pachete si calculam media numarului de noduri vizitate conform fiecarui protocol:
In [45]:
nr=50
lentr1=0
lentr2=0
for k in range(nr):
traj1=forward(Q)
traj2=forward(S)
lentr1+=len(traj1)-1
lentr2+=len(traj2)-1
print 'Media timpului petrecut in retea de un pachet forwardat cf CARQ este: ', float(lentr1)/nr
print 'Media timpului petrecut in retea de un pachet forwardat cf NCARQ este: ', float(lentr2)/nr
Media teoretica era respectiv 54.50, respectiv 58.82.
Prin urmare atat teoretic cat si experimental protocolul CARQ este mai performant.
In [46]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom.css", "r").read()
return HTML(styles)
css_styling()
Out[46]: