Afin de faciliter la reproduction de ces courbes, commençons par fixer la graine et afficher les premières valeurs pseudo-aléatoires.
set.seed(1:6)
library(ggplot2)
runif(10)
## [1] 0.26551 0.37212 0.57285 0.90821 0.20168 0.89839 0.94468 0.66080
## [9] 0.62911 0.06179
version
## _
## platform x86_64-pc-linux-gnu
## arch x86_64
## os linux-gnu
## system x86_64, linux-gnu
## status
## major 3
## minor 0.2
## year 2013
## month 09
## day 25
## svn rev 63987
## language R
## version.string R version 3.0.2 (2013-09-25)
## nickname Frisbee Sailing
Intuitivement, c'est symétrique entre les deux villes et donc, ça doit bien se répartir entre les deux villes. Chacune d'elle doit avoir à peu près la moitié de la population totale. Ceci étant dit, il n'est pas évident pour moi que ça converge. À tout instant, on pourrait ne pas avoir de “chance” et, même si c'est “improbable”, tirer un très grand nombre de fois de suite la même ville, ce qui déséquilibrerait la situation.
Après, je pense que si on part de situations où il y a une population totale plus importante, probablement qu'on converge plus rapidement vers l'équilibre.
Enfin, probablement que si on part d'une situation déséquilibrée, on doit conserver ce déséquilibre…
Écrivons la fonction calculant la dynamique de la population de nos deux villes.
city_evolution <- function(n = 500, p1 = 1, p2 = 1) {
n1 <- c(p1)
n2 <- c(p2)
for (t in 1:(n - 1)) {
if (runif(1) < (n1[t]/(n1[t] + n2[t]))) {
n1[t + 1] = n1[t] + 1
n2[t + 1] = n2[t]
} else {
n2[t + 1] = n2[t] + 1
n1[t + 1] = n1[t]
}
}
data.frame(t = (1:n), n1 = n1, n2 = n2)
}
Puis affichons une dizaine de trajectoires.
df = data.frame()
for (trial in 1:10) {
d = city_evolution()
d$trial = trial
df = rbind(df, d)
}
ggplot(data = df, aes(x = t, y = n1/(n1 + n2), color = factor(trial))) + geom_line()
On notera le factor(trial) qui permet d'avoir une échelle de couleurs discrète puisque les trajectoires n'ont rien à voir les unes avec les autres.
limit_sample <- function(trials = 400, p1 = 1, p2 = 1, n = 500) {
limit = c()
for (trial in 1:trials) {
d = city_evolution(n, p1, p2)
limit[trial] = d[n, ]$n1/(d[n, ]$n1 + d[n, ]$n2)
}
limit
}
d <- data.frame(lim = limit_sample(trials = 1000, p1 = 1, p2 = 1), p1 = 1, p2 = 1)
ggplot(data = d, aes(x = lim)) + stat_bin(aes(y = ..count../sum(..count..)),
colour = "darkgreen", fill = "white") + ylab("Frequency") + xlim(0, 1)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Ce n'est pas vraiment plat mais aucun endroit n'a vraiment l'air d'être franchement favorisé par rapport à un autre. La fonction de répartition qui y ressemble le plus est donc plutôt la loi uniforme. Pour tester cette hypothèse, j'ai augmenté le nombre d'échantillons:
d <- data.frame(lim = limit_sample(trials = 50000, p1 = 1, p2 = 1), p1 = 1,
p2 = 1)
ggplot(data = d, aes(x = lim)) + stat_bin(aes(y = ..count../sum(..count..)),
colour = "darkgreen", fill = "white") + ylab("Frequency") + xlim(0, 1)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Aaah, cette fois-ci, c'est bien plus flagrant.
Commençons par étudier l'influence de la taille de la population initiale
df2 = data.frame()
for (p2 in c(1, 10, 100)) {
d = data.frame(lim = limit_sample(trials = 1000, p1 = p2, p2 = p2), p1 = p2,
p2 = p2)
df2 = rbind(df2, d)
}
ggplot(data = df2, aes(x = lim)) + stat_bin(aes(y = ..count../sum(..count..)),
colour = "darkgreen", fill = "white") + ylab("Frequency") + xlim(0, 1) +
facet_wrap(~p2, scales = "free_y")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Comme on pouvait s'y attendre, plus la population initiale est importante, plus il devient difficile de s'écarter de la proportion initiale. La forme est intéressante. On croit reconnaître une jolie courbe en cloche mais à la différence d'une Gaussienne, le support est fini. Peut-être que ça tend vers une gaussienne ? Ou bien vers un dirac plus simplement ? Et quand on parle de convergence, de toutes façons, ça ne doit pas être trivial puisqu'on a plusieurs limites: la taille de la population, la taille de la population initiale, …
Bref, regardons maintenant l'influence du déséquilibre de la population initiale
df = data.frame()
for (p2 in c(1, 10, 100)) {
d = data.frame(lim = limit_sample(trials = 1000, p1 = 1, p2 = p2), p1 = 1,
p2 = p2)
df = rbind(df, d)
}
ggplot(data = df, aes(x = lim)) + stat_bin(aes(y = ..count../sum(..count..)),
colour = "darkgreen", fill = "white") + ylab("Frequency") + xlim(0, 1) +
facet_wrap(~p2, scales = "free_y")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Notons qu'on combine deux effets là, le déséquilibre et la taille de la population initiale mais on voit bien comment l'état initial influence la suite.
Et bien, je n'avais pas imaginé que l'on puisse converger vers une loi uniforme… Le fait que ça converge bien à chaque fois est intuitif mais pas évident. J'imagine qu'on doit parler de convergence presque sûrement (ou quelque chose comme ça), encore faut-il le définir bien comme il faut. Après, le reste l'influence de la population de départ était assez logique encore que la forme de la distribution n'est pas triviale.
Quelles autres hypothèses pourrait-on tester ?
Sinon, c'est quand même bien lent ces calculs de trajectoires. Je me demande comment accélérer tout ça…