DM: Sujet 2 comparaison de méthode et influence de la forme de la chaine de Markov
Delise Antoine
9 mars 2017
Partie 1: Les chaines de test
Tout d’abord, contruisons des fonctions afins de générer différentes formes de chaines: (Celà demande n’est qu’une simple écriture de matrice je ne pense qu’il soit nécessaire de commenter plus)
Line
Mline = function(N) {
M = matrix(data=0, nrow = N, ncol = N)
for(i in 1:(N-1)) {
M[i, i+1] = 1
M[i+1, i] = 1
M[i,] = M[i,] / sum(M[i,])
}
return(M)
}Circle
MCircle = function(N) {
M = matrix(data=0, nrow = N, ncol = N)
for(i in 1:N) {
M[i, 1+(i%%N)] = 1
M[i, ((i-2)%%(N))+1] = 1
M[i,] = M[i,] / sum(M[i,])
}
return(M)
}Sucette
MSucette = function(N) {
#afin de pas s'embêter on créer une sucette de 2N en reprenant le code de line N et circle N
#après réflexion on aurait aussi pu créer une line N et joindre une extrémité à un noeud intérieur, on ne cherche pas la performance de la génération donc je laisse comme celà.
N = N / 2
M = matrix(data=0, nrow = N*2, ncol = N*2)
for(i in 1:N) {
M[i, 1+(i%%N)] = 1
M[i, ((i-2)%%(N))+1] = 1
#M[i,] = M[i,] / sum(M[i,])
}
M[N,N+1] = 1
for(i in (N+1):(2*N-1)) {
M[i, i+1] = 1
M[i, i-1] = 1
#M[i,] = M[i,] / sum(M[i,])
}
M[2*N,2*N-1] = 1
for(i in 1:(2*N)) {
M[i,] = M[i,] / sum(M[i,])
}
return(M)
}Partie 2: L’état stationnaire
Par itération
Pour commencer nous allons calculer l’état stationnaire de cette chaine de markov de la manière la plus simple en itérant la matrice de transition sur une vecteur de probabilité, celà revient à une simple multiplication (on pourra voir que cette méthode converge bien assez rapidement sur des matrices résonnable donc il n’est pas nécessaire de chercher à accélérer le calcul par des méthode de trigonalisation et de puissance).
Les résultats de ces test seront mieux exposés et commentés dans la troisième partie
#
N = 12
M = MSucette(12)
P = matrix(data=1, nrow = 1, ncol = N)
P[1,] = P[1,] / sum(P[1,])
P## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 for(i in 1:1000) {
P = P %*% M
}
P## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.125
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.04166667Version fonctionnelle
Afin de faire des comparaisons plus simplement dans la dernière partie on créer une version fonctionnelle de ce code
#
testMul = function(M, K = 10000) {
N = length(M[1,]) # On lit la taille sur la première ligne, à l'origine à cause de de la sucette 2N, et puis ça retire un paramêtre
P = matrix(data=1, nrow = 1, ncol = N)
P[1,] = P[1,] / sum(P[1,])
P
for(i in 1:K) {
P = P %*% M
}
return(P)
}Par vecteur de probabilité
La deuxième méthode, elle aussi assez simple à réalisé consiste à faire se promener un pion sur la chaine de markov et de calculer ses fréquences d’apparitions
M = Mline(6)
N = length(M[1,])
X = 1
P = matrix(data=0, nrow = 1, ncol = N)
for(i in 1:100000) {
P[X] = P[X] + 1
X = sample(1:N, size=1, prob = M[X,])
}
P = P / sum(P)
P## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.10152 0.20228 0.20022 0.19833 0.19826 0.09939la version fonctionnelle:
testIT = function(M, K = 10000) {
N = length(M[1,])
X = 1
P = matrix(data=0, nrow = 1, ncol = N)
for(i in 1:K) {
P[X] = P[X] + 1
X = sample(1:N, size=1, prob = M[X,])
}
P = P / sum(P)
return(P)
}Par les valeurs propres
Finnalement la méthode la plus exacte mais aussi un peu plus technique, on calcule l’état stationnaire grâce à ses valeurs propres.
On prends la valeur propre la plus grande de la transposé et divise par la somme: On sait que 1 est une valeur propre (car matrice d’une chaine de markov) et que de plus toute les valeurs propre sont comprises entre -1 et 1, on ne prends pas directement 1 car en raison des impressisions cette dernière peut être effectivement écrite comme étant 1.000000000 dans la console mais être en réalité légrement supérieur ou inférieur à 1, d’où le choix du maximum qui fix le problème. De plus la somme des valeurs du vecteur propre associé à 1 vaut 1 ou -1 on divise donc afin d’avoir une valeur positive.
eigen(t(M))$vectors[,which((eigen(t(M))$values) == max(eigen(t(M))$values))]/(sum(eigen(t(M))$vectors[,which((eigen(t(M))$values) == max(eigen(t(M))$values))]))## [1] 0.1 0.2 0.2 0.2 0.2 0.1et la version fonctionnelle
testEig = function(M) {
return(eigen(t(M))$vectors[,which((eigen(t(M))$values) == max(eigen(t(M))$values))]/(sum(eigen(t(M))$vectors[,which((eigen(t(M))$values) == max(eigen(t(M))$values))])))
}Partie 3 Petite analyse des résultat
Puisque la méthode 3 donne la valeurs exacte (aux incertitudes de calcul près…) ous allons l’utiliser pour commenter un peu les états stationaires.
La ligne
Intuition:
Je pensais que la répartitions serait en forme de cloche ou du moins de bosse (fréquence faible aux extrémités et élevé au centre) car bien que les probabiltés soient les mêmes le noeud qui précède a moins de chance que quelqu’un arrive de l’extrémité puisque de même il y a moins de chance que quelqu’un arrive sur l’extrémité, et de proche en proche on aurait le milieu qui serait le moins inpacté.
testEig(Mline(42))## [1] 0.01219512 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024
## [7] 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024
## [13] 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024
## [19] 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024
## [25] 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024
## [31] 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024
## [37] 0.02439024 0.02439024 0.02439024 0.02439024 0.02439024 0.01219512 testEig(Mline(20))## [1] 0.02631579 0.05263158 0.05263158 0.05263158 0.05263158 0.05263158
## [7] 0.05263158 0.05263158 0.05263158 0.05263158 0.05263158 0.05263158
## [13] 0.05263158 0.05263158 0.05263158 0.05263158 0.05263158 0.05263158
## [19] 0.05263158 0.02631579 testEig(Mline(10))## [1] 0.05555556 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111
## [7] 0.11111111 0.11111111 0.11111111 0.05555556 testEig(Mline(5))## [1] 0.125 0.250 0.250 0.250 0.125On peu voir que mon intuition fût complètement fausse les extrémités ont une probabilité exactement la moitié des autres et ces derniers ont une répartition uniforme.
Le cercle
Intuition:
La probabilité est uniforme puisque le contexte est le même en tout noeuds, il n’est donc pas possible qu’un noeud aie une probabilité plus forte ou faible.
testEig(MCircle(42))## [1] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [7] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [13] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [19] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [25] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [31] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [37] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 testEig(MCircle(20))## [1] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
## [15] 0.05 0.05 0.05 0.05 0.05 0.05 testEig(MCircle(10))## [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 testEig(MCircle(5))## [1] 0.2 0.2 0.2 0.2 0.2Mon intuition a été bonne, il n’y a pas grand chose à commenter ici.
La sucette
Intuition:
Je pense que la probanilité est plus faible près de l’extrémité et plus forte près de l’intersection.
testEig(MSucette(42))## [1] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [7] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [13] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [19] 0.02380952 0.02380952 0.03571429 0.02380952 0.02380952 0.02380952
## [25] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [31] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952
## [37] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.01190476 testEig(MSucette(20))## [1] 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.075 0.050
## [12] 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.025 testEig(MSucette(10))## [1] 0.10 0.10 0.10 0.10 0.15 0.10 0.10 0.10 0.10 0.05 testEig(MSucette(6))## [1] 0.16666667 0.16666667 0.25000000 0.16666667 0.16666667 0.08333333Mon intuition était faussé de la même manière que pour la ligne et on a donc une répartition uniforme partout sauf à l’extrémité qui à une probabilité qui vaut la moitié des autres et l’intersection avec le baton qui a une probabilité 50% plus forte.
Partie 4: Convergence
Nous allons maintenant nous consentrer sur la vitesse de convergence et l’efficacité de ces algorythmes pour ces différentes formes. Pour celà nous considéreront la distance euclidienne comme mesure de comparaison et la méthode 3 comme valeur cible.
Par itération:
Ligne
library(ggplot2)
M = Mline(30)
cible = testEig(M)
a = proc.time()[3]
res = lapply(rep(seq(100, 10000, by=100), each=10), function(K) testIT(M,K))
b = proc.time()[3]
qplot(rep(seq(10, 10000, by=100), each=10),as.numeric(lapply(res, function(PP) dist(rbind(cible, PP)))), geom='smooth', span =0.5)## `geom_smooth()` using method = 'gam'
b-a## elapsed
## 31.527Cercle
library(ggplot2)
M = MCircle(30)
cible = testEig(M)
a = proc.time()[3]
res = lapply(rep(seq(100, 10000, by=100), each=10), function(K) testIT(M,K))
b = proc.time()[3]
qplot(rep(seq(10, 10000, by=100), each=10),as.numeric(lapply(res, function(PP) dist(rbind(cible, PP)))), geom='smooth', span =0.5)## `geom_smooth()` using method = 'gam'
b-a## elapsed
## 31.344Sucette
library(ggplot2)
M = MSucette(30)
cible = testEig(M)
a = proc.time()[3]
res = lapply(rep(seq(100, 10000, by=100), each=10), function(K) testIT(M,K))
b = proc.time()[3]
qplot(rep(seq(10, 10000, by=100), each=10),as.numeric(lapply(res, function(PP) dist(rbind(cible, PP)))), geom='smooth', span =0.5)## `geom_smooth()` using method = 'gam'
b-a## elapsed
## 31.781On peut voir que la distance moyenne à la cible pour la sucette diminue bien plus doucement que pour les deux autres. Par ailleurs la ligne converge un peu moins rapidement que le cercle. De plus la distance converge vers 0.05, calculer au delà de 10000 ne semble pas très efficace.
Par multiplication:
Ligne
library(ggplot2)
M = Mline(30)
cible = testEig(M)
a = proc.time()[3]
res = lapply(rep(seq(1, 100, by=1), each=10), function(K) testMul(M,K))
b = proc.time()[3]
qplot(rep(seq(1, 100, by=1), each=10),as.numeric(lapply(res, function(PP) dist(rbind(cible, PP)))), geom='smooth', span =0.5)## `geom_smooth()` using method = 'gam'
b-a## elapsed
## 0.35Cercle
library(ggplot2)
M = MCircle(30)
cible = testEig(M)
a = proc.time()[3]
res = lapply(rep(seq(1, 100, by=1), each=10), function(K) testMul(M,K))
b = proc.time()[3]
qplot(rep(seq(1, 100, by=1), each=10),as.numeric(lapply(res, function(PP) dist(rbind(cible, PP)))), geom='smooth', span =0.5)## `geom_smooth()` using method = 'gam'
b-a## elapsed
## 0.229Sucette
library(ggplot2)
M = MSucette(30)
cible = testEig(M)
a = proc.time()[3]
res = lapply(rep(seq(1, 100, by=1), each=10), function(K) testMul(M,K))
b = proc.time()[3]
qplot(rep(seq(1, 100, by=1), each=10),as.numeric(lapply(res, function(PP) dist(rbind(cible, PP)))), geom='smooth', span =0.5)## `geom_smooth()` using method = 'gam'
b-a## elapsed
## 0.234Pour le cercle les valeurs sont celles initial et ainsi la valeur exacte est trouvé instantanément, celà montre qu’avoir une intuition des probabilités de l’état stationnaire peu considérablement accélerer la convergence. De plus la valeur final pour la ligne est deux fois plus petite pour la ligne que pour la sucette qui converge un peu moins vite. Néanmoins les valeurs continu de se raprocher rapidement de 0, ce qui veux dire qu’il peut être interressant des nombre d’itération encore plus grands pour les tests. Enfin en seulement 100 itération cette méthode donne une distance 10 fois plus petite que la méthode précédente en 10000 itération. De plus il y a une différence colossale de temps de calcul nécessaire pour atteindre une même distance, on voit que la méthode par multiplication est beauoup plus efficace que celle par itéation. Enfin si les probabilité variait énormément il faudrait effectuer beaucoup plus de calcul avec les itération pour s’assurer d’être passé partout.
Le temps
Afin de comparer les temps on mesure aussi avec la troisième méthode ####Ligne
M = Mline(30)
a = proc.time()[3]
res = testMul(M,100)
b = proc.time()[3]
print("mul")## [1] "mul"b-a## elapsed
## 0.002a = proc.time()[3]
cible = testEig(M)
b = proc.time()[3]
print("eig")## [1] "eig"b-a## elapsed
## 0.005print("dist")## [1] "dist"dist(rbind(cible, res))## cible
## 0.004815672Cercle
library(ggplot2)
M = MCircle(30)
a = proc.time()[3]
res = testMul(M,100)
b = proc.time()[3]
print("mul")## [1] "mul"b-a## elapsed
## 0.002a = proc.time()[3]
cible = testEig(M)
b = proc.time()[3]
print("eig")## [1] "eig"b-a## elapsed
## 0.004print("dist")## [1] "dist"dist(rbind(cible, res))## cible
## 4.361579e-16Sucette
M = MSucette(30)
a = proc.time()[3]
res = testMul(M,100)
b = proc.time()[3]
print("mul")## [1] "mul"b-a## elapsed
## 0.002a = proc.time()[3]
cible = testEig(M)
b = proc.time()[3]
print("eig")## [1] "eig"b-a## elapsed
## 0.005print("dist")## [1] "dist"dist(rbind(cible, res))## cible
## 0.00580905La méthode par itération est considérablement plus lente que les autres pour avoir des résultats bien inférieur, à tel point qu’il n’est pas utile de mesure son temps d’execution d’un test pour le comparer aux autres. Par contre le temps mis par la méthode des multiplication de matrice est nettement inférieur à celui mis par la méthode des valeurs propres (facteur 5 ici) pour une distance extrèmement faible.
Conclusion :
Ceci montre que les deux méthodes potentiellement utile sont: Par multiplicationde, très rapide et converge bien mais néanmoins sera toujours qu’une approximation. Par les valeurs propre, notablement plus lente, de plus il n’est éventuellement pas possible de réaliser les calculs nécessaires si les matrices deviennent trop grosse.
Enfin sur un système ou la matrice de transition n’est pas réalisable le calcule de la méthode deux peux être recréer en faisant se déplacer de nombreuses instances réparties sur les noeuds