Probabilités et statistiques
Travaux pratiques avec Matlab
Stefan Le Coz1
1. Ce document est librement inspiré de polycopiés dont je disposais lors de son écriture (2005). À l’époque, je n’avais pas pris la peine de faire de bibliographie et je n’ai pas gardé trace des document dont je me suis inspiré. Malgré les apparences, il n’y a donc aucune prétention à l’originalité quant à la présentation qui est faite dans ce document, qui est plus une compilation des différentes sources dont je disposais à l’époque qu’un travail véritablement original.
Table des matières
1 Introduction à Matlab 2
1.1 Remarques préliminaires concernant l’apprentissage de Matlab . . 2
1.2 Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Calculs matriciels élémentaires . . . . . . . . . . . . . . . . . . . . 5
1.4 Calculs matriciels plus élaborés . . . . . . . . . . . . . . . . . . . 7
1.5 Graphiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.6 Entrées et sorties . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7 Opérations logiques, boucles et exécutions conditionnelles . . . . . 8
1.8 Fonctions et fichiers .m . . . . . . . . . . . . . . . . . . . . . . . . 9
1.9 Récursivité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.10 Fonctions en ligne . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Statistique descriptive 12
3 Simulation de variables aléatoires 14
3.1 Loi uniforme et simulation . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Simulation de lois discrètes . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Simulation de lois par la méthode du rejet . . . . . . . . . . . . . 16
3.4 Simulation des lois gaussiennes . . . . . . . . . . . . . . . . . . . 16
3.4.1 Méthode de Box-Müller . . . . . . . . . . . . . . . . . . . 16
3.4.2 Méthode polaire . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4.3 Quelques précisions . . . . . . . . . . . . . . . . . . . . . . 17
3.5 Simulation de lois par leur fonction de répartition . . . . . . . . . 18
3.6 Étude de la planche de Galton . . . . . . . . . . . . . . . . . . . . 18
4 Statistique inférentielle 20
4.1 Estimations et intervalles de confiance . . . . . . . . . . . . . . . 20
4.1.1 Estimations ponctuelle . . . . . . . . . . . . . . . . . . . . 21
4.1.2 Intervalles de confiance . . . . . . . . . . . . . . . . . . . . 21
4.1.3 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.2 Test du ?2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.2.1 Ajustement par la loi binomiale . . . . . . . . . . . . . . . 22
4.2.2 Ajustement par la loi de Poisson . . . . . . . . . . . . . . . 22
4.2.3 Ajustement par la loi normale . . . . . . . . . . . . . . . . 23
1
Chapitre 1
Introduction à Matlab
Un langage de programmation, tel Matlab, ne s’apprend pas tout à fait comme une théorie mathématique. En particulier, il ne faut pas espérer l’apprendre de manière très linéaire.
Ce chapitre se veut une aide pour commencer à apprivoiser Matlab. Il introduit, par petites touches, ce qui est utile pour programmer.
S’il ne faut retenir qu’une seule chose de cette introduction, c’est la remarque 1.2 (À l’aide!).
Matlab est un logiciel commercial de calcul matriciel développé par la société MathWorks. Son nom est la contraction de "Matrix Laboratory". Il consiste essentiellement en un interpréteur de commandes, écrites dans un langage de programmation spécifique appelé langage Matlab. Les commandes Matlab sont saisies et interprétées ligne à ligne dans une fenêtre (console). Comme nous le verrons plus loin, elles peuvent également être regroupées dans un fichier dont le nom se termine par .m.
Les variables Matlab sont toutes des tableaux, "arrays" en anglais, définies au moment de leur affectation. Il n’y a donc pas besoin de les déclarer. Un nombre complexe (ou réel, ou entier) est un tableau de taille 1 × 1, un vecteur ligne est un tableau de taille 1 × n et un vecteur colonne est un tableau de taille n × 1. Le langage Matlab a été conçu pour faciliter les opérations sur les matrices.
Le langage Matlab permet de manipuler des données de différents types dont certains sont imbriqués : entiers, nombres réels, nombres complexes, caractères, booléens. Ces données peuvent être assemblées en tableaux de différents types imbriqués : vecteurs, matrices, tenseurs (matrices à plus de 2 dimensions). Les chaînes de caractères sont des vecteurs de caractères.
2
En Matlab, les nombres réels sont représentés en virgule flottante avec 52 chiffres significatifs en base 2, soit un peu moins de 16 en base 10 (regarder help eps). Les opérations élémentaires sur les nombres réels sont :
+ ? ? / \ ˆ
qui représentent respectivement l’addition, la soustraction, la multiplication, la division à gauche et la division à droite, et enfin l’élévation à une puissance. Ces opérations s’appliquent également aux matrices, et dans ce cas, \ et / sont différentes.
Comme nous allons le voir, Matlab différentie minuscules et majuscules dans les noms de variables. L’affectation est notée =.
Matlab effectue ses calculs dans l’ensemble des tableaux à plusieurs dimensions, et donc en particulier les matrices, carrées ou non, et les vecteurs ligne et colonne. Bien entendu les opérations sur ces tableaux ne sont faites que lorsqu’elles ont un sens. En général, Matlab donne un sens assez intuitif aux opérations entre matrices. Cela dit, son langage contient de nombreuses spécificités que nous allons apprendre à utiliser.
Voici un exemple de code Matlab, à saisir ligne par ligne dans la fenêtre de commandes. Exemple 1.1 a=0
On crée une variable réelle nommée a, initialisée à 0. Matlab affiche sa valeur une fois que l’on a appuyé sur la touche Entr\’{e}e. En suffixant par un point-virgule, on évite l’affichage de la valeur de a.
a=0;
Pour connaître le contenu d’une variable, il suffit d’invoquer son nom.
a
Matlab différencie majuscules et minuscules. Ainsi, on peut créer la variable A, différente de a.
A=1
Vous pouvez rappeler les commandes précédement exécutées au moyen des touches fléchées de votre clavier (haut et bas).
Voici un calcul compliqué à base des variables a et A précédentes.
A*a+cos(a)/(1+sqrt(1+A^2))
La fonction sqrt donne la racine carrée ("square root"). La variable spéciale ans contient la dernière réponse de Matlab qui n’a pas été affectée à une variable. On peut lister les variables actuellement définies avec la commande whos. On voit que pour Matlab, les variables a et A sont des matrices 1 × 1.
whos
Pour une liste plus succincte, on peut utiliser la commande who. On peut détruire une variable au moyen de la commande clear. clear a % destruction de la variable a
Vérifions que la variable a n’existe plus.
who
On peut aussi détruire toutes les variables à l’aide de la commande clear. Pour obtenir de l’aide (en anglais!) sur une commande, on peut utiliser la commande help. help who
La commande helpwin permet d’obtenir une liste des commandes Matlab classée par thème. Nous pouvons être amené à rencontrer certaine valeurs particulières. Tester les commandes suivantes.
| 1/0 | % Infini positif Inf |
| -1/0 | % Infini n\’{e}gatif -Inf |
| 0/0 | % Not a Number (NaN) |
| 0*Inf | % Not a Number (NaN) |
| 1/Inf | % Donne bien z\’{e}ro |
| pi | % Donne bien le nombre pi |
help pi % Donne sa d\’{e}finition pour Matlab i % Racine carr\’{e}e (complexe) de -1. help i % Sa d\’{e}finition pour Matlab.
Remarquons qu’il est quand même possible d’utiliser i comme variable
Remarque 1.2 [À l’aide!] La commande help permet d’obtenir de l’aide sur les commandes Matlab, help help donne de l’aide sur l’aide! Enfin, la commande lookfor permet de lister les commandes Matlab par mots clés.
Remarque 1.3 [Affichages] Si X est un tableau :
— la commande disp(X) affiche son contenu sans afficher son nom. S’il est vide, rien n’est affiché. La variable X vide est symbolisée par [ ], on peut tester si X est vide grâce à la commande isempty(X).
— la commande display(X) affichera le nom de la variable (X ici) ainsi que son contenu, même s’il est vide. Cette commande est utilisée automatiquement par Matlab lorsque qu’une expression ne se termine pas par un pointvirgule.
Exercice 1.4 Tester et interpréter les commandes suivantes.
A= [2,4,8;3,9,27] A(2,3) b=A(2,:) c=A(:,3) B=A(1:2,1:2) A U=[l:20]
V=[l:2:20] V-U(l:2:20) W=[pi/2:1:pi]
Les commandes comme sqrt qui prennent un ou plusieurs arguments (ou paramètres) entre parenthèses constituent des fonctions. Nous verrons plus loin comment en créer de nouvelles.
Les tableaux suivants indiquent différents calculs possibles sur les matrices.
| Commande | Sens |
| A=[3,4;1,0] | Matrice 2 × 2 |
| B=(A>0) | B = matrice tq Bij=1 si Aij>0 et 0 sinon |
| V=[5;5] | Vecteur colonne de dimension 2, s’écrit aussi [5,5]’ |
| length (V) | Renvoie la longueur du vecteur V |
| size(A) | La fonction size renvoie la taille de A |
Figure 1.1 – Quelques exemples.
| A’ | transposée de A |
| A*B | Multiplication matricielle de A et B |
| A+B | Addition matricielle de A et B |
| A^99 | Puissance matricielle |
| expm(A) | Exponentielle matricielle |
| inv(A) | Matrice inverse |
| sqrtm(A) | Racine carrée matricielle (cf. help sqrtm) |
| logm(A) | Logarithme matriciel (cf. help logm) |
| A*V | Image du vecteur colonne V par la matrice carrée A |
| V’*A | Devinez |
| A\V | Solution du système linéaire AX=V (par pivot, et pas par inv(A)) |
| help slash | Aide explicative sur la division matricielle précédente |
Figure 1.2 – Opérations matricielles classiques.
| A=ones(n,n) | Matrice de taille n × n dont tous les éléments valent 1 |
| B=zeros(n,n) | Matrice nulle de taille n × n |
| C=eye(n,n) | Matrice identité de taille n × n. Que donne eye(3,7)? |
| [A,B,C] | juxtaposition horizontale des matrices A, B et |
| [A;B;C] | juxtaposition verticale des matrices A, B et C |
| N=0*J | Astuce pour obtenir une matrice N nulle de même taille que J sans connaître la taille de J. |
| sparse | Création d’une matrice «creuse» (cf help sparse) |
Figure 1.3 – Création de matrices.
| C=A.*B | C(i,j)=A(i,j)*B(i,j) |
| C=A./B | C(i,j)=A(i,j)/B(i,j) |
| C=A.3 | C(i,j)=A(i,j)*3 |
| C=3+A | C(i,j)=3+A(i,j) |
| C=3*A | C(i,j)=3*A(i,j) |
| C=A./3 | C(i,j)=A(i,j)/3 |
| C=cos(A) | C(i,j)=cos(A(i,j)) |
| C=log(A) | C(i,j)=log(A(i,j)). Ne pas confondre avec logm(A)! |
| G=sqrt (A) | C(i,j)=sqrt(A(i,j)). Ne pas confondre avec sqrtm(A)! |
| C=exp(A) | C(i,j)=exp(A(i,j)). Ne pas confondre avec expm(A)! |
| C=abs(A) | C(i,j)=|A(i,j)| |
Figure 1.4 – Opérations entrée par entrée.
L’aide sur les opérations élémentaires s’obtient par la commande help ops et celle sur les opérations élémentaires matricielles par help elmat.
L’ajout d’un point devant un opérateur arithmétique indique à Matlbab que les opérations sur la matrice se font composantes par composantes. Vous pouvez faire help arith pour de l’aide sur les opérations arithmétiques et help @ pour de l’aide sur les opérateurs en général et les caractères spéciaux.
Comme on vient de le voir, les sous matrices s’obtiennent en spécifiant des intervalles d’indices. Si i,j, k sont des entiers relatifs, alors :
— i:j est identique à i,i + 1, ,j,
— i:k:j est identique à i,i + k,i + 2k, ,j.
L’intervalle vide est représenté par la matrice vide [ ]. Une ligne ou une colonne entière peut être obtenue en utilisant le caractère : seul. Enfin, une matrice d’entiers E peut servir à spécifier les indices d’une matrice A, en écrivant A(E). Il ne faut pas confondre les expressions de la forme [ ], qui permettent de fabriquer des matrices, avec celles de la forme M( ), qui permettent de considérer une sous-matrice de la matrice M.
Le tableau 1.4 présente des opérations matricielle plus élaborées que les précédentes. Remarquons que les opérateurs de type sum servent très souvent lorsqu’il s’agit d’éviter les boucles.
| det(M) | Déterminant |
| rank(M) | Rang |
| trace(M) | Trace |
| D=eig(M) | Renvoie le vecteur colonne des valeurs propres de M |
| [P,D]=eig(M) | Diagonalise M=P*D*P?1 avec D diagonale |
| sum(V) | Somme des éléments du vecteur V |
| sum(M) | Somme des colonnes de la matrice M (renvoie un vecteur ligne) |
| sum(M,2) | Somme les éléments de la matrice M selon la dimension 2 |
| sum(sum(M) | Somme totale des éléments de la matrice M |
| cumsum(V) | Sommes cumulatives des entrées du vecteur V |
| cumsum(M) | Matrice des sommes cumulatives des colonnes de M |
| cumprod(V) | Produits cumulatifs des entrées du vecteur V |
| cumprod(M) | Matrice des produits cumulatifs des colonnes de M |
| max(M) | Renvoie un vect. ligne = au max sur chaque col . de M |
| max(max(M)) | Maximum des entrée de la matrice M. |
Figure 1.5 – Calculs matriciels élaborés
Pour l’aide sur les opérations matricielles élémentaires et sur les fonctions matricielles, on peut utiliser helpwin elmat ainsi que helpwin matfun
Le principe général des représentations graphiques est de se ramener à des calculs sur des matrices ou des vecteurs. Ainsi la repésentation d’une fonction de R dans R commencera par la création d’un vecteur d’abscisses, en général régulièrement espacées, auxquelles on applique la fonction pour créer le vecteur des ordonnées. La commande hold on permet de superposer des tracés successifs sur une même figure.
La commande input permet de demander à l’utilisateur de saisir des valeurs de variables. La commande pause permet de stopper l’exécution de Matlab pendant un temps déterminé, cf. help pause. La forme spéciale pause off désactive les pauses tandis que pause on les réactive. La commande save permet de sauvegarder le contenu des variables en cours ainsi que leur nom dans un fichier, dont
| mean(Y) | Moyenne arithmétique |
| median(Y) | Médiane |
| std(Y) | écart type normalisé en N ? 1 |
| plot(Y) | Tracé du vecteur Y. Abscisses? Ordonnées? |
| plot(Y,’r-’) | Autre tracé avec couleur et type de ligne |
| plot (X,Y,’r-’) | Vous comprenez |
| hist(Y) | Histogramme (10 classes par défaut, cf. help hist) |
| [Eff,Cl]=hist(Y,50) | Histogramme à 50 classes de Matlab |
| figure | Création d’une nouvelle fenêtre graphique |
| title (’Nom’) | Titre de la figure |
| xlabel(’Nom’) | Titre des abscisses |
| ylabel(’Nom’) | Titre des ordonnées |
| imagesc(A); | Trace la matrice, couleurs = valeurs relatives |
| hold off | On ne va pas utiliser la même fenêtre graphique |
| contour(A) | Trace les lignes de niveau |
Figure 1.6 – Fonctions graphiques.
le nom est par défaut . Ce fichier peut être lu par la commande load, qui restaure donc toutes les variables.
Pour Matlab, tout nombre peut être considéré comme une valeur logique, en identifiant les nombres non nuls à vrai (noté 1) et le nombre zéro à faux (noté 0). La commande boolean permet d’obtenir le booléen associé à un nombre vu comme une valeur logique. Les principales opérations sur les valeurs logiques sont : — la négation logique ~ , — le "ou" logique | , — le "et" logique &.
On peut également utiliser respectivement les fonctions or, not et and. Il y en a d’autres, faites donc un help ops. On peut créer des valeurs logiques à partir de nombres avec les opérations de comparaisons, qui activent automatiquement la nature logique de leurs arguments :
— "a égal à b" s’écrit a == b
— "a différent de b" s’écrit a ~= b
— "a supérieur strictement à b" s’écrit a > b
— "a supérieur ou égal à b" s’écrit a >= b
— "a inférieur strictement à b" s’écrit a < b
— "a inférieur ou égal à b" s’écrit : a <= b
On peut également utiliser respectivement les fonctions eq, ne, gt, ge, lt et leq. Le code Matlab R=randn(5,5);M=(R>=2.1) crée une matrice M de même taille que R qui contient des 0 là où R est inférieure à 2.1 et des 1 là où R est supérieure ou égale à 2.1. C’est donc une matrice "booléenne". La commande any(V) renvoie 1 si au moins l’un des éléments du vecteur V est non nul, et 0 sinon. La fonction all(V) renvoie 1 si tous les éléments du vecteur V sont non nuls, et 0 sinon. Enfin la fonction find(V) renvoie les indices correspondants aux éléments non nul du vecteur V. Ainsi, V(find(V>2)) renvoie les valeurs supérieures à 2 du vecteur V.
Les commandes Matlab isreal, ischar, issparse, isnan, isinf, isfinite permettent de tester la nature des variables. Ainsi, isreal(a) renverra 1 si la variable a est un réel et 0 sinon.
Le langage Matlab comprend les boucles du type for et while ainsi que les structures d’exécutions conditionnelles du type if. L’argument d’un if et d’un while est de type logique, ce qui n’est pas le cas de celui d’un for, qui est de type matriciel. La commande break permet de sortir de la boucle while ou for la plus interne.
Les commandes Matlab peuvent toutes être considérées comme des fonctions, c’est-à-dire des entitées nommées, qui prennent des paramètres éventuels (arguments) et qui renvoient des résultats (valeurs de retour). Beaucoup de fonctions Matlab, comme mean par exemple, sont en réalité écrites en Matlab, et le code Matlab correspondant est stocké dans un fichier dont le nom se termine par .m. Pour mean, il s’agit de mean.m. Pour ajouter de nouvelles fonctions à Matlab, il nous suffit d’écrire de nouveaux fichiers de ce type. Vous aurez compris qu’un fichier .m ne contient qu’une seule fonction, qui a le même nom que le fichier, au suffixe .m près.
Voici un code qui définit une fonction stat, qui prend comme paramètre un vecteur x et qui renvoie sa moyenne et son écart type. Ce code devra être stocké dans le fichier nommé stat.m pour que Matlab fasse le lien avec la fonction stat. Le commentaire de la deuxième ligne constitue l’aide qui est affichée lorsque l’utilisateur tape help stat.
— function [moyenne,ecartype]=stat(x)
— %STAT Renvoie la moyenne et l’\’{e}cart type du vecteur pass\’{e} en argument."
— n=length(x);
— moyenne=sum(x)/n;
— ecartype=sqrt(sum((x-moyenne)/n);
— return % En fait, inutile en fin de fonction
Le fichier stat.m devra être placé dans un répertoire que Matlab scrutera. En général, Matlab cherche automatiquement dans le répertoire en cours. Pour afficher le répertoire en cours, utilisez la fonction pwd, pour lister son contenu, utilisez la fonction ls et pour changer de répertoire courant, utilisez la fonction cd.
La commande Matlab type permet de lister le contenu du fichier .m d’une fonction. Ainsi, type mean va vous montrer le code source Matlab de la fonction mean. Un certain nombre de fonctions Matlab ne correspondent à aucun fichier .m, elles sont "internes" à Matlab pour une plus grande efficacité et l’on parle de fonctions "built-in". C’est par exemple le cas de la fonction type elle-même! Cette séparation entre fonctions internes et externes se retrouve dans la plupart des interpréteurs de langages. La commande exist permet de connaître le type d’une commande ou d’une variable, cf. help exist.
Sur la plupart des systèmes informatiques, les fichiers sont organisés en une structure arborescente de répertoires qui contiennent des fichiers et des sousrépertoires. La commande Matlab which donne le chemin du fichier .m associé à une fonction Matlab, lorsqu’il existe. Par exemple, which mean -ALL affichera tous les fichiers mean.m avec leurs chemins. Seul le premier sera utilisé par Matlab! La commande what liste les fichiers .m du répertoire en cours. Faites help what et help which pour plus d’information.
Matlab a besoin de savoir où chercher les fichiers .m dans l’arborescence des fichiers du système. Il recherchera automatiquement dans une liste de réperoires appelée "PATH". La commande path vous liste les répertoires du "PATH". Les commandes addpath et rmpath permettent d’ajouter et d’enlever des répertoires de la liste "PATH". Ces modifications ne sont pas permanentes et sont perdues lorsque vous quittez Matlab.
Par défaut, la liste "PATH" contient les répertoires de Matlab, qui contiennent eux-même les fichiers .m qui constituent Matlab. Afin d’accélérer la lecture des fichiers .m, Matlab maintient en permanence une base de données des fichiers .m avec le répertoire associé qui les contient. La commande rehash permet de mettre à jour cette base de données en relisant la liste "PATH".
Il faut enfin savoir que Matlab cherche d’abord les fichiers .m des fonctions dans le répertoire courant ("working directory"). Répétons-le, la commande pwd affiche le répertoire courant ("print working directory"), la commande cd change de répertoire courant ("change directory"), et les commande ls ou dir listent le contenu du répertoire courant.
En général, vous n’aurez pas à modifier la liste "PATH". Sa modificaton permanente dépend du système utilisé (Unix, MS-Windows, ).
Le langage Matlab autorise la récursivité, comme le montre l’exemple archiclassique suivant :
function fn=fact(n);
% fact cacule de fa\c{c}on r\’{e}cursive (et particuli\‘{e}rement inefficace) % la factorielle de n. if (n<=0), fn=1; else fn=n*fact(n-1); end
Exercice 1.5 En utilisant la fonction prod, écrire une fontion calculant n! de façon non récursive (et beaucoup plus efficace).
Pour des fonctions très légères, on peut procéder autrement que par la création d’un fichier .m : la directive inline permet de créer des fonction nouvelles à la volée. Voici un exemple :
On définit une nouvelle fonction trisin qui prend 3 paramètres a, b et c trisin=inline(’sin(a)*sin(b)*sin(c)’); On teste notre nouvelle fonction :
trisin(pi/2,pi/4,pi/2)
Matlab détermine automatiquement le nombre de paramètres de la fonction. Bien entendu, une fois que vous quitterez Matlab, cette fonction sera perdue car elle n’est pas sauvegardée dans un fichier.
Exercice 1.6 Apprenez à utiliser la fonction fplot pour tracer les graphes de fonctions définies avec inline.
Chapitre 2
Statistique descriptive
Exercice 2.1 On se donne les séries statistiques suivantes. Écrire un programme donnant pour chacune d’entre elles une représentation graphique des données (diagramme en baton ou histogramme selon les cas), tracant la fonction de répartition, et indiquant les caractéristiques de positions (mode, médiane, moyenne, quartiles) et de dispertion (étendue, variance, écart-type).
| Ménages | Effectifs | Fréquences |
| 1 personne | 5 845 140 | 0,2713 |
| 2 personnes | 6 366 948 | 0,2957 |
| 3 personnes | 3 821 700 | 0,1774 |
| 4 personnes | 3 371 484 | 0,1565 |
| 5 personnes | 1 439 144 | 0,0668 |
| 6 personnes ou plus | 695 736 | 0,0323 |
| Somme | 21 542 152 | 1,0000 |
Figure 2.1 – Ménages suivant le nombre de personnes du ménage en 1990, en
France.
Exercice 2.2 [Séries statistiques à deux variables] On se donne la série statistique double suivante.
Écrire un programme donnant une représentation graphique des données (nuage de points), traçant la droite de régression (en utilisant la méthode des moindres carrés) et indiquant le coefficient de corrélation entre les deux variables.
12
CHAPITRE 2. STATISTIQUE DESCRIPTIVE
| Classes d’âge | Effectifs | Classes d’âge | Effectifs | Classes d’âge | Effectifs |
| [15,17[ | 800 | [29,31[ | 51200 | [41,43[ | 800 |
| [17,19[ | 7 000 | [31,33[ | 32 000 | [43,45[ | 200 |
| [19,21[ | 14 350 | [33,35[ | 23 200 | [45,47[ | 80 |
| [21,23[ | 28 000 | [35,37[ | 16 000 | [47,49[ | 40 |
| [23,25[ | 44 000 | [37,391 | 8 000 | ||
| [25,271 | 54 800 | [39,41[ | 2 400 | ||
| [27,29[ | 56000 |
Figure 2.2 – Naissances en France en 1994 en fonction de l’âge de la mère (en millier).
| Firmes | CA | Firmes | CA | Firmes | CA |
| General Motors | 860 | Nestlé | 231 | PSA | 166 |
| Ford Motor | 713 | FIAT | 226 | BASF | 159 |
| Exxon | 563 | Sony | 216 | ABB | 158 |
| Royal Dutch Shell | 526 | Honda Motor | 215 | Pepsico | 158 |
| Toyota Motor | 441 | Elf Aquitaine | 208 | PEMEX | 156 |
| Hitachi | 412 | NEC | 205 | Nippon Steel | 156 |
| Matsushita | 377 | Daewoo | 197 | Mitsubishi Heavy I. | 155 |
| General Electric | 359 | Du Pont de Nemours | 194 | Amoco | 150 |
| Daimler Benz | 356 | RWE | 191 | Bayer | 148 |
| IBM | 355 | Japan Tobacco | 190 | Nippon Oil | 144 |
| Mobil | 331 | Philips | 186 | BMW | 144 |
| Nissan Motor | 317 | Mitsubishi Motors | 185 | Hewlett-Packard | 139 |
| Philip Morris | 298 | Texaco | 185 | Kansai Electric P. | 138 |
| Chrysler | 290 | EDF | 183 | Total | 137 |
| Siemens | 289 | Renault | 179 | Conagra | 131 |
| British Petroleum | 281 | Fujitsu | 177 | Motorola | 123 |
| IRÏ | 279 | Mitsubishi Electric | 176 | Petroleos de Venezuela | 123 |
| Volskwagen | 273 | ENI | 174 | Bat industries | 122 |
| Tokyo Electric P. | 271 | Chevron | 172 | Boeing | 122 |
| Toshiba | 260 | Hoeschst | 170 | Mazda Motor | 120 |
| Unilever | 251 | Procter & Gamble | 168 | Thyssen | 119 |
| VEBA | 243 | Alcatel Alsthom | 168 |
Figure 2.3 – Chiffres d’affaires (en milliards de francs), en 1994, des principaux groupes industriels mondiaux.
| Xi | 29 | 4 | 69 | 65 | 98 | 55 | 40 | 19 | 62 | 73 | 37 | 0 | 41 | 75 | 79 |
| Yi | 26 | 4 | 65 | 66 | 90 | 56 | 43 | 26 | 66 | 76 | 35 | 2 | 45 | 68 | 86 |
Chapitre 3
Simulation de variables aléatoires
Considérons une suite de chiffres x1, x2, ,xn de l’intervalle [0,1] ou d’entiers entre 0 et 9. Quand peut-on dire qu’elle est constituée de réalisations indépendantes d’une loi uniforme sur [0,1] (resp. sur {0, ,9})? Et bien dans tout les cas bien sûr. Cependant, l’intuition nous indique que ce sont les propriétés "statistiques" des suites qui sont importantes, et pas celles d’une suite donnée Pour certains, une suite qui peut être générée par un algorithme n’est pas aléatoire. De ce point de vue, une suite finie n’est jamais aléatoire!
Pour nous, simuler une loi de probabilité L consistera à écrire un programme informatique pouvant générer des suites finies dont on considérera que ce sont des réalisations indépendantes de loi L. Bien entendu, comme nous ferons appel le plus souvent à des algorithmes déterministes, les suites générées ne seront pas du tout aléatoires et encore moins indépendantes. On parlera de suites pseudoaléatoires. Cependant, elles se montrent suffisantes dans la pratique et l’on pourra vérifier expérimentalement qu’elles ont de "bonnes" propriétés en utilisant les tests statistiques.
La loi uniforme est la première loi que l’on simule. Diverses méthodes permettent ensuite de simuler à partir de la loi uniforme, un grand nombre de lois, plus ou moins quelconques. Dans la pratique, on met en œuvre des algorithmes rapides de génération de suites pseudo-aléatoires.
Matlab permet de simuler une loi uniforme (resp. normale) via la fonction rand (resp. randn).
| rand | renvoie un nombre selon la loi uniforme entre 0 et 1 |
| randn | renvoie un nombre selon la loi normale |
| rand(4,4) | Une matrice 4x4 dont les entrées sont des réalisations de loi uniforme sur [0,1] |
| randn(4,4) | Idem mais avec la loi normale |
| rand(1000,1) | 1000 réalisations pseudo-uniformes et pseudo-indépendantes |
Figure 3.1 – Uilisation de rand et randn
14
Les générateurs de nombres aléatoires sont très importants dans les applications. Ils permettent d’une part de simuler des phénomènes dont la modélisation fait appel à des lois de probabilités et d’autre part de sécuriser les échanges de données, dans les réseaux téléinformatiques par exemple. Dans la pratique, les générateurs aléatoires sont basés sur des algorithmes déterministes ou encore sur une source d’"entropie" extérieure difficilement prédictible liée au fonctionnement de la machinerie informatique (par exemple le temps séparant deux interruptions du processeur), ces deux approches pouvant être combinées dans un même processus.
Simuler la loi uniforme consiste à produire une suite de nombres qui peuvent être considérés comme autant de réalisations indépendantes de variables aléatoires uniformes sur [0,1]. Bien entendu, un algorithme est toujours déterministe lorsqu’on le connaît, et il vaut donc mieux parler de générateur pseudo-aléatoire de la loi uniforme. Mathématiquement, les n sorties successives d’un tel générateur seront considérées comme la donnée de U1(?), ,Un(?) pour un ? ? ? où les Ui sont des v.a.r. i.i.d. Ui : (?,F,P) ? [0,1] de loi uniforme.
La mémoire des ordinateurs étant finie, il est impossible de stocker des suites infinies. En conséquence, les nombres réels sont toujours approximés, avec une certaine précision. Soit donc p ? N? qui va jouer le rôle de précision. Un réel xn dans [0,1] s’écrira yn/p où yn est un entier dans {0, ,p}. Un algorithme très usité pour générer les yn consiste à utiliser la relation de récurrence xn ? yn/p où yn+1 ? ayn[m], a et m sont des entiers et m est assez grand.
En général, la période de (yn)n est plus petite que m. Plus précisément, si a et m sont premiers entre eux, et si 0 < y0< m, alors yn n’est jamais nul et la période de (yn)n n’est rien d’autre que l’ordre de a dans Zm. On rappelle que d’après le (petit) théorème de Fermat-Euler, l’ordre de a divise l’indicateur d’Euler de a (qui vaut m ? 1 quand m est premier).
Matlab 4 utilise un algorithme de ce type pour implémenter la fonction rand mais les versions plus récentes de Matlab utilisent un algorithme différent et plus performant. Dans certaines bibliothèques de calcul scientifique, on utilise (a,m) = (1313,259) ou encore (a,m) = (75,231 ? 1).
Soit p1, ,pn des nombres dans [0,1] tels que p1+···+pn = 1. En partitionnant l’intervalle [0,1] en morceaux adjacents de longueurs p1, ,pn, il vient :
Théorème 3.1 (Simulation d’une loi discrète à support fini) SoitU une v.a.r. de loi uniforme sur [0,1] etP := p1?x1 + ··· + pn?xn une loi discrète où lesxi sont tous différents. Alors la v.a.r.
x1I{U<p1} + x2I{p16U6p1+p2} + ··· + xnI{p1+···pn?16U61}
suit la loi P.
Exercice 3.2
(1) Écrire une fonction implémentant le résultat du Théorème 3.1.
(2) Simuler, sans utiliser d’instruction conditionnelle, la loi uniforme discrète sur {1, ,10}.
(3) Idem pour la loi binomiale B(n,p).
(4) Montrer que la méthode issue du Théorème 3.1 s’étend sans difficulté aux lois discrètes à support dénombrable. Donnez un programme d’illustration.
(5) Simuler la loi de Poisson.
(6) Illustrer graphiquement la convergence de la loi binomiale B(n, p) vers la loi de Poisson P(?) lorsque np ? ?.
(7) Illustrer graphiquement la convergence de la loi binomiale B(n, p) vers la loi normale N(np, p(1 ? p)) (théorème de Moïvre-Laplace).
On attribue parfois cette méthode à J. von Neumann. Soit une loi L de densité f sur R, que l’on désire simuler. Supposons que f est continue et à support compact [a,b]. Le graphe de f est donc contenu dans un rectangle [a,b] × [0,M]. On considère alors des vecteurs aléatoires indépendants ((Xn,Yn))n de loi uniforme sur [a,b]×[0,M]. On définit ensuite T comme le plus petit n tel que f(Xn) 6 Yn. Alors on montre que T est fini presque sûrement et que XT a pour loi L.
Exercice 3.3 Simuler la loi uniforme sur le disque unité dans R2 par la méthode du rejet.
On utilise le résultat suivant
Proposition 3.4 SoientU1etU2deux v.a.r. indépendantes de loi uniforme U([0, 1]).
Posons
,
AlorsX1etX2sont deux v.a.r. indépendantes de loi normale centrée réduite N(0, 1).
Exercice 3.5
(1) En utilisant la Proposition 3.4, écrire une fonction de paramètre n simulant n réalisations de la loi normale centrée réduite N(0, 1).
(2) Écrire une fonction de paramètres n, m et ? simulant n réalisations d’une loi normale N(m,?).
(3) Sachant qu’une v.a. suivant une loi du ?2 à ? degrés de liberté peut être vu comme la somme de ? v.a. iid de loi normale N(0, 1), simuler n réalisations d’une v.a. de loi ?2(?). Faire varier ?. Qu’obtient-on pour ? grand (? > 30).
Proposition 3.6 Soient (U1,U2) un couple de loi uniforme sur le disque unité. Posons
,
AlorsX1etX2sont deux v.a.r. indépendantes de loi normale centrée réduite N(0, 1).
Exercice 3.7 En utilisant la Proposition 3.6, écrire une fonction de paramètre n simulant n réalisations de la loi normale centrée réduite N(0, 1) (Indication : simuler d’abord la loi uniforme sur le disque unité par la méthode du rejet).
Ces méthodes sont à la fois très classiques et très lentes, si bien que dans la pratique, on utilise les méthodes du type "Marsaglia", qui semblent être de loin les plus performantes. Il est préférable d’utiliser la fonction Matlab randn qui est déjà toute programmée et qui fait appel à une méthode plus efficace. En effet, la méthode de simulation de Box-Müller est coûteuse en temps de calcul car? elle nécessite l’évaluation des fonctions log, et cos. On peut se passer du cos en faisant appel à la méthode polaire, qui est celle utilisée dans Matlab 4 pour implémenter la fonction randn. Elle demeure cependant assez lente, en raison à?
la fois des 21% de rejet et de l’évaluation des fonctions log et .
Les versions 5 et 6 de Matlab utilisent plutôt une méthode due à Marsaglia, appelée "algorithme du Ziggurat". Elle consiste à utiliser la méthode du rejet conjointement à des tables donnant une approximation de la densité gaussienne par des rectangles. Ces tables étant précalculées, la génération de nombres aléatoires gaussiens est très rapide, quasiment aussi rapide que celle des nombres aléatoires uniformes.
Soit L une loi sur R de fonction de répartition F. Si G désigne l’inverse continue à gauche de F, définie pour tout y dans [0,1] par
G(p) := inf{x : F(x) > p},
alors, pour toute v.a. uniforme U sur [0,1], la v.a. G(U) est de loi L. Par exemple, pour la loi de Cauchy, de densité (?(1+x2))?1 sur R, on a F(x) = arctan(x)/?+1/2 et G(p) = tan(?(p ? 1/2)).
Exercice 3.8 [Simulation de la loi exponentielle] En utilisant le résultat suivant, simuler 100 réalisations de la loi exponentielle.
Proposition 3.9 SoitU une v.a.r. de loi uniforme sur [0,1] et? > 0 un réel. Alors la v.a.r. ???1 log(U) suit la loi exponentielle de paramètre?.
Cette méthode s’applique bien entendu aux lois de Weibull, dont les lois exponentielles ne constituent qu’un cas très particulier. Solution 1
lambda=0.5;clf;hold on; title(’Simulation d’une loi exponentielle’) ylabel(’Densit\’{e}’);xlabel(’Valeurs’); [E,C]=histo(-log(rand(100,1))/lambda,100,0,1) plot(C,lambda*exp(-lambda*C),’k-*’) legend(’Empirique’,’Th\’{e}orique’)
La planche de Galton est une planche verticale sur laquelle sont plantées n lignes de clous disposés en pyramide : 1 clou sur la première ligne, 2 clous sur la deuxième ligne, , n clous sur la n-ième ligne. Des billes sont lâchées l’une après l’autre à partir du sommet. A chaque clou rencontré, la bille a
— une probabilité p = 1/2 de passer à droite (succès), — une probabilité 1 ? p = 1/2 de passer à gauche (échec).
Sous la planche, les billes tombent dans n + 1 compartiments, numérotés de 0 à n de gauche à droite, dans lesquels elles s’empilent.
Soit X la variable aléatoire qui à chaque bille associe le numéro du compartiment atteint. Ainsi, l’événement (X = k) est réalisé quand la bille passe k fois à droite des clous (k succès). Il est évident que la v.a. X suit la loi binomiale B(n,p).
Exercice 3.10
Figure 3.2 – Planche de Galton à 7 lignes
1. Écrire un programme de simulation de la planche de Galton. Ce programmeconsiste à
(a) visualiser l’évolution des effectifs expérimentaux de billes tombées danschaque compartiment (lors de la chute des billes l’une après l’autre),
(b) donner les effectifs théoriques des billes tombées dans chaque compartiment,
(c) comparer graphiquement les piles de billes tombées dans les compartiments avec les piles de billes calculées théoriquement.
2. Expliquer théoriquement la symétrie de la répartition des effectifs dans lescompartiments.
3. Changer la valeur de p et expliquer les résultats.
Pour visualiser un joli programme d’illustration de la planche de Galton, vous pouvez consulter le site suivant. http
Remarque 3.11 Il existe d’autres méthodes pour simuler de façon "concrète" la loi binomiale, notamment la machine de Bittering. On peut également réaliser d’autres expériences pratiques où les probabilités permettent de mesurer certaines
grandeurs (aiguille de Buffon, cercle répétiteur, )
Chapitre 4
Statistique inférentielle
Soit X une variable aléatoire sur une population dont on ne peut faire une étude complète. On suppose que la variable aléatoire X est de moyenne et de variance finie, et on note
E[X] = µ; V [X] = ?2.
Notons que ces paramètres sont inconnues a priori. On se propose d’une part, de donner une estimation ponctuelle de ces paramètres; d’autre part, de donner un intervalle de confiance de µ au niveau de confiance 1 ? ?. Pour cela, nous nous appuyerons sur les deux résultats théoriques suivants.
Théorème 4.1 (Loi faible des grands nombres) Soit (Xn)n?Nune suite de variables aléatoires indépendantes de même loi et d’espéranceµ finie. Alors
en probabilité lorsquen ? +?.
Remarque 4.2 On dit qu’une suite (Sn)n de v.a. converge en probabilité vers µ lorsque
?? > 0 lim P(|Sn ? µ| > ?) = 0.
n?+?
Théorème 4.3 (Théorème central-limite) Soit (Xn)n?Nune suite de variables aléatoires indépendantes de même loi, d’espéranceµ et de variance?2finies. Posons
?n ? N? et .
Alors la suite (Zn)nconverge en loi vers la loi normale N(0,1).
Remarque 4.4 On dit qu’une suite (Zn)n de v.a. converge en loi vers la loi normale lorsque
?a, b ? R tel que a < b, on a
20
Le problème qui se pose est de donner, à partir d’un n-échantillon (x1, ,xn), une estimation des valeurs caractéristiques d’une v.a. définie sur une population. Le tableau suivant donne une première approche du problème.
| valeur théorique | valeur empirique | |
| Moyenne | µ | |
| Variance | ?2 |
Figure 4.1 – Deux estimateurs
Mais ces estimateurs ne sont pas tous les "meilleurs" (en certain sens que nous ne chercherons pas à préciser). Le tableau suivant donne les "meilleurs" estimateurs des paramètres précédents.
| valeur théorique | "meilleur" estimateur | |
| Moyenne | µ | me |
| Variance | ?2 |
Figure 4.2 – "Meilleurs" estimateurs
À tout n-échantillon (x1, ,xn) de moyenne empirique me on associe l’intervalle de confiance au niveau de confiance 1 ? ? de la moyenne µ suivant :
,
où t? est le réel défini par , avec ? la fonction de répartition de la loi normale N(0, 1). Noter que si la variance ?2 est inconnue on la remplace par son "meilleur" estimateur . Remarquons que l’inégalité de Bienaymé-
Tchebychev permet également de trouver des intervalles de confiance, mais de façon bien moins précise que le théorème central-limite.
Remarque 4.5 Le réel t? s’obtient par le code Matlab sqrt(2)*erf(inv(1-a)).
Les objectifs de cette simulation sont
1. montrer expérimentalement le caractère relatif des intervalles de confiances;
2. comprendre la notion de confiance qu’on peut leur attacher.
Exercice 4.6 [Comparaison du pourcentage d’intervalles de confiance contenant effectivement µ avec le seuil de confiance.]
1. Simuler m échantillons de tailles n d’une variable aléatoire suivant une loi de votre choix.
2. Pour un niveau de confiance 1?? fixé, déterminer m intervalles de confiance de µ à partir des m échantillons.
3. Calculer le pourcentage d’intervalles contenant réellement µ.
4. Donner une illustration graphique de ces résultats.
Exercice 4.7 Un prélèvement de 100 échantillons de même effectif 40 pièces a montré l’existence d’un certain nombre de pièces défectueuses. Les résultats obtenus sont :
| Nb de pièces défectueuses | Effectifs observés |
| 0 | 28 |
| 1 | 40 |
| 2 | 21 |
| 3 | 7 |
| 4 ou plus | 4 |
1. Tracer l’histogramme de la répartition de l’échantillon.
2. Quel type de loi cela inspire-t-il pour ajuster le phénomène?
3. Déterminer les paramètres de la loi supposée et tester l’ajustement.
Exercice 4.8 On observe les nombres de vols d’objets dans un immeuble de 4 étages. Les données sont les suivantes.
1. Est-il raisonnable de supposer que les vols suivent une loi de Poisson?
| Étage de l’immeuble | Nombre d’objets volés |
| 0 | 545 |
| 1 | 325 |
| 2 | 110 |
| 3 | 15 |
| 4 | 5 |
2. Peut-on supposer que le nombre d’objets volés suit une loi de Poisson deparamètre ? = 1?
3. Si la réponse est oui pour la question 1., donner une estimation expérimentaledu paramètre ??
Exercice 4.9 A l’école X, il est connu que grâce aux différentes notes (dont on fait une moyenne pondérée) on détermine pour chaque étudiant sa note finale. Sur l’ensemble de la promotion, la moyenne est de 10 sur 20 (précisément); mais cela cache bien des disparités :
| Note finale | Nombre d’étudiants |
| Moins de 5 | 2 |
| De 5 à moins de 7.5 | 21 |
| De 7.5 à moins de 9.75 | 36 |
| De 9.75 à moins de 10.25 | 30 |
| De 10.25 à moins de 12.5 | 37 |
| De 12.5 à moins de 15 | 23 |
| 15 et plus | 1 |
1. Tracer l’histogramme de répartition des étudiants. Quel type de loi celainspire-t-il pour ajuster le phénomène?
2. Droite de Henry : Le principe de la droite de Henry repose sur l’hypothèse, supposée vérifiée, de normalité d’un phénomène. S’il est acquis que la v.a. X suit une loi normale N(m,?), la v.a. T = (X ? m)/? suit la loi N(0,1). Ceci s’apparente à une équation de droite de pente 1/? et d’ordonnée à l’origine ?m/?. Lorsqu’on n’est pas sûr de la normalité du phénomène, il faut représenter graphiquement les points (xi,ti) et constater s’ils sont ou non alignés. Les valeurs de ti sont obtenues à partir des valeurs de ?(ti) assimilées aux valeurs de fréquences cumulées croissantes.
(a) Compléter le tableau suivant. Pour déterminer les ti ?à partir de pi =
?(ti) on utilise la fonction inverse qu’on a déjà vu ti = 2erfinf(2pi?1)
| Note | Nb étudiants | Effectifs cumulés | Fréquences cumulées | ti |
| Moins de 5 | 2 | |||
| De 5 à 7.5 | 21 | |||
| De 7.5 à 9.75 | 36 | |||
| De 9.75 à 10.25 | 30 | |||
| De 10.25 à 12.5 | 37 | |||
| De 12.5 à 15 | 23 | |||
| 15 et plus | 1 |
(b) Tracer le nuage de points (xi,ti) : (5,t1),(6.25,t2), ,(15,t6) (on prend les milieux des classes sauf pour les extrémités).
(c) Déterminer la droite de Henry qui n’est autre que la droite de régressionde la série (xi,tj) (Pour cela, on peut utiliser la commande Matlab polyfit; voir help). En déduire les caractéristiques de la loi normale supposée : m = E(X) et ? = ?(X) (graphiquement ou (et) par calcul).
3. Tester à l’aide du ?2 la validité de cette hypothèse (loi normale) en utilisant la moyenne et l’écart type obtenus ci-dessus. N’oublier pas de déterminer le degré de liberté de la loi ?2 et de regrouper les classes de faible effectif (moins de 5).
