DM 2 : Comparaisons de méthodes pour calculer l’état stationnaire

Méthode de travail

J’ai fait mon DM en avec d’autres personnes de la classe, il se peut donc que certaines fonctions se ressemblent.

Initialisation

  #Nombre d'etat dans les graphes
  N = 10;
  set.seed(42);

On prend un petit nombre d’états pour pouvoir éprouver les algorithmes sur le nombre de répétition.

Construction des graphes

  ligne = matrix(c(0),nrow = N, ncol = N);
  for(i in 2:N){
    ligne[i,i-1] = 1
    ligne[i-1, i] = 1
  }
  for(i in 1:N){
    ligne[i,] = ligne[i,]/sum(ligne[i,])
  }
  
  anneau = matrix(c(0), nrow = N, ncol = N);
  for(i in 2:N){
    anneau[i,i-1] = 1
    anneau[i-1, i] = 1
  }
  anneau[1,N] = 1
  anneau[N, 1] =1
   for(i in 1:N){
    anneau[i,] = anneau[i,]/sum(anneau[i,])
  }
  
  #Pour la sucette on reprend un anneau puis on le "coupe" et "recolle" aléatoirement.
  sucette = matrix(c(0), nrow = N, ncol = N);
  for(i in 2:N){
    sucette[i,i-1] = 1
    sucette[i-1, i] = 1
  }
  sucette[1,N] = 1
  sucette[N, 1] =1
  
  node1 = sample(1:N, 1);
  node2 = sample(1:N, 1);
  while(node1==node2 && sucette[node1,node2] == 1){ node2 = sample(1:N, 1);}
  sucette[node1,((node1+1)%%N)] = 0;
  sucette[((node1+1)%%N), node1] = 0;
  sucette[node1,node2] = 1;
  sucette[node2,node1] = 1;
  for(i in 1:N){
    sucette[i,] = sucette[i,]/sum(sucette[i,])
  }

Marche aléatoire

  #Marche aleatoire
  marcheAleatoire = function(nbNoeud, matrice, nombreDePas = 1000000){
    compteur  = rep(0, nbNoeud);
    etat = sample(1:nbNoeud,1);
    for(i in 1:nombreDePas){
      compteur[etat] = compteur[etat] +1;
      etat=sample(1:nbNoeud, 1, prob= matrice[etat,]);
    }
    compteur = compteur/nombreDePas
    compteur
  }

Ligne

  plot(marcheAleatoire(N,ligne))

On voit que pour une ligne, c’est plutot uniformément répartie (légérement supérieur à 1/10), sauf pour les deux extrémitées.

Anneau

  plot(marcheAleatoire(N,anneau))

On voit que pour un anneau, c’est très aléatoire. Augmentons le nombre de répétitions pour voir si on a atteind l’état stationnaire.

Anneau

  plot(marcheAleatoire(N,anneau,10000000))

On a donc quelque chose qui semble aléatoire. On se sait pas trop comment interpreter.

Sucette

  plot(marcheAleatoire(N,sucette))

Ici on a quelque chose qui semble stable et linéaire. Cependant vu que la sucette est générée aléatoirement ont sait pas non plus comment elle est et donc comment interpreter.

Vecteurs de probabilités

  #vecteur de proba arbitraire
  vecteurProba = function(nbNoeud, matrice, nombreDePas = 100000){
    proba = rep(1/nbNoeud, nbNoeud)
    for(i in 1:nombreDePas){
      proba = proba %*% matrice;
    }
    proba
  }

Ligne

  vecteurProba(N, ligne)
##            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.05555556 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
##           [,7]      [,8]      [,9]      [,10]
## [1,] 0.1111111 0.1111111 0.1111111 0.05555556
  plot(vecteurProba(N,ligne)[1,]);

On est en accord avec les résultats précedents. Les valeurs sont plus stables que pour l’algorithme précédent.

Anneau

  plot(vecteurProba(N,anneau)[1,]);

On voit enfin ce que j’espèrais voir sur le premier algorithme. Ici on est dans un état stationnaire.

Sucette

  plot(vecteurProba(N,sucette)[1,]);

On retrouve ici le même genre de graphique que pour l’algo précédent.

Veleur propre

  #Code récupéré sur http://mescal.imag.fr/membres/arnaud.legrand/teaching/2015/RICM4_EP.php 
  #pour la modelisation d'un cache web
  
  #Valeur propre ligne
  E=head(eigen(t(ligne)))
  head(E$values)
## [1] -1.0000000  1.0000000  0.9396926 -0.9396926 -0.7660444  0.7660444
  #Ok ici c'est la premiere qui vaut 1 qui nous interesse soit ici la deuxieme
  ssp=E$vectors[,2];
  ssp = Re(ssp/(sum(ssp)))
  options(digits=7);
  ssp
##  [1] 0.05555556 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111
##  [7] 0.11111111 0.11111111 0.11111111 0.05555556
  #Valeur propre anneau
  E=head(eigen(t(anneau)))
  head(E$values)
## [1]  1.000000  0.809017  0.809017  0.309017  0.309017 -0.309017
  #Ok ici c'est la premiere qui vaut 1 qui nous interesse soit ici la premiere
  ssp=E$vectors[,1];
  ssp = Re(ssp/(sum(ssp)))
  options(digits=7);
  ssp
##  [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
  #Valeur propre sucette
  E=head(eigen(t(sucette)))
  head(E$values)
## [1]  1.0000000 -0.9863613  0.9458172 -0.8794738  0.7891405 -0.6772816
  #Ok ici c'est la premiere qui vaut 1 qui nous interesse soit ici la deuxieme
  ssp=E$vectors[,2];
  ssp = Re(ssp/(sum(ssp)))
  options(digits=7);
  ssp
##  [1]  1.698414e+15 -3.350499e+15  3.212778e+15 -2.987421e+15  2.680574e+15
##  [6] -2.300609e+15  1.857889e+15 -1.364490e+15  8.338719e+14 -2.805078e+14
  valeurPropre = function(matrice){
    E=head(eigen(t(matrice)))
    #On cherche l'indice
    indice = which.max(E$values)
    ssp=E$vectors[,indice];
    ssp = Re(ssp/(sum(ssp)))
    options(digits=7);
    return(ssp)
    
  }
  plot(valeurPropre(ligne))

Ici on a donc les valeurs théoriques qui sont bien celles observées la dernière foi.

  plot(valeurPropre(anneau))

De même que précédement.

  plot(valeurPropre(sucette))

De même que précedement.

Analyse

Pour les résultats, on voit que la première méthode ne marche pas toujours (ou alors il faut un grand nombre de répétitions). Pour départager les deux autres nous allons augmenter le nombre d’état.

  N = 1000;

Ligne

  t1 = Sys.time();
  plot(vecteurProba(N,ligne)[1,]);

  t2= Sys.time();
  tvp = difftime(t2,t1)
  tvp
## Time difference of 8.599496 mins
  t1 = Sys.time();
  plot(valeurPropre(ligne));

  t2= Sys.time();
  tvpropre = difftime(t2,t1)
  tvpropre
## Time difference of 9.093838 secs

L’utilisation des valeurs propres et bien plus rapide ! ##Anneau

  t1 = Sys.time();
  plot(vecteurProba(N,anneau)[1,]);

  t2= Sys.time();
  tvp = difftime(t2,t1)
  tvp
## Time difference of 8.68155 mins
  t1 = Sys.time();
  plot(valeurPropre(anneau));

  t2= Sys.time();
  tvpropre = difftime(t2,t1)
  tvpropre
## Time difference of 1.777416 secs

On voit ici que je m’étais trompé sur les anneaux. c’est pas une ligne mais une parabolle. Maintenant que je le vois ca parait logique. Je suis quand même étonné que ce soit centré (ca voudrait dire qu’on part pile poile du milieu).

Ici aussi, les valeurs propres sont plus efficaces.

Sucette

  t1 = Sys.time();
  plot(vecteurProba(N,sucette)[1,]);

  t2= Sys.time();
  tvp = difftime(t2,t1)
  tvp
## Time difference of 8.530877 mins
  t1 = Sys.time();
  plot(valeurPropre(sucette));

  t2= Sys.time();
  tvpropre = difftime(t2,t1)
  tvpropre
## Time difference of 8.575259 secs

Ce dernier test prouve que les calcul par valeur propre est bien plus efficace que les autres.