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.04166667

Version 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.09939

la 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.1

et 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.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.

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.2

Mon 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.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.

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.527

Cercle

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

Sucette

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.

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'

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

Line

Circle

Sucette

Partie 2: L’état stationnaire

Par itération

Version fonctionnelle

Par vecteur de probabilité

la version fonctionnelle:

Par les valeurs propres

et la version fonctionnelle

Partie 3 Petite analyse des résultat

La ligne

Intuition:

Le cercle

Intuition:

La sucette

Intuition:

Partie 4: Convergence

Par itération:

Ligne

Cercle

Sucette

Par multiplication:

Ligne