set.seed(42)

Q0: Décrire votre intuition

Mon intuition initiale était que les proportions de chaque type d’individus évoluait très lentement et que tout type de dérive était possible.

Je me disais néanmoins qu’il était peut-être possible, vu le caractère synchrone de l’évolution, d’avoir des évolutions non continues avec des genres de cycles (une certaine proportion pour les itérations paires et une autre proportion pour les itérations impaires).

Enfin, une extinction d’un des types d’individu me semblait possible mais mon intuition était que c’était très peu probable, surtout si la population est grande.

Codons gaiement

Je suis reparti de ce que nous avions écrit en TD la dernière fois.

new_population = function(P = 10000, bb = 0.1, mm = 0.3) {
  Pere = sample(size=P, x=c(0,1,2), replace=T, prob=c(mm,1-(bb+mm),bb));
  Mere = sample(size=P, x=c(0,1,2), replace=T, prob=c(mm,1-(bb+mm),bb));
  
  Enfant_P=ifelse(Pere==0,0,
                  ifelse(Pere==1,sample(size = P,x=c(0,1),replace = T),1));
  Enfant_M=ifelse(Mere==0,0,
                  ifelse(Mere==1,sample(size = P,x=c(0,1),replace = T),1));
  Enfant = Enfant_M + Enfant_P;
  list(BB=sum(Enfant==2),MM=sum(Enfant==0))
}
new_population(); 
## $BB
## [1] 1605
## 
## $MM
## [1] 3555

Je n’avais alors plus qu’à rajouter une fonction qui génère une trajectoire (et mette le tout dans une “dataframe” car c’est plus propre comme ça):

population_evolution = function(P=20, Imax=20, bb=0.1, mm=0.3) {
  df=data.frame(Timestep=0,BB=bb,MM=mm); 
  for(i in 1:Imax) {
    res=new_population(P,bb,mm);
    bb=res$BB/P;mm=res$MM/P;
    df=rbind(df,data.frame(Timestep=i,BB=bb,MM=mm));
  }
  df
}
population_evolution(P=20);
##    Timestep   BB   MM
## 1         0 0.10 0.30
## 2         1 0.15 0.40
## 3         2 0.25 0.20
## 4         3 0.25 0.30
## 5         4 0.35 0.20
## 6         5 0.30 0.15
## 7         6 0.25 0.15
## 8         7 0.35 0.20
## 9         8 0.35 0.05
## 10        9 0.60 0.15
## 11       10 0.45 0.05
## 12       11 0.40 0.05
## 13       12 0.40 0.20
## 14       13 0.30 0.10
## 15       14 0.45 0.05
## 16       15 0.40 0.10
## 17       16 0.30 0.05
## 18       17 0.40 0.10
## 19       18 0.50 0.15
## 20       19 0.35 0.10
## 21       20 0.35 0.15

Puis la même chose pour générer plusieurs trajectoires:

population_evolutions = function(P=20, Imax=20, N=10, bb=0.1, mm=0.3) {
  df=data.frame()
  for(i in 1:N) {
    traj=population_evolution(P,Imax,bb,mm);
    traj$Idx=i;
    df=rbind(df,traj);
  }
  df
}
population_evolutions(P=20,Imax=4,N=4);
##    Timestep   BB   MM Idx
## 1         0 0.10 0.30   1
## 2         1 0.15 0.30   1
## 3         2 0.15 0.35   1
## 4         3 0.10 0.60   1
## 5         4 0.05 0.60   1
## 6         0 0.10 0.30   2
## 7         1 0.15 0.25   2
## 8         2 0.20 0.25   2
## 9         3 0.15 0.45   2
## 10        4 0.15 0.55   2
## 11        0 0.10 0.30   3
## 12        1 0.05 0.40   3
## 13        2 0.05 0.45   3
## 14        3 0.15 0.60   3
## 15        4 0.10 0.50   3
## 16        0 0.10 0.30   4
## 17        1 0.05 0.45   4
## 18        2 0.05 0.50   4
## 19        3 0.05 0.60   4
## 20        4 0.00 0.60   4

Q1 Cas d’une petite population

On s’intéresse dans cette question au cas d’une population de taille réduite sur un horizon réduit.

set.seed(2)

Premiers pas

Commençons par générer puis visualiser une seule trajectoire:

pop = population_evolution(P = 20,Imax = 20, bb = 4/20, mm = 12/20);

Alors bien, sûr, on peut utiliser le R de base.

xrange = c(min(pop$Timestep),max(pop$Timestep))
plot(x = NULL, y = NULL, xlim = xrange,ylim = c(0,1), type="o", xlab="Génération i", ylab="Proportion d'individus de chaque type");
lines(pop$Timestep, pop$BB, type = "o",col = "blue");
lines(pop$Timestep, pop$MM, type = "o",col = "brown");

Certains m’ont même fait des versions superbes pour bien illustrer les trois types:

I=max(pop$Timestep);
plot(c(0,I),c(0,1),type="n",xlab="Generation", ylab="Proportion d'individus de chaque type")
polygon(c(c(0,I),c(I,0)),c(c(0,0),c(1,1)), col="burlywood")
polygon(c(0:I,rev(0:I)),c(pop$BB,rep(0,I+1)), col="skyblue")
polygon(c(0:I,rev(0:I)),c(1-pop$MM,rep(1,I+1)), col="tan4")

Ceci dit, dès qu’on essaye de dessiner plusieurs trajectoires, on va être un peu embêtés et plus y comprendre grand chose. En particulier, on ne saura plus quelle courbe bleue correspond à quelle courbe marron.

Finalement, l’état de la population, a deux coordonnées, pourquoi ne pas représenter les choses dans cet espace là, comme le suggérait l’ennoncé:

plot(pop$BB,pop$MM, xlim=c(0,1), ylim=c(0,1))
lines(pop$BB,pop$MM)

Bon, c’est un peu moche aussi et ça ne permet pas de savoir dans quel sens on va mais par contre, je devrais être capable de dessiner plusieurs trajectoires. Essayons avec ggplot2.

library(ggplot2)
ggplot(pop,aes(x=BB,y=MM,color=Timestep)) + geom_path() + 
  xlim(0,1) +  ylim(0,1) + theme_bw() + coord_fixed()