Simulations de chaînes de Markov

Li-Thiao-Té Sébastien

1 Chaînes à nombre d’états fini

On suppose dans cette section que les chaînes de Markov sont finies à espace d’états [1, n]. Pour les illustrations, on utilise les matrices suivantes :

P1 = [0 1 0;
      0 0 1;
      1 0 0] // matrice circulante
P2 = [ 0.9  0.1  0;
       0    0.9  0.1;
       0    0    1]
P3 = [0.4  0.6  0.0;
      0.2  0.5  0.3;
      0.1  0.7  0.2;]
// météo    1 neige   2 nuages   3 soleil

1.1 Tirages et trajectoires

La loi est fournie sous la forme d’un vecteur ligne proba. On tire un nombre r au hasard selon une loi uniforme et on compare la fonction de répartition à r :

function [x] = tirer(proba)
   // renvoie un nombre tiré au hasard selon la proba
   r = rand(1,1,'uniform');
   res = find([cumsum(proba)] > r);
   x = res(1,1)
endfunction;

La trajectoire de la chaine de Markov s’obtient en itérant les tirages aléatoires :

function [traj] = trajectoire(P,depart,nb_pas)
   // trajectoire d'une chaine de Markov
   traj = zeros(nb_pas,1)
   traj(1,1) = depart
   for i=2:nb_pas,
      traj(i,1) = tirer(P(traj(i-1,1),:))
   end;
endfunction;

Fig. 1 – Trajectoires de P1, P2, P3 2

1.2 Critère d’irréductibilité

Matrice potentielle On peut calculer simplement les sommes partielles de la matrice potentielle :

function [u] = potentielle(P,k)
   // matrice potentielle
   [n,m] = size(P);
   u = eye(n,m);
   temp = u;
   for i = 1:k,
      temp = temp * P;
      u = u + temp;
   end;
endfunction;

Une chaîne de Markov (finie) est irréductible ssi ∀x, y U(x, y) = +∞. Malheureusement, ce critère n’est pas effectif car on ne peut calculer que des sommes partielles de la matrice potentielle en général. Par contre, on peut s’intéresser à la matrice d’adjacence P̃ définie par :

P̃(x, y) = (P (x, y) > 0)

La matrice d’adjacence est une matrice de booléens et P̃(x, y) indique si il y a une arête entre x et y ou non. Alors P̃^{k}(x, y) indique si il y a un chemin entre x et y de longueur k. On définit alors une matrice Ũ :

Ũ(x, y) = ∑_{k} P̃^{k}(x, y)

Cette matrice est calculable car on peut se ramener aux chemins de longueur inférieure à n. En Scilab, les opérations matricielles ne sont pas définies sur les booléens, on utilise un calcul avec renormalisation à chaque étape :

function [x] = irreductible(P)
   [n,m] = size(P);
   pk = eye(n,m);
   u = pk;
   p = (P>0)*1;
   for i = 1:n,
      pk = ((pk * p)>0)*1;
      u = u + pk;
   end;
   x = prod((u > 0)*1);
endfunction;
-->irreductible(P1)
 ans  =    1.
-->irreductible(P2)
 ans  =    0.
-->irreductible(P3)
 ans  =    1.

1.3 Calcul de la période

On rappelle que la période est définie comme le pgcd des I(x) avec :

I(x) = {n > 0, P^{n}(x, x) > 0}

La période est donc le pgcd des longueurs des cycles dans le graphe correspondant à la chaîne de Markov. On peut se restreindre aux circuits dans le graphe (cycles entre des états distincts). Il suffit donc de considérer le pgcd de l’ensemble Ĩ = ⋃ Ĩ(x) avec :

Ĩ_{x} = {n ∈ [1, n], P̃^{n}(x, x) = T rue}

function [g] = periode(P)
   [n,m] = size(P);
   toutes_l = zeros(1,n); // initialisé à 0 car pgcd divise 0
   p = (P>0)*1;
   pk = eye(n,m);
   for k=1:n
      pk = ((pk * p)>0)*1;
      if (sum(eye(pk).*pk)>0) then
      toutes_l(1,k) = k;
      end;
   end
   g = gcd(uint32(toutes_l));
endfunction

1.4 Mesure invariante

Loi des Grands Nombres On utilise le corollaire 19.14.1 du théorème ergodique : si X est une chaîne de Markov récurrente positive de probabilité invariante π, alors on a :

\frac{1}{n} ∑_{k=0}^{n−1} P^{k}(y, x) ⟶ π(x)

On y reconnait une somme partielle de la matrice potentielle pour tout y. Afin de s’affranchir du choix arbitraire de y, on fait une moyenne sur tous les états de départ possibles :

function [mu] = proba(P,k)
   [n,m] = size(P);
   u = potentielle(P,k-1);
   mu = sum(u,'r') / (k * n)
endfunction

Théorème de Perron-Froebenius On écrit P^{n} ≃ 1.π, on peut donc obtenir π en faisant une moyenne sur les colonnes :

function [mu] = probap(P,k)
   [n,m] = size(P);
   mu = sum(P^k,'r') / n
endfunction;

Comparaison des vitesses de convergence On visualise la convergence des fréquences en fonction du nombre d’itérations :

function [] = plot_convergence(P,fin)
   [n,m] = size(P);
   temp = zeros(fin,2*n);
   u = eye(n,m);
   pn = u;
   for i=1:fin,
      pn = pn * P;
      u = u + pn;
      temp(i,:) = [(sum(u,'r') / (i * n)) (sum(pn,'r')/n)];
   end;
   msg = "LGN(" + strcat(string(1:n),")@LGN(") + ")
        @PF(" + strcat(string(1:n),")@PF(") + ")";
   plot2d(1:fin,temp,leg=msg);
   xtitle(["Convergence vers la mesure invariante" "Comparaison
    entre la Loi des Grands Nombres et le Théorème de
    Perron-Froebenius"],"trajectoire","");
endfunction;

Fig. 2 – Comparaison pour la chaîne de la météo.

Calcul exact Pour un nombre fini d’états, la mesure invariante est un simple vecteur, solution du système linéaire μ × P = μ. Ce système est équivalent à (P^{t} − I)μ^{t} = 0, mais (P^{t} − I) n’est pas inversible car le vecteur 1 est un vecteur propre de P^{t} associé à la valeur propre 1. On remplace donc l’une des contraintes par ∑ μ_{i} = 1.

function [mu] = probas(P)
   [n,m] = size(P);
   A = P' - eye(m,n);
   A(1,:) = ones(1,n);
   mu = (A\( [1;zeros(n-1,1)] ))';
endfunction;

2 Chaînes à nombre d’états dénombrable

2.1 Marches aléatoires

Exemples de marches aléatoires On s’est restreint aux marches aléatoires sur ℤ pour des raisons de simplicité. Pour simuler des champs de Gibbs par exemple, il faut sensiblement plus de travail :

P0 = [0.5 0 0.5] // marche aléatoire symétrique
P1 = [0.55 0 0.45]
P2 = [0.45 0 0.55]

Tirages aléatoires On utilise ici l’homogénéité spatiale des marches aléatoires. On se donne une probabilité sur [−n, n] codée par un vecteur ligne de taille 2n+1. On fait les tirages aléatoires comme précédemment, mais on enlève n + 1 au résultat afin d’obtenir un "offset" ou décalage :

function [x] = tirer(proba)
   // renvoie un nombre tiré au hasard selon la proba
   [n,m] = size(proba);
   r = rand(1,1,'uniform');
   res = find([cumsum(proba)] > r);
   x = res(1)-floor(m/2)-1
endfunction;
function [traj] = trajectoire(P,depart,nb_pas)
   // trajectoire d'une marche aléatoire
   traj = zeros(nb_pas,1)
   traj(1,1) = depart
   for i=2:nb_pas,
      traj(i,1) = traj(i-1,1) + tirer(P(1,:))
   end;
endfunction;

2.2 Convergence en loi

D’après le théorème ergodique, les probabilités empiriques pour une chaîne irréductible récurrente tendent vers :

– la probabilité invariante si la chaîne est récurrente positive

– 0 si la chaîne est récurrente nulle

Fig. 3 – Marche aléatoire symétrique sur ℤ

Fig. 4 – Distance à l’origine (écart-type) de la marche aléatoire symétrique sur ℤ. Elle est en \sqrt{n}.

A l’aide des fonctions suivantes, on visualise la convergence des fréquences empiriques sur une trajectoire :

function [x,y] = distribution(v,a,b)

// a et b sont des entiers
x = (a:b)';
y = zeros(x);
for i=1:(b-a+1),
  y(i,1) = length(find(v==x(i,1)));
end
n=length(find((v>=a) *1 .* (v<=b)));
y = y/n;
endfunction;
function [] = plot_convergence(P,depart,pas)
// convergence en loi d'une marche aléatoire
// affiche les lois successives (restreintes à [a b])
// évaluées aux instants donnés par le vecteur pas
l = length(pas);
tr = trajectoire(P, depart, pas(l));
a = min(tr);
b = max(tr);
x = (a:b)';
y = zeros(b-a+1,l);
for i=1:l,
   [u,v] = distribution(tr(1:pas(i),1),a,b);
   y(:,i) = v;
end
msg = "n = " + strcat(string(pas),"@n = ");
plot2d(x,y,leg=msg);
xtitle("Convergence en loi de la distribution empirique")
endfunction;

2.3 Système ALOHA

La chaîne de Markov du système ALOHA n’est malheureusement pas homogène en espace, on a donc une fois de plus redéfini une fonction de tirage aléatoire. Une fonction proba_aloha(i,a,nu) se charge de calculer le vecteur des probabilités des états suivants (pris par rapport à l’état courant entre -1 et le nombre maximal de nouveaux messages). Puis, on utilise la fonction tirer de la première section pour tirer l’état suivant.

function [proba] = proba_aloha(i,a,nu)
// renvoie l'état suivant à partir de i
[n,m] = size(a);
if i == 0 then
   proba = [ 0, a(1) + a(2), 0 , a(3:$) ];
else
   b = binomial(nu,i);
   proba = [ b(2)*a(1), (1-b(2))*a(1) + b(1)*a(2), (1-b(1))*a(2), a(3:$) ];
end
endfunction;
function [traj] = trajectoire(anu,depart,nb_pas)
// trajectoire d'un ALOHA instable

Fig. 5 – La marche aléatoire symétrique sur ℤ est récurrente nulle.

a = anu(1:($-1));
nu = anu($);
traj = zeros(nb_pas,1)
traj(1,1) = depart
for i=2:nb_pas,
   traj(i,1) = traj(i-1,1) + tirer(proba_aloha(traj(i-1,1), a, nu))-2
end;
endfunction;

Stabilisation Pour avoir une chaîne récurrente positive, on module la valeur de ν suivant le nombre de messages. Une condition suffisante de stabilité est :

λ < a_{0} exp \{\frac{a_{1}}{a_{0}} − 1\}

et l’on peut choisir alors :

ν(k) = \frac{a_{0} − a_{1}}{ka_{0} − a_{1}}

function [traj] = trajectoire(anu,depart,nb_pas)
// trajectoire d'un ALOHA stable
a = anu;
traj = zeros(nb_pas,1)
traj(1,1) = depart
for i=2:nb_pas,
   nu = (a(1)-a(2)) / (traj(i-1,1)*a(1) - a(2));
   traj(i,1) = traj(i-1,1) + tirer(proba_aloha(traj(i-1,1), a, nu))-2

end;
endfunction;

Fig. 6 – Trajectoire d’un ALOHA instable.

2.4 Files d’attentes

Pour la simulation des files d’attentes, on profite du générateur de nombres aléatoires de Scilab qui connaît déja les loi de Poisson et géométrique :

function [traj] = trajectoire(anu,depart,nb_pas)
// trajectoire d'une file d'attente
lambda = anu(1);
p = anu(2);
traj = zeros(nb_pas,1)
traj(1,1) = depart
for i=2:nb_pas,
   traj(i,1) = max(traj(i-1,1) - grand(1,1,'geom',p),0)
      + grand(1,1,'poi',lambda)
end;
endfunction;

Il suffit de vérifier que l’on traite en moyenne davantage de personnes que d’arrivées pour avoir la stabilité.

Fig. 7 – Convergence en loi d’un ALOHA instable. Fig. 8 – Trajectoire d’un ALOHA stable.

Fig. 9 – Convergence en loi d’un ALOHA stable. Fig. 10 – Trajectoire d’une file stable.

Fig. 11 – Convergence en loi d’une file stable. 13