Cours sur la statistique avec le logiciel R
Cours sur la statistique avec le logiciel R
Chapitre 1 Introduction
L’objectif de ce cours est de mettre en évidence les liens entre toutes les notions en statistique rencontrées dans les années antérieures par les étudiants : statistiques des-criptives, variables aléatoires et statistique mathématique. Le logiciel R est utilisé pour illustrer les applications des outils de statistique.
1.1 Plan de cours
– Statistiques descriptives
– Graphiques sous R : personnalisation des graphes
– Modèle linéaire : régression simple et multiple, ANOVA sous R
– Estimations ponctuelles et par intervalles de confiance
– Tests paramétriques
– Tests non paramétriques
– Outils R pour les valeurs extrêmes
1.2 Vocabulaire
Une statistique
La statistique / Les statistiques
- un nombre calculé à partir des observations
- un domaine
- le résultat de l’application de méthode
- la collecte des données
- le traitement (descriptive)
- l’interprétation (exploratoire)
- la prévision et la décision (décisionnelle)
1.3 Statistique et probabilité
Les probabilités sont l’outil de base du statisticien, car le “hasard” intervient à plusieurs niveaux en statistique : la répartition des données, le bruit, etc..
Introduction
1.4 Les avantages de R
– S-PLUS
– méthodes récentes
– multi-plateforme
– gratuit
…
4. Fonction de répartition empirique
…
ecdf(x)
0.0 0.2 0.4 0.6 0.8 1.0 Fn(x) |
5. Régression linéaire
yi = a + bxi + "i
On considère un jeu de donnés cars fourni par R.
speed dist | ||
1 | 4 | 2 |
2 | 4 | 10 |
1.5 Quelques exemples | ||
3 | 7 | 4 |
4 | 7 | 22 |
5 | 8 | 16 |
6 | 9 | 10 |
7 | 10 | 18 |
8 | 10 | 26 |
9 | 10 | 34 |
- 11 17
- 11 28
- 12 14
...
Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared: 0.6511,Adjusted R-squared: 0.6438
F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12
…
Chapitre 2 Statistique descriptive et première session de R
Nous présentons dans ce chapitre les objets et les commandes élémentaires du logiciel R. Dans la suite, sont écrites à gauche des commandes à taper, et à droite des commentaires sur ces commandes ou des questions relatives.
Le symbole > apparaît automatiquement en début de chaque ligne de commandes.
Le symbole + apparaît en début de ligne si la précédente est incomplète.
Le symbole # permet d’insérer un commentaire.
Une petite astuce très utile lorsque vous tapez des commandes directement dans la console : en utilisant les flèches Haut et Bas du clavier, vous pouvez naviguer dans l’histo-rique des commandes tapées précédemment, que vous pouvez alors facilement réexécuter ou modifier.
2.1 Les commandes élémentaires à connaître
help()
help.start() L’aide au format “html” (web) q() Quitter R
example(plot)
demo()
2.2 Premiers pas
x=c(1,4,9) | La fonction c() concatène des scalaires ou des vecteurs. |
y=c(x,2,3) | |
c(1:4) | Le premier terme est 10, le dernier est 100 et le pas est 10. |
seq(10,100,10) | |
x=rep(0,10) | On crée un vecteur constitué de 0 répété 10 fois. |
rep(y,10) | On simule 20 v.a. i.i.d. suivant la loi normale standard. |
x=rnorm(20) | |
y=rexp(20) | On simule 20 v.a. i.i.d. suivant la loi exponentielle d’espérance 1. |
median(x) | |
mean(x) | |
var(x) | |
sd(x) | |
summary(x) | |
sum(x) | |
length(x) | Pour tracer la première “courbe” |
plot(x) | |
lines(x) | Pour ajouter une ligne |
points(y) | Pour ajouter un nuage de points |
hist(x) | |
boxplot(x) | |
barplot(x) |
Exercice
Pour une loi normale d’espérance 5 et de variance 2 :
- Faire une représentation graphique de sa fonction de répartition et de sa densité sur
[0; 10].
- Calculer la probabilité des évènements : X 0, X 5, 1 < X 3 et X > 10.
- Calculer entre quelles valeurs 95% des tirages de X sont compris.
2.3 Langage de programmation
- Les différents objets : variable, vecteur, matrice, liste, table, fonction
a. list : Les listes sont très pratiques pour retourner des objets complexes notamment en sortie de fonctions.
data=c(1,2,3,4)
name="ABC"
maliste=list(donnees=data,nom=name) list permet de combiner les objets de différents types et différents tailles dans une même structure.
b. data.frame : Les data.frame sont très utiles pour stocker des données.
before=c(1,2,3,4)
after=c(4,5,6,7)
type=c("A","B","C","D")
x=data.frame(before,after,type) data.frame permet de combiner des objets “homogènes” dans une même structure.
typeof(x$type)
x=data.frame(before,after,I(type)) I permet de garder le type initial.
- les données avec labels : factor factor(c(1,3,2,2,1),levels=1:3)
- matrice
matrix(1:12,nrow=3,ncol=4)
matrix(1:12,nrow=3,byrow=TRUE) Par défaut la matrice est remplie par colonne. cbind(1:3,4:6,7:9)
rbind(1:3,4:6,7:9)
e. fonction masomme=function(x,y){s=x+y s}
2. Importer et enregistrer les données : fichiers plats, read.table, write.table, save y=data.frame(a=I("a"),b=pi) write.table(y,file="y.csv",sep=",",col.names=T)
y2=read.table("y.csv")
y2=read.table("y.csv",sep=’,’,header=T) write.table(x,file="x.txt",sep=" ",col.names=T) write.table(x,file="x.txt",sep=" ",col.names=T,row.names=F) save(x,y,file="xy.RData")
load("xy.RData")
Exercice
- Créez deux vecteurs de dimensions quelconques et insérez le second vecteur entre les 2ème et 3ème éléments du premier vecteur.
- Saisissez deux matrices 3 3 A et B et calculez la somme et le produit de ces matrices et enfin inversez la matrice A si possible.
- Soient x un vecteur avec n “levels” et y un vecteur numérique de dimension n, que se passe-t-il si l’on tape y[x] ?
2.4 Statistique descriptive
- variable : quantitative, qualitative, modalité, effectif, tableaux de contingence Exemple : chickwts –masses de poulets en fonction de leur alimentation Nourriture Traduction
casein caséine horsebean fève
linseed | graine de lin |
meatmeal | farine animale |
soybean | soja |
sunflower | tournesol |
chickwts Données fournies par R
unique(chickwts$feed) Modalité de la variable “feed”
w=cut(chickwts$weight,3) Découper et transformer un vecteur en facteur.
table(w,chickwts$feed) Table de contingence
summary(chickwts)
plot(weight feed,data=chickwts)
- paramètre : position (min, max, moyenne, médiane, quantile), dispersion (étendue, variance, écart-type)
- lois de probabilités
nom de loi | nom en R | paramètres | ||
beta | beta | shape1 | shape2 | ncp |
binomial | binom | size | prob | |
Cauchy | cauchy | location | scale | |
chi-squared | chisq | df | ncp | |
exponential | exp | rate | ||
F | f | df1 | df2 | ncp |
gamma | gamma | shape | scale | |
geometric | geom | prob | ||
hypergeometric | hyper | m | n | k |
log-normal | lnorm | meanlog | sdlog | |
logistic | logis | location | scale | |
negative binomial | nbinom | size | prob | |
normal | norm | mean | sd | |
Poisson | pois | lambda | ||
Student | t | df | ncp | |
uniform | unif | min | max | |
Weibull | weibull | shape | scale | |
Wicoxon | wilcox | m | n |
Exercice
1. Tirez au hasard 100 nombres xi dans l’intervalle [0; 1] avec une probabilité uniforme puis donnez le nombre de xi supérieurs à 0:5.
- Tirez 100 couples de points (x; y) aléatoirement dans le carré [0; 1] [0; 1]. Calculez le centre de gravité du nuage de points, représentez le nuage de points obtenus dans une fenêtre graphique, puis y ajouter en rouge le centre de gravité.
- Un bulbe fleurit avec une probabilité égale à 0:7. On plante 100 bulbes dans un pré carré de 11 mètres de côté en respectant un espace d’environ 1 mètre entre chaque bulbe. Simulez la floraison de ces 100 bulbes. Combien de bulbes ont fleuri ? Sauvegardez les données dans un unique fichier texte.
2.5 Représentations graphiques
- Un graphique standard
plot(x,y,main="Titre de la figure",sub="sous-titre", xlab="abscisse",ylab="ordonnée")
Titre de la figure
0.0 0.5 1.0 1.5 2.0
abscisse sous-titre
2. Il est possible d’ajouter une légende, du texte, une ligne ... legend(-2.5,0.5,"sin(x)",col=4,pch=1,lty=1)
…
Histogram of xHistogram of x
En utilisant les données disponibles chickwts dans R, tracer la figure suivante.
Croissance de poulets
…
Chapitre 3 Modèle linéaire
3.1 Régression linéaire simple
4 postulats du modèle pour tous i = 1; : : : ; n,
– IE ("i) = 0
– Var ("i) = 2
– "i sont i.i.d. de loi gaussienne.
– Xi sont déterministes.
Les estimateurs
…
Loi du 2 àndegrés de liberté
Soient X1; : : : ; Xn n v. a. indép. de loi N (0; 1), alors
S = X12 + + Xn2 suit une loi du 2 à n degrés de liberté, notée 2(n).
Loi de Student à n degrés de liberté
La loi de Student à n degrés de liberté, notée T (n), est la loi du quotient
Modèle linéaire
où N suit une loi N (0; 1) et S suit une loi 2(n),NetSétant deux v. a. indép..
Loi de Fisher à n1 et n2 degrés de liberté
Soient S1 et S2 deux v. a. indép. de loi respectives 2(n1)et2(n2). Alors le quotient
- =S1=n1 S2=n2
suit une loi de Fisher à n1 et n2 degrés de liberté, notée F (n1; n2).
Théorème de Cochran
Soient E1 et E2 deux sous-espaces vectoriels orthogonaux de E = IRd de dimensions respectives k1 et k2 et soit Y un vecteur aléatoire de IRd de loi normale centrée isotrope de variance 2. Alors PE1 (Y ) et PE2 (Y ) sont deux v. a. gaussienne centrées indépendantes et kPE1 (Y )k2 (resp. kPE2 (Y )k2) est une loi 2 2(k1) (resp. 2 2(k2)).
Exercice 1
Montrer les propriétés des estimateurs suivantes.
…
Exercice 2
Pn ^222
Montrer que SCR = i=1(Yi Yi) est une v. a. de loi (n 2).
Exercice 3
Montrer que sous l’hypothèse nulle : = 0, t suit la loi T (n 2).
Exercice 4
Montrer que sous l’hypothèse nulle : = 0, F suit la loi F (1; n 2).
3.2 Régression linéaire multiple
…
Les 4 formules fondamentales sont issues de la minimisation en de la somme des carrés résiduelle (SCR), somme qui peut s’écrire matriciellement sous la forme
SCR ( ) = kY X k2 = (Y X )0(Y X ):
3.3 Autres utilisations du test de Student
Il peut être intéressant de tester la position des paramètres par rapport à des valeurs particulières, ce qu’autorise le test de Student, comme l’indique le tableau ci dessous.
…
La notation désigne le coefficient j de la variable X(j) ou le terme constant du modèle linéaire estimé par la MMCO.
3.4 Étude de cas : CHAMPA
ln Y = + 1 ln X(1) + 2 ln X(2) + 3 ln X(3)
Y : la consommation de champagne en volume par personne
X(1) : le revenu par personne
X(2) : le prix du champagne en euro constants
X(3) : le prix des liqueurs et apéritifs en euro constants
Y | X1 | X2 | X3 | |
1 | 42.53 | 57.21 | 76.60 | 73.60 |
2 | 38.73 | 59.11 | 80.65 | 72.91 |
3 | 40.00 | 61.47 | 86.76 | 66.97 |
4 | 45.39 | 64.03 | 85.35 | 63.20 |
5 | 51.74 | 67.60 | 84.11 | 55.08 |
6 | 65.39 | 71.71 | 81.74 | 59.20 |
7 | 72.38 | 75.50 | 80.95 | 65.60 |
8 | 59.04 | 76.16 | 91.73 | 59.54 |
9 | 61.26 | 78.01 | 96.46 | 56.80 |
10 | 75.55 | 81.91 | 95.41 | 61.25 |
11 | 82.53 | 86.64 | 96.16 | 73.02 |
12 | 90.47 | 93.04 | 98.40 | 88.91 |
Modèle linéaire
13 100.00 100.00 100.00 100.00
14 110.47 105.36 99.30 111.42
15 127.93 109.76 97.84 123.31
16 139.04 115.03 95.96 136.91
17 143.80 120.53 104.27 150.62
Dans une première étape, tester le modèle log-linéaire suivant.
ln Y = + 1 ln X(1) + 2 ln X(2) + 3 ln X(3)
- Le modèle est il globalement significatif ?
- Quelle est la signification économique des coefficients ? Les coefficients ont-ils un signe conforme à vos attentes ?
- Une variable ne doit-elle pas être éliminée ? Laquelle ? Pourquoi ?
Dans une seconde étape, tester le modèle log-linéaire suivant.
ln Y = + 1 ln X(1) + 2 ln X(2)
- Ce modèle est il statistiquement satisfaisant ?
- Testez l’hypothèse 1 > 1.
- Testez l’hypothèse « la baisse des prix du champagne entraîne un progrès de la consommation en volume ». Peut on dire la même chose de la consommation en valeur ?
Réalisez une régression ascendante.
- Présentez de manière synthétique les résultats obtenus.
- La régression ascendante donne t il les mêmes résultats que la régression descen-dante ? Qu’en concluez-vous ?
Chapitre 4 Les lois à queue épaisse
4.1 Lois -stables
Il existe plusieurs définitions équivalentes de lois stables ainsi que plusieurs paramé-trisations. Nous présentons ici trois définitions équivalentes.
…
Le nombre c dans (1) vérifie c = a + b pour un nombre 0 < 6 2. La loi de X est dite -stable. Si D = 0 on dit que X est strictement -stable. Remarquons que si X est de loi gaussienne N ( ; 2), alors les termes de gauche dans (1) sont respectivement N (a ; (a )2) et N (b ; (b )2), et le terme de droite est N (c + D; (c )2). Si on prend c2 = a2 + b2 et D = (a + b c) , l’égalité (1) est établie, c’est-à-dire que les lois normales sont 2-stables.
La difficulté technique dans l’étude de lois -stables est que, sauf dans les cas particu-liers ( = 0:5; 1 et 2), il n’y a pas de forme explicite pour la densité. Les seules informations utilisables pour les v.a. -stables sont leur fonction caractéristique.
Définition 4.1.2 La fonction caractéristique d’un v.a. -stable (0 < < 2) dans Rd
…
Les lois à queue épaisse
Cette définition montre que la loi stable dans Rd est spécifiée par un nombre entre 0 et 2, indice de stabilité, une mesure finie sur Sd1 dite mesure spectrale et un vecteur dans Rd.
Dans le cas unidimensionnel la sphère unité ne contient que deux points, i.e. S = f 1; 1g. La mesure spectrale se réduit à deux valeurs ( 1) et (1). La fonction carac-téristique peut être écrite comme suit
…
où = (1) + ( 1) et = ( (1) ( 1))= . La loi stable unidimensionnelle est déterminée par quatre paramètres : 2 (0; 2], 2 [ 1; 1], > 0 et 2 R1. Le paramètre contrôle l’asymétrie. La loi est dite biaisée totalement vers la droite si = 1 et biaisée totalement vers la gauche si = 1. Si = 0 la densité est symétrique par rapport à . Les paramètres et sont les paramètres d’échelle et de position.
Nous utilisons la notation S1( ; ; ; ) pour la loi -stable unidimensionnelle. Le fait que X a la loi S1( ; ; ; ) est noté par l’écriture “X S1( ; ; ; )”. La Figure présente l’influence des paramètres sur la forme de la densité d’une loi stable.
…
Figure 4.1 – Influence des paramètres sur la forme de la loi stable. Les densités de la loi
S1(1:5; ( ; 1); 0) = 1; 0; 1 (a) et de la loi S1( ; (0; ); ) (b).
Propriété Si X a la loi -stable avec 0 < < 2 alors
EkXkp < 1 pour tout 0 < p < ;
EkXkp = 1 pour tout p > :
La troisième définition montre que les lois -stables sont la seule limite possible non-triviale des sommes normalisées des v.a. i.i.d..
Définition 4.1.3 Un v.a. X dans Rd a une loi stable s’il possède un domaine d’attraction, i.e. s’il existe une suite de v.a. i.i.d. Y1; Y2; : : : dans Rd, une suite de nombres positifs fbng et une suite de vecteurs réels fang telles que où ) représente la convergence faible.
Si X est une variable aléatoire gaussienne et si les Yi sont des variables aléatoires i.i.d. de variance finie, alors est le théorème central limite ordinaire. En général, bn = n1= L(n) où L(x) est une fonction à variation lente.
Si X S1( ; ; ; ), alors on a
…
La Figure présente le log-log graphe de la queue de lois stables.
Log−log plot
…
Figure 4.2 – Log-log graphe de PfX > xg où X S1( ; 0; 1; 0), = 0:5; 1; 1:5 et 2.
4.2 Lois de valeurs extrêmes
Soit (Xn; n 1) une suite de variables aléatoires i.i.d.. Notons leurs maximum Mn = max(Xi; : : : ; Xn). Supposons qu’il existe une suite ((an; bn); n 1) telle que an > 0 et la suite ((Mn bn)=an; n 1) converge en loi vers une limite non triviale. Alors à une
Les lois à queue épaisse
translation et un changement d’échelle près la fonction de répartition de la limite est de la forme suivante :
…
Exemples de convergence du maximum renormalisé
– Si X1 Unif[0; 1], alors la suite (n(Mn 1); n 1) converge en loi vers une loi de Weibull.
– Si X1 Exp(1), alors la suite (Mn log(n); n 1) converge en loi vers une loi de Gumbel.
Définition 4.2.1 La loi L est dite max-stable si pour tout n 2, (X1; : : : ; Xn) étant des variables aléatoires indépendantes de loi L, il existe an > 0 et bn 2 R, tels que
(max(X1; : : : ; Xn) bn)=an suit la loi L.
L’exercice suivant permet de vérifier que les lois de Weibull, Gumbel et Fréchet sont max-stables.
Exercice
Soit (Xn; n 1) une suite de variables aléatoires indépendantes, de même loi que X. On pose Mn = max(Xi; : : : ; Xn). Montrer que si X suit la loi de
– Weibull de paramètre , alors Mn a même loi que n 1= X ;
– Gumbel, alors Mn a même loi que X + log(n) ;
– Fréchet de paramètre , alors Mn a même loi que n1= X.
4.3 Estimateur par paquet
On divise l’échantillon 1; 2; : : : ; N en n groupes disjoints Gm;1; : : : ; Gm;n, chacun contient m éléments. En pratique on choisit d’abord m et on pose n = [N=m] ([a] représente la partie entière d’un nombre a positif). Quand N tend vers l’infini on a nm = [N=m]m N. De plus on demande n; m ! 1 quand N ! 1.
On estime d’abord le paramètre . Posons
Mm;i(1) = maxfk k j 2 Gm;ig; i = 1; : : : ; n; | (6) |
c’est-à-dire Mm;i(1) est la plus grande norme dans le groupe Gm;i. Notons m;i tel que | |
k m;ik = Mm;i(1); | (7) |
et Mm;i(2) = maxfk k j 2 Gm;inf m;igg; i = 1; : : : ; n; | (8) |
4.3 Estimateur par paquet
Mm;i(2) est alors la deuxième grande norme dans le même groupe.
On a pour chaque i
bm; | bm | ! ) c( 1 | ; 2 | ) quand m ! 1; | (9) | ||
Mm;i(1) | Mm;i(2) | 1= | 1= | ||||
où bm = m1= L(m) et c = (S)1= . Ici i= Pj, et 1; 2; : : : sont des variables aléatoires
j=1
…
où m;i est défini par 7. Les e.a. m;1; : : : ; m;n sont i.i.d. à valeurs dans S. On peut prouver que pour chaque i,
Sous les conditions convenables sur m et n, on a aussi la normalité asymptotique de l’estimateur de .