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)
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)
}
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)
}
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)
}
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.04166667
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)
}
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.09939
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)
}
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.1
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))])))
}
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.
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.125
On 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.
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.2
Mon intuition a été bonne, il n’y a pas grand chose à commenter ici.
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.08333333
Mon 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.
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.
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.527
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.344
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.781
On 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.
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'