Plan

1. Objets et Opérations

2. Les fonctions
Notes de Cours sur le logiciel R
3. Les graphiques

Anne PHILIPPE 4. Graphiques avec ggplot2

22 janvier 2018
5. Structures de contrôle et Itérations
Université de Nantes,
Laboratoire de Mathématiques Jean Leray
email : [email protected]
6. Autour des lois de probabilités

7. Outils graphiques en statistique

8. Inférence statistique
1 2

9. Séries Chronologiques

Installation Documentations

• Documents sur le logiciel R :

Le logiciel R est un freeware disponible sur le site http://www.math.sciences.univ-nantes.fr/⇠philippe/R_freeware.html

http://cran.r-project.org/
• Site consacré aux graphiques
Il existe des versions
addictedtor.free.fr/graphiques/
• Windows • Collection spécifique UseR chez Springer
• MacOS X • Plus de 80 livres,
• Linux ... par exemple
• Introductory Statistics With R
Outils disponibles : • Bayesian Computation With R
• Applied Statistical Genetics With R :
• un langage de programmation orienté objet
• Generalized Additive Models : An Introduction with R
• des fonctions de "base" • Extending the Linear Model With R
• des librairies/packages complémentaires (1800 sur le site CRAN) • Time Series Analysis And Its Applications : With R Examples

3 4
Au démarrage Sous linux

> apparaît
automatiquement en
début de chaque ligne
de commandes
+ apparaît en début de
ligne si la ligne
précédente est
incomplète

5 6

Utiliser l’aide Éditeur

Sous MacOS et Windows, un éditeur de texte intégré au logiciel R

> help ( " p l o t " )
> ? plot
> help . search ( " p l o t " )
> ?? p l o t

Les démos :
# pour o b t e n i r l a l i s t e d e s demos
> demo ( )
> demo ( g r a p h i c s )

Les exemples :
La fonction example exécute les exemples
généralement inclus à la fin des fichiers • Ctrl R exécute la ligne sur laquelle se trouve le curseur ou les lignes
d’aide. d’un bloc sélectionné.

> example (FUN)

• source("nom-du-fichier.R") pour exécuter le code contenu dans
le fichier nom-du-fichier.R

7 8
Librairies

• Toutes les librairies ne sont pas chargées au lancement du logiciel

• library() retourne la liste des librairies installées.
Objets et Opérations
• library(LIB) charge la librairie LIB
• library(help = LIB) retourne la liste des fonctions de la librairie
LIB
• search(), searchpaths() retourne la liste des librairies chargées.

Opérations élémentaires Objets

1. Opérations élémentaires sur les scalaires : ⇤, , +, /, ˆ

>2+4
6 Les objets de base sont

2. Opérations avec affectation (avec ou sans affichage) • vecteurs, matrices

x=2+4 • data.frames, listes
x
6
( x=2+4) # a v e c a f f i c h a g e du r é s u l t a t Quelques fonctions génériques :
6
• ls() retourne la liste des objets de la session.
3. Les principaux types sont
• rm(a) supprime l’objet a
• entier , réel, complexe
• caractère
• logique : TRUE, FALSE, NA (not available)

10 11
Fonctions is/as Créer des vecteurs

• is.xxx(obj) teste si obj est un objet de type xxx

• as.xxx(obj) contraint si possible obj au type d’objet xxx où xxx • la fonction c( ) concatène des scalaires ou des vecteurs :
représente un type d’objet (complex, real, vector matrix etc...) > x=c ( 1 , 4 , 9 )
> y=c ( x , 2 , 3 )
> x=3 > y
> as . complex ( x ) [1] 1 4 9 2 3
> is . real (x)
[ 1 ] 3+0 i
[ 1 ] TRUE
> as . c h a r a c t e r ( x ) • Suites arithmétiques de raison 1 ou -1 : c(a:b).
> i s . complex ( x )
[ 1 ] "3"
[ 1 ] FALSE
> c (1:4) > c (4:1)
# a<b r a i s o n 1 # a>b r a i s o n 1
Remarque : [1] 1 2 3 4 [1] 4 3 2 1
Conversion de TRUE / FALSE en valeur numérique :
# a b n ’ e s t p a s un e n t i e r
> as . i n t e g e r (T) > c (1.4:7)
[1] 1 [ 1 ] 1.4 2.4 3.4 4.4 5.4 6.4
> as . i n t e g e r ( F )
[1] 0

12 13

Créer des matrices

• Généralisation : seq(a,b,t) où a est premier terme, le dernier  b
et la raison t
seq ( from , t o ) la raison est 1 Les matrices sont créées avec la fonction matrix() à partir d’un vecteur.
seq ( from , to , by= ) on f i x e l a r a i s o n
On doit fixer le nombre de colonnes ncol et/ou le nombre de lignes nrow.
seq ( from , to , l e n g t h . o u t= ) on f i x e l e nb de t e r m e s
> x = m a t r i x ( c ( 2 , 3 , 5 , 7 , 1 1 , 1 3 ) , n c o l =2)
par exemple
> seq ( 1 , 4 , by =0.1) Par défaut la matrice est remplie colonne par colonne. Pour remplir ligne
[ 1 ] 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 .... par ligne, on ajoute l’argument byrow=T
[26] 3.5 3.6 3.7 3.8 3.9 4.0
> y = m a t r i x ( c ( 2 , 3 , 5 , 7 , 1 1 , 1 3 ) , n c o l =2 , byrow=T)
> x > y
• x=rep(y ,n) pour créer un vecteur constitué de l’élément y répété
[ ,1] [ ,2] [ ,1] [ ,2]
n fois. (y peut être un scalaire ou un vecteur) par exemple [1 ,] 2 7 [1 ,] 2 3
> rep ( 1 , 4 ) [2 ,] 3 11 [2 ,] 5 7
[1] 1 1 1 1 [3 ,] 5 13 [3 ,] 11 13

14 15
 Quelques matrices particulières : diagonale, Toeplitz

Attention : si la dimension du vecteur n’est pas égale au produit (ncol

> diag ( 1 : 4 )
⇥ nrow) alors l’opération effectuée est la suivante : [ ,1] [ ,2] [ ,3] [ ,4]
> toeplitz (1:4)
[ ,1] [ ,2] [ ,3] [ ,4]
> m a t r i x ( c ( 1 : 3 ) , n c o l =2 , nrow=3) [1 ,] 1 0 0 0
[1 ,] 1 2 3 4
[ ,1] [ ,2] [2 ,] 0 2 0 0
[2 ,] 2 1 2 3
[1 ,] 1 1 [3 ,] 0 0 3 0
[3 ,] 3 2 1 2
[2 ,] 2 2 [4 ,] 0 0 0 4
[4 ,] 4 3 2 1
[3 ,] 3 3 >

> m a t r i x ( c ( 1 : 3 ) , n c o l =2) diag

[ ,1] [ ,2]
La fonction diag retourne une matrice diagonale lorsque le paramètre
[1 ,] 1 3
[2 ,] 2 1 d’entrée est un vecteur.
Si le paramètre d’entrée est une matrice, alors elle retourne un vecteur
constitué de la diagonale de la matrice

16 17

Concaténer des vecteurs/matrices Extraire des éléments d’un vecteur ou d’une matrice
> mat=m a t r i x ( v e c t , n c o l =3 , nrow=3)
> v e c t=c ( 1 . 5 : 9 . 5 )
[ ,1] [ ,2] [ ,3]
> vect
• rbind [1 ,] 1.5 4.5 7.5
[ 1 ] 1.5 2.5 3.5 4.5 5.5
[2 ,] 2.5 5.5 8.5
6.5 7.5 8.5 9.5
[3 ,] 3.5 6.5 9.5
• cbind
Colonne/ligne d’une matrice
Extraire un élément
> x =1:10
> mat [ , 1 ] > mat [ 3 , ]
> y=x ^2 > vect [ 1 ] > mat [ 2 , 1 ]
[ 1 ] 1.5 2.5 3.5
> rbind (x , y ) [ 1 ] 1.5 [ 1 ] 2.5
[ 1 ] 3.5 6.5 9.5
[ ,1] [ ,2] [ ,3] [ ,4] [ ,5] [ ,6] [ ,7] [ ,8] [ ,9] [ ,10]
x 1 2 3 4 5 6 7 8 9 10
y 1 4 9 16 25 36 49 64 81 100 Extraire un bloc ou plusieurs coordonnées
> cbind ( x , y ) > mat [ 2 : 3 ,1:2] > vect [ c (1 ,3 ,7)]
x y [ ,1] [ ,2] [ 1 ] 1.5 3.5 7.5
[1 ,] 1 1 [1 ,] 2.5 5.5
[2 ,] 2 4 [2 ,] 3.5 6.5
[3 ,] 3 9
[4 ,] 4 16 Attention : vect[-j] retourne le vecteur vect sans la j ème coordonnée
[5 ,] 5 25
[6 ,] 6 36 > vect[ c (1 ,3 ,7)] retourne 2.5 4.5 5.5 6.5 8.5 9.5
etc
18 19
Opérations sur les Matrices/Vecteurs

• Les opérations + * - / entre 2 vecteurs ou matrices de même Attention

dimension sont des opérations terme à terme.
Si les vecteurs ne sont pas de même longueur,
> A le plus court est complété automatiquement .
[ ,1] [ ,2]
[1 ,] 2 1 > x =c ( 1 : 5 )
[2 ,] 4 9 > x x : 1 2 3 4 5
> x
> B [1] 1 2 3 4 5 y : 1 2 1 2 1
[1] 1 2 3 4 5
[ ,1] [ ,2] > y =c ( 1 , 2 )
> y
[1 ,] 0 2 > y x+y : 2 4 4 6 6
[1] 0 0 0 1 1
[2 ,] 1 1 [1] 1 2
> A⇤B > x + y
> x⇤y
[ ,1] [ ,2] [1] 2 4 4 6 6
[1] 0 0 0 4 5
[1 ,] 0 2
[2 ,] 4 9
>

20 21

Action d’une fonction sur un vecteur ou une matrice

Soit FUN une fonction définie sur les scalaires qui retourne un scalaire.
Quelques opérations particulières sur les matrices
Par exemple
> a=m a t r i x ( 1 , n c o l =2 , nrow=2)
> a sqrt square root
[ ,1] [ ,2] abs absolute value
[1 ,] 1 1 sin cos tan trigonometric functions ( radians )
[2 ,] 1 1 exp log e x p o n e n t i a l and n a t u r a l l o g a r i t h m
> a+3 #m a t r i c e + s c a l a i r e log10 common l o g a r i t h m
[ ,1] [ ,2] gamma lgamma gamma f u n c t i o n and i t s n a t u r a l l o g
[1 ,] 4 4
[2 ,] 4 4 Lorsque le paramètre d’entrée est un vecteur (respectivement une
> a+c ( 1 : 2 ) #m a t r i c e + v e c t e u r matrice), la fonction FUN est appliquée sur chacune des composantes.
[ ,1] [ ,2]
L’objet retourné est un vecteur (respectivement une matrice).
[1 ,] 2 2
[2 ,] 3 3 s
Exemple
Si A = (ai,j ) est une matrice, alors exp(A) retourne une matrice
constituée des éléments e ai,j .
22 23
Quelques fonctions sur les matrices Objets booléens et instructions logiques

• Le produit matriciel est obtenu avec % ⇤ % • Les opérations logiques : < , > , <= , >= , != [différent],
• Calcul des valeurs/vecteurs propres :eigen == [égal] retournent TRUE ou FALSE.
• Calcul du déterminant : det • La comparaison entre deux vecteurs est une comparaison terme à
• t(A) retourne la transposée de la matrice A terme.
• décomposition de Choleski :chol (X) retourne R telle que • Si les vecteurs ne sont pas de même longueur, le plus court est
X = R 0 R où R est une matrice triangulaire supérieure et R 0 est la complété automatiquement.
transposée de R.
> a= 1 : 5 ; b =2.5
• décomposition svd : svd(X) retourne (U,D,V) telles que X = UDV 0 > a<b
où U et V sont orthogonales et D est diagonale. [ 1 ] TRUE TRUE FALSE FALSE FALSE

solve Il est possible de définir plusieurs conditions à remplir avec les opérateurs

• solve(A) retourne l’inverse de la matrice A • ET : &

• solve(A,b) retourne x tel que Ax = b • OU : |

24 25

Effet de la précision sur la comparaison de réels

Pour extraire les éléments d’un vecteur vect, on peut utiliser des
instructions logiques.

• Soit I un vecteur de booléens de même longueur que vect :

vect[I] retourne les coordonnées vect[j] telles que p 2
I [j] = TRUE . Est ce que 2 = 2?

Applications > ( s q r t ( 2 ) ^ 2 == 2 )
[ 1 ] FALSE
• extraire les composantes >8
vect[vect>8] : vect>8 est un vecteur de TRUE et FALSE, on Une solution
extrait les composantes affectées à TRUE.
> a l l . equal ( sqrt (2)^2 ,2)
• extraire les composantes >8 ou <2 [ 1 ] TRUE
vect[(vect>8) | (vect<2)] #ou
>isTRUE ( a l l . e q u a l ( s q r t ( 2 ) ^ 2 , 2 ) )
• extraire les composantes >8 et <10
vect[(vect>8) & (vect<10)]

26 27
Fonction which
Ces fonctions retournent un scalaire :
P Q
• sum() (somme i xi ), prod() (produit i xi ), mean() (moyenne
Soit vec un vecteur logique. La fonction which(vec) retourne les 1
P n
i=1 xi )
indices des coordonnées du vecteur vec qui prennent la valeur TRUE n
• max(), min()
> x =(1:10)^2
> x • length() (longueur du vecteur),
[1] 1 4 9 16 25 36 49 64 81 100
• dim(), ncol(), nrow() (dimension de la matrice/nombre de
> which ( x== 2 5 )
[1] 5 lignes / nombre de colonnes.)
> which ( x > 2 1 )
[1] 5 6 7 8 9 10 Ces fonctions retournent un vecteur :
Pn
• cumsum() (sommes cumulées (x1 , x1 + x2 , . . . , i=1 xi ), cumprod()
Exemple
(produits cumulés),
Les commandes x[x>1] et x[which(x>1)] retournent le même vecteur.
• sort (tri), order, unique
Cas particulier remarque : sort(x) = x[order(x)]
• which.max(x) retourne which(x==max(x)) • fft() (transformé de Fourier)

28 29

Définition des data.frames Opérations sur les dataframes

1. Pour visualiser les premières lignes head()

C’est une matrice dont toutes les colonnes ne sont pas nécessairement du 2. Pour definir ou visualiser le nom des lignes row.names
même type : scalaire, booléen, caractère. Par exemple 3. Pour definir ou visualiser le nom des colonnes names
> d a t a 1= data . frame ( x1 =1 , x2 = 1 : 1 0 , a= l e t t e r s [ 1 : 1 0 ] ) 4. La dimension de l’objet est donnée par dim
x1 x2 a
> names ( d a t a 1 )
1 1 1 a
[ 1 ] " x1 " " x2 " " a "
2 1 2 b
> names ( d a t a 1 )< c ( " c1 " , " c2 " , " c3 " )
3 1 3 c
> head ( d a ta 1 , 3 )
4 1 4 d
c1 c2 c3
5 1 5 e
1 1 1 a
6 1 6 f
2 1 2 b
7 1 7 g
3 1 3 c
8 1 8 h
> dim ( d a t a 1 )
9 1 9 i
[ 1 ] 10 3
10 1 10 j
>row . names ( d a t a 1 ) < l e t t e r s [ 1 : 1 0 ]
#l e v e c t e u r s l e t t e r s c o n t i e n t l e s l e t t r e s de l ’ a l p h a b e t
Par défaut les lignes sont numérotées 1,2 etc. > head ( d a ta 1 , 2 )
c1 c2 c3
a 1 1 a
30 31
b 1 2 b
Opérations sur les dataframes Opérations sur les dataframes

Les opérations entre des dataframes sont opérations terme à terme

Pour extraire un élément ou un bloc, la syntaxe est la même que pour les
comme pour les matrices.
matrices.
A = data . frame ( x = 1 : 3 , y =2:4)
B = data . frame ( xx =1 , yy =1:3) Pour extraire une colonne les deux syntaxes suivantes peuvent être
C= data . frame ( x = 1 : 3 , y=r e p ( " a " , 3 ) ) utilisées
> A$ x
> A > B > A+B > C [1] 1 2 3
x y xx yy x y x y > A[ ,1]
1 1 2 1 1 1 1 2 3 1 1 a [1] 1 2 3
2 2 3 2 1 2 2 3 5 2 2 a
3 3 4 3 1 3 3 4 7 3 3 a Pour concaténer des dataframes ayant le même nombre de lignes
data . frame (A , B)
> A+C
Warning m e s s a g e : x y xx yy
x y
I n Ops . f a c t o r ( l e f t , r i g h t ) : 1 1 2 1 1
1 2 NA
c e c i n ’ e s t par p e r t i n e n t pour des 2 2 3 1 2
2 4 NA
variables facteurs 3 3 4 1 3
3 6 NA

32 33

Définition d’une liste Opérations sur les listes

C’est une structure qui regroupe des objets (pas nécessairement de même
• Pour visualiser la liste des composantes d’une liste
type). On crée les listes avec la fonction list
>names ( r d n )
Exemple [ 1 ] " s e r i e " " t a i l l e " " type "
On construit une liste appelée rnd qui contient 3 objets : > summary ( r d n )
L e n g t h C l a s s Mode
• un vecteur dans serie serie 100 none numeric
taille 1 none numeric
• un scalaire dans taille type 1 none character
• une chaîne de caractères dans type
• Pour atteindre les composantes d’une liste
La syntaxe est la suivante >r d n $ t a i l l e OU >r n d [ [ 2 ] ]
>r d n= l i s t ( s e r i e =c ( 1 : 1 0 0 ) , t a i l l e =100 , t y p e=" a r i t h m " ) [ 1 ] 100 [ 1 ] 100

attention
attention
Si la liste a été créée sans spécifier de noms aux variables, il n’y a pas de
Une liste peut être créée sans donner des noms aux variables c’est a dire
nom par défaut et la seule la première syntaxe est utilisable.
rdn=list(c(1:100),100,"arithm") .
34 35
 Importer/exporter des données

• Pour extraire les objets d’une liste 1. Importer une suite : x=scan("data.dat") : pour créer un vecteur
>a t t a c h ( r d n ) à partir de données stockées dans un fichier, ici data.dat.
" s e r i e " " t a i l l e " " type "
2. Importer un tableau :
• supprimer les objets créés à la ligne précédente : x=read.table("data.dat")
x=read.table("data.dat", header=TRUE)
>detach ( r d n )
L’instruction header=TRUE permet de préciser que la première ligne
du fichier contient le nom des colonnes du tableau.
3. Exporter : write, write.table

36 37

Importer/exporter des objets R

1. Sauvegarder des objets R dans un fichier

> x =1:10
> y= l i s t ( a=1 , b=TRUE , c=" e x e m p l e " )
> save ( x , y , f i l e=’ sav . rda ’ )

2. Lire un fichier qui contient des objets R

> load ( " sav . rda " )
> x Les fonctions
[1] 1 2 3 4 5 6 7 8 9 10

Attention si un objet R appelé x (ou y) existait avant l’appel de la

fonction load, il a été remplacé par celui contenu dans le fichier
sav.rda
3. saveRDS peut aussi être utiliser si on sauvegarde une unique liste. Le
fichier est lu avec readRDS. On peut changer le nom de la liste à la
lecture du fichier
> b = readRDS ( ’ s a v . r d s ’ )
> b
> a= l i s t ( x =1 , y =3) $x 38
> saveRDS ( a , ’ s a v . r d s ’ ) [1] 1
$y
[1] 3
Structure générale pour créer des fonctions Exemple

La fonction suivante retourne le résultat de n lancers d’une pièce.

PF = f u n c t i o n ( n , p r o b a . p i l e )
{
• La structure générale d’une fonction est
#nb a l é a t o i r e s u i v a n t U n i f ( 0 , 1 )
>FUN=f u n c t i o n ( l i s t e _d e s_p a r a m è t r e s ) u=r u n i f ( n )
{ p f =(u<p r o b a . p i l e )
commandes p f = as . i n t e g e r ( p f )
r e t u r n ( o b j e t s_r e t o u r n é s ) return ( pf )
} }

• Les accolades { et } définissent le début et la fin de la fonction. La sortie de la fonction est

• La dernière instruction return contient le ou les objets retournés PF ( 1 0 , 1 / 2 )
par la fonction. [1] 0 0 1 0 1 0 1 1 0 0
# a v e c a f f e c t a t i o n de l a s o r t i e d a n s l e v e c t e u r x
• Exécuter la fonction : FUN(...) x = PF ( 1 0 , 1 / 2 )
x
[1] 0 1 1 1 0 0 0 0 0 0

39 40

Renvois multi-arguments Paramètres par défaut

La fonction return interdit les sorties avec plusieurs arguments : il faut On peut affecter des valeurs par défaut aux paramètres d’entrée d’une
les regrouper dans un seul objet sous la forme d’une liste. fonction.
PF = f u n c t i o n ( n , p r o b a . p i l e )
Modification de la fonction PF
{
u=r u n i f ( n ) #nb a l é a t o i r e s u i v a n t U n i f ( 0 , 1 ) par défaut on suppose que la pièce est équilibrée.
p f =(u<p r o b a . p i l e )
PF = f u n c t i o n ( n , p r o b a . p i l e =1/ 2 )
p f = as . i n t e g e r ( p f )
{
f =mean ( p f )
u=r u n i f ( n ) #nb a l é a t o i r e s u i v a n t U n i f ( 0 , 1 )
r e t u r n ( l i s t ( e c h a n t i l l o n = pf , f r e q u e n c e = f ) )
p f =(u<p r o b a . p i l e )
}
p f = as . i n t e g e r ( p f )
return ( pf )
Exécution de la fonction : }

PF ( 4 , 1 / 2 ) l= PF ( 4 , 1 / 2 )
$echantillon l$echantillon Les commandes
[1] 0 0 1 0 [1] 0 0 1 1 PF ( 1 0 ) OU PF ( 1 0 , 1 / 2 )
$frequence l$f
[ 1 ] 0.25 [ 1 ] 0.5
retournent le même résultat
41 42
Paramètres d’entrée Remarque sur les valeurs par défaut

Il y a trois façons de spécifier les paramètres d’entrée d’une fonction Modification de la fonction PF :
on inverse l’ordre des paramètres d’entrée :
• par la position : les paramètres d’entrée sont affectés aux premiers
PF = f u n c t i o n ( p r o b a . p i l e =1/ 2 , n )
arguments de la fonction. PF(3,1/2) : les paramètres d’entrée sont {
n=3 et proba.pile=1/2 ...
}
• par le nom : il s’agit du moyen le plus sûr, les noms des arguments
sont précisés de manière explicite. On peut alors écrire PF ( 1 0 )
PF(proba.pile=1/2 ,n=3), l’ordre n’est plus prioritaire E r r e u r d a n s . I n t e r n a l ( r u n i f ( n , min , max ) ) : ’ n ’ e s t manquant
PF ( n=10)
• avec des valeurs par défaut : ces valeurs par défaut seront utilisées si
[1] 0 0 1 0 1 0 1 0 1 0
les paramètres d’entrée ne sont pas spécifiés. On peut alors écrire
PF(3) ou PF(n=3) : les paramètres d’entrée sont n=3 et la valeur
IL est donc préférable de placer les paramètres sans valeur par défaut en
par défaut pour proba.pile c’est à dire 1/2.
premier dans la déclaration des variables d’entrée.

43 44

Pour les fonctions de deux variables la fonction Vectorize

• f une fonction de deux variables f : (x, y ) 7! f (x, y ) Soit f une fonction dont le paramètre d’entrée x est un scalaire.
• x et y deux vecteurs de même dimension. Vectorize transforme la fonction f en une fonction vectorielle c’est à
La commande f(x,y) retourne le vecteur constitué des éléments f (xi , yi ) dire une fonction qui évalue la fonction f en chaque point d’un vecteur
. d’entrée.

Si x et y ne sont pas de même dimension, celui de plus petite dimension Soit x = (x1 , ..., xn ), on veut évaluer f aux points xi .
est répété. # on t r a n s f o r m e l a f o n c t i o n { f } en une f o n c t i o n v e c t o r i e l l e df .
> df = Vectorize ( f , ’ x ’ )
Tableaux croisés >y=d f ( x )

La fonction outer retourne une matrice de la forme Autres programmations possibles

M(i, j) = fun(xi , yj ) 1. avec une boucle for

> y = rep (0 , n )
x =1:5 > for ( i in 1:n) y [ i ] = f (x [ i ])
y =1:5
M=o u t e r ( x , y , ’ f u n ’ ) 2. avec la fonction sapply
> y = sapply (x , ’f ’)
fun peut aussi être une opération élémentaire +/-*
45 46
Illustration sur des fonctions de 2 variables suite
> f=f u n c t i o n ( x , y ) s i n ( x+y ^2)
> f (1 ,1)
[ 1 ] 0.9092974
> # > f (1 , y )
> x = 1:3 [1] 0.9092974 0.9589243 0.5440211
> y= 1 : 3 > # identique à
> # on c a l c u l e f ( x [ i ] , y [ i ] ) > f ( rep ( 3 , 1 ) , y )
> f (x , y) [1] 0.7568025 0 . 6 5 6 9 8 6 6 0.5365729
[1] 0.9092974 0.2794155 0.5365729 >
> z =1:2 > #c a l c u l du t a b l e a u c r o i s é f (x[ i ] , y[ j ])
> f (x , z) > df = Vectorize ( f , ’ x ’ )
[1] 0.9092974 0.2794155 0.7568025 > df ( x , z )
Message d ’ a v i s : [ ,1] [ ,2] [ ,3]
I n x + y ^2 : [1 ,] 0.9092974 0.1411200 0.7568025
l a t a i l l e d ’ un o b j e t p l u s l o n g n ’ e s t p a s m u l t i p l e de l a t a i l l e d ’ un o b j e t p l u s c o[ 2u r, t] 0.9589243 0.2794155 0.6569866
> # le calcul effectué est > # les vecteurs colonnes sont f (1 ,1:2) f (2 ,1:2) f (3 ,1:2)
> #f ( x [ 1 ] , z [ 1 ] ) f ( x [ 2 ] , z [ 2 ] ) f ( x [ 3 ] , z [ 1 ] )
> #i d e n t i q u e à
> f (x , c(z , z [1]))
[1] 0.9092974 0.2794155 0.7568025
47 48

Les fonctions usuelles plot(), lines(), points()

• plot est la fonction centrale

• Le fonctions points ou lines sont utilisées pour superposer des
courbes ou des nuages de points.
Premier exemple : représenter des vecteurs plot(y) / plot(x,y)
x= seq ( 4 ,4 ,.1) y=l o g ( x^2+1/ x ^2)

Les graphiques plot(y)

1 2 3 4
y
0 20 40 60 80

Index

plot(x,y,pch=3)
1 2 3 4
y

-4 -2 0 2 4

49
Quelques arguments de la fonction plot Graphique en 3D

type=p
x < seq ( 10 , 1 0 , l e n g t h= 3 0 )
y < x

4
f < f u n c t i o n ( x , y ) { r < s q r t ( x^2+y ^ 2 ) ; 10 ⇤ s i n ( r ) / r }

2
0
• pour fixer les limites des axes 0 5 10 15 20 25 30 z < outer ( x , y , f )
p e r s p ( x , y , z , t h e t a = 3 0 , p h i = 3 0 , expand = 0 . 5 , c o l = " l i g h t b l u e " )
Index

• ylim=c(ay,by) et xlim=c(ax,bx) type=l

• par défaut les bornes sont optimisées

4
sur la première courbe tracée

2
0
0 5 10 15 20 25 30

• type="p" (points) ou "l" Index

(ligne) : pour tracer une ligne ou un type=s

nuage de points.

4
x

2
0
• pch : type de points 0 5 10 15 20 25 30

z
Index

• lty : type de lignes.

y
type=h
x

• col : couleur

4
x

2
0
0 5 10 15 20 25 30

Index

50 51

Représentation graphique d’une matrice Superposition de courbes

Maunga Whau Volcano
600

superposer des courbes

500

commentaire

5
x=rnorm ( 2 0 )
400

y=r e x p ( 2 0 )
#nuage de p o i n t s
300
y

4
plot (x , y)
200

#a j o u t e r un nuage de p o i n t s

3
100

p o i n t s ( x + . 1 , y + . 1 , pch =2)

y
100 200 300 400 500 600 700 800
#a j o u t e r une l i g n e
l i n e s ( s o r t ( x ) , y , l t y =2)

2
x

#a j o u t e r une l i g n e h o r i z o n t a l e
a b l i n e ( h=3)

1
#t e x t e + f r a m e t i t l e
z

text (1 ,5 , " commentaire " )

t i t l e ( " superposer des courbes " )

0
-1 0 1 2
y

x x

52 53
légende legend Autour de la fonction plot
legend ( 1 , 1 . 9 , c ( " x^2 1" , " s i n " , " c o s " ) , c o l = c ( 3 , 4 , 6 ) ,
l t y = c (2 , 1, 1 ) , pch = c ( 1, 3 , 4))

Le graphique produit par la fonction plot(x) dépend de la classe de

x^2-1 l’objet x.
sin les emplacements
cos
methods ( p l o t )
1.5

prédéfinis :
[ 1 ] " p l o t . data . frame " " plot . default "
[3] " plot . density " " plot . factor"
1.0

top

1.0
topright, inset = .02
topleft, inset = .05
(x,y)
(x,y)
(x,y) [ 5 ] " plot . formula " " plot . function "
[ 7 ] " plot . histogram " " plot . lm "
0.5

0.5
x^2 - 1

[ 9 ] " p l o t . mlm" " plot . mts "

left center right [ 1 1 ] " p l o t . new" " plot . POSIXct "
0.0

0.0
y
(x,y) (x,y) (x,y)

[ 1 3 ] " p l o t . POSIXlt " " plot . table "

[15] " plot . ts " " plot . window "

-0.5
-0.5

[ 1 7 ] " p l o t . xy "
bottomleft bottom bottomright

-1.0
(x,y) (x,y) (x,y)
-1.0

0.0 0.2 0.4 0.6 0.8 1.0

-3 -2 -1 0 1 2 3

x 54 55

Illustrations Représentation graphique d’une matrice ou dataframe

• sur une fonction (par ex sin) Les outils graphiques matplot et pairs sont adaptés aux matrices dont
p l o t ( s i n , x l i m=c( p i , p i ) ) les colonnes correspondent à des variables.
1.0

Exemple
300

• sur un tableau
0.5
sin (x)

200

>data ( i r i s )
0.0

x=r p o i s ( 1 0 0 0 , 1 )
>i r i s
100
-0.5

y=t a b l e ( x )
-1.0

y
0

-3 -2 -1 0 1 2 3 0 1 2 3 4 5 Sepal . Length Sepal . Height P e t a l . Length P e t a l . Height Species

x 0 1 2 3 4 6
x x 1 5.1 3.5 1.4 0.2 setosa
374 372 162 71 20 1
2 4.9 3.0 1.4 0.2 setosa
plot (y) Histogram of r density.default(x = r)
3 4.7 3.2 1.3 0.2 setosa
0.4

4
• sur un histogramme ou une
20

...
0.3
15
Frequency

densité
Density

54 5.5 2.3 4.0 1.3 v e r s i c o l o r

0.2
10

r=rnorm ( 1 0 0 , 1 ) 55 6.5 2.8 4.6 1.5 v e r s i c o l o r

0.1
5

z=h i s t ( r , p l o t=F ) ...

0.0
0

plot ( z ) -1 0 1 2 3 -2 -1 0 1 2 3 4 147 6.3 2.5 5.0 1.9 virginica

w=d e n s i t y ( r ) r N = 100 Bandwidth = 0.3564 148 6.5 3.0 5.2 2.0 virginica
p l o t (w)
56 57
Fonction matplot(matrice ...) Utilisation des biplots , pairs(matrice , ....)

Cette fonction représente sur un même graphique les colonnes d’une Cette fonction représente tous les nuages de points possibles entre les
matrice ou d’une data.frame. différentes colonnes.

Anderson’s Iris Data – 3 species

Anderson's Iris Data -- 3 species
2.0 2.5 3.0 3.5 4.0 0.5 1.0 1.5 2.0 2.5

4.5 5.5 6.5 7.5

8

Sepal.Length

2.0 2.5 3.0 3.5 4.0

6

Sepal.Width

1 2 3 4 5 6 7
iris

Petal.Length

0.5 1.0 1.5 2.0 2.5

2

Petal.Width

4.5 5.5 6.5 7.5 1 2 3 4 5 6 7

0

58 59
0 50 100 150

Construction d’un polygone autre exemple

1. On commence par fixer les axes des abscisses et des ordonnés à Distance Between Brownian Motions

-5 0 5 10 15 20
l’aide d’un graphique vide.

Distance
p l o t ( c ( 1 , 9 ) , c ( 0 , 2 ) , t y p e="n" )

2. Avec la fonction polygon, on trace le polygone défini pas ses

sommets
0 20 40 60 80 100

Time

polygon ( c ( 4 . 5 , 5 . 5 , 5 . 5 , 4 . 5 ) ,c (0 ,0 ,1 ,1))
n=100
3. Arguments supplémentaires z . h a u t=c ( 0 , cumsum ( rnorm ( n ) ) )
z . b a s= c ( 0 , cumsum ( rnorm ( n ) ) )
p o l y g o n ( 1 : 9 , c ( 2 , 1 , 2 , 1 ,NA, 2 , 1 , 2 , 1 ) , c o l=c ( " r e d " , " b l u e " ) ,
b o r d e r=c ( " g r e e n " , " y e l l o w " , lwd =3 , l t y =c ( " d a s h e d " , " s o l i d " ) )
xx < c ( 0 : n , n : 0 )
yy < c ( z . bas , r e v ( z . h a u t ) )
etape 1 etape 2 etape 3
#g r a p h i q u e v i d e p o u r f i x e r l e s d i m e n s i o n s
2.0

2.0

2.0

plot ( xx , yy , t y p e="n" , x l a b=" Time " , y l a b=" D i s t a n c e " )

1.5

1.5

1.5

#t r a c e r l e p o l y g o n e
c(0, 2)

c(0, 2)

c(0, 2)
1.0

1.0

1.0

p o l y g o n ( xx , yy , c o l=" g r a y " , b o r d e r = " r e d " )

0.5

0.5

0.5

t i t l e ( " D i s t a n c e Between B r o w n i a n M o t i o n s " )

0.0

0.0

0.0

2 4 6 8 2 4 6 8 2 4 6 8

c(1, 9) c(1, 9) c(1, 9)

60 61
 • Pour sauvegarder un graphique :
• utilisation de la fonction dev.print.
• pour obtenir un fichier postscript :
dev.print(postscript, file="essai.ps")
• pour obtenir un fichier pdf :
dev.print(pdf, file="essai.pdf")

• utilisation des menus (sous widows ou mac seulement)

Graphiques avec ggplot2
• La fenêtre graphique peut être fractionnée en utilisant
• par(mfrow=c(n,m)), on obtient alors n ⇥ m graphiques sur une
même page organisés sur n lignes et m colonnes
• split.screen(m,n)
• screen(i), screen(i,FALSE) pour sélectionner la sous fenêtre
• erase.screen()
• close.screen(all = TRUE)

62

Description La fonction qplot ou quickplot

La librairie ggplot2 fournit un outil de visualisation des données à l’aide x= seq ( 4 , 4 , . 1 )

y=l o g ( x^2+1/ x ^2)
de jolis graphiques qplot (x , y)
q p l o t ( x , y geom=" l i n e " )
1. La fonction pour produire rapidement des graphiques simples est q p l o t ( x , y , geom=" l i n e " , c o l o u r=I ( " g r e e n " ) )
qplot (elle remplace la fonction plot)
2. Pour des graphiques plus complexes la fonction a utilisée est ggplot
• Toutes les données doievent être au format data.frame
• les variables sont organisées en colonne.

Référence :

Représentation par des points (par défaut si deux vecteurs en entrée), par
une ligne et par une ligne en modifiant la couleur noire par défaut.
63 64
Option de la fonction qplot hitogramme avec qplot

On conserve les options "xlab", "ylab", "xlim" et "ylim". S’il y a un seul vecteur en entrée le graphique par défaut est un
histogramme. On controle la largeur des classes avec binwidth. Par
Pour changer la couleur ou le type de point
défaut, l’histogramme est dessiné avec 30 classes.
y = rnorm ( 2 0 )
x= 1 : 2 0 x=rnorm ( 1 0 0 )
q p l o t ( x , y , s h a p e = I ( 1 : 2 0 ) , c o l o u r=I ( rainbow ( 2 0 ) ) q p l o t ( x , c o l o u r=I ( " p u r p l e " ) , f i l l =I ( " p i n k 3 " ) , b i n w i d t h =.05)
qplot (x , y , c o l o u r=I ( rainbow ( 2 0 ) ) , s i z e=I ( ( 1 : 2 0 ) / 2 ) ) q p l o t ( x , c o l o u r=I ( " p u r p l e " ) , f i l l =I ( " p i n k " ) , b i n w i d t h =.5)
q p l o t x , c o l o u r=I ( " p u r p l e " ) , f i l l =I ( " p i n k 1 " ) , b i n w i d t h =1)

Estimateur à noyau de la densité

q p l o t ( y , geom = c ( " d e n s i t y " ) , c o l=I ( " c y a n 4 " ) )

65 66

Couplage avec une fonction statistique de lissage La fonction ggplot

• La fonction ggplot est utilisée pour construire l’objet plot initial.

1. Le lissage : par défaut loes • elle est presque toujours suivi de + pour ajouter un composant au
tracé.
q p l o t ( x , y , geom = c ( " p o i n t " , " smooth " ) )
• La fonction aes liste les variables utilisées dans la construction du
2. regression linéaire ou polynomiale graphique x et y ainsi que les variables qui contrôlent l’esthétique du
q p l o t ( x , y , geom = c ( " p o i n t " , " smooth " ) , method=" lm " )
graphique colour (col dans plot) / shape (pch dans plot) / size (cex
q p l o t ( x , y , geom = c ( " p o i n t " , " smooth " ) , method=" lm " , dans plot)
f o r m u l a=y~p o l y ( x , 2 ) )
Il existe differentes façons d’appeler la fonction ggplot
1. si toutes les calques utilisent les mêmes données dans l’objet
data.frame df et les mêmes arguments de style. graphique
g g p l o t ( df , a e s ( x , y , <a u t r e s v a r i a b l e s >))

2. si toutes les tracés utilisent les mêmes données

ggplot ( df )

3. si les calculs utilisent différents jeux de données

67 ggplot () 68
ggplot sur les données data(iris) Autre librairie : (GGally)

• Equivalent du graphique avec matplot

data ( i r i s )
Une version sophistiquée de pairs
i r i s =data . frame ( x = 1 : 1 5 0 , i r i s ) ggpairs ( i r i s , aes ( colour = Species , alpha = 0.4))
ggplot ( i r i s ) +
geom_l i n e ( a e s ( x , S e p a l . L e n g t h ) , c o l o u r = I (1)) +
geom_l i n e ( a e s ( x , S e p a l . Width ) , c o l o u r = I (2)) +
geom_l i n e ( a e s ( x , P e t a l . L e n g t h ) , c o l o u r = I (3)) +
geom_l i n e ( a e s ( x , P e t a l . Width ) , c o l o u r = I (4)) +
ylab (" IRIS ") +
x l a b ( "" )

• la couleur dépend de la variable Species

ggplot ( i r i s ) +
geom_l i n e ( a e s ( x , S e p a l . Length , c o l o u r = S p e c i e s ) )

69 70

Instructions conditionnelles

La syntaxe
if (condition) {instructions} permet de calculer les instructions
uniquement si la condition est vraie.
Structures de contrôle et if (condition) { A } else {B} calcule les instructions A si la
Itérations condition est vraie et les instructions B sinon.
Par exemple,
if (x>0) y=x*log(x) else y=0
Remarque : Si les instructions se limitent à un seul calcul comme dans
cet exemple on peut utiliser la fonction ifelse
y=ifelse(x>0,x*log(x),0)

71
Opération non vectorielle Itérations

1. Avec la commande if (condition) ... else ... la condition

On utilise les boucles pour exécuter plusieurs fois une instruction ou un
ne peut pas être vectorielle
bloc d’instructions
Par exemple
> x= 1 : 7 Les trois types de boucle sont
> i f ( x >2) p r i n t ( "A" ) e l s e p r i n t ( "B" )
[ 1 ] "B" • for :

Dans la condition, x a une longueur supérieure à 1. Seul le premier for(var in seq) {commandes}
élément est utilisé : • while :
(x>2) correspond à (x[1] >2) while(cond) { commandes}
2. ifelse permet d’appliquer une instruction conditionnelle sur • repeat :
chacune des coordonnées d’un vecteur.
repeat { commandes ; if (cond) break }
> i f e l s e ( x <2 , "A" , "B" )
[ 1 ] "A" "B" "B" "B" "B" "B" "B"
1. Dans une boucle for, le nombre d’itérations est fixé.
Alternative :
2. La durée d’execution des boucles while/repeat peut être infinie !
> g = f u n c t i o n ( x ) i f ( x >2) p r i n t ( "A" ) e l s e p r i n t ( "B" )
>vg = V e c t o r i z e ( g , ’ x ’ )
> vg ( 1 : 1 0 ) 72 73

Exemple Exemple

On veut simuler n variables aléatoires suivant la loi de X1 + ... + Xp où

Les 3 programmations suivantes retournent le même résultat. les Xi sont iid suivant la loi uniforme sur [0, 1].
1. Avec l’instruction for On stocke l’échantillon dans Y .
for ( i in 1:10) # i n i t i a l i s a t i o n de Y a v e c d e s 0
{ 3. Avec l’instruction repeat Y = rep (0 , n )
commandes for ( i in 0:n )
} i= 1 { Y [ i ] = sum ( r u n i f ( p , 0 , 1 ) ) }
repeat Une autre programmation :
2. Avec l’instruction while. { Une autre façon de remplir Y # i n i t i a l i s a t i o n de Y a v e c d e s 0
commandes par concaténation : Y = rep (0 , n )
i= 1 i = i +1 for ( i in 1:p)
w h i l e ( i <= 10 ) i f ( i >10) b r e a k # Y e s t un v e c t e u r v i d e
{ Y = Y + runif (n ,0 ,1) }
{ } Y = NULL
commandes for ( i in 0:n )
i = i +1 {
} Y = c (Y , sum ( r u n i f ( p , 0 , 1 ) ) )
}

74 75
Exemple Exemple

On dispose des températures moyennes mensuelles relevées à Nottingham

On simule des variables aléatoires suivant la loi de Bernoulli B(1/2) pendant 10 ans de 1920 à 1939.
jusqu’à l’obtention du premier 1. > nottem
Jan Feb Mar Apr May Jun J u l Aug Sep Oct Nov
Le nombre de variables (noté N dans le code ci dessous) simulées suivant Dec
la loi de Bernoulli suit une loi géométrique de paramètre 1/2. 1920 4 0 . 6 4 0 . 8 4 4 . 4 4 6 . 7 5 4 . 1 5 8 . 5 5 7 . 7 5 6 . 4 5 4 . 3 5 0 . 5 4 2 . 9 3 9 . 8
1921 4 4 . 2 3 9 . 8 4 5 . 1 4 7 . 0 5 4 . 1 5 8 . 7 6 6 . 3 5 9 . 9 5 7 . 0 5 4 . 2 3 9 . 7 4 2 . 8
etc
En utilisant while En utilisant repeat > nottem=m a t r i x ( nottem , n c o l =12 , byrow=T)

x=rbinom ( 1 , 1 , . 5 ) On souhaite calculer un profil moyen annuel et le stocker dans le vecteur

N=1 N=0
w h i l e ( x != 1 ) repeat {
temp.
{ x=rbinom ( 1 , 1 , . 5 ) temp = r e p ( 0 , 1 2 )
x=rbinom ( 1 , 1 , . 5 ) N=N+1 f o r ( i i n 1 : 1 2 ) temp [ i ] = mean ( nottem [ , i ] )
N=N+1 i f ( x==1) b r e a k
} } On peut aussi initialiser le vecteur temp comme le vecteur vide, puis on le
remplit en concaténant les résultats
temp = NULL
76 for ( i in 1:12) temp = c ( temp , mean ( nottem [ , i ] ) ) 77

La fonction apply Exemple d’utilisation de la fonction apply

apply(matrice , MARGIN, FUN, ARG.COMMUN)

La fonction apply() permet d’appliquer la même fonction FUN sur

toutes les lignes (MARGIN=1) ou les colonnes (MARGIN=2) d’une matrice
Retour à l’exemple des températures
MAT
Calcul du profil annuel :
• Chaque ligne ( (MARGIN=1) ou les colonnes (MARGIN=2) ) est
> temp = a p p l y ( nottem , 2 , mean ) # moyenne s u r l e s c o l o n n e s
affectée au premier paramètre de la fonction FUN [ 1 ] 39.695 39.190 42.195 46.290 52.560 ....
• La syntaxe de FUN est FUN( x, ARG.COMMUN)
où Calcul des moyennes annuelles
• x est un vecteur > temp . a n n u e l l e = a p p l y ( nottem , 1 , mean ) # moyenne s u r l e s lignes
• ARG.COMMUN représente éventuellement des paramètres
supplémentaires qui sont communs à toutes les exécutions.

Cette fonction remplace une boucle sur le nombre de colonnes

ou de lignes

78 79
Supprimer une boucle en utilisant apply replicate

Objectif : Estimer la distribution de la médiane empirique

La fonction replicate évalue n fois une même expression
• On simule N échantillons de taille n iid suivant la loi gaussienne
standard. [utile l’expression implique des nombres aléatoires ! !].
• Pour chaque échantillon, on calcule et on stocke la valeur de la
médiane empirique dans le vecteur med. La syntaxe : replicate(n, expression)
1. avec boucle Exemples
n =100 > r e p l i c a t e (10 ,1+1)
N=500 [1] 2 2 2 2 2 2 2 2 2 2
med = 1 :N
f o r ( i i n 1 :N) med [ i ] = q u a n t i l e ( rnorm ( n , 0 , 1 ) , p r o b s =1/ 2 ) > x=rnorm ( 1 0 )
Histogram of med

2. Sans boucle > r e p l i c a t e ( 5 , mean ( x ) )

[1] 0.5073058 0.5073058 0.5073058 0.5073058 0.5073058

3
A l e a = m a t r i x ( rnorm ( n⇤N, 0 , 1 ) , n c o l= N)
med = a p p l y ( A l e a , 2 , q u a n t i l e , p r o b s =1/ 2 )

2
Density
> r e p l i c a t e ( 5 , mean ( rnorm ( 1 0 ) ) )
[1] 0.22137912 0.09330663 0.12511439 0.02613061 0.19182371

1
0
−0.2 0.0 0.2 0.4
med

80 81

Suite de l’exemple sur la médiane empirique sapply,lapply

Ces fonctions calculent la même fonction sur tous les éléments d’un
vecteur ou d’une liste.
La syntaxe :
La fonction replicate permet de programmer sans boucle le code
suivant lapply(X,FUN, ARG.COMMUN)
n =100
N=500 • La fonction lapply applique une fonction FUN à tous les éléments
med = 1 :N du vecteur ou de la liste X.
f o r ( i i n 1 :N) med [ i ] = q u a n t i l e ( rnorm ( n , 0 , 1 ) , p r o b s =1/ 2 )
• Les valeurs de X sont affectées au premier argument de la fonction
FUN.
La syntaxe est la suivante
• Si la fonction FUN a plusieurs paramètres d’entrée, ils sont spécifiés
med = r e p l i c a t e ( 5 0 0 , q u a n t i l e ( rnorm ( 1 0 0 , 0 , 1 ) , p r o b s =1/ 2 ) )
après le nom de la fonction : ARG.COMMUN
• Cette fonction retourne le résultat sous la forme de listes.
• sapply est une fonction similaire à lapply mais le résultat est
retourné sous forme de vecteurs, si possible.

82 83
Exemple d’utilisation de lapply Supprimer une boucle à l’aide de la fonction sapply

cars est une liste constituée de deux vecteurs speed et dist Il est très souvent possible de supprimer les boucles for en utilisant
On calcule la moyenne et les quantiles des deux composantes de la liste. lapply ou sapply.
> cars
speed d i s t Soit x = (x1 , ...., xt ) une série représentant l’évolution du prix d’une
1 4 2
action sur une durée de t jours.
2 4 10
3 7 4 On veut calculer Mi le prix maximum sur la période 1 à i pour
>l a p p l y ( c a r s , mean )
i = 1, ..., t.
$speed $dist
[ 1 ] 15.4 [ 1 ] 42.98
1. avec une boucle
>l a p p l y ( c a r s , q u a n t i l e , p r o b s = ( 0 : 4 ) / 4 ) M = 1: t
$speed for ( i in 1: t )
0% 25% 50% 75% 100% {
4 12 15 19 25 M[ i ] = max ( x [ 1 : i ] )
$dist }
0% 25% 50% 75% 100%
2 26 36 56 120 2. sans boucle
FUN . max = f u n c t i o n ( n , y ) max ( y [ 1 : n ] )
84 85

M = s a p p l y ( 1 : t , FUN . max , y = x )

La fonction tapply Evaluer le temps de calcul

Le type factor est un objet vectoriel qui permet de spécifier une

classification discrète (nombre fini de groupes).
La fonction tapply applique une fonction FUN sur les sous groupes d’un
vecteur X définis par une variable de type factor GRP
proc.time évalue combien de temps réel et CPU (en secondes) le
tapply(X,GRP,FUN,...) processus en cours d’exécution a déjà pris.
> note > ptm < p r o c . time ( )
[ 1 ] 10.83676 11.63757 11.07312 13.79699 10.84186 > f o r ( i i n 1 : 5 0 ) mean ( r u n i f ( 1 0 0 0 ) )
10.72562 13.58680 14.85070 13.15659 14.36744 > p r o c . time ( ) ptm
> xgr utilisateur système écoulé
[ 1 ] " a " "b" " a " "b" " a " " a " "b" " a " "b" "b" 0.016 0.001 0.027
> g r=f a c t o r ( x g r )
> gr
[1] a b a b a a b a b b
Levels : a b
> t a p p l y ( n o t e , gr , mean )
a b
11.66561 13.30908
86 87

La fonction aggregate permet aussi d’évaluer une fonction sur des sous
Comparaison de trois programmations for/replicate/apply

( n= 10^3 , M= 10^4)
user system e l a p s e d
wfor 0.680 0.051 0.760
wrep 0 . 6 1 8 0.018 0.632
#b o u c l e wapp 1 . 1 3 2 0.157 1.286
> x =1:M
> f o r ( i i n 1 :M) ( n= 10^2 , M= 10^5)
x [ i ] = mean ( r u n i f ( n ) ) user system e l a p s e d Autour des lois de probabilités
#r e p l i c a t e wfor 3.160 0.239 3.451
> x=r e p l i c a t e (M, mean ( r u n i f ( n ) ) ) wrep 2 . 5 2 2 0.066 2.611
#a p p l y wapp 3 . 0 3 1 0.253 3.456
> Z=m a t r i x ( r u n i f ( n⇤M) , n c o l=M)
> x=a p p l y ( Z , 2 , mean ) ( n= 10 , M= 10^6)
user system e l a p s e d
wfor 32.077 0.419 43.655
wrep 3 7 . 2 1 9 0.481 41.147
wapp 2 9 . 9 1 1 0.392 30.357

88

Généralités Quelques lois disponibles

Soit X une variable aléatoire de loi PX

(P 1. Lois discrètes
x2A P(X = x) loi discrète
PX (A) = P(X 2 A) = R • Loi binomiale (n,p) binom
f (x)d.x loi continue
A • Loi hypergéométrique (N,n,k) hyper
• Loi de Poisson (a) pois
Pour les lois classiques, des fonctions existent pour • Loi géométrique (p) geom
• Loi à support fini {(ai , pi ), i = 1...m} sample
• calculer 8
<P(X = x) pour les lois discrètes
2. Lois continues
• la densité • Loi Gaussienne (m, 2 ) norm
:f (x) pour les lois continues
• Loi uniforme sur [a, b] unif
• la fonction de répartition F (x) = P(X  x)
• Loi de Student à ⌫ degrés de liberté t
• les quantiles F (u) = inf{x : F (x) u} 1
• Loi du 2 à ⌫ degrés de liberté chisq
• simuler des nombres aléatoires suivant la même loi que X .

1. Si F est une bijection alors F =F 1

89 90
 Exemple : loi gaussienne de moyenne 0 et de variance 1

• dnorm(x,0,1) : densité au point x

• pnorm(x,0,1) : fonction de répartition au point x
• qnorm(↵,0,1) : quantile d’ordre ↵
Si www représente le nom d’une des lois alors
• rnorm(n,0,1) : échantillon de taille n
• dwww(x,...) calcule la densité de la loi www au point x Gaussian N(0,1) distribution
cdf
density

• pwww(x,... ) : calcule la fonction de répartition au point x

0.4

1.0
0.8
x=seq ( 5 , 5 , . 0 1 )

0.3

0.6
dnorm(x)

pnorm(x)
• qwww(↵,... ) : calcule le quantile d’ordre ↵ u=seq ( 0 , 1 , . 0 1 )

0.2

0.4
0.1

0.2
• rwww(n,...) retourne un échantillon de taille n

0.0

0.0
p l o t ( x , dnorm ( x , 0 , 1 ) , t y p e=" l " ) -4 -2 0 2 4 -4 -2 0 2 4

x x

les . . . représentent les paramètres spécifiques à chaque loi.

p l o t ( x , pnorm ( x , 0 , 1 ) , t y p e=" l " ) quantile function random values

0.0 0.5 1.0 1.5 2.0

4
p l o t ( u , qnorm ( u , 0 , 1 ) , t y p e=" l " )

rnorm(10)
qnorm(u)

0
-2
p l o t ( rnorm ( 1 0 , 0 , 1 ) )

-4
0.0 0.2 0.4 0.6 0.8 1.0 2 4 6 8 10

u Index

91 92

Simulation d’une loi discrète finie Simuler le résultat d’un dé à 6 faces

Soit Z une variable aléatoire à valeurs dans {x1 , ....., xk } telle que On simule suivant la loi uniforme sur l’ensemble {1, ..., 6}.
On simulen = 100 réalisations de façon indépendante
p(Z = xi ) = probi
# Les r é a l i s a t i o n s
pour tout i = 1...k > DE
[1] 3 2 5 3 4 5 4 4 2 4 2 5 2 6 5 6 4 3 5 4 5 1OO
3 lancers
5
• La fonction sample est un générateur de nombres aléatoires pour les [24] 3 5 4 4 5 5 4 4 4 3 5 5 2 5 3 5 4 3 1 3 1 3 4
[47] 2 5 6 2 1 5 1 6 2 1 2 3 5 1 3 4 2 3 1 2 6 1 6

20
lois discrète à support fini.
[70] 6 4 1 3 4 5 6 1 6 3 5 4 1 5 2 2 2 3 3 2 2 6 6
• La syntaxe est sample(x, n, replace= TRUE, prob)

15
[93] 6 4 2 5 5 1 5 4

table(DE)
où # c a l c u l du nombre de 1 , . . . , 6 d a n s

10
#l ’ é c h a n t i l l o n s i m u l é :
• x contient les valeurs prises par Z
> t a b l e (DE)

5
• prob contient les valeurs des probabilités.
DE
lorsque l’option prob n’est pas spécifié, par défaut c’est la loi

0
1 2 3 4 5 6 1 2 3 4 5 6

uniforme Lorsque replace=TRUE indique que le tirage est avec 12 17 17 19 23 12 DE

remise. Autrement dit on simule une réalisation de n variables # r e p r é s e n t a t i o n de l a r é p a r t i t i o n

aléatoires indépendantes et de même loi > p l o t ( t a b l e (DE ) )

93 94
Construire une fonction càdlàg constante par morceaux Exemple : loi binomiale

La fonction stepfun retourne un objet de type function qui défini une On trace la fonction de répartition de la loi binomiale B(10, 1/2)
fonction constante par morceaux et continue à droite.
points . discont = 0:10
La syntaxe est F . p o i n t s . d i s c o n t = pbinom ( p o i n t s . d i s c o n t , 1 0 , 0 . 3 )
F=stepfun(x, z)
F = stepfun ( points . discont , c (0 ,F . points . discont ) )
où p l o t ( F , v e r t i c a l =FALSE , c o l =4 , main=" f c t r e p l o i b i n o m i a l e " , y l a b="F" )

• le vecteur x contient les points de discontinuité de la fonction,

• le vecteur z est de la forme c(a, y) où a est la valeur prise par F

1.0
● ● ● ● ●

avant le premier point de discontinuité et y les valeurs de la fonction

●

0.8
aux points de discontinuités x. ●

0.6
f(x)
Pour représenter graphique la fonction F on utilise plot ou lines (pour

0.4
●

superposer) : avec l’argument vertical=FALSE

0.2
●

> plot (F , v e r t i c a l =FALSE ) ou l i n e s ( F , v e r t i c a l =FALSE )

●

0.0
0 2 4 6 8 10

x
Application : les fonctions de répartition des lois discrètes
95 96

Les données Island

islands data: area (sq. miles) Africa Antarctica Asia

11506 5500 16988
Asia Australia Axel Heiberg Baffin
Africa
North America 2968 16 184
South America
Antarctica Banks Borneo Britain
Europe
Australia 23 280 84
Greenland
New Guinea Celebes Celon Cuba
Borneo 73 25 43
Madagascar
Baffin Devon Ellesmere Europe
Sumatra
Honshu 21 82 3745
Britain
Victoria Greenland Hainan Hispaniola
Ellesmere

Outils graphiques en statistique

Celebes 840 13 30
New Zealand (S)
Java Hokkaido Honshu Iceland
New Zealand (N)
Newfoundland 30 89 40
Cuba
Luzon Ireland Java Kyushu
Iceland
Mindanao 33 49 14
Ireland
Novaya Zemlya Luzon Madagascar Melville
Hokkaido
Hispaniola 42 227 16
Sakhalin
Moluccas Mindanao Moluccas New Britain
Tasmania
Celon 36 29 15
Banks
Devon New Guinea New Zealand (N) New Zealand (S)
Tierra del Fuego
Southampton 306 44 58
Melville
Axel Heiberg Newfoundland North America Novaya Zemlya
Spitsbergen
New Britain 43 9390 32
Taiwan
Kyushu Prince of Wales Sakhalin South America
Timor 13 29 6795
Prince of Wales
Hainan Southampton Spitsbergen Sumatra
Vancouver
16 15 183
0 5000 10000 15000 Taiwan Tasmania Tierra del Fuego
14 26 19
Timor Vancouver Victoria
13 12 82
97
Statistique descriptive sur les données islands Questions :

> mean ( x )
boxplot
[ 1 ] 3.734162 data : sqrt(islands)

> var ( x )

120
[ 1 ] 9.134685

100
> quantile (x , c (.25 ,.5 ,.75))

80
25% 50% 75%

60
Comment représenter graphiquement l’échantillon d’une variable
1.794741 2.990066 4.326417

40
Q1 mediane Q3

20
• continue

0
>b o x p l o t ( x ) • discrète
• qualitative ?
Construction du boxplot :
Les points représentent les valeurs
Comment représenter un tableau de données numériques ?
aberrantes. Elles sont définies comme
les observations en dehors de l’intervalle

I = [Q1 1.5(Q3 Q1 ), Q3 +1.5(Q3 Q1 )].

98 99

Histogramme Estimation de la densité pour des lois continues

On compare sur des données simulées suivant la loi gaussienne deux

estimateurs de la densité (hist / density) avec la densité théorique de la
h i s t ( s q r t ( i s l a n d s ) , b r e a k s = 12) loi
h i s t ( s q r t ( i s l a n d s ) , b r e a k s = c ( 4 ⇤ 0 : 5 , 10 ⇤ 3 : 5 , 7 0 , 1 0 0 , 1 4 0 ) )
p a r ( mfrow=c ( 1 , 2 ) )
Histogram of dat
x=rnorm ( 1 0 0 0 0 )
equidistant breaks NON - equidistant breaks
y=d e n s i t y ( x )

0.4
0.15

plot (y)
0.00 0.02 0.04 0.06

Density
19
0.10

0.2
Density

Density

11 h i s t ( x , p r o b a=T)
0.05

5
lines (y)

0.0
32
1 0 0 2 3 2
0.00

-3 -2 -1 0 1 2 3 4
0 20 40 60 80 100 140 0 20 40 60 80 100 140 z=seq ( min ( x ) , max ( x ) , 0 . 0 1 )
dat
sqrt(islands) sqrt(islands)
#s u p e r p o s e d e n s i t é t h é o r i q u e
WRONG histogram
density.default(x = dat) l i n e s ( z , dnorm ( z , 0 , 1 ) , l t y =2)
les options
0.4
15
Frequency

Remarque sur la fonction hist

10

Density

L’argument breaks fixe le nombre ou les bornes des classes. On peut

0.2
5

aussi utiliser nclass pour fixer le nombre de classes. proba =T : l’aire en dessous de la
0

0.0

courbe est égale à 1. Par défaut proba

0 20 40 60 80 100 140

sqrt(islands) -4 -2 0 2 4

N = 1000 Bandwidth = 0.2179

=F : l’aire est égale au nombre
100
d’observations. 101
Variables qualitatives Représentation graphique d’une variable quantitative
t a b=t a b l e ( c e n t r a l . p a r k . c l o u d )
On relève pendant un mois la méteo à Cental Park. La variable est à p i e ( tab )
barplot ( tab )
valeurs dans {clear, partly.cloudy, cloudy }.
> c e n t r a l . park . cloud Comment ordonner les modalités ?
[ 1 ] p ar tl y . cloudy p a r t ly . cloudy p a r t ly . cloudy clear en général, on les classe suivant l’ordre croissant des effectifs.
[ 5 ] p ar tl y . cloudy p a r t ly . cloudy clear cloudy
barplot ( sort ( tab ))
[ 9 ] p ar tl y . cloudy c l e a r cloudy p a r t l y . cloudy
[13] cloudy cloudy clear p a r tl y . cloudy
[17] p a r t l y . cloudy c l e a r clear clear pie(tab)

Pour résumer l’information : on calcule la répartition (ou la fréquence)

0 4 8
clear
des modalités dans l’échantillon.
clear partly.cloudy cloudy
> table ( c e n t r a l . park . cloud )
c e n t r a l . park . cloud
partly.cloudy
c l e a r par tly . cloudy cloudy
modalités ordonnés par effectif croissant
11 11 9
cloudy

0 4 8
102 cloudy clear partly.cloudy
103

Remarque sur le choix de la plage des ordonnées (ylim) Variables discrètes

On relève le salaire de
3 personnes On relève N le nombre d’appels téléphoniques reçus par un central
pendant un mois.
40

John
> sales
John J a c k P e t e r
30

Jack
> N [ 1 ] 6 16 3 6 8 10 9 9 10
45 44 46 etc
20

> t = t a b l e (N)
10

Représentation graphique

15
Peter
N
pie ( sales )
0

2 3 4 5 6 7 8

nb d'appels

10
John Jack Peter

barplot ( s a l e s ) 9 10 11 12 13 14 15 16 17

5
2 4 3
Modification de l’axe des ordonnés pour améliorer la lisibilité du

0
7 16 15 25 18 23 18 25 jours

graphique : 9 6 9 4 2
>
47

Représentation de la distribution et comparaison avec la densité de la loi

46

de Poisson de paramètre la moyenne empirique de l’échantillon (points en

45

b a r p l o t ( s a l e s , y l i m=c ( 4 0 , 4 7 ) , xpd=F )
rouge).
44
43

0.12
John Jack Peter

J= 2 : 1 7

table(N)/length(N)
104 105

0.08
p l o t ( t a b l e (N) / l e n g t h (N) , y l a b=" " )

0.04
p o i n t s ( J , d p o i s ( J , mean (N ) ) )

0.00
2 3 4 5 6 7 8 9 11 13 15 17

N
La fonction barplot pour des variables quantitatives Comparaison de 2 variables quantitatives indépendantes

On compare le nombre d’appels par mois de deux standards

0.00 0.02 0.04 0.06 0.08 0.10 0.12

fréquence
>b a r p l o t ( t / l e n g t h (N) )
comparaison des quantiles

10 15 20 25

10 12 14 16
2 4 6 8 10 12 14 16

nd d appels par jour

15
Attention : Les labels des barres ne coincident pas avec la coordonnée

N2
10

8
sur l’axe des abscisses.

6
5
5

4
0
Illustration : on veut ajouter sur graphique précédent la densité de la loi
1 2 2 5 8 12 17 5 10 15

de Poisson de paramètre égal à la moyenne empirique de l’échantillon

b o x p l o t (N, N2 , c o l=c ( " g r e e n " , " o r a n g e " ) )
(points en rouge), on obtient t=t a b l e ( N2 )
x t = as . i n t e g e r ( names ( t ) ) + . 3
J= 2 : 1 7
# GRAPHIQUE CORRIGÉ p l o t ( t a b l e (N ) , c o l=" g r e e n " , lwd =2 , x l a b=" " , y l a b=" " )
# GRAPHIQUE FAUX
b r = b a r p l o t ( t a b l e (N) / l e n g t h (N) ) l i n e s ( x t , t , c o l=" o r a n g e " , t y p e="h" , lwd =2)
b a r p l o t ( t a b l e (N) / l e n g t h (N) )
p o i n t s ( br , d p o i s ( J , mean (N ) ) )
p o i n t s ( J , d p o i s ( J , mean (N ) ) )
q q p l o t (N, N2 )
GRAPHIQUE CORRIGÉ

GRAPHIQUE FAUX abline (0 ,1)

0.00 0.02 0.04 0.06 0.08 0.10 0.12
0.00 0.02 0.04 0.06 0.08 0.10 0.12

fréquence

Pour les variables aléatoires continues, on compare les estimations de la

fréquence

2 4 6 8 10 12 14 16
2 4 6 8 10 12

nd d appels par jour

14 16

densité : histogrammes, estimateur à noyau.

106 107
nd d appels par jour

Lien entre des variables qualitatives dépendantes Autres graphiques pour résumer la table de contingence

prev grade
Données : pour 122 étudiants, on 1 B+ B+
2 A- A-
dispose de la note obtenue cette 3 B+ A-

14

14
20
A A
année (A+ ....F) et celle de l”année 4 F F A- A-

12

12
5 F F B+ B+
précédente B B

15
6 A B

10

10
B- B-
7 A A C+ C+
C C

8
8 C D D D

10
F F

6
table(grades)
> table(grades)

4
5
A A- B+ B B- C+ C D F
grade

2
prev A A- B+ B B- C+ C D F
A 15 3 1 4 0 0 3 2 0

0
A

A- 3 1 1 0 0 0 0 0 0 A B+ B- C F A A- B B- C F A B+ B- C F
grade

B+ 0 2 2 1 2 0 0 1 1 grade prev prev

B+ A-

B 0 1 1 4 3 1 3 0 2
B- 0 1 0 2 0 0 1 0 0
C C+B- B

C+ 1 1 0 0 0 0 1 0 0 > b a r p l o t ( tab , x l a b=" g r a d e " , l e g e n d . t e x t=T , c o l=rainbow ( 1 0 ) , b e s i d e=T)

C 1 0 0 1 1 3 5 9 7 > b a r p l o t ( t ( t a b ) , x l a b=" p r e v " , c o l=rainbow ( 1 0 ) )
F D

prev D 0 0 0 1 0 0 4 3 1 > b a r p l o t ( t ( t a b ) , x l a b=" p r e v " , l e g e n d . t e x t=T , c o l=rainbow ( 1 0 ) , b e s i d e=T)

F 1 0 0 1 1 1 3 4 11
> plot(table(grades))

108 109
Utilisations des couleurs Utilisation des couleurs, de la taille et du type de marques

Données : pour 150 enfants, on relève les informations suivantes : age / On représente la taille en fonction du poids.
poids / taille /sexe
Pour visualiser le lien entre les

60

60
> kid.weights
variables : on utilise la fonction

50

50
age weight height gender 150

pairs sur les 3 premières colonnes.

height

height
1 58 38 38 M

40

40
2 103 87 43 M
On peut ajouter en option la couleur

30

30
3 87 50 48 M

20

20
4 138 98 61 M des points définies à partir de la 0

10

10
5 82 47 47 F
variable gender
20 40 60 80 100 120 140 20 40 60 80 100 120 140

weight weight

20 60 100 140 20 60 100 140

80 120

80 120
age age

La couleur des points est définie par Le type de points est définie par la
40

40
0

0
la variable gender variable gender
60 100 140

60 100 140
weight weight

La taille des points est Le dégradé de couleurs est défini à

20

20

partir de la variable age.

50

50
height height
proportionnelle à l’age des individus.
30

30
10

10
0 40 80 120 10 30 50 0 40 80 120 10 30 50
plot ( weight , height , c e x =ag e /max ( ag e ) ⇤ 2 , c o l=as . i n t e g e r ( g e n d e r ) )
p l o t ( w e i g h t , h e i g h t , c o l = rainbow ( 1 5 0 ) [ ag e ] , pch=as . i n t e g e r ( g e n d e r ) )
p a i r s (X [ , 1 : 3 ] ) p a i r s (X [ , 1 : 3 ] , c o l=X [ , 4 ] ) 110 111

histogramme : le choix du nombre de classes

1000 Normal Random Variates

0.4

0.4
1. On simule un échantillon

0.3

0.3
Density

Density
suivant la loi gaussienne de

0.2

0.2
0.1

0.1
taille 1000

0.0

0.0
2. On trace l’histogramme pour
-4 -2 0 2 4 -4 -2 0 2 4

x x

Inférence statistique différentes valeurs du nombre

de classes.

0.4

0.4
3. On compare l’histogramme
0.3

0.3
Density

Density
0.2

0.2
avec la densité théorique de loi
0.1

0.1
gaussienne (courbe en rouge)
0.0

0.0
-4 -2 0 2 4 -4 -2 0 2 4

x x

hist(x, nclass=*** ,proba=T) ou hist(x,proba = T) par

défaut le nombre de classes est optimisé pour des échantillons gaussiens

112
Précision sur la fonction density Illustration sur un mélange gaussien

Soit X1 , · · · , Xn n variables aléatoires i.i.d. suivant la loi de densité f . On teste la fonction density sur des données simulées 2 suivant un
La fonction density calcule l’estimateur de f suivant mélange de lois gaussiennes :
n ✓ ◆ 1 1 3 1
ˆ 1X 1 x Xk f (x) = p e
1
2 (x 5)2
+ p e
1
2 (x+5)
2

fn (x) = kern 4 2⇡ 4 2⇡
n bwn bwn
k=1
same bandwidth, bandwidths = 0.1 bandwidths = 1
6 different kernels

• le choix du noyau kern

0.4

0.15
0.10

0.3
Density

Density

0.10
gaussian epanechnikov rectangular

0.2
1.0

1.0

1.0

0.05
0.8

0.8

0.8

0.08

0.1
0.6

0.6

0.6
Density

Density

Density

0.00
0.0
0.4

0.4

0.4
−5 0 5 −10 −5 0 5 10
0.2

0.2

0.2

0.06
Density
0.0

0.0

0.0
N = 50 Bandwidth = 0.1 N = 50 Bandwidth = 1
−1.0 0.0 0.5 1.0 −1.0 0.0 0.5 1.0 −1.0 0.0 0.5 1.0

bandwidths = 2 bandwidths = 3

0.04
triangular biweight cosine

0.08
1.0

1.0

1.0
0.8

0.8

0.8

0.06
0.08
0.6

0.6

0.6

0.02
Density

Density

Density

Density

Density
0.4

0.4

0.4

0.04
0.04
0.2

0.2

0.2

0.02
0.0

0.0

0.0

0.00
−1.0 0.0 0.5 1.0 −1.0 0.0 0.5 1.0 −1.0 0.0 0.5 1.0

0.00

0.00
−10 −5 0 5 10
−10 −5 0 5 10 −15 −5 0 5 10 15

• la dimension de la fenêtre bw. Par défaut le paramètre est optimisé N = 50 Bandwidth = 2.096 N = 50 Bandwidth = 2 N = 50 Bandwidth = 3

2. Simulation des données w=rbinom(50,1,1/4)

pour un échantillon gaussien sample = w*rnorm(50,5) +(1-w) * rnorm(50,-5)
113 114

Fonction de répartition empirique (ecdf) Estimation de la moyenne : comportement asymptotique.

La fonction de répartition empirique est définie par Soit X1 , . . . , Xn une suite de variables aléatoires iid, si E |X1 | < 1 alors
n
1X n
F̂n (t) = 1] 1t] (Xi )
n!1
! F (t) 1X
n p.s. Sn = Xi ! E (X1 )
i=1 n
i=1
C’est un estimateur de la fonction de répartition.
evolution de Sn en fonction de n : cas gausien
n=5000

1.5
1.0

1.0
#l o i gaussienne

0.5
0.8

x=rnorm ( n , 0 , 1 )

0.0
> x= rnorm ( 1 0 0 )
y=cumsum ( x ) / ( 1 : n ) 0 1000 2000 3000 4000 5000
c(0:(n − 1))/n

0.6

>Fn=e c d f ( x )
p l o t ( y , t y p e= ’ ’ l ’ ’ ) Index

>p l o t ( Fn , c o l=" g r e e n " )

0.4

evolution de Sn en fonction de n : cas cauchy

> z=seq ( min ( x ) , max ( x ) , 0 . 0 1 )
# l o i de

0 1 2 3 4
> l i n e s ( z , pnorm ( z ) , c o l =1 , lwd =2)
0.2

Cauchy ou S t u d e n t ( 1 d d l )

yc
x c=r t ( n , 1 )
0.0

y c=cumsum ( x c ) / ( 1 : n )

-2
−3 −2 −1 0 1 2 3
p l o t ( yc , t y p e= ’ ’ l ’ ’ ) 0 1000 2000 3000 4000 5000

sort(x) Index

115 116
Les fonctions pour les tests statistiques classiques : Exemple : Test de Student t.test()

X1 , . . . , Xn iid N (1, 1) et Y1 , . . . , Ym iid E(1)

Test H0 : E (X ) = E (Y ) vs H1 : E (X ) 6= E (Y )
t . t e s t ( x , mu=5 , a l t=" two . s i d e d " ) #s t u d e n t
t . t e s t ( x , y , a l t=" l e s s " , c o n f =.95) > x = rnorm ( 1 0 0 , 1 , 1 )
> y = rexp (200 ,1)
var . text ( x , y ) # comparaison v a r i a n c e > t . test (x , y)
cor . t e s t ( x , y ) # non c o r r e l a t i o n
chisq . test (x , y) #i n d é p e n d a n c e Welch Two Sample t t e s t
Box . t e s t ( z , l a g = 1 ) #non c o r r e l a t i o n data : x and y
t = 0.2178 , d f = 1 7 8 . 4 4 6 , p v a l u e = 0 . 8 2 7 8
shapiro . test (x) #n o r m a l i t é a l t e r n a t i v e h y p o t h e s i s : t r u e d i f f e r e n c e i n means i s n o t e q u a l t o 0
k s . t e s t ( x , " pnorm " ) #n o r m a l i t é K S 95 p e r c e n t c o n f i d e n c e i n t e r v a l :
ks . t e s t ( x , y ) # même d i s t r i b u t i o n 0.2648092 0 . 2 1 2 1 6 0 8
sample e s t i m a t e s : mean o f x : 0 . 9 5 4 4 1 2 7 mean o f y : 0 . 9 8 0 7 3 6 9

117 118

Test d’ajustement Validité de l’algorithme

On construit un générateur de nombres aléatoires suivant la loi normale Pour différentes valeurs de n = 1, 5, 10
en utilisant le Théorème Central Limite appliqué à des variables aléatoires
1. On teste la normalité à l’aide de deux procédures de test, le test de
iid suivant la loi uniforme sur (0, 1).
Kolmogorov et celui de Shapiro.
2. On compare la fonction de répartition empirique et la fonction de
Méthode de simulation non exacte ....
répartition théorique.
U = 1, . . . , Un iid suivant la loi uniforme 3. On compare l’estimation de la densité par un histogramme et la
r n densité théorique.
n 1 1X
(Ūn ) ) X ⇠ N(0, 1) Ūn = Ui X = simU ( 1 0 0 0 , n )
12 2 n
i=1
t e s t 1=k s . t e s t (X , " pnorm " , 0 , 1 ) $p . v a l u e
La générateur s’écrit t e s t 2 = s h a p i r o . t e s t (X) $p . v a l u e
simU< f u n c t i o n ( t a i l l e , s i z e )
p l o t ( s o r t (X ) , ( 1 : s i z e ) / s i z e , t y p e=" s " , c o l=" o r a n g e " )
{
l i n e s ( seq ( 3 , 3 , . 1 ) , pnorm ( seq ( 3 , 3 , . 1 ) ) , c o l=" r e d " )
y = m a t r i x ( r u n i f ( t a i l l e ⇤ s i z e ) , n c o l=s i z e )
t i t l e ( p a s t e ( "n=" , n , " e t p v a l u e ( k s )=" , f l o o r ( c ( t e s t 1 ) ⇤ 1 0 0 ) / 1 0 0 ) )
( a p p l y ( y , 1 , mean) 1/ 2 ) ⇤ s q r t ( t a i l l e / 1 2 )
h i s t (X , c o l=" o r a n g e " , x l i m=c ( 3 , 3 ) , p r o b a=T , main=" " )
}
l i n e s ( seq ( 3 , 3 , . 1 ) , dnorm ( seq ( 3 , 3 , . 1 ) ) , c o l=" r e d " )
t i t l e ( p a s t e ( " p v a l u e ( s h a p i r o )=" , f l o o r ( c ( t e s t 2 ) ⇤ 1 0 0 ) / 1 0 0 ) )
119 120
Les résultats Régression linéaire

Le modèle le plus simple

n= 1 et pvalue(ks)= 0 n= 5 et pvalue(ks)= 0.01 n= 10 et pvalue(ks)= 0.81
y = ax + b + ✏
(1:size)/size

(1:size)/size

(1:size)/size
0.8

0.8

0.8
Pour réaliser une régression linéaire, par la méthode des moindres carrés,
0.0

0.0

0.0
on utilise la fonction lm .
−1.5 −0.5 0.5 1.0 1.5 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 Sur les données cars
sort(X) sort(X) sort(X) > lm ( c a r s $ d i s t ~ c a r s $ s p e e d )

120
Call :

100
lm ( f o r m u l a = c a r s $ d i s t ~ c a r s $ s p e e d )

80
pvalue(shapiro)= 0 pvalue(shapiro)= 0.03 pvalue(shapiro)= 0.18 Coefficients :

cars$dist

60
( Intercept ) cars $speed
0.4
Density

Density

Density

40
0.3
17.579 3.932

20
0.00

0.0

0.0
> p l o t ( c a r s $ speed , c a r s $ d i s t )

0
> a b l i n e ( r e g , c o l =3 , lwd = 2 ) 5 10 15 20 25

−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 cars$speed

X X X 121 122

Visualisation : plot.lm Extension

Residuals vs Leverage Residuals vs Fitted

3

49 23 49
23
40

0.5

35
2
Standardized residuals

20
1

Residuals

0
0
-1

-20

39
-2

Cook's distance

0.00 0.02 0.04 0.06 0.08 0.10 0 20 40 60 80

Leverage Fitted values

lm(cars$dist ~ cars$speed) lm(cars$dist ~ cars$speed)

Normal Q-Q Scale-Location

• Régression multiple lm(v ⇠ v1 + v2 + v3)
49
23

• Régression linéaire généralisée glm

3

49
23
1.5

35

35
2

Standardized residuals
Standardized residuals

• etc
1.0
1
0

• Méthodes non paramétriques

0.5
-1

0.0
-2

-2 -1 0 1 2 0 20 40 60 80

Theoretical Quantiles Fitted values

lm(cars$dist ~ cars$speed) lm(cars$dist ~ cars$speed)

Remarque : Si les données sont disponibles sous la forme d’une liste on

peut utiliser la syntaxe
> f i t = lm ( d i s t ~ s p e e d , c a r s )

123 124
Régression non paramétrique Polynômes Locaux

>data ( c a r s )
>a t t a c h ( c a r s )
>p l o t ( s p e e d , d i s t )
> l i n e s ( ksmooth ( s p e e d , d i s t , " n o r m a l " , b a n d w i d t h =2) , c o l =2)
> l i n e s ( ksmooth ( s p e e d , d i s t , " n o r m a l " , b a n d w i d t h =5) , c o l =3)

120
> l i n e s ( ksmooth ( s p e e d , d i s t , " n o r m a l " , b a n d w i d t h =10) , c o l =4)

100
>data ( c a r s )

80
>c a r s . l o = l o e s s ( d i s t ~ s p e e d , c a r s )

dist

60
>p = p r e d i c t ( c a r s . l o )
20 40 60 80 100 120

>p l o t ( c a r s )

40
> l i n e s ( seq ( 5 , 3 0 , 1 ) , p$ f i t , c o l =3)

20
0
dist

5 10 15 20 25

speed
0

5 10 15 20 25

speed

125 126

L’objet “série chronologique” est une liste qui contient

• les valeurs observées,

Séries Chronologiques • la fréquence des observations,
• la date de la première observation
• la date de la dernière, etc...

127
Exemple
Pour créer une série chronologique, on utilise la fonction ts
> data ( USAccDeaths ) > b r u i t . b l a n c=t s ( rnorm ( 1 0 0 ) , f r e q u e n c y = 1 , s t a r t = c ( 1 ) , end=c ( 1 0 0 ) )
> USAccDeaths >p l o t ( b r u i t . b l a n c )
>a c f ( b r u i t . b l a n c , t y p e=" c o r r e l a t i o n " )
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec #ou t y p e =" c o v a r i a n c e " ou t y p e =" p a r t i a l "
1973 9007 8106 8928 9137 10017 10826 11317 10744 9713 9938 9161 8927
1974 7750 6981 8038 8422 8714 9512 10120 9823 8743 9129 8710 8680
1975 8162 7306 8124 7870 9387 9556 10093 9620 8285 8466 8160 8034
Series Bruit.Blanc
1976 7717 7461 7767 7925 8623 8945 10078 9179 8037 8488 7874 8647
1977 7792 6957 7726 8106 8890 9299 10625 9302 8314 8850 8265 8796

1.0
1978 7836 6892 7791 8192 9115 9434 10484 9827 9110 9070 8633 9240

0.8
> p l o t ( USAccDeaths )

0.6
Bruit.Blanc

ACF

0.4
7000 8000 9000 10000 11000

0.2
-1
USAccDeaths

0.0
-2

-0.2
0 20 40 60 80 100 0 5 10 15 20

Time Lag
1973 1974 1975 1976 1977 1978 1979

Time

128 129

Filtre Étude préliminaire de la série USAccDeaths

periodogram = function ( t r a j )
{
n = length ( t r a j )
• Premier exemple de filtre : (I Lp )d , on utilise la fonction diff f r e q = 1 : ( n %/% 2 ) / ( n )
> d i f f ( x , l a g=p , d i f f e r e n c e s =d ) p e r i o d o g r a m = Mod( f f t ( t r a j ) ) [ 2 : ( n %/% 2 + 1 ) ] ^
2/ ( 2 ⇤ p i ⇤ n )
plot ( freq , periodogram , type = " l " )
• La fonction filter permet d ’appliquer des filtres linéaires :
}
y= f i l t e r ( x , s i d e s= 1 , method = ‘ ‘ c o n v o l u t i o n ’ ’ , f i l t e r =c ( 2 , 3 ) )
Series USAccDeaths Series USAccDeaths
# y [ i ] = 2 ⇤ x [ i ] + 3 ⇤ x [ i 1]

0.2 0.6 1.0

y= f i l t e r ( x , s i d e s= 2 , method = ‘ ‘ c o n v o l u t i o n ’ ’ , f i l t e r =c ( 2 , 3 , 4 ) )

−0.4 0.0 0.4

Partial ACF
ACF
# y [ i ] = 2 ⇤ x [ i 1] + 3 ⇤ x [ i ] +4] ⇤ x [ i +1]

−0.4
y= f i l t e r ( x , method = ’ ’ r e c u r c i v e ’ ’ , f i l t e r =c ( 2 , 3 ) ) 0.0 0.5 1.0 1.5 0.2 0.6 1.0 1.4

#y [ i ] = x [ i ] + 2 ⇤ y [ i 1] + 3 ⇤ y [ i 2] >a c f ( USAccDeaths ) Lag Lag

>p a c f ( USAccDeaths )
>p e r i o d o g r a m ( USAccDeaths )
Ces fonctions permettent en particulier de créer un générateur de

periodogram

1500000
processus ARMA.

0
0.0 0.1 0.2 0.3 0.4 0.5

freq

130 131
modélisation et prévision AR

La fonction ar permet d’estimer les paramètres d’un modèle AR

>a r ( l h , a i c = FALSE , o r d e r . max = 4 ) # on f i x e p=4
Ap (L)X [t] = e[t] Call :
a r ( x = l h , a i c = FALSE , o r d e r . max = 4 )
Si l’ordre p n’est pas précisé, le meilleur modèle AR pour le critère AIC Coefficients :
est sélectionné. 1 2 3 4
0.6767 0.0571 0.2941 0.1028
>data ( l h ) O r d e r s e l e c t e d 4 s i g m a ^2 e s t i m a t e d as 0 . 1 9 8 3
>a r ( l h )
Call : ar ( x = lh )
ATTENTION x doit être une série chronologique (x=ts(x))
Coefficients :
1 2 3
0.6534 0.0636 0.2269
Order s e l e c t e d 3 s i g m a ^2 e s t i m a t e d as 0.1959

132 133

Autour des modèles ARMA Modèle ARMA : Estimation

• La fonction ar permet d’estimer les paramètres d’un processus AR.

• Pour les modèles ARMA d’ordre (p,q)

• ARMAacf pour le calcul des covariances théorique d’un modèle X [t] = a[1]X [t 1]+...+a[p]X [t p]+e[t]+b[1]e[t 1]+...+b[q]e[t q]
ARMA
on utilise la fonction arima, la syntaxe est
• ARMAtoMA pour le développement en MA infinie d’un modèle ARMA
out=arima(x,order=c(p,0,q))
• arima.sim pour simuler des trajectoires d’un modèle ARMA ou
la sortie out est une liste contenant :
ARIMA
out$coef : estimation des coefficients,
out$resid : estimation des résidus e[t]
ATTENTION x doit être une série chronologique (x=ts(x))

134 135
modélisation et prévision SARIMA Exemple

Plus généralement, la fonction arima permet d’estimer les paramètres data ( USAccDeaths )
d’un modèle SARIMA a = c ( USAccDeaths )
USAcc= t s ( a [ 1 : 6 0 ] , f r e q u e n c y =12 , s t a r t=c ( 1 9 7 3 , 1 ) )
Ap (L)↵P (Ls )Y [t] = Q (L
s
)Bq (L)e[t] avec Y [t] = (I L)d (I Ls )D X [t]
f i t = a r i m a ( USAcc , o r d e r=c ( 0 , 1 , 1 ) , s e a s o n a l= l i s t ( o r d e r=c ( 0 , 1 , 1 ) ) )
Call :
la syntaxe est la suivante :
a r i m a ( x = USAcc , o r d e r = c ( 0 , 1 , 1 ) ,
o u t=a r i m a ( x , o r d e r=c ( p , d , q ) , s e a s o n a l = l i s t ( order = c (0 ,1 , 1)))
s e a s o n a l= l i s t ( o r d e r=c (P , D,Q) , p e r i o d=s ) ) Coefficients :
ma1 sma1
la sortie est une liste contenant : 0.4343 0.4419
Approx s t a n d a r d e r r o r s :
out$coef : estimation des coefficients, ma1 sma1
out$aic : critère AIC, 0.1368 0.0122
out$resid : estimation des résidus e[t] s i g m a ^2 e s t i m a t e d 1 1 4 2 7 6 : l o g l i k e l i h o o d = 341.73 , a i c = 687.46
option : include.mean=F ou T

136 137

Prévision pour des modèles SARIMA Plus généralement : la fonction predict

Pour prévoir à l’horizon h, on utilise la fonction predict

>o u t=a r i m a 0 ( . . . )
> p=p r e d i c t ( out , h )
p$ p r e d #c o n t i e n t l e s p r é v i s i o n s
p$ s e #e r r e u r s de p r é v i s i o n ( é c a r t t y p e )
> methods ( p r e d i c t )
[ 1 ] " p r e d i c t . ar " " predict . arima "
[3] "predict . loess " " predict . ppr "
Exemple :
11000

[ 5 ] " p r e d i c t . smooth . s p l i n e " " predict . smooth . s p l i n e . f i t "

[ 7 ] " p r e d i c t . glm " " predict . lm "
10000

p = p r e d i c t ( f i t , n . ah ea d =12)
[ 9 ] " p r e d i c t . mlm"
USAccDeaths

p l o t ( USAccDeaths , c o l=" o r a n g e " )

9000

l i n e s ( p$ p r e d , c o l=" r e d 3 " )
l i n e s ( p$ p r e d +1.96 ⇤p$ se , c o l=" b l u e 2 " )
8000

l i n e s ( p$ p r e d 1.96 ⇤p$ se , c o l=" b l u e 2 " )

7000

1973 1974 1975 1976 1977 1978 1979

Time

138 139
Lissage exponentiel La librairie forecast

La méthode de lissage de Holt & Winter est disponible sous R. Les

fonctions s’appelle HoltWinters et predict.HoltWinters.
Cette librarie développée par Rob J Hyndman contient des méthodes
Exemple : sur la série data(co2)
pour l’analyse et la prévision de séries temporelles univariée
data ( co2 )
m = H o l t W i n t e r s ( co2 ) #l i s s a g e 1. le lissage exponentiel, Holt Winter, ...
p = p r e d i c t (m, 5 0 , p r e d i c t i o n . i n t e r v a l = TRUE) # prevision
p l o t (m, p ) 2. Modèles BATS (Exponential smoothing state space model with
Holt−Winters filtering
Box-Cox transformation, ARMA errors, Trend and Seasonal
components)

370
cette fonction peut aussi être utilisée 3. modélisation automatique ARIMA.

360
pour 4. modélisation automatique SARIMA .

350
Observed / Fitted
• des modèles multiplicatifs
Par défaut cette librarie utilise le traitement parallèle pour accélérer les
340
• le lissage simple si l’on impose 330
calculs dans la sélection automatique des modèles .
les coefficients.
320

1960 1970 1980 1990 2000

Time

140 141