Calculs statistiques en R

Calculs statistiques en R Yousri Henchiri 14-11-2021 Table des matières 1 Statistiques descriptives 1 1.1 Fonctions retournant une seule statistique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Traitement des données manquantes et argument na.rm . . . . . . . . . . . . . . . . . . . . . 3 1.3 Fonctions pouvant retourner plusieurs statistiques . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Fonctions retournant un vecteur de statistiques . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 Calcul de fréquences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Distributions de probabilité 9 2.1 Fonction de densité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Fonction de répartition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Fonction quantile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3 Génération de nombres pseudo-aléatoires 14 3.1 Fonction sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Germe de la génération pseudo-aléatoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 R est un environnement spécialisé dans les calculs statistiques. Voyons comment réaliser de production de statistiques descriptives, de manipulation des distributions de probabilité, de génération de nombres pseudo-aléatoires. 1 Statistiques descriptives Les fonctions permettant de calculer des statistiques descriptives en R sont très nombreuses. Les principales sont présentées ici, mais il en existe d’autres. Le tableau suivant classe les fonctions mentionnées ci-dessous selon le type de mesure calculée (mesure de position, de tendance centrale, de dispersion ou de fréquences) et selon le type de fonctionnement de la fonction (calcul vectoriel, données combinées en une ou plusieurs valeurs). 1 1.1 Fonctions retournant une seule statistique Certaines fonctions de calcul de statistiques descriptives retournent en sortie une seule valeur. C’est le cas des fonctions suivantes : mesures de position : min , max ; mesures de tendance centrale : (https://stat.ethz.ch/R-manual/R-patched/library/base/html/mean.html) , (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/median.html) ; mesure de dispersion : (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/sd.html) (écart-type), (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/IQR.html) (écart interquartile), (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/mad.html) (écart absolu médian). Utilisons le jeu de données cars du package datasets pour présenter quelques exemples. Ce jeu de données contient 50 observations de 2 variables numériques. str(cars) ## 'data.frame': 50 obs. of 2 variables: ## $ speed: num 4 4 7 7 8 9 10 10 10 11 ... ## $ dist : num 2 10 4 22 16 10 18 26 34 17 ... # Moyenne des observations de la variable dist mean(cars$dist) ## [1] 42.98 Si l’objet en entrée a plus d’une dimension, la sortie est tout de même de longueur 1. Donc tous les éléments contenus dans l’objet sont mis en commun pour faire le calcul. max(cars) ## [1] 120 Fonctions which.max et which.min Les fonctions min et max retournent respectivement la valeur la plus petite et la valeur la plus grande parmi les éléments d’un objet. Les fonctions which.max et which.min retournent pour leur part la position dans l’objet du premier maximum ou minimum. which.min(cars$speed) ## [1] 1 Dans l’exemple précédent, il y a en fait deux observations qui prennent la valeur minimum de min(cars$speed) . La commande suivante permet de trouver la position de toutes les observations prenant la valeur minimale. which(cars$speed == min(cars$speed)) ## [1] 1 2 2 1.2 Traitement des données manquantes et argument na.rm Les fonctions min , max , mean , median et sd , ainsi que quelques autres fonctions vues dans ces notes, ont un argument en commun nommé na.rm . Cet argument sert à indiquer à la fonction comment agir en présence de données manquantes ( NA ). Par défaut, na.rm prend la valeur FALSE pour ces fonctions. Cette valeur signifie que les données manquantes ne doivent pas être retirées avant d’effectuer le calcul. Cependant, en présence de données manquantes, ces fonctions ne sont pas en mesure de calculer des statistiques. Par exemple, supposons que nous voulions calculer la médiane des données dans le vecteur suivant. x <- c(3, 6, NA, 8, 11, 15, 23) Si nous ne retirons pas la donnée manquante, nous obtenons le résultat suivant. median(x) ## [1] NA Pour calculer plutôt la médiane des observations non manquantes, il faut donner la valeur TRUE à l’argument na.rm comme suit. median(x, na.rm = TRUE) ## [1] 9.5 Notons que la fonction na.omit permet de retirer les observations manquantes d’un objet R. Si l’objet est un vecteur, les éléments contenant NA sont retirés. na.omit(x) ## [1] 3 6 8 11 15 23 ## attr(,"na.action") ## [1] 3 ## attr(,"class") ## [1] "omit" La fonction na.omit ajoute deux attributs à l’objet, dont un pour identifier les observations retirées. Remarquons que les deux commandes suivantes retournent le même résultat. median(x, na.rm = TRUE) ## [1] 9.5 median(na.omit(x)) ## [1] 9.5 Si la fonction na.omit reçoit en entrée une matrice ou un data frame, elle retire toutes les lignes contenant au moins un NA , comme dans cet exemple : ex_jeu <- data.frame(x, y = c(2, NA, 8, 9, 6, NA, 2)); ex_jeu ## x y ## 1 3 2 ## 2 6 NA ## 3 NA 8 ## 4 8 9 ## 5 11 6 ## 6 15 NA ## 7 23 2 3 na.omit(ex_jeu) ## x y ## 1 3 2 ## 4 8 9 ## 5 11 6 ## 7 23 2 1.3 Fonctions pouvant retourner plusieurs statistiques D’autres fonctions peuvent retourner plus d’une statistique : mesures de position : (https://stat.ethz.ch/R-manual/R-patched/library/base/html/range.html) , (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/quantile.html) ; résumé comprenant plusieurs mesures : (https://stat.ethz.ch/R-manual/R-patched/library/base/html/summary.html) ; variances, covariances et corrélations : var , (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/cor.html) , cor . Fonctions range et quantile La fonction range retourne à la fois le minimum est le maximum, comme dans cet exemple : range(cars$speed) ## [1] 4 25 Une façon simple d’obtenir l’étendue d’observations à partir de la sortie de la fonction range est de procéder comme suit : diff(range(cars$speed)) ## [1] 21 La fonction (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/quantile.html) calcule des quantiles empiriques. Par défaut, elle retourne le minimum, le maximum et les quartiles, comme dans cet exemple : quantile(cars$speed) ## 0% 25% 50% 75% 100% ## 4 12 15 19 25 L’argument probs permet de demander n’importe quels quantiles. Dans l’exemple suivant, les premiers et neuvièmes déciles sont demandés. quantile(cars$speed, probs = c(0.1, 0.9)) ## 10% 90% ## 8.9 23.1 Fonction summary La fonction summary retourne les statistiques suivantes selon l’entrée qu’elle reçoit : vecteur numérique : minimum, premier quartile, médiane, moyenne, troisième quartile, maximum ; facteur : fréquences des modalités (comme la fonction table vue plus loin) ; matrice ou data frame : la fonction summary est appliquée séparément à chacune des colonnes. Utilisons le jeu de données Puromycin du package datasets pour présenter quelques exemples d’utilisation de la fonction summary . Ce jeu de données contient 23 observations de 3 variables, dont deux numériques et une catégorique, stockée sous forme de facteur. str(Puromycin) 4 ## 'data.frame': 23 obs. of 3 variables: ## $ conc : num 0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 ... ## $ rate : num 76 47 97 107 123 139 159 152 191 201 ... ## $ state: Factor w/ 2 levels "treated","untreated": 1 1 1 1 1 1 1 1 1 1 ... ## - attr(*, "reference")= chr "A1.3, p. 269" # Vecteur numérique en entrée : summary(Puromycin$rate) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 47.0 91.5 124.0 126.8 158.5 207.0 # Facteur en entrée : summary(Puromycin$state) ## treated untreated ## 12 11 # Data frame entier en entrée : summary(Puromycin) ## conc rate state ## Min. :0.0200 Min. : 47.0 treated :12 ## 1st Qu.:0.0600 1st Qu.: 91.5 untreated:11 ## Median :0.1100 Median :124.0 ## Mean :0.3122 Mean :126.8 ## 3rd Qu.:0.5600 3rd Qu.:158.5 ## Max. :1.1000 Max. :207.0 Fonctions var , cov et cor La fonction var peut prendre en entrée un vecteur ou un objet à deux dimensions. Si elle reçoit en entrée un vecteur, elle calcule la variance empirique de toutes les valeurs, comme dans cet exemple : var(cars$speed) ## [1] 27.95918 Cependant, si elle reçoit en entrée une matrice ou un data frame de valeurs numériques, elle considère que chaque colonne contient les observations d’une variable aléatoire. Elle va calculer une matrice de variances-covariances, comme dans cet exemple : var(cars) ## speed dist ## speed 27.95918 109.9469 ## dist 109.94694 664.0608 La fonction cov fait exactement le même calcul par défaut. cov(cars) ## speed dist ## speed 27.95918 109.9469 ## dist 109.94694 664.0608 Elle peut cependant calculer des covariances de Kendall ou de Spearman (toutes deux des statistiques non paramétriques basées sur les rangs des observations) au lieu de la covariance classique de Pearson. La fonction (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/cor.html) calcule des corrélations plutôt que des covariances. Elle aussi peut utiliser les définitions de Pearson (par défaut), Kendall et Spearman. Voici un exemple de calcul de matrice de corrélations de Spearman. 5 cor(cars, method = "spearman") ## speed dist ## speed 1.0000000 0.8303568 ## dist 0.8303568 1.0000000 1.4 Fonctions retournant un vecteur de statistiques Certaines fonctions, telles que les suivantes, retournent autant de statistiques qu’il y a d’éléments dans l’objet donné en entrée. mesures de position : cummin , cummax , pmin , pmax ; rangs : (https://stat.ethz.ch/R-manual/R-patched/library/base/html/rank.html) . Fonctions cummin et cummax Les fonctions (https://stat.ethz.ch/R-manual/R-patched/library/base/html/cumsum.html) et cummax calculent les minimums et les maximums cumulatifs. Comme nous pouvons le constater dans l’exemple suivant, la valeur en position i du vecteur retourné par une de ces deux fonctions est la valeur minimale ou maximale dans le sous-vecteur x[1:i] . cummin(x = c(-2, 4, -3, 4, 7, -6, 0)) ## [1] -2 -2 -3 -3 -3 -6 -6 Fonctions pmin et pmax Les fonctions (https://stat.ethz.ch/R-manual/R-patched/library/base/html/Extremes.html) et pmax calculent le minimum et le maximum par position, entre autant de vecteurs que désiré, comme dans l’exemple suivant. pmax( c(-2, 4, -3, 4, 7, -6, 0), c( 1, 2, 3, 4, 5, 6, 7), c( 5, 0, -2, 4, 5, 3, 3) ) ## [1] 5 4 3 4 7 6 7 Ces fonctions sont utiles pour remplacer des valeurs par un seuil. Par exemple, l’instruction suivante permet de remplacer par zéro toute valeur négative contenue dans le vecteur en entrée. pmax(c(-2, 4, -3, 4, 7, -6, 0), 0) ## [1] 0 4 0 4 7 0 0 Fonction rank Certains tests statistiques non paramétriques utilisent des statistiques basées sur les rangs des observations. Voici un exemple d’obtention de ces rangs avec la fonction rank . rank(c(-2, 4, -3, 4, 7, -6, 0)) ## [1] 3.0 5.5 2.0 5.5 7.0 1.0 4.0 Par défaut, en cas d’égalité, le rang moyen est utilisé. Pour changer cette option, il faut modifier la valeur de l’argument ties.method . Dans l’exemple suivant, le rang minimum est retourné en cas d’égalité. rank(c(-2, 4, -3, 4, 7, -6, 0), ties.method = "min") ## [1] 3 5 2 5 7 1 4 6 1.5 Calcul de fréquences Les fonctions (https://stat.ethz.ch/R-manual/R-patched/library/base/html/table.html) , (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/xtabs.html) et (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/ftable.html) permettent de calculer des fréquences. Voici un petit jeu de données pour illustrer l’utilisation de ces fonctions. Il contient des observations concernant 7 individus fictifs : la couleur de leurs yeux, la couleur de leurs cheveux et leur genre. sondage <- data.frame( yeux = c("brun", "brun", "bleu", "brun", "vert", "brun", "bleu" ), cheveux = c("brun", "noir", "blond", "brun", "brun", "blond", "brun" ), genre = c("féminin", "masculin", "féminin", "féminin", "masculin", "féminin", "masculin") ) sondage ## yeux cheveux genre ## 1 brun brun féminin ## 2 brun noir masculin ## 3 bleu blond féminin ## 4 brun brun féminin ## 5 vert brun masculin ## 6 brun blond féminin ## 7 bleu brun masculin Fonctions table La fonction table permet de compter le nombre d’occurrences de chacune des modalités d’une variable catégorique dans des données. Demandons, par exemple, à table de compter le nombre d’individus dans les données sondage classés dans chacune des catégories de couleurs de cheveux. table(sondage$cheveux) ## ## blond brun noir ## 2 4 1 La fonction table produit un tableau de fréquences à une variable si elle reçoit les observations d’une seule variable. Elle peut aussi produire des tableaux de fréquences croisées à deux variables ou plus. # Exemple de tableau de fréquences à deux variables (avec variables nommées) table(yeux = sondage$yeux, cheveux = sondage$cheveux) ## cheveux ## yeux blond brun noir ## bleu 1 1 0 ## brun 1 2 1 ## vert 0 1 0 # Exemple de tableau de fréquences à trois variables (data frame en entrée à table) t3 <- table(sondage) t3 ##,, genre = féminin ## ## cheveux ## yeux blond brun noir ## bleu 1 0 0 ## brun 1 2 0 ## vert 0 0 0 ## 7 ##,, genre = masculin ## ## cheveux ## yeux blond brun noir ## bleu 0 1 0 ## brun 0 0 1 ## vert 0 1 0 Fonctions ftable La fonction ftable retourne un tableau de fréquences sous la forme d’une table « plate » (en anglais flat , d’où le f dans le nom de la fonction) dans le cas d’un croisement de 3 variables ou plus, plutôt que sous la forme d’un array comme le fait la fonction table . Elle accepte les mêmes types d’entrées que table (série d’objets atomiques à une dimension ou objet récursif dont les éléments sont interprétables en facteurs) et peut aussi recevoir une sortie de la fonction table , comme dans l’exemple suivant. ftable(t3) ## genre féminin masculin ## yeux cheveux ## bleu blond 1 0 ## brun 0 1 ## noir 0 0 ## brun blond 1 0 ## brun 2 0 ## noir 0 1 ## vert blond 0 0 ## brun 0 1 ## noir 0 0 Fonctions xtabs La fonction xtabs fait le même calcul que les fonctions précédentes, mais elle prend en entrée une formule. Le tableau de fréquences à deux variables créé précédemment peut être réobtenu de la façon suivante avec xtabs . xtabs(~ yeux + cheveux, data = sondage) ## cheveux ## yeux blond brun noir ## bleu 1 1 0 ## brun 1 2 1 ## vert 0 1 0 La fonction xtabs est utile lorsque les données que nous avons en main contiennent déjà des fréquences, car il est possible d’inclure une variable réponse contenant des dénombrements dans la formule que nous lui fournissons en entrée. Par exemple, xtabs permet de facilement retrouver le tableau de fréquences marginales croisées entre les variables yeux et cheveux à partir du tableau de fréquences à trois variables produit précédemment mis sous forme de data frame, qui a l’allure suivante. t3_df <- as.data.frame(t3) t3_df ## yeux cheveux genre Freq ## 1 bleu blond féminin 1 ## 2 brun blond féminin 1 ## 3 vert blond féminin 0 ## 4 bleu brun féminin 0 8 ## 5 brun brun féminin 2 ## 6 vert brun féminin 0 ## 7 bleu noir féminin 0 ## 8 brun noir féminin 0 ## 9 vert noir féminin 0 ## 10 bleu blond masculin 0 ## 11 brun blond masculin 0 ## 12 vert blond masculin 0 ## 13 bleu brun masculin 1 ## 14 brun brun masculin 0 ## 15 vert brun masculin 1 ## 16 bleu noir masculin 0 ## 17 brun noir masculin 1 ## 18 vert noir masculin 0 Il suffit de procéder comme suit : t2 <- xtabs(Freq ~ yeux + cheveux, data = t3_df) t2 ## cheveux ## yeux blond brun noir ## bleu 1 1 0 ## brun 1 2 1 ## vert 0 1 0 ## **2 Distributions de probabilité** Le package stats de l’installation de base de R comprend, pour plusieurs distributions de probabilité, des fonctions de calcul de : la fonction de densité (forme d* où * change selon la distribution), la fonction de répartition (forme p* ) et la fonction quantile (forme q* ). La fiche d’aide ouverte par la commande help(Distributions) énumère toutes les distributions de probabilité offertes dans le package stats . Il existe aussi des fonctions relatives à d’autres distributions de probabilité dans des packages sur le CRAN (voir https://CRAN.R-project.org/view=Distributions pour découvrir ce qui est offert). 2.1 Fonction de densité Les fonctions R implémentant des fonctions de densité ont un nom qui débute par le lettre d pour density . Dans le cas d’une variable aléatoire discrète, la fonction de densité est plus justement appelée fonction de masse. Il s’agit alors d’une probabilité, pour une variable aléatoire suivant une certaine distribution, de prendre une certaine valeur. Exemple : distribution binomiale Soit X une variable aléatoire représentant le nombre de 6 obtenus lors de 5 lancés d’un dé. Cette variable aléatoire suit une distribution binomiale de paramètres n = 5 et p = 1 / 6, donc X ∼ Bin ( n = 5 , p = 1 / 6). Calculons P ( X = 2), soit la probabilité que la variable aléatoire X prenne la valeur 2. 9 dbinom(x = 2, size = 5, prob = 1/6) ## [1] 0.160751 La fonction dbinom implémente donc la fonction de masse d’une distribution binomiale. Les arguments size et prob sont les paramètres n et p de la distribution, selon la notation utilisée ci-dessus. Les fonctions de la famille d* peuvent calculer plusieurs valeurs de densité par un seul appel de la fonction, car celles-ci acceptent des valeurs d’arguments de longueur supérieure à 1. Ces fonctions travaillent donc de façon vectorielle, comme presque toutes les fonctions de calcul en R. Voici un exemple de code R permettant de représenter graphiquement la densité Bin ( n = 5 , p = 1 / 6) complète, en mettant en évidence la valeur calculée ci-dessus, soit P ( X = 2). barplot( height = dbinom(0:5, size = 5, prob = 1/6), names.arg = 0:5, main = "Densité Binomiale(5,1/6)", xlab = "x", ylab = "dbinom(x, size = 5, prob = 1/6)" ) barplot( height = c(0, 0, dbinom(2, size = 5, prob = 1/6), 0, 0, 0), col = "blue", add = TRUE ) Densité Binomiale(5,1/6) 0 1 2 3 4 5 x Exemple : distribution normale standard Supposons maintenant que X est une variable aléatoire continue de distribution normale standard, donc X ∼ N ( µ = 0 , σ [2] = 1). La valeur de la fonction de densité pour cette distribution en la valeur x = 1, souvent notée fX (1), vaut : dnorm(x = 1) ## [1] 0.2419707 10 Voici une représentation graphique de la densité complète dans laquelle la valeur calculée ci-dessus est mise en évidence. curve(expr = dnorm, xlim = c(-3, 3), main = "Densité N(0,1)") segments( x0 = c(-4, 1), y0 = dnorm(1), x1 = 1, y1 = c(dnorm(1), -1), lty = 2, col = "blue" ) Densité N(0,1) 0.3 0.4 Col2 Col3 0.3 0.4 2 3 Ici, nous n’avons pas eu besoin de fournir des valeurs aux arguments de la fonction dnorm relatifs aux paramètres de la distribution, parce que nous avons utilisés leurs valeurs par défaut. Ces paramètres, pour la densité N ( µ, σ [2] ), sont représentés par les arguments mean = µ et sd = σ (remarquez que l’argument de la fonction R représente l’écart-type, pas la variance). Tout comme le premier argument, nommé x , les arguments des fonctions de la famille d* représentant des paramètres de la distribution acceptent aussi en entrée plus d’une valeur. Voici un exemple, qui permet de calculer la densité en x = 1 pour X ∼ N ( µ = − 2 , σ [2] = 1), X ∼ N ( µ = 0 , σ [2] = 2 . 25) et X ∼ N ( µ = 1 , σ [2] = 4) en un seul appel à la fonction dnorm . dnorm(x = 1, mean = c(-2, 0, 1), sd = c(1, 1.5, 2)) ## [1] 0.004431848 0.212965337 0.199471140 2.2 Fonction de répartition La fonction de répartition d’une variable aléatoire X est définie par FX ( x ) = P ( X ≤ x ). Il s’agit donc toujours d’une probabilité, d’où le p au début des noms des fonctions R implémentant des fonctions de répartition. Exemple : distribution normale standard Prenons encore comme exemple la distribution normale standard. Nous avons donc X ∼ N (0 , 1). Calculons la valeur de la fonction de répartition de cette variable aléatoire en x = 1. pnorm(q = 1) ## [1] 0.8413447 Il s’agit de la valeur de la probabilité P ( X ≤ 1). 11 Le premier argument des fonctions de la famille p* ne se nomme pas x , il se nomme plutôt q . Cette lettre réfère au mot quantile et souligne le lien entre les fonctions de répartition et les fonctions quantiles. Les arguments suivants des fonctions de la famille p* permettent de spécifier les valeurs des paramètres de la distribution, comme pour les fonctions de la famille p* . Voici une illustration graphique du lien entre la fonction de densité et la fonction de répartition en utilisant la distribution normale standard. La valeur de la fonction de répartition en un point quelconque x (dans le graphique x = 1) est égale à l’aire sous la courbe de densité entre −∞ et x . par_default <- par(mfrow = c(1, 2)) curve(expr = dnorm, xlim = c(-3, 3), main = "Densité N(0,1)") x <- seq(from = -3, to = 1, length = 100) polygon(x = c(x, 1, -3), y = c(dnorm(x), 0, 0), col = "blue") curve( expr = pnorm, xlim = c(-3, 3), main = "Fonction de répartition N(0,1)", xname = "q" ) segments( x0 = c(-4, 1), y0 = pnorm(1), x1 = 1, y1 = c(pnorm(1), -1), lty = 2, col = "blue" ) par(par_default) Densité N(0,1) −3 −2 −1 0 1 2 3 x Fonction de répartition N(0,1) 1.0 Col2 Col3 −3 −2 −1 0.0 0.2 0.4 0.6 0.8 −3 −2 −1 0.0 0.2 0.4 0.6 0.8 0 1 2 3 q 2.3 Fonction quantile La fonction quantile est l’inverse généralisé de la fonction de répartition. Les fonctions R implémentant des fonctions quantile ont un nom qui débute par le lettre q pour quantile . Exemple : distribution normale standard Voyons de quoi à l’air la fonction quantile de la distribution normale standard. Premièrement, calculons la valeur de la fonction quantile en un point, disons en p = 0 . 8413447. 12 qnorm(p = 0.8413447) ## [1] 0.9999998 Il s’agit de la valeur x pour laquelle P ( X ≤ x ) = 0 . 8413447, où X ∼ N (0 , 1). Le premier argument d’une fonction de la famille q* se nomme p . Cette notation peut nous aider à nous rappeler que cet argument représente une probabilité et accepte donc seulement des valeurs entre 0 et 1. Encore une fois, les arguments suivants des fonctions de la famille q* permettent de spécifier les valeurs des paramètres de la distribution. Le graphique suivant illustre le lien entre la fonction de répartition et la fonction quantile. Le graphique de la fonction quantile est simplement obtenu en inversant les axes du graphique de la fonction de répartition. par_default <- par(mfrow = c(1,2), pty = "s") curve( expr = pnorm, xlim = c(-3, 3), n = 1000, ylim = c(0, 1), asp = 6, main = "Fonction de répartition N(0,1)", xname = "q" ) segments( x0 = c(-4, 1), y0 = pnorm(1), x1 = 1, y1 = c(pnorm(1), -1), lty = 2, col = "blue" ) curve( expr = qnorm, xlim = c(0, 1), n = 1000, ylim = c(-3, 3), asp = 1/6, main = "Fonction quantile N(0,1)", xname = "p" ) segments( x0 = c(-0.2, pnorm(1)), y0 = qnorm(pnorm(1)), x1 = pnorm(1), y1 = c(qnorm(pnorm(1)), -4), lty = 2, col = "blue" ) par(par_default) Fonction de répartition N(0,1) Fonction quantile N(0,1) 2 3 Col2 Col3 0.0 0.2 0.4 0.6 0. −3 −2 −1 0 0.0 0.2 0.4 0.6 0. −3 −2 −1 0 1. p 1.0 Col2 Col3 −3 −2 −1 0 1 0.0 0.2 0.4 0.6 0. q −3 −2 −1 0 1 0.0 0.2 0.4 0.6 0. q 2 3 13 3 Génération de nombres pseudo-aléatoires Les fonctions r* permettent de générer pseudo-aléatoirement des observations selon une certaine distribution désignée par * . La lettre r au début de leur nom signifie random . Par exemple, voici la représentation graphique de 3 échantillons générés aléatoirement selon 3 distributions différentes : des distributions normale, uniforme continue et khi-deux. Pour chaque échantillon, nous traçons l’histogramme des observations simulées pour représenter leur densité empirique. Nous superposons à cet histogramme la courbe de densité de la distribution théorique à partir de laquelle les observations ont été générées. par_default <- par(mfrow = c(1, 3)) # Densité normale d'espérance 5 et de variance 4 x1 <- rnorm(1000, mean = 5, sd = 2) hist(x = x1, freq = FALSE, ylab = "Proportion", main = "Densité empirique de x1") curve(expr = dnorm(x, mean = 5, sd = 2), add = TRUE, col = "blue") title(main = "courbe bleue = densité visée N(5,4)", line = 0.5) # Densité uniforme continue entre 2 et 8 x2 <- runif(1000, min = 2, max = 8) hist(x = x2, freq = FALSE, ylab = "Proportion", main = "Densité empirique de x2") curve(expr = dunif(x, min = 2, max = 8), add = TRUE, col = "blue") title(main = "courbe bleue = densité visée U(2,8)", line = 0.5) # Densité chi-carré à 4 degrés de liberté x3 <- rchisq(1000, df = 4) hist(x = x3, freq = FALSE, ylab = "Proportion", main = "Densité empirique de x3") curve(expr = dchisq(x, df = 4), add = TRUE, col = "blue") title(main = expression(paste("courbe bleue = densité visée ", chi[4]ˆ2)), line = 0.5) par(par_default) 14 Densité empirique de x1 courbe bleue = densité visée N(5,4) −2 0 2 4 6 8 10 x1 2 3 4 5 6 7 8 x2 Densité empirique de x2 courbe bleue = densité visée U(2,8) Densité empirique de x3 courbe bleue = densité visée χ24 0 5 10 15 x3 Comme nous pouvons le constater sur ces graphiques, la distribution empirique des observations générées avec ces fonctions se rapproche vraiment de la distribution théorique demandée. En générant un nombre encore plus grand d’observations (ici nous en avons généré 1000 pour chaque distribution), la densité empirique se rapprocherait encore plus de la densité théorique. 3.1 Fonction sample La fonction (https://stat.ethz.ch/R-manual/R-patched/library/base/html/sample.html) permet de tirer un échantillon aléatoire parmi un ensemble d’éléments fourni à son argument x sous forme de vecteur. sample(x = c("Luc", "Kim", "Paul", "Ève")) ## [1] "Kim" "Luc" "Paul" "Ève" Ce vecteur peut contenir des données de n’importe quel type. sample(x = 1:6) ## [1] 4 5 3 2 6 1 Par défaut, sample effectue un tirage sans remise et sélectionne autant d’éléments que l’ensemble de départ en contient. La fonction effectue donc une permutation aléatoire des éléments de l’ensemble de départ. L’argument size permet cependant de contrôler la taille de l’échantillon tiré. sample(x = 1:6, size = 5) ## [1] 4 5 2 1 6 L’argument replace permet quant à lui de spécifier si le tirage doit être effectué avec ( replace = TRUE ) ou sans remise ( replace = FALSE , valeur par défaut). sample(x = 1:6, size = 10, replace = TRUE) ## [1] 5 1 6 5 3 1 6 6 3 4 Il est aussi possible d’attribuer des probabilité de sélection à chacun des éléments grâce à l’argument prob . Par défaut, tous les éléments ont une probabilité égale d’être tiré. 15 sample(x = 1:6, size = 20, replace = TRUE, prob = c(1/2, rep(1/10, times = 5))) ## [1] 5 3 3 3 1 1 1 1 6 1 6 1 6 3 1 1 4 3 2 6 La fonction (https://stat.ethz.ch/R-manual/R-patched/library/base/html/sample.html) est très similaire à la fonction sample , mais elle prend comme premier argument un seul entier, n , et tire aléatoirement size entiers entre 1 et n . sample.int(n = 6, size = 4) ## [1] 5 2 6 3 Si nous souhaitons sélectionner aléatoirement des observations (lignes) dans un jeu de données, il est d’usage de sélectionner d’abord des entiers compris entre 1 et le nombre total d’observations, puis d’extraire du jeu de donnée les observations sur les lignes portant les numéros sélectionnés. Voici un exemple, utilisant célèbres données iris, incluses dans l’installation de base de R dans le data frame nommé (https://stat.ethz.ch/R-manual/R-patched/library/datasets/html/iris.html) (du package datasets ). str(iris) ## 'data.frame': 150 obs. of 5 variables: ## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ... ## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ... ## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ... ## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ... ## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ... index_ech <- sample.int(n = nrow(iris), size = 5, replace = FALSE) index_ech ## [1] 9 53 84 118 80 iris_ech <- iris iris_ech ## Sepal.Length Sepal.Width Petal.Length Petal.Width Species ## 9 4.4 2.9 1.4 0.2 setosa ## 53 6.9 3.1 4.9 1.5 versicolor ## 84 6.0 2.7 5.1 1.6 versicolor ## 118 7.7 3.8 6.7 2.2 virginica ## 80 5.7 2.6 3.5 1.0 versicolor 3.2 Germe de la génération pseudo-aléatoire Les nombres générés avec les fonctions r* et les échantillons tirés avec sample sont qualifiés de pseudoaléatoires, car ils proviennent d’un algorithme déterministe qui tente de reproduire le hasard. Un tel algorithme est nommé en anglais random number generator (RNG) ou pseudo-random number generator . La fiche d’aide ouverte par la commande (https://stat.ethz.ch/R-manual/R-patched/library/base/html/Random.html) contient de l’information sur la génération de nombres pseudoaléatoires en R. Plusieurs RNG sont implémentés en R. Celui utilisé par défaut, nommé Mersenne Twister, a été choisi parce qu’il est réputé être bon. Sans entrer dans les détails du fonctionnement des RNG implémentés dans le package base de R, il faut savoir qu’ils travaillent tous à partir d’une séquence de nombres appelée germe (en anglais seed ). En contrôlant ce germe, il est possible de générer de nouveau, autant de fois que désiré, les mêmes valeurs. Par défaut, R contrôle le germe des RNG de façon automatique. Chaque fois qu’une commande faisant intervenir un RNG est évaluée, R crée un nouveau germe à partir, notamment, de l’heure à laquelle la commande est soumise. Ainsi, deux générations pseudo-aléatoires consécutives ne produisent en général pas le même résultat. 16 sample(x = letters, size = 5) ## [1] "c" "v" "a" "u" "x" sample(x = letters, size = 5) ## [1] "t" "k" "m" "f" "p" En tout temps, il est possible de connaître le germe du RNG en R. Il est stocké dans un objet nommé .Random.seed . Cet objet est un vecteur numérique de longueur 626 pour le RNG Mersenne Twister. str(.Random.seed) ## int [1:626] 10403 70 -738736558 -2074313779 -2024393337 -394776466 -210556856 122496.. sample(x = letters, size = 5) ## [1] "p" "b" "g" "n" "c" str(.Random.seed) ## int [1:626] 10403 75 -738736558 -2074313779 -2024393337 -394776466 -210556856 122496.. sample(x = letters, size = 5) ## [1] "t" "n" "f" "e" "d" str(.Random.seed) ## int [1:626] 10403 82 -738736558 -2074313779 -2024393337 -394776466 -210556856 122496.. Nous constatons qu’au moins un élément de ce vecteur (le deuxième élément) change à chaque fois que nous appelons la fonction sample . La fonction (https://stat.ethz.ch/R-manual/R-patched/library/base/html/Random.html) permet de fixer le germe du RNG à partir d’une seule valeur entière. Par exemple, la commande suivante : set.seed(753) À n’importe quel moment, dans n’importe qu’elle session R, nous obtenons le même échantillon set.seed(753) sample(letters, size = 5) ## [1] "v" "r" "g" "i" "t" si nous soumettons la commande set.seed(753) avant la commande sample(letters, size = 5) . L’utilité de fixer le germe d’une génération pseudo-aléatoire est de rendre son résultat reproductible. 17