R Notebook Tutorial: Continuous Random Variables
Introduction
Le but de ce R Notebook (Tutorial) est, entre autres, de vous montrer comment travailler avec des histogrammes
et pourquoi il faut parfois transformer des données avant de les visualiser. En définitive, il vous mènera à
travers une analyse du concept de densité d’une variable continue en illustrant comment on construit des
fonctions de densité et de répartition sur R.
1) Ouvrir R et définir le répertoire de travail.
setwd("path/to/your/created/directory")
2) Créer un nouveau script et le nommer tutorialcrv.R. Travailler à partir de ce script.
3) Télécharger les fichiers [Link] et [Link] dans le
répertoire de travail.
4) Ouvrir les deux fichiers avec un éditeur de texte. Qu’est-ce que vous voyez? En quoi sont-ils différents?
Notez que les fichiers ont des dimensions différentes - ce sont des fichiers pour deux exercices différents.
5) Importer le contenu des deux fichiers sur R.
Forbes = [Link]("[Link]", header = T)
attach(Forbes) #permet d'ajouter les noms des variables de 'Forbes' au système R.
CEO = [Link]("[Link]", header = T)
attach(CEO) #permet d'ajouter les noms des variables de 'CEO' au système R.
1) Motivation pour l’histogramme
Parfois, les valeurs prises par un caractère sont si nombreuses qu’un tableau de fréquences ou une fonction de
densité deviennent peu lisibles.
1) Afin d’illustrer ce problème, nous allons utiliser un exemple où nous analyserons les fortunes des gens
les plus riches (selon le magazine “Forbes”) d’une année non-spécifiée.
# Insérez les commandes suivantes dans votre script R:
f_a = [Link](table(Forbes$Age))
barplot(f_a, xlab = "Age", ylab = "Proportions", main = "Barplot: Forbes Age",
ylim = c(0, 0.04))
box()
1
Barplot: Forbes Age
0.04
0.03
Proportions
0.02
0.01
0.00
27 33 38 43 48 53 58 63 68 73 78 83 88 93
Age
2) Comme solution, un regroupement des valeurs sous la forme d’un histogramme est conseillé. Les
commandes en dessous vous montrent comment construire un histogramme sur R.
par(mfrow = c(1, 2))
barplot(f_a, xlab = "Age", ylab = "Proportions", main = "Barplot: Forbes Age",
ylim = c(0, 0.04))
box()
hist(Forbes$Age, probability = TRUE, xlab = "Age", main = "Histogram: Forbes Age",
ylim = c(0, 0.04))
box()
2
Barplot: Forbes Age Histogram: Forbes Age
0.04
0.04
0.03
0.03
Proportions
Density
0.02
0.02
0.01
0.01
0.00
0.00
27 40 52 64 76 88 40 60 80 100
Age Age
3) Ensuite, regardez comment l’aspect de l’histogramme peut changer en variant le paramètre “nclass”:
par(mfrow = c(2, 2))
hist(Forbes$Age, nclass = 4, probability = TRUE, xlab = "Age",
main = "nclass=4", ylim = c(0, 0.04))
box()
hist(Forbes$Age, nclass = 8, probability = TRUE, xlab = "Age",
main = "nclass=8", ylim = c(0, 0.04))
box()
hist(Forbes$Age, nclass = 16, probability = TRUE, xlab = "Age",
main = "nclass=16", ylim = c(0, 0.04))
box()
barplot(f_a, xlab = "Age", ylab = "Proportions", main = "nclass=automatic",
ylim = c(0, 0.04))
box()
3
nclass=4 nclass=8
0.04
0.04
0.03
0.03
Density
Density
0.02
0.02
0.01
0.01
0.00
0.00
20 40 60 80 100 20 40 60 80 100
Age Age
nclass=16 nclass=automatic
0.04
0.04
0.03
0.03
Proportions
Density
0.02
0.02
0.01
0.01
0.00
0.00
40 60 80 100 27 38 48 58 68 78 88
Age Age
4) En jouant un peu avec les paramètres, vous découvrirez vite:
i) Quel est le role du paramètre “probability” dans la fonction “hist”,
ii) À quoi sert la commande “box()”.
à mettre un cadre, carré à l'histogramme
4
2) Transformation logarithmique des données
Parfois, il est judicieux de transformer au préalable les données qu’on dispose afin qu’elle soient plus facilement
interprétables par la suite.
1) Avant de se pencher sur un exemple concret, nous allons commencer par analyser les données “CEO” et
la variabilité de la rémunération d’un CEO en Amérique d’une année non-spécifiée.
plot(CEO[, "Compensation"], ylab = "Compensation", main = "Plot: Compensation")
Plot: Compensation
8e+07
Compensation
4e+07
0e+00
0 5000 10000 15000 20000 25000
Index
dim(CEO)
## [1] 27643 3
names(CEO)
## [1] "Name" "Compensation" "Firm"
CEO[1:10, ]
## Name Compensation
## 1 1 Jay Sugarman 90192027
## 2 2 Steven P Jobs 74750000
## 3 3 Larry C Glasscock 46212719
## 4 4 Eugene N Melnyk 41917908
## 5 5 L George Klaus 39723164
## 6 6 Mario J Gabelli 38702110
## 7 7 Warren J Spector 38255431
## 8 8 James E Cayne 33925412
## 9 9 H Lawrence Culp Jr. 31083042
5
## 10 10 Richard E Dauch 30957693
## Firm
## 1 iStar Financial Inc.
## 2 Apple Computer Inc.
## 3 Anthem Inc.
## 4 Biovail Corporation
## 5 Epicor Software Corporation
## 6 Gabelli Asset Management Inc.
## 7 The Bear Stearns Companies Inc.
## 8 The Bear Stearns Companies Inc.
## 9 Danaher Corporation
## 10 American Axle & Manufact Holdings Inc.
range(Compensation)
## [1] 0 90192027
2) La fonction “range” nous montre que certains dirigeants ont une compensation égale a zéro. Voici le
nombre de dirigeants dans cette situation:
nrow(subset(CEO, CEO[, "Compensation"] == 0))
## [1] 434
# Alternativement, comme nous avons ajouté les noms des
# variables de 'CEO' au système R auparavant
# ('attach(CEO)'), vous pouvez aussi accèder à
# l'information demandée comme suit:
nrow(subset(CEO, Compensation == 0))
## [1] 434
# Où bien vous l'analysez de la façon suivane:
length(Compensation[Compensation == 0])
## [1] 434
# Raisonement économique: Ces CEO, ont-elles/ils vraiment
# gagné zéro US-Dollars?
3) Essayons maintenant de présenter la variable représentant la rémunération (Compensation) comme si
elle était de type discrète avec un “barplot”. Vous verrez immédiatement pourquoi cela n’est pas la
meilleure solution pour présenter ces données. . .
barplot([Link](table(Compensation)), ylab = "Proportions",
main = "Barplot: Compensations", ylim = c(0, 0.02))
box()
6
Barplot: Compensations
0.020
Proportions
0.010
0.000
0 134810 201308 271875 370668 529000 888624 2901990
4) Est-ce qu’un histogramme à travers la fonction “hist” est un bonne solution dans cette situation ? Nous
pouvons constater que ce n’est pas la meilleure solution non plus. . .
hist(Compensation, ylim = c(0, 30000))
box()
Histogram of Compensation
30000
20000
Frequency
10000
0
0e+00 2e+07 4e+07 6e+07 8e+07
Compensation
5) À ce stade, nous pouvons déjà arriver à certaines conclusions:
7
• La variable “Compensation” pourrait profiter d’une modélisation en tant que variable continue, et non
pas en tant que variable discrète,
• La distribution des salaires (Compensation) est fortement asymétrique,
• Le jeu de données que nous sommes en train d’analyser a une taille conséquente.
Combien d’individus dans notre population ont des compensation au-delà de 2e+7 (i.e. 2*10ˆ7, autrement
dit: USD 20 Mio).
# Seulement 33 CEOs (0.11% des individus représentés dans
# nos données), gagnent plus que USD 20 Mio:
Compensation[Compensation > 2 * 10ˆ7]
## [1] 90192027 74750000 46212719 41917908 39723164 38702110
## [7] 38255431 33925412 31083042 30957693 30674065 29556378
## [13] 28853282 27783140 26812148 26331461 26159146 25344265
## [19] 25096042 23852064 22863045 22634271 22303437 22217160
## [25] 22142704 21530638 21400579 20921299 20569000 20199625
## [31] 20123253 20070000 20070000
length(Compensation[Compensation > 2 * 10ˆ7])
## [1] 33
length(Compensation[Compensation > 2 * 10ˆ7])/length(Compensation)
## [1] 0.001193792
# Les rémunérations de ces 33 personnes s'étendent
# d'approx. USD 20 Mio. à approx. USD 90 Mio:
range(Compensation[Compensation > 2 * 10ˆ7])
## [1] 20070000 90192027
# Concernant notre situation: cette minorité de 33
# personnes est problématique car elle occupera une grande
# partie de l'espace visuel (à cause de la plage de valeurs
# allant de USD 20 Mio à USD 90 Mio) de n'importe quel
# graphique que l'on puisse faire.
# Par contre, le reste des individus de notre population
# (représentant plus de 99% de notre population de CEOs)
# seront situés dans les premiers 20 Mio. de Dollars. No
# graph will resist - peu importe les graphiques que l'on
# fasse, si nous ne prenons pas en considération cette
# asymétrie (skewness),ils ne seront jamais très
# esthétiques.
6) Solution possible: Transformer les données de façon à ce que les (grandes) différences entre les
rémunérations ne soient plus gonflées, mais plutôt aplaties.
Une façon d’arriver à cela est de logarithmiser les données (“log10”, comme vous l’avez vu-e-s au Gymnase /
à l’École Professionnelle). Voici un exemple pour une séquence de valeurs “x” entre 1 et 1000:
x = seq(1, 10ˆ3, 10)
plot(x, log10(x), type = "l", xlab = "example values (x)", ylab = "log10(x)",
main = "Log10 transformation")
box()
# Notez le saut sur l'ordonnée (y-axis) quand on monte par
8
# exemple de x=10 à x=110:
abline(v = 10, lty = 2)
abline(v = 110, lty = 2)
abline(h = log10(10), lty = 2, lwd = 2, col = "red")
abline(h = log10(110), lty = 2, lwd = 2, col = "red")
# Notez le saut sur l'ordonnée (y-axis) quand on monte par
# exemple de x=810 à x=910:
abline(v = 810, lty = 2)
abline(v = 910, lty = 2)
abline(h = log10(810), lty = 2, lwd = 2, col = "red")
abline(h = log10(910), lty = 2, lwd = 2, col = "red")
Log10 transformation
3.0
2.5
2.0
log10(x)
1.5
1.0
0.5
0.0
0 200 400 600 800 1000
example values (x)
7) Deux élément sont à souligner:
• Un des grands avantages de l’utilisation des “log10” est que malgré la transformation, les chiffres restent
interprétables. Comme vous l’avez appris au Gymnase / à l’École Professionnelle, le log10 de 1000 est
égal à 3, parce qu’il représente la puissance qui nous mène de 10 à 1000 (10ˆ3 = 1000).
• Le log10 de 0 n’existe pas, car aucune puissance peut mener 10 à devenir 0.
En sachant cela, nous allons commencer par nettoyer la variable “Compensation” des valeurs égales à zéro.
Puis, nous allons représenter nos données logarithmisées à l’aide d’un historgramme:
9
compensation_nette = subset(CEO$Compensation, Compensation >
0)
compensation_nette_et_logrithmisee = log10(compensation_nette)
# Juste pour qu'on ait une variable avec un nom plus court
# pour continuer ;-):
complog = compensation_nette_et_logrithmisee
par(mfrow = c(1, 2))
hist(complog, probability = TRUE, ylim = c(0, 1.5))
# nombre de colonnes défini par R.
box()
hist(complog, probability = TRUE, ylim = c(0, 1.5), n = 50)
# avec plus de colonnes visuelles.
box()
Histogram of complog Histogram of complog
1.5
1.5
1.0
1.0
Density
Density
0.5
0.5
0.0
0.0
0 2 4 6 8 0 2 4 6 8
complog complog
3) Densité
En augmentant le nombre d’intervalles dans un histogramme, nous pouvons voir apparaître une fonction
continue. . .
1) Analysez ce processus en créant des histogrammes avec un nombre croissant de colonnes pour le log10
de la variable “Compensation” de la base de données “CEO”.
Présentez la fonction de densité (probability distibution function, i.e. pdf) du log10 des rémunérations et
superposez-la aux histogrammes que vous avez crées auparavant.
10
Faites attention à l’interprétation de l’ordonnée (y-axis) d’une fonction de densité:
[[Link]
par(mfrow = c(2, 2))
hist(complog, probability = TRUE, 15, col = "green", main = "Density complog (15 bins)",
ylim = c(0, 1.5))
lines(density(complog), lwd = 3, col = "red")
box()
hist(complog, probability = TRUE, 30, col = "green", main = "Density complog (30 bins)",
ylim = c(0, 1.5))
lines(density(complog), lwd = 3, col = "red")
box()
hist(complog, probability = TRUE, 60, col = "green", main = "Density complog (60 bins)",
ylim = c(0, 1.5))
lines(density(complog), lwd = 3, col = "red")
box()
hist(complog, probability = TRUE, 120, col = "green", main = "Density complog (120 bins)",
ylim = c(0, 1.5))
lines(density(complog), lwd = 3, col = "red")
box()
11
1.5
Density complog (15 bins) Density complog (30 bins)
1.5
1.0
1.0
Density
Density
0.5
0.5
0.0
0.0
0 2 4 6 8 0 2 4 6 8
complog complog
Density complog (60 bins) Density complog (120 bins)
1.5
1.5
1.0
1.0
Density
Density
0.5
0.5
0.0
0.0
0 2 4 6 8 0 2 4 6 8
complog complog
2) Créez la fonction de répartition (cumulative distibution function, i.e. cdf) du log10 de la variable
“Compensation” et répondez aux questions suivantes:
i) De combien est P(compensation < USD 150’000)?
ii) De combien est P(compensation < USD 360’000)?
Flog = ecdf(complog) #Flog est modélisée pour être une fonction R!
xx = seq(-1, 10, 0.01)
plot(xx, Flog(xx), type = "l", xlab = "x=complog", ylab = "F(x)",
main = "CDF: complog")
12
CDF: complog
1.0
0.8
0.6
F(x)
0.4
0.2
0.0
0 2 4 6 8 10
x=complog
cdf_analyse = ecdf(Compensation)
cdf_analyse(150000) #réponse à la question i).
## [1] 0.1658286
cdf_analyse(360000) #réponse à la question ii).
## [1] 0.5578266
3) Représentez la fonction de densité (PDF) et la fonction de répartition (CDF) du log10 des rémunérations
l’une à côté de l’autre - et voilà, vous avez représenté la distribution de la variable log10(Compensation)
avec succès :-)!
par(mfrow = c(1, 2))
plot(density(complog), main = "PDF: complog", ylim = c(0, 1.2))
plot(xx, Flog(xx), type = "l", xlab = "x=complog", ylab = "F(x)",
main = "CDF: complog", ylim = c(0, 1.2))
13
PDF: complog CDF: complog
1.2
1.2
1.0
1.0
0.8
0.8
Density
F(x)
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0 2 4 6 8 0.0 0 2 4 6 8 10
N = 27209 Bandwidth = 0.04511 x=complog
14