--- title: "Corrigé de la séance de TD RICM4 PS 2" author: "Arnaud Legrand" date: "16 septembre 2016" output: html_document --- # Simulation d'une famille de 5 enfants ## Première version (à la java/C/...) ```{r} N = 100000; Famille = rep(x = 0,N); for(i in 1:N) { for (j in 1:5) { Famille[i] = Famille[i] + sample(x = c(0,1),size = 1,replace = T); } } nb = 0 for(i in 1:N) { if(Famille[i]==2) {nb = nb +1;} } nb/N ``` Ça marche bien mais on peut remarquer que c'est un peu lent (et long à écrire... ;). ## Deuxième version qui exploite les vecteurs de R ```{r} N = 100000; Famille = rep(x = 0,N); for(i in 1:N) { Famille[i] = sum(sample(x = c(0,1),size = 5,replace = T)); } nb = length(Famille[Famille==2]) nb/N ``` On a viré deux boucles de façon assez différente... C'est un poil plus court mais on ne voit pas de changement de vitesse notable. ## Troisième version bien plus proche de la formulation mathématique Si on devait modéliser ce problème, on pourrait noter $X_i$ la variable aléatoire valant 1 si le $i$ème enfant est un garçon et 0 sinon. Ainsi, le nombre de garçons de la famille s'exprime $N_G=X_1+X_2+X_3+X_4+X_5$. Ce codage avec un 0 ou un 1 n'est donc pas complètement arbitraire. Il a une sémantique que j'utilise pour recombiner mes variables en d'autres variables dont la sémantique est non ambigüe. La probabilité d'avoir un garçon s'écrit alors $P(N_G==2)$. ```{r} N = 100000; X1= sample(x = c(0,1),size = N,replace = T); X2= sample(x = c(0,1),size = N,replace = T); X3= sample(x = c(0,1),size = N,replace = T); X4= sample(x = c(0,1),size = N,replace = T); X5= sample(x = c(0,1),size = N,replace = T); NG = X1 + X2 + X3 + X4 +X5; nb = length(NG[NG==2]) nb = sum(NG==2) nb/N ``` Non seulement c'est court mais en plus, ça va vite!!! :) Alors oui, bien sûr, ce copier coller pour X1 à X5, c'est un peu moche et on pourrait remplacer par une boucle mais ça, par contre, ça ne change pas fondamentalement les choses... # Les yeux noirs... Sachant que les yeux bleus sont un caractère récessif, on s'intéresse au "jeu" suivant. On se donne une population $P^0$ de taille $N$. Pour générer la population $P^1$: on prend deux personnes au hasard dans $P^0$; ces deux personnes engendrent un enfant qui va dans $P^1$; on recommence jusqu'à ce que $P^1$ soit composé de N personnes. Au tout départ, la populatio $P^0$ comporte une certaine fraction d'individus ayant les allèles $ \{B,B\}$, $\{B,M\}$ et $\{M,M\}$. La question que l'on se pose est comment ces proportions évoluent-elles entre $P^0$ et $P^1$. On pourra éventuellement générer $P^2$, $P^3$ de la même façon. A-t-on une extinction des yeux bleus ? Des cycles ? ... Pour modéliser ce problème, on considère comme variable aléatoire le nombre d'alèles bleus pour le père et pour la mère. Ces variables aléatoires se recombinent alors aisément pour définir le nombre d'alèles bleus de l'enfant. ```{r} # J'ai enrobé tout ça avec une fonction après coup pour que ça soit plus pratique... new_population = function(N = 10000, a = 0.1, b = 0.3) { Pere = sample(size=N, x=c(0,1,2), replace=T, prob=c(a,1-(a+b),b)); Mere = sample(size=N, x=c(0,1,2), replace=T, prob=c(a,1-(a+b),b)); Enfant_P=ifelse(Pere==0,0, ifelse(Pere==1,sample(size = N,x=c(0,1),replace = T),1)); Enfant_M=ifelse(Mere==0,0, ifelse(Mere==1,sample(size = N,x=c(0,1),replace = T),1)); Enfant = Enfant_M + Enfant_P; Enfant } a=0.1; b=0.3; Enfant = new_population(a=a,b=b); ``` On peut alors regarder la nouvelle distribution (quand $N$ est grand, c'est relativement stable) ```{r,fig.width=5} hist(Enfant,probability = T) ``` Bon, un histogramme en fait, c'est pas vraiment fait pour ça. Donc, voici mieux (google et stackoverflow sont vos amis) et je compare en plus avec la distribution initiale: ```{r,fig.width=5} counts <- table(Enfant) bar = barplot(counts / sum(counts), main="", xlab="Nombre d'alèles bleus") points(x = bar, y=c(a,1-(a+b),b)) ``` Bon, ben, ça a l'air d'avoir "significativement" changé: ```{r} counts/sum(counts) ``` ```{r,fig.width=5} bb = counts[["0"]]/sum(counts) mm = counts[["2"]]/sum(counts) Enfant2 = new_population(a=bb,b=mm) counts2 <- table(Enfant2) bar = barplot(counts2 / sum(counts2), main="", xlab="Nombre d'alèles bleus") points(x = bar, y=c(bb,1-(bb+mm),mm)) ``` Ah ben, cette fois-ci, ça n'a pas l'air d'avoir vraiment changé. Faudrait tenter pour d'autres valeurs...