\documentclass[a4paper,10pt]{article}
\usepackage[french]{babel}
\usepackage[T1]{fontenc}
\usepackage{fancybox}
\usepackage[fleqn]{amsmath}
\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{bbm}
\usepackage{graphicx}

\newcommand{\E}[2][]{\ensuremath{\mathbb{E}_{#1} \left[#2 \right]}}
\newcommand{\Proba}[2][]{\ensuremath{\mathbb{P}_{#1} \left[#2 \right]}}
\newcommand{\ud}{\mathrm{d}}

%\setlength{\textwidth}{14cm}

%opening
\title{Simulations de chaînes de Markov}
\author{Li-Thiao-Té Sébastien}
\date{}

\begin{document}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{Chaînes à nombre d'états fini}

\section{Chaînes à nombre d'états fini}

\paragraph*{}
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 :
\begin{verbatim}
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
\end{verbatim}

\subsection{Tirages et trajectoires}
\paragraph*{}
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$ :
\begin{verbatim}
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;
\end{verbatim}
\paragraph*{}
La trajectoire de la chaine de Markov s'obtient en itérant les tirages aléatoires :
\begin{verbatim}
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;
\end{verbatim}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.7\textwidth]{trajectoire1.eps}
\includegraphics[width=.7\textwidth]{trajectoire2.eps}
\includegraphics[width=.75\textwidth]{trajectoire3.eps}
\caption{Trajectoires de P1, P2, P3 }
\end{center}
\end{figure}

\subsection{Critère d'irréductibilité}
\paragraph*{Matrice potentielle}
On peut calculer simplement les sommes partielles de la matrice potentielle :
\begin{verbatim}
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;
\end{verbatim}
\paragraph*{}
Une chaîne de Markov (finie) est irréductible ssi $\forall x,y \quad U(x,y) = + \infty$. 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 $\tilde{P}$ définie par : \[\tilde{P}(x,y) = (P(x,y) > 0)\]
\paragraph*{}
La matrice d'adjacence est une matrice de booléens et $\tilde{P}(x,y)$ indique si il y a une arête entre $x$ et $y$ ou non. Alors $\tilde{P}^k(x,y)$ indique si il y a un chemin entre $x$ et $y$ de longueur $k$. On définit alors une matrice $\tilde{U}$ :
\[\tilde{U}(x,y) = \sum_k \tilde{P}^k(x,y) \]
\paragraph*{}
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 :
\begin{verbatim}
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.
 \end{verbatim}

\subsection{Calcul de la période}
\paragraph*{}
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 $\tilde{I} = \bigcup \, \tilde{I(x)}$ avec :
\[\tilde{I}_{x} = \{ n \in [1,n], \tilde{P}^n (x,x) = True\} \]
\begin{verbatim}
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
\end{verbatim}

\subsection{Mesure invariante}
\paragraph*{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 $\pi$, alors on a :
\[ \frac{1}{n} \sum_{k=0}^{n-1} P^k(y,x) \longrightarrow \pi(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 :
\begin{verbatim}
function [mu] = proba(P,k)
   [n,m] = size(P);
   u = potentielle(P,k-1);
   mu = sum(u,'r') / (k * n)
endfunction
\end{verbatim}
\paragraph*{Théorème de Perron-Froebenius}
On écrit $P^n \simeq 1.\pi$, on peut donc obtenir $\pi$ en faisant une moyenne sur les colonnes :
\begin{verbatim}
function [mu] = probap(P,k)
   [n,m] = size(P);
   mu = sum(P^k,'r') / n
endfunction;
\end{verbatim}
\paragraph*{Comparaison des vitesses de convergence}
On visualise la convergence des fréquences en fonction du nombre d'itérations :
\begin{verbatim}
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;
\end{verbatim}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=.95\textwidth]{cv3.eps}
\caption{Comparaison pour la chaîne de la météo.}
\end{center}
\end{figure}
\paragraph*{Calcul exact} Pour un nombre fini d'états, la mesure invariante est un simple vecteur, solution du système linéaire $\mu \times P = \mu$. Ce système est équivalent à $(P^t - I) \mu^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 $\sum \mu_i = 1$.
\begin{verbatim}
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;
\end{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{Chaînes à nombre d'états dénombrable}
\section{Chaînes à nombre d'états dénombrable}

\subsection{Marches aléatoires}
\paragraph*{Exemples de marches aléatoires}
On s'est restreint aux marches aléatoires sur $\mathbb{Z}$ pour des raisons de simplicité. Pour simuler des champs de Gibbs par exemple, il faut sensiblement plus de travail :
\begin{verbatim}
P0 = [0.5 0 0.5] // marche aléatoire symétrique
P1 = [0.55 0 0.45]
P2 = [0.45 0 0.55]
\end{verbatim}

\paragraph*{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 :
\begin{verbatim}
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;
\end{verbatim}
\begin{verbatim}
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;
\end{verbatim}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.9\textwidth]{marche1.eps}
\caption{Marche aléatoire symétrique sur $\mathbb{Z}$}
\end{center}
\end{figure}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.9\textwidth]{marche2.eps}
\caption{Distance à l'origine (écart-type) de la marche aléatoire symétrique sur $\mathbb{Z}$. Elle est en $\sqrt{n}$.}
\end{center}
\end{figure}

\subsection{Convergence en loi}
\paragraph*{}
D'après le théorème ergodique, les probabilités empiriques pour une chaîne irréductible récurrente tendent vers :
\begin{itemize}
\item la probabilité invariante si la chaîne est récurrente positive
\item 0 si la chaîne est récurrente nulle
\end{itemize}
\paragraph*{}
A l'aide des fonctions suivantes, on visualise la convergence des fréquences empiriques sur une trajectoire :
\begin{verbatim}
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;
\end{verbatim}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.9\textwidth]{distrib1.eps}
\caption{La marche aléatoire symétrique sur $\mathbb{Z}$ est récurrente nulle.}
\end{center}
\end{figure}

\subsection{Système ALOHA}
\paragraph*{}
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 \verb|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 \verb|tirer| de la première section pour tirer l'état suivant.
\begin{verbatim}
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
   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;
\end{verbatim}
\paragraph*{Stabilisation}
Pour avoir une chaîne récurrente positive, on module la valeur de $\nu$ suivant le nombre de messages. Une condition suffisante de stabilité est :
\[ \lambda < a_0  \exp \left\{ \frac{a_1}{a_0} -1 \right\} \]
et l'on peut choisir alors :
\[ \nu(k) = \frac{a_0 - a_1}{k a_0 - a_1} \]
\begin{verbatim}
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;
\end{verbatim}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=1.\textwidth]{aloha1.eps}
\caption{Trajectoire d'un ALOHA instable.}
\end{center}
\end{figure}
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=1.\textwidth]{aloha2.eps}
\caption{Convergence en loi d'un ALOHA instable.}
\end{center}
\end{figure}
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=1.\textwidth]{aloha3.eps}
\caption{Trajectoire d'un ALOHA stable.}
\end{center}
\end{figure}
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=1.\textwidth]{aloha4.eps}
\caption{Convergence en loi d'un ALOHA stable.}
\end{center}
\end{figure}


\subsection{Files d'attentes}
\paragraph*{}
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 :
\begin{verbatim}
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;
\end{verbatim}

\paragraph*{}
Il suffit de vérifier que l'on traite en moyenne davantage de personnes que d'arrivées pour avoir la stabilité.
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.65\textwidth]{file1.eps}
\caption{Trajectoire d'une file stable.}
\end{center}
\end{figure}
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.65\textwidth]{file2.eps}
\caption{Convergence en loi d'une file stable.}
\end{center}
\end{figure}



\end{document}
