TD1 : Analyse de données avec R (RICM4, 2013)

1. Régulation des naissances

Premiers pas

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
            }
        }
    }
    c(girls, boys)
}
new_population(1000)
## [1] 738 759

Au fait, une façon d'écrire ça, “à la R”:

new_population2 <- function(n = 100) {
    x <- floor(2 * runif(n))
    y <- floor(2 * runif(n))
    c(sum(1 - x) + sum((1 - x) * (1 - y)), sum(x) + sum((1 - x) * y))
}
new_population2(1000)
## [1] 693 775

Étonamment, c'est au moins 10 fois plus rapide mais également assez illisible donc tant que posssible… évitez ce style. ;)

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 ?…

p <- new_population2(2000)
p[2]/(p[1] + p[2])
## [1] 0.4856

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…

ratio_exp_dist <- function(nb_couples = 100, n = 200) {
    girls = c()
    boys = c()
    for (i in 1:200) {
        p = new_population2(nb_couples)
        girls <- c(girls, p[1])
        boys <- c(boys, p[2])
    }
    ratio = boys/(girls + boys)
}
hist(ratio_exp_dist(nb_couples = 2000))

plot of chunk unnamed-chunk-4

À 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:

pop_generate_dist <- function(nb_couples = 100, n = 200) {
    girls = c()
    boys = c()
    for (i in 1:200) {
        p = new_population2(nb_couples)
        girls <- c(girls, p[1])
        boys <- c(boys, p[2])
    }
    list(girls = girls, boys = boys)
}
pop <- pop_generate_dist(nb_couples = 2000)
plot(pop$girls, pop$boys)

plot of chunk unnamed-chunk-5

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 ?

test_pop = 2000
p <- new_population2(test_pop)
(p[1] + p[2])/(2 * test_pop)
## [1] 0.7518

Bon, et bien à vous de montrer que le taux de renouvellement est bien de ¾…

Une infinité de filles!!!

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
            }
        }
    }
    c(girls, boys)
}
new_population_repeat(1000)
## [1] 1004 1000

Bon et donc mêmes questions:

  1. vers quoi converge le pourcentage de garçons ?

    test_pop = 2000
    p <- new_population_repeat(test_pop)
    (p[2])/(p[1] + p[2])
    
    ## [1] 0.5044
    
  2. Et cette fois, la taille de la population ? Ça augmente, ça diminue, ça reste stable ?

    test_pop = 2000
    p <- new_population_repeat(test_pop)
    (p[1] + p[2])/(2 * test_pop)
    
    ## [1] 0.9992
    

2. Analyse statistique

Proportion de voyelles

Dans le texte du premier paragraphe de Candide (en Francais), nous avons observé une population des voyelles.

Il y avait en Westphalie, dans le chateau de M. le baron de Thunder-ten-tronckh, un jeune garcon a qui la nature avait donne les moeurs les plus douces. Sa physionomie annoncait son ame. Il avait le jugement assez droit, avec l'esprit le plus simple; c'est, je crois, pour cette raison qu'on le nommait Candide. Les anciens domestiques de la maison soupconnaient qu'il etait fils de la soeur de monsieur le baron et d'un bon et honnete gentilhomme du voisinage, que cette demoiselle ne voulut jamais epouser parce qu'il n'avait pu prouver que soixante et onze quartiers, et que le reste de son arbre genealogique avait ete perdu par l'injure du temps.

Pareil pour le texte en anglais. On souhaite comparer les fréquences d'apparition de ces voyelles entre les deux langues ainsi qu'à celles des langues anglaises et françaises.

df <- data.frame(voyelle = c("a", "e", "i", "o", "u", "y"), nombre = c(44, 81, 
    37, 35, 35, 2))
# Hint: utilisez plutôt read.csv pour les imports de données...
df$freq = df$nombre/sum(df$nombre)
df_en <- data.frame(voyelle = c("a", "e", "i", "o", "u", "y"), nombre = c(40, 
    69, 25, 34, 12, 8))
df_freq <- data.frame(voyelle = c("a", "e", "i", "o", "u", "y"), freq_th = c(0.1795, 
    0.3739, 0.1677, 0.1227, 0.1489, 0.00073))

df_en <- data.frame(voyelle = c("a", "e", "i", "o", "u", "y"), nombre = c(40, 
    69, 25, 34, 12, 8))
df_en$freq = df_en$nombre/sum(df_en$nombre)
df_en_freq <- data.frame(voyelle = c("a", "e", "i", "o", "u", "y"), freq_th = c(0.196, 
    0.315, 0.176, 0.193, 0.07, 0.05))

Bon, maintenant que les données sont dans des data-frames, on pourrait utiliser les fonctions classiques de R pour faire les dessins, comme ça:

par(mfrow = c(1, 2))
plot(df$voyelle, df$freq, main = "Français")
plot(df_en$voyelle, df_en$freq, main = "Anglais")

plot of chunk unnamed-chunk-11

Bon, mais c'est assez moche finalement et configurer les légendes et tout devient vite cauchemardesque. Je vous invite donc plutôt à utiliser ggplot2 qui demande un autre type de gymnastique mais au moins vous ouvrira des horizons insoupçonnés… ;)

library(ggplot2)
df <- merge(df, df_freq)
df_en <- merge(df_en, df_en_freq)
df$language = "Français"
df_en$language = "Anglais"
df <- rbind(df, df_en)
ggplot(data = df, aes(x = voyelle, y = freq, fill = voyelle)) + geom_bar(stat = "identity") + 
    facet_wrap(~language) + theme_bw() + geom_point(aes(y = freq_th))

plot of chunk unnamed-chunk-12

Ou encore (ne me demandez pas pourquoi le factor(1), je me suis juste inspiré de la documentation!):

ggplot(data = df, aes(x = factor(1), y = freq, fill = voyelle)) + geom_bar(stat = "identity", 
    width = 1) + facet_wrap(~language) + theme_bw() + coord_polar(theta = "y")

plot of chunk unnamed-chunk-13

Mmmh, bof. Vraiment pas clair que les camemberts soient la meilleure repréentation pour ce genre de comparaisons…

S'il vous arrivait que les données ne soient pas vraiment structurées comme il faut. Jetez à tout hasard un coup d'oeil au paquet reshape.

Longueur des mots

On etudie la variable statistique \( X \) qui représente la longueur des mots dans le même texte en Francais. On commence dont par séparer chacun des mots (merci sed):

cat dat/chap1-candide.txt | sed -e "s/['\.;,]/ /g" -e 's/  */ /g' | sed -e 's/ /\n/g' > dat/toto.txt

avant de faire des statistiques dessus

lst <- read.table("./dat/toto.txt")
X <- nchar(as.character(lst$V1))

Pour calculer les fréquences cumulées, vous pourriez bien sûr utiliser des fonctions comme mean et cumsum mais il y a bien sûr déjà ce qu'il faut…

par(mfrow = c(1, 2))
hist(X)
plot(ecdf(X))

plot of chunk unnamed-chunk-15

par(mfrow = c(1, 1))
summary(X)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    2.00    4.00    4.37    6.00   19.00

Vu qu'on s'intéresse à des variables discrètes ici et qu'il y a peu de valeurs différentes, il vaut mieux s'assurer que les bins des histogrammes tombent bien comme il faut…

hist(X, breaks = max(X))

plot of chunk unnamed-chunk-16

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 évitez les boucles autant que possible. :)

Convincing...