--- title: "RICM4/PS 2014: Régulation des naissances" author: "Arnaud Legrand" date: "07/10/2014" output: html_document: toc: true theme: cosmo highlight: tango --- # Versions simples... ## La compacte naturelle Écrivons une première version naïve, toute simple mais qui renvoie les résultats proprement avec une dataframe. ```{r} new_population <- function(nb_couples = 100) { girls = 0 boys = 0 for(i in 1:nb_couples) { if(runif(1) > 0.5) { boys <- boys + 1 } else { girls <- girls + 1 if(runif(1) > 0.5) { boys <- boys + 1 } else { girls <- girls + 1 } } } data.frame(girls,boys) } new_population(1000) ``` ## La détaillée On peut aussi faire une version où on a le détail de la population. ```{r} new_population_detailed <- function(nb_couples = 100) { df=data.frame(); for(i in 1:nb_couples) { girls = 0 boys = 0 if(runif(1) > 0.5) { boys <- boys + 1 } else { girls <- girls + 1 if(runif(1) > 0.5) { boys <- boys + 1 } else { girls <- girls + 1 } } df=rbind(df,data.frame(girls,boys)) } df } pop=new_population_detailed(1000) ``` Dans ce cas, on a accès au détail: ```{r} head(pop); ``` Et si on veut le nombre total d'enfants de chaque sexe: ```{r} sum(pop$boys) sum(pop$girls) ``` ## La illisible mais efficace Au fait, une façon d'écrire ça, "à la R" pourrait être ```{r} new_population2 <- function(n = 100) { x <- floor(2*runif(n)) y <- floor(2*runif(n)) data.frame(girls=sum(1-x) + sum((1-x)*(1-y)), boys = sum(x) + sum((1-x)*y)) } new_population2(1000) ``` Étonamment, c'est au moins 10 fois plus rapide mais également assez illisible donc à manier avec précaution. # Tout ça pour quoi ? Bref, dans tous les cas, vous vous rendez vite compte que quand n=100, il peut y avoir beaucoup plus de filles que de garçons ou l'inverse. Pour n=1000, l'écart se resserre. Et donc, le taux de masculinité tend vers ?... ```{r} p <- new_population2(2000) p$boys/(p$boys+p$girls) ``` Surprenant, non ? Mmmh, mais ça change à chaque fois donc comment être sûr. Visualisons celà en regardant par exemple la distribution expérimentale de la proportion de garçons... ```{r fig.width=7, fig.height=4} pop_generate_dist <- function (nb_couples=100,n=200) { df= data.frame() for(i in 1:n) { df = rbind(df,new_population2(nb_couples)) } df$ratio = df$boys/(df$girls+df$boys) df } hist(pop_generate_dist(nb_couples=2000)$ratio) ``` À partir de là, on peut voir plein de choses bizarres quand il y a peu de couples (nb_couples faible) ou qu'on ne fait que peu d'échantillons (n faible). Allez, encore une variation sur le thème: ```{r} pop <- pop_generate_dist(nb_couples=2000) plot(pop$girls,pop$boys) ``` Mmmh, mais que peut bien vouloir dire une telle représentation ? Qu'il y a une correlation entre le nombre de filles et le nombre de garçons ? Est-ce vraiment ce qui nous intéresse ? Bref, le choix de la représentation graphique est vraiment primordial dans l'analyse que l'on fait et dans le message que l'on souhaite passer. Et au fait, la taille de la population ? Ça augmente, ça diminue, ça reste stable ? ```{r} population_init=2000 p <- new_population2(population_init) (p$boys+p$girls)/(2*population_init) ``` Bon, et bien à vous de montrer que le taux de renouvellement est bien de 3/4... # Variations sur le thème ## Une infinité de filles!!! ```{r} new_population_repeat <- function(nb_couples = 100) { girls = 0 boys = 0 for(i in 1:nb_couples) { repeat { x <- runif(1) if(x>0.5) { boys <- boys + 1 break } else { girls <- girls + 1 } } } data.frame(girls,boys) } new_population_repeat(1000) ``` Bon et donc mêmes questions: 1. vers quoi converge le pourcentage de garçons ? ```{r} test_pop = 2000 p <- new_population_repeat(test_pop) (p$boys)/(p$boys+p$girls) ``` 2. Et cette fois, la taille de la population ? Ça augmente, ça diminue, ça reste stable ? ```{r} (p$boys+p$girls)/(2*test_pop) ``` ## Et voilà ce qui arrive quand on laisse trop le choix... ```{r} new_population_repeat <- function(nb_couples = 100, proba_no_child=.01, proba_girl=c(.5,.45,.4,.35,.3,.25,.2,.15,.1)) { df=data.frame(); for(i in 1:nb_couples) { girls = 0 boys = 0 x <- runif(1) if(x>proba_no_child) { for(p in proba_girl) { x <- runif(1) if(x<1-p) { boys <- boys + 1; break } else { girls <- girls + 1 } } } df=rbind(df,data.frame(girls,boys)) } df } pop = new_population_repeat(nb_couples = 1000) ``` # Moralité * R dispose déjà de toutes les fonctionnalités dont vous avez besoin et les choix par défaut sont souvent les bons (pour autant qu'on sache un minimum dans quels cas ils s'appliquent). Donc, ne vous amusez jamais à recoder des fonctionnalités statistiques ou graphique et encore une fois évitez les boucles autant que possible. :) * Quand vous faites une simulation, utilisez des data.frame pour faire des choses propres. * Évitez de faire les statistiques au fur et à mesure de la simulation. Il est généralement plus souple de générer plein de scénarios et de faire les statistiques après coup dessus. ![Convincing...](http://imgs.xkcd.com/comics/convincing.png)