J’ai fait mon DM en avec d’autres personnes de la classe, il se peut donc que certaines fonctions se ressemblent.
#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.
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 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
}
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.
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.
plot(marcheAleatoire(N,anneau,10000000))
On a donc quelque chose qui semble aléatoire. On se sait pas trop comment interpreter.
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.
#vecteur de proba arbitraire
vecteurProba = function(nbNoeud, matrice, nombreDePas = 100000){
proba = rep(1/nbNoeud, nbNoeud)
for(i in 1:nombreDePas){
proba = proba %*% matrice;
}
proba
}
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.
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.
plot(vecteurProba(N,sucette)[1,]);
On retrouve ici le même genre de graphique que pour l’algo précédent.
#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.
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;
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.
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.