Introduction a Matlab et Octave`
B. Torrésani
Université de Provence
Master Mathématiques et Applications
Spécialité GSI
Année Universitaire 2009-10, premier semestre
2
| 1 |
Quelques lignes d’introduction.
Matlab est un logiciel (commercial) de calcul numérique, de visualisation et de programmation très performant et convivial2. Le nom de Matlab vient de MATrix LABoratory, les éléments de données de base manipulés par Matlab étant des matrices (pouvant bien évidemment se réduire à des vecteurs et des scalaires) qui ne nécessitent ni dimensionnement ni déclaration de type. Contrairement aux langages de programmation classiques (scalaires et à compiler), les opérateurs et fonctions Matlab permettent de manipuler directement et interactivement ces données matricielles, rendant ainsi Matlab particulièrement efficace en calcul numérique, analyse et visualisation de données en particulier.
Mais Matlab est aussi un environnement de développement (”progiciel”) à part entière : son langage d’assez haut niveau, doté notamment de structures de contrôles, fonctions d’entrée-sortie et de visualisation 2D et 3D, éditeur/debugger, outils de construction d’interface utilisateur graphique (GUI) permet à l’utilisateur d’élaborer ses propres fonctions ainsi que de véritables programmes (”M-files”) appelés scripts vu le caractère interprété de ce langage.
Matlab est disponible sur tous les systèmes d’exploitation standards (Windows, Unix/Linux, MacOS X ) et son architecture est relativement ouverte. Le champ d’utilisation de Matlab peut être étendu aux systèmes non linéaires et aux problèmes associés de simulation avec le produit complémentaire SIMULINK. Les capacités de Matlab peuvent en outre être enrichies par des fonctions spécialisées regroupées au sein de dizaines de ”toolboxes” (boˆ?tes à outils qui sont des collections de ”M-files”) couvrant des domaines très variés tels que :
– Traitement de signaux (et du son en particulier)
– Traitement d’image, cartographie
– Analyse de données
– Statistiques
– Finance et mathématiques financières
– Mathématiques symboliques (accès au noyau Maple V)
– Analyse numérique (accès aux routines NAG)
Une interface de programmation applicative (API) rend finalement possible l’interaction entre Matlab et les environnements de développement classiques (exécution de routines C ou Fortran depuis Matlab , ou accès aux fonctions Matlab depuis des programmes C ou Fortran). Matlab permet en outre de déployer de véritables applications à l’aide des outils de conversion optionnels suivants :
– Matlab ? code C/C++, avec le Matlab Compiler
– Matlab ? Excel add-ins, avec le Matlab Excel Builder
– Matlab ? objets COM Windows, avec le Matlab COM Builder
Toutes ces caractéristiques font aujourd’hui de Matlab un standard incontournable en milieu académique, dans les différents domaines de l’ingénieur et la recherche scientifique.
GNU Octave , et autres alternatives à Matlab : Matlab est un logiciel commercial qui couˆte relativement cher, même au tarif académique. Cependant, il existe des logiciels libres et open-source analogues voire compatibles avec Matlab , donc gratuits, également multi-plateformes :
– GNU Octave : logiciel offrant la meilleure compatibilité par rapport à Matlab (qualifiable de ”clone Matlab ”, surtout depuis la version Octave 2.9/3.0 et avec les packages du repository Octave
-Forge).
– FreeMat : logiciel libre compatible avec Matlab et Octave , déjà très abouti, avec IDE comprenant : editor/debugger, history, workspace tool, path tool, file browser, 2D/3D graphics
– Scilab : logiciel libre ”analogue” à Matlab et Octave en terme de fonctionnalités, très abouti, plus jeune que Octave mais beaucoup moins compatible avec Matlab (fonctions et syntaxe différentes )
3
et ne permettant donc pas une réutilisation aisée de scripts.
– R : un logiciel libre, davantage orienté vers les statistiques.
– NumPy : Numerical Python, une extension de Python tournée vers le calcul scientifique –
Dans ce cours, on se focalisera sur l’utilisation de Matlab et Octave , qui sont en très grande partie compatibles. Les TP étant effectués avec Matlab , les instructions données sont des instruction Matlab; ceci dit, elles sont également opérationnelles sous Octave .
Octave est téléchargeable librement sur le site
http
En utilisant la commande mkdir, commencer par créer un répertoire, par exemple MonMatlab (ou autre), dans lequel vous stockerez votre travail. Placez vous dans ce répertoire (en utilisant la commande cd).
Téléchargez dans ce répertoire les fichiers se trouvant sur http http 2.wav
En ligne de commande (dans une fenêtre “terminal”), taper matlab & (rappel : le symbole “&” signifie que le logiciel est lancé en tâche de fond). Normalement, l’interface de matlab s’affiche. Examiner les différentes composantes de celle-ci : fenêtre de commande, espace de travail, répertoire courant, historique des commandes,
| Exercice 1.1 Dans la fenêtre de commande, taper >> a = 1+1 puis >> b = 1+2; La première ligne calcule 1+1, l’affecte dans la variable a, et affiche le résultat. La seconde calcule 1+2, l’affecte dans la variable b, et n’affiche pas le résultat à cause de la présence du “;” en fin de ligne. Pour afficher de nouveau le résultat, taper par exemple >> a Noter que les variables ainsi créées sont également visibles dans l’onglet “Workspace”. Noter aussi que lorsque le résultat du calcul n’est pas affecté à une variable, il est affecté à une variable par défaut notée ans. |
La commande qui sauve : help, accessible soit dans la fenêtre de commande
>> help doc
soit sous format HTML, grâce à la commande doc, soit directement à partir du menu de l’interface Matlab (en haut à droite). Il est aussi possible de rechercher une fonction à partir d’un mot clé, grâce à la commande lookfor. Par exemple >> lookfor(’variance’)
Exercice 1.2 Se documenter sur les variables who, whos, et clear, et les essayer.
Quitter Matlab :
>> exit
1.1.2 Gestion des répertoires, fichiers,
Relancer Matlab depuis votre répertoire de travail. Dans la fenêtre de commande, tester les instructions pwd et ls. Les informations obtenues via ces commandes sont également accessibles dans l’onglet “Current Directory”.
1.2. VECTEURS, MATRICES,
Matlab permet de sauvegarder des variables dans des fichiers de données, dont le nom porte l’extension “.mat”, et de charger de tels fichiers. Se documenter sur les instructions save et load. Pour vérifier, interpréter puis effectuer les commandes suivantes :
>> A = 1 :10;
>> B = 3;
>> who >> save
>> ls
>> save A
>> clear all
>> load
>> who >> load
>> ls
Ceci vous permettra de sauvegarder des résultats d’une séance à l’autre.
Il existe également des fonctions de lecture/écriture de plus bas niveau, plus proche du langage C (fopen, fclose, fprintf, fscanf, fwrite, ). Il est également possible de lire (et écrire) des fichiers plus spécifiques, par exemple des fichiers son aux formats .wav ou .au (fonctions wavread, auread, ou des fichiers image au format .jpg (fonction readjpg, ). On y reviendra plus tard.
1.2 Vecteurs, matrices,
L’élément de base sous Matlab est la matrice (le vecteur, ligne ou colonne, étant considéré comme une matrice particulière). On peut accéder aux dimensions des matrices en utilisant l’instruction size, par exemple
>> A = rand(2,3);
>> tmp = size(A);
>> tmp(1)
>> tmp(2)
Les règles de multiplication sont par défaut les règles de multiplication matricielle.
Pour générer “manuellement” une matrice, il faut savoir que le séparateur de colonnes est la virgule (ou l’espace), et le séparateur de lignes est le “point virgule”. Pour vérifier :
>> M = [1,2,3;4,5,6;7,8,9]
>> M = [1 2 3;4 5 6;7 8 9]
Dans les deux expressions suivantes
>> x = 1:100
>> x = 1:10:100
on génère deux vecteurs (ligne) constitués des 100 premiers entiers positifs dans le premier cas, et des entiers de 1 à 100 par pas de 10 dans le second cas.
L’opération “ ’ ” représente la conjugaison Hermitienne des matrices (la transposition suivie de la conjugaison complexe). La transposition des matrices s’obtient en associant la conjugaison Hermitienne à la conjugaison complexe. Pour vérifier :
>> A=[1,i;2i,2;3i,3]
>> AH = A’
>> AT = conj(A’)
Exercice 1.3 Après en avoir vérifié la syntaxe à l’aide du help, tester les fonctions ones, zeros, eye, rand et randn.
| Exercice 1.4 1. Générer deux matrices “multipliables” A et B, et les multiplier. 2. Générer deux matrices “non-multipliables” A et B, et essayer de les multiplier. |
Matlab implémente également un autre type d’opération sur les matrices (et les vecteurs) : les opérations “point par point”, dites “pointées”. Ces opérations ne sont possibles qu’entre matrices de même taille. Tester cette opération, par exemple
>> A = [1,2,3;4,5,6];
>> B = [4,5,6;7,8,9];
>> A.*B
>> A./B ou encore
>> x = 1:5
>> x.^ 2
Matlab a été originellement créé pour l’algèbre linéaire (matlab signifie « matrix laboratory »). Ainsi, l’objet de base est la matrice, et de nombreuses opérations matricielles sont déjà implémentées et optimisées.
| Exercice 1.5 1. Se documenter sur les instructions inv et det. Générer deux matrices carrées de même taille A et B, vérifier qu’elles sont inversibles. Pour vérifier : >> A/B >> A*inv(B) >> A\B >> inv(A)*B
2. Générer deux vecteurs de même taille (ligne ou colonne), et effectuer leur produit scalaire (ou leur produit Hermitien s’ils sont complexes). On pourra procéder de deux fac¸ons : – en utilisant l’instruction sum et une multiplication pointée – en utilisant une transposition et un produit matriciel Dans les deux cas, ¸ca ne prend qu’une ligne. 3. Utiliser les opérations matricielles pour résoudre le système linéaire suivant (après avoir testé l’existence de solutions) : 4. Se documenter sur les commande null et eig. Utiliser ces deux commandes pour résoudre le système ci-dessus en utilisant la diagonalisation des matrices. |
La plupart des fonctions usuelles (trigonométriques, exponentielles, ) existent également sous forme vectorielle et matricielle. Par exemple, tester la suite d’expressions
>> xxx = linspace(1,10,100);
>> C = cos(2*pi*xxx);
On pourra utiliser l’instruction
>> plot(C)
pour visualiser le résultat.
De même, pour les matrices, l’expression exp(A) renvoit la matrice dont les éléments sont les exponentielles des éléments de matrice de A : c’est une opération pointée,
>> A = eye(2)
>> B = exp(A)
1.3. GRAPHIQUES EL´ EMENTAIRES´
à ne pas confondre avec l’exponentielle matricielle des matrices carrées
pour laquelle est définie la fonction expm (se documenter)
>> B = expm(A)
Matlab permet également de manipuler des tableaux à plus d’une ou deux entrées, comme par exemple des tableaux “cubiques”, ou plus. Par exemple >> A = randn(2,2,2)
Matlab possède d’intéressantes possibilités graphiques. Il est possible d’ouvrir une (nouvelle) fenêtre graphique grâce à l’instruction figure. Les fenêtres sont numérotées. On se place dans une fenêtre donnée (par exemple la quatrième) en utilisant de nouveau l’instruction figure (dans ce cas, en effectuant figure(4). Il est également possible de tracer plusieurs graphiques dans la même figure grâce à l’instruction subplot (se documenter).
La fonction de base est la fonction plot. Par exemple, on pourra (après les avoir interprêtées) essayer les instructions
>> x = linspace(0,10,100);
>> f = exp(-x/10).* cos(2*pi*x);
>> g = f + randn(size(x))/3;
>> plot(x,g)
>> hold on
>> plot(x,f,’r’)
>> xlabel(’Abscisses’) >> ylabel(’Ordonnées’)
>> title(’Ceci est un titre’)
On a utilisé ici l’instruction hold on, qui permet de conserver un tracé et en superposer un autre (par défaut, un nouveau tracé efface le précédent, ce qui correspond à l’instruction hold off), et l’argument ’r’ de la fonction plot, qui permet de tracer en rouge au lieu du bleu (couleur par défaut). D’autres couleurs sont possibles, se documenter. Se documenter sur les fonctions xlabel et ylabel. Le résultat de ces quelques lignes se trouve en Fig. 1.1 (figure de gauche).
Il existe également de multiples fonctions graphiques avec lesquelles l’on peut jouer pour améliorer l’aspect du graphe généré. Par exemple, en rempla¸cant les trois dernières lignes ci-dessus par
>> xlabel(’Abscisses’),’fontsize’,20)
>> ylabel(’Ordonnées’,’fontsize’,20)
>> title(’Ceci est un titre’,’fontsize’,24)
on modifie les tailles des caractères dans les légendes et le titre (voir tracé de droite). Il est aussi possible d’effectuer ces modifications directement sur la fenêtre graphique, en l’éditant.
L’instruction line permet de tracer une ligne dans la fenêtre graphique. La syntaxe est la suivante : line([x1,x2, xn],[y1,y2, yn]) trace un polygone à n côtés, dont les sommets ont pour coordonnées (x1,y1), Par exemple, les lignes
>> x = [1,1.5,2,1.5,1];
>> y = [1,0,1,0.5,1];
>> line(x,y)
génèrent un “chevron” passant par les points de coordonnées (1,1), (1.5,0), (2,1), (1.5,0.5) et (1,1) (voir Figure 1.2
Fig. 1.1 – Exemple de représentation graphique
Fig. 1.2 – Un “chevron” tracé avec l’instruction line.
1.3. GRAPHIQUES EL´ EMENTAIRES´
Fig. 1.3 – Différentes visualisations graphiques d’une fonction de deux variables
Exercice 1.6 Pour vérifier, tracer un losange centré sur un point donné.
– Les instructions image, imagesc permettent de représenter une matrice sous forme d’une image. Se documenter, et imager une matrice aléatoire
– Il existe des instructions permettant de faire du graphique tridimensionnel (plot3, mesh, surf, ), on les verra plus tard.
Les quelques lignes suivantes permettent de donner différentes représentations graphiques d’une fonction de 2 variables (la fonction peaks, un classique de Matlab )
>> Z = peaks(64);
>> imagesc(Z); axis xy (représentation sous forme d’image)
>> contour(Z,32); (32 lignes de niveau)
>> surf(Z);
>> surfl(Z); shading interp; colormap(’copper’);
Les résultats se troucvent dans la Fig. 1.3 (ou` on a également joué avec la colormap pour la dernière image).
Une fois un graphique tracé, il est toujours possible de le modifier en utilisant les outils graphiques de Matlab . On peut par exemple modifier les couleurs des courbes, les fontes des labels, titres et légendes et leur taille,
Il est également possible de sauvegarder un graphique sous divers formats. Le meilleur est a priori le format “.fig”, qui permet d’éditer le graphique lors de sessions futures. Il est également possible de sauvegarder aux formats ps (postscript), eps (postscript encapsulé, c’est à dire sans information de page), jpg (format JPEG, compression avec perte),
Un script est un fichier texte, contenant une série d’instructions Matlab . Ces instructions sont exécutées lorsque le script est exécuté. Les variables crées lors de l’exécution du script restent en mémoire après exécution. L’extension ’standard’ d’un fichier de script est ’.m’.
Une fonction diffère d’un script par l’existence de variables d’entrée et de sortie. Les variables introduites dans la fonction, autres que les variables de sortie, sont effacées après exécution de la fonction. La syntaxe est la suivante :
function [variables sortie] = mafonction(variables entrée) Cette fonction sera sauvée dans le fichier “mafonction.m”.
Par exemple, la fonction ci-dessousprend comme variable d’entrée un vecteur de longueur quelconque, et retourne sa moyenne et sa variance
function [m,v] = moyennevariance(X) % function : moyennevariance % usage : [m,v] = moyennevariance(X) % Variable d’entrée :
% X : vecteur
% Variables de sortie :
% m : moyenne (scalaire) % v : variance (scalaire) longueur = length(X); m = sum(X)/longueur; v = sum((X-m).^2)/(longueur-1);
Dans cette fonction, les lignes précédées d’un signe % sont des lignes de commentaires. L’intérêt d’insérer de telles lignes dans une fonction est double : d’une part, elles permettent de se souvenir des significations des variables, et d’autre part, ce commentaire apparaˆ?t lorsque l’on effectue >> help(moyennevariance)
L’interface utilisateur Matlab contient un éditeur de texte permettant d’écrire des scripts et des fonctions.
| Exercice 1.7 Ecrire une fonction de visualisation de la transformée de Fourier discrète (TFD) N?1 Xˆk = X Xne?2i?kn/N n=0 d’un vecteur de taille fixée N. On utilisera pour cela la fonction fft, qui calcule la TFD. La fonction aura pour argument d’entrée un vecteur, calculera sa transformée de Fourier, et tracera le module (fonction abs) de sa transformée de Fourier. Elle retournera la transformée de Fourier. |
1.5. DEVOIR
Les exercices suivants sont à faire individuellement, et à rendre sous forme d’un fichier “archive”, à envoyer à l’adresse . Pour cela, sous unix/linux, créer un répertoire nommé nom prenom tp1, et copier tous les fichiers *.* dans ce répertoire. Ensuite, dans le répertoire supérieur, effectuer tar -zcvf nom prenom nom prenom tp1/*. Sous Windows, on pourra utiliser un logiciel de compactage (winzip, winrar,7zip, ) pour créer l’archive.
Exercice 1.8 Ecrire une fonction HistConv.m, prenant pour argument un entier positif N, qui génère N nombres pseudo-aléatoires uniformément distribués, trace l’histogramme correspondant et affiche la moyenne et la variance de la population ainsi créée. On pourra utiliser la fonction lookfor pour trouver les fonctions Matlab nécessaires.
| Exercice 1.9 Ecrire une fonction calculecorr.m qui prend comme variables d’entrée deux vecteurs de même dimension, calcule le coefficient de corrélation entre ces vecteurs et l’affiche dans la fenêtre de commande. On rappelle que le coefficient de corrélation est donné par , ou` X est la moyenne de X, ?X est l’écart-type de X, et de même pour Y . On pourra tester cette fonction en utilisant des vecteurs aléatoires (randn). |
| Exercice 1.10 Ecrire une fonction resolsystlin.m, prenant en entrée un vecteur à 4 composantes y = (y1,y2,y3,y4), et retournant en sortie la solution de l’équation ? x1 + 2x2 + 3x3 + x4 = y1 ??? 4x1 + 5x2 + 6x3 + x4 = y2 7x1 + 8x2 + 10x3 + x4 = y3??? 9x1 + 11x2 + 12x3 + x4 = y4 |
Exercice 1.11 Ecrire une fonction tracepolygone.m prenant en entrée un entier positif N, et dessinant dans la fenêtre graphique un polygone régulier à N côtés.
| Exercice 1.12 Ecrire une fonction proj2.m prenant en entrée une matrice A réelle symétrique N ×N et un vecteur V à N lignes, et effectuant les opérations suivantes : – Diagonalisation de A, en utilisant la fonction eig. – Sélection des deux valeurs propres les plus grandes, en utilisant la fonction de tri sort – Projection du vecteur V sur le plan engendré par les deux vecteurs propres correspondants. Les variables de sortie seront la projection de V sur le plan, ainsi que les composantes de V sur les deux vecteurs propres. |
12 1. INTRODUCTION
| 2 |
Pour commencer, deux définitions :
– Algorithme : ensemble de règles précises, qui permettent d’obtenir un résultat à partir d’opérations
élémentaires.
– Algorithme numérique : Les données et les résultats sont des nombres.
Cela va sans dire, mais mieux encore en le disant : un programme ne se con¸coit pas sans une analyse préalable du problème posé, c’est à dire une analyse des points suivants :
– Quel est le problème posé? quel est le résultat attendu?
– Quelles sont les données de départ?
– Quelles sont les données que doit retourner le programme? les représentations graphiques?
– Quel est l’algorithme de calcul?
– Comment “optimiser” le calcul? Comment le décomposer en étapes indépendantes?
Une fois cette analyse faite et écrite, la phase de codage proprement dite peut commencer.
Un programme Matlab ne diffère pas fondamentalement d’un programme dans un langage classique.
La seule différence importante tient dans le fait que la majorité des opérations par défaut dans Matlab
étant vectorisées, il est généralement préférable d’exploiter cette caractéristique, par exemple en évitant des boucles dans les programmes.
Les connecteurs de relation (qui retournent soit 0 soit 1 suivant que l’expression correspondante est vraie ou fausse), sont résumés dans le tableau suivant.
| Connecteur | Signification |
>> a == b
retourne 1 si le contenu de la variable a est égal au contenu de la variable b, et 0 sinon.
Exercice 2.1 Ecrire une fonction monmax.m, prenant en entrée deux nombres, et retournant le plus grand des deux (sans utiliser la fonction max, sinon c’est de la triche).
| Exercice 2.2 Interpréter puis tester (ou l’inverse ) les deux lignes suivantes >> x = randn(10,1); >> y = (abs(x)>=1) On essaiera aussi la fonction find : >> z = find(abs(x)>=1) |
13
Quant aux connecteurs logiques (opérant en binaire), les plus importants se trouvent dans le tableau suivant
| Connecteur | Signification |
| & | et |
| | | ou |
| xor | ou exclusif |
| ? | non |
Exercice 2.3 tester ces instructions, après en avoir vérifié la syntaxe à l’aide du help.
Exercice 2.4 Construire les tables des trois connecteurs &, | et xor.
| Exercice 2.5 Interpréter les commandes : >> x=3; >> y=2; >> (x==3)&(y==2) >> (x==3)|(y==0) Dans chaque cas, quel sera le résultat? |
Les boucles sont l’un des moyens de contrôle de l’exécution du programme. La syntaxe est identique à ce qu’elle est dans la plupart des langages.
– Boucles for : on peut créer une boucle en utilisant for end.
>> k2 = zeros(1,5);
>> for k = 1 :5
>> k2(k) = k^2;
>> end
Il est important de signaler que cette boucle est bien moins efficace que la simple ligne
>> k2 = (1 :5).^2
On peut aussi réaliser des boucles for imbriquées.
– Boucle while : On peut créer une boucle en utilisant while end.
>> k = 1;
>> k2 = zeros(1,5);
>> while tmp <= 25
>> tmp = k^2;
>> k2(k) = tmp;
>> k = k+1;
>> end
Exercice 2.6 Ecrire deux fonctions fact1.m et fact2.m calculant la factorielle d’un entier, en utilisant une boucle for puis une boucle while. On pourra utiliser la fonction factorial pour vérifier le résultat.
La construction switch-case permet d’énumérer un certain nombre de cas de figures. Par exemple,
>> x = ceil(10*rand);
>> switch x
>> case {1,2}
>> disp(’Probabilite = 20%’) >> case {3,4,5}
>> disp(’Probabilite = 30%’) >> otherwise
>> disp(’Probabilite = 50%’)
>> end
2.1. QUELQUES EL´ EMENTS DE PROGRAMMATION SOUS´ MATLAB ET OCTAVE 15
Exercice 2.7 Exécuter et interpréter ce script (qui utilise la fonction ceil qui arrondit un réel à l’entier le plus proche, et la fonction disp, qui affiche du texte à l’écran).
Remarque : l’exemple précédent permet de voir une utilisation de chaˆ?ne de caractères : l’expression ’Probabilité = 20%’ est une chaˆ?ne de caractères.
La suite d’instructions if elseif else permet de choisir entre plusieurs options. La syntaxe est là encore similaire à la syntaxe “classique” : par exemple
>> x = randn(1);
>> if abs(x) <= 1
>> y = 0;
>> elseif abs(x) > 5
>> y = x^ 2;
>> else
>> y = x;
>> end
Exercice 2.8 Que fait cette suite d’instructions?
| Exercice 2.9 En utilisant la fonction rand, – Générer des points uniformément distribués dans un carré de côté 2 et de centre 0. – Tester la distance à l’origine du point généré. – Effectuer ces opérations à l’intérieur d’une boucle, stocker dans un tableau à 2 colonnes les coordonnées des points dont la distance à l’origine est inférieure à 1, et tracer ces derniers. – Pour compléter, on pourra également tracer les histogrammes des coordonnées Cartésiennes des points ainsi générés, ainsi que de leurs coordonnées polaires. |
Remarque : Dans l’exercice précédent, on pourra remarquer que l’on a généré des points pseudoaléatoires, distribués sur le disque de rayon 1 avec une distribution uniforme, c’est à dire une densité de probabilités de la forme
si px2 + y2 ? 1
0 sinon
En notant R et ? les variables aléatoires correspondant aux coordonnées radiale et angulaire des points générés, il est facile de calculer leur densité de probabilités, et de les comparer aux résultats obtenus.
Matlab offre de nombreux moyens de communication entre un programme et l’utilisateur, par exemple :
– Comme on l’a déjà vu, on peut afficher un message, une valeur à l’écran avec l’instruction disp :
>> disp(’Ceci est un test’)
– On peut entrer une valeur avec l’instruction input :
>> x = input(’Valeur de x = ’)
Matlab attend alors qu’un nombre soit entré sur le clavier.
– Il est possible d’afficher un message d’erreur, en utilisant error.
Il faut remarquer qu’il existe des fonctions plus sophistiquées, telles que inputdlg, ou errordlg, sur lesquelles on peut se documenter grâce à l’aide. Ces fonctions ouvrent une fenêtre de dialogue.
Matlab permet également d’effectuer des opérations d’entrée-sortie avec des fichiers (et la fenêtre de commande), en s’inspirant de la syntaxe de la programmation en langage C. Par exemple, la commande
>> x = pi;
>> fprintf(’Le resultat du calcul est x=% f\n’,x)
affiche la valeur de x à la fin de la phrase. Le symbole %f est le format dans lequel sera affichée la valeur de x (ici, float), et le symbole \ n est le “retour charriot”.
Les entrées-sorties avec des fichiers s’effectuent avec une syntaxe similaire, en n’oubliant toutefois pas d’ouvrir le fichier avant de lire ou écrire (instruction fopen), et de le refermer ensuite (instruction close). fopen affecte un identifiant au fichier ouvert, identifiant qu’il faut préciser pour lire ou écrire. Se documenter sur ces instructions, au besoin en utilisant l’exercice qui suit.
| Exercice 2.10 En vous aidant du help, interpréter puis effectuer les lignes de code suivantes : >> fp1 = fopen(’’,’w’);
>> x = pi; >> fprintf(fp1,’%f,x); >> fclose(fp1); >> fp2 = fopen(’’,’r’); >> y = fscanf(fp2,’%f’); >> fclose(fp2); >> fprintf(’y = %e\n’,y); |
Exercice 2.11 Compléter votre fonction factorielle en lui faisant vérifier que la variable d’entrée est un entier positif, et la faisant sortir avec un message d’erreur si tel n’est pas le cas. Pour tester que l’argument est un entier, on pourra utiliser la fonction floor. Pour tester le signe, utiliser sign.
On s’intéresse à un marcheur, se promenant sur l’axe réel, et se dépla¸cant dans une direction ou l’autre par pas entiers. On suppose qu’à chaque déplacement n, il a une probabilité p d’aller vers la droite (?n = 1), et une probabilité q = 1 ? p de se diriger vers la gauche (?n = ?1). Sa position à l’instant N est donc
.
| Exercice 2.12 – Ecrire une fonction hasardbinaire.m, prenant en entrée une probabilité p ? [0,1], et retournant en sortie un nombre pseudoaléatoire valant 1 avec probabilité p et ?1 avec probabilité 1 ? p. – Ecrire une fonction marchehasard1.m, prenant comme variables d’entrée le nombre de pas de temps N considéré et la probabilité p, et retournant comme variable de sortie la position XN, et trac¸ant optionnellement la trajectoire du marcheur, c’est à dire Xn en fonction de n. – Dans une fonction testmarchehasard1.m, faire une boucle appelant marchehasard1 un grand nombre M de fois (avec ) pour les mêmes valeurs de p et N, et tracer l’histogramme (en utilisant hist) des valeurs de X ainsi obtenues. Il sera utile de faire plusieurs simulations pour pouvoir interpréter le résultat. On pourra interpréter le résultat en utilisant le théorème de la limite centrale. |
Le but est ici de simuler numériquement une suite de N nombres pseudo-aléatoires, prenant des valeurs x1, xK, avec probabilités p1, pK fixées.
2.4. DEVOIR 17
Fig. 2.1 – Deux trajectoires de marche au hasard 1D : position du marcheur en fonction du temps (10 000 pas). Haut : p = 0.5; bas : p = 0.48.
| Exercice 2.13 Générer pour cela une fonction VAdiscrete.m, prenant comme variables d’entrée l’ensemble des valeurs possibles X = {x1, xK} et le vecteur de probabilités correspondant P = {p1, pK}, ainsi que le nombre N de valeurs pseudo-aléatoires désirées. La variable de sortie sera la suite de N valeurs ainsi générées. On s’appuiera sur le générateur de nombres pseudo-aléatoires rand, qui fournit des nombres pseudo-aléatoires U uniformément distribués entre 0 et 1, et on remarquera que pk = P(U ? [pk?1,pk?1 + pk[ (avec la convention p0 = 0). La fonction effectuera les opérations suivantes : – Vérification que la somme des probabilités vaut 1 – Vérification que X et P sont de même taille – Génération des nombres pseudo-aléatoires. – Calcul des fréquences des nombres obtenus, et comparaison (par exemple, par l’erreur relative) avec les probabilités théoriques. Pour ce qui concerne le dernier point, on pourra utiliser l’instruction >> sum(Z == x)/length(Z) ou` Z est le nombre généré, et x une valeur donnée à condition de la comprendre. |
Les exercices suivants sont à faire par binôme, et à rendre sous forme d’un fichier “archive”, à en-
voyer à l’adresse . Pour cela, sous unix/linux, créer un répertoire nommé nom prenom tp2, et copier tous les fichiers *.* dans ce répertoire. Ensuite, dans le répertoire supérieur, effectuer tar -zcvf nom prenom nom prenom tp2/*. Alternativement, en utilisant le gestionnaire de fenêtres (généralement gnome ou kde), cliquer “droit” sur le répertoire et choisir un format d’archivage et de compactage (de préférence, .zip, ou (ou .tgz, qui est le même format).
Sous Windows, on pourra utiliser un logiciel de compactage (winzip, winrar,7zip, ) pour créer l’archive.
Exercice 2.14 Implémenter une fonction matmult.m prenant en entrée deux matrices A et B telles que le produit matriciel AB soit bien défini, testant leurs dimensions (et retournant une erreur si AB n’existe pas) et effectuant le produit matriciel en utilisant des boucles.
Exercice 2.15 Implémenter une fonction permettant de comparer les performances de la fonction précédente et de la multiplication matricielle de Matlab . Cette fonction génèrera des matrices carrées aléatoires de tailles croissantes, les multipliera en utilisant les deux méthodes, et enregistrera le temps de calcul dans les deux cas (en utilisant l’instruction cputime). Les temps de calcul seront stockés dans deux tableaux, et tracés dans une fenêtre graphique.
Exercice 2.16 (optionnel) Finaliser la fonction testmarchehasard1.m cidessus, en superposant à l’histogramme obtenu une Gaussienne de moyenne et variance égales aux valeurs prédites par le théorème de la limite centrale. Alternativement, utiliser la fonction subplot pour partitionner la fenêtre graphique en deux sous-fenêtres, tracer l’histogramme à gauche, et la densité théorique à droite.
Exercice 2.17 On s’intéresse maintenant à une marche aléatoire dans le plan. Partant du point de coordonnées X0 = (0,0), on génère une suite de points dans le plan de la fa¸con suivante : Xn+1 = Xn+ (u,v) ou` (u,v) est un vecteur pseudo-aléatoire valant (1,0), (?1,0), (0,1) et (0,?1) avec probabilités égales à 1/4.
Ecrire une fonction marchehasard2.m, qui génère une telle trajectoire. On
écrira tout d’abord une fonction générant les vecteurs (u,v) ci-dessus. Cette fonction sera appelée par marchehasard2 à l’intérieur d’une boucle.
Exercice 2.18 (optionnel) On pourra également compléter cette fonction en ajoutant d’autres affichages graphiques, par exemple en utilisant la fonction plot3. On peut également utiliser la fonction subplot pour partitionner la fenêtre graphique en deux sous-fenêtres, et tracer la trajectoire en 2D à gauche, et la trajectoire 3D à droite ou faire plein d’autres choses
Fig. 2.2 – Trajectoire de marche au hasard 2D : traces de pas du marcheur dans le plan (20 000 pas).
| 3 |
L’objectif de cette troisième séance est d’une part d’approfondir les aspects de Matlab et Octave liés à l’algèbre linéaire, et de voir (ensuite) une application à la régression linéaire.
Matlab est un langage spécialement adapté à l’algèbre linéaire. Ainsi de nombreuses opérations (telles que par exemple la multiplication matricielle) sont implémentées, et sont beaucoup plus efficaces que si on les reprogramme en utilisant des boucles. Par exemple, la fonction suivante testmatvectmult.m compare deux versions de la multiplication matrice-vecteur (les instructions tic et toc permettent d’avoir une estimation du temps de calcul écoulé entre le tic et le toc).
function [W1,W2] = testmatvectmult(A,V)
% testmatvectmult : teste la multiplication matrice-vecteur
% usage : [W1,W2] = testmatvectmult(A,V)
% Arguments : matrice A et vecteur V, multipliables
[m A,n A] = size(A); [m V,n V] = size(V); if ?(n V==1) error(’V doit etre un vecteur (colonne)’);
end
if n A?=m V
error(’matrice et vecteur non multipliables’);
end
tic
W1 = zeros(mA,1); for m = 1 :mA for k = 1 :n A
W1(m) = W1(m) + A(m,k)*V(k); end
end fprintf(’Multiplication avec boucle : ’) toc
tic
W2 = A*V;
fprintf(’Multiplication vectorielle : ’) toc
Exercice 3.1 Télécharger cette fonction sur le site du cours, et la tester sur des matrices de tailles 10, 100, 200,
19
Il existe un certain nombre d’algorithmes permettant de décomposer une matrice en produit de matrices plus simples. On peut par exemple citer
– la décomposition QR (fonction qr), qui décompose une matrice en un produit d’une matrice triangulaire supérieure R et d’une matrice unitaire Q,
– la décomposition LU (fonction lu), qui décompose en un produit d’une matrice triangulaire supérieure U et d’une matrice L, produit d’une matrice triangulaire inférieure et d’une matrice de permutation (ne pas confondre ces dernières avec les parties triangulaires supérieure et inférieure obtenues avec tril et triu),
– la décomposition de Cholesky (fonction chol), qui factorise une matrice définie positive A en un produit A = LL? ou` L est triangulaire inférieure.
| Exercice 3.2 – Tester ces fonctions sur de petites matrices. – Ecrire une fonction QRresol.m prenant comme variable d’entrée une matrice carrée A ? M(N) et un vecteur y ? M(N,1), et utilisant la décomposition QR pour – calculer le déterminant de A et donc tester si A est inversible – si A est inversible, résoudre le système Ax = y (qui devient donc Rx = Q?1y). On rappelle que le déterminant d’une matrice unitaire est égal à 1, et que son inverse est égal à sa conjuguée hermitienne (complexe conjuguée de la transposée). On aura besoin d’expliciter l’inverse d’une matrice triangulaire supérieure, et de le coder (en utilisant une boucle). |
L’instruction det permet de calculer le déterminant d’une matrice carrée. Il est aussi instructif de calculer ce déterminant directement. On peut pour celà utiliser l’expression en fonction des cofacteurs. On rappelle que les cofacteurs d’une matrice A = {ak`} ? MN sont définis par
ck`(A) = (?1)k+` det(A(k,`)) ,
ou` A(k,`) est la matrice réduite, obtenue à partir de A en lui enlevant la k-ième ligne et la `-ième colonne. Le déterminant s’exprime alors sous la forme
det(A) = Xak0`ck0`(A) = Xak`0ck`0(A) ,
` k
ou` `0 (resp. k0) est une colonne (resp. ligne) fixée.
Exercice 3.3 Ecrire une fonction cofact.m, prenant comme variables d’entrée une matrice carrée A (à N lignes et N colonnes), et deux entiers 1 ? k,` ? N, et retournant le cofacteur correspondant ck`(A).
Exercice 3.4 Ecrire une fonction mondet.m, calculant le déterminant d’une matrice carrée en utilisant la fonction cofact.m ci-dessus.
Exercice 3.5 Ecrire une fonction cofact2.m et une fonction mondet2.m s’appelant l’une l’autre. Comparer le résultat obtenu avec det et mondet sur des matrices de tailles 3,4,5,6, Attention, mondet2 devient vite inefficace lorsque N grandit
L’instruction rref analyse l’espace image d’une matrice, qu’elle retourne sous forme réduite. Par exemple, la séquence d’instructions
>> A = magic(4);% Génère une matrice particulière
>> [R,jb] = rref(A);% réduit la matrice
>> Rang = length(jb);
>> Base = R(jb);
3.1. GEN´ ERALIT´ ES´
Fig. 3.1 – Trois trajectoires solutions du système X0(t) = AX(t), pour trois conditions initiales différentes, et deux matrices A = A1 (gauche) et A = A2 (droite). Les valeurs du temps sont prises sur une échelle logarithmique.
donne une évaluation du rang de A, et une matrice (ici appelée Base) dont les colonnes forment une base de l’image de A. Une base orthonormale de l’image peut aussi être directement obtenue grace à orth, et l’instruction rank donne le rang.
L’instruction null retourne une base orthonormale du noyau d’une matrice carrée.
On a déjà rencontré les fonctions concernant la diagonalisation et l’inversion des matrices. L’inverse d’une matrice carrée est obtenue via la fonction inv. Se documenter sur cette instruction (doc inv). Lorsque la matrice est proche d’être singulière, un message d’erreur s’affiche dans la fenêtre de commande. Voir aussi l’instruction cond. Les fonction mldivide et rldivide retournent des inverses à gauche et à droite d’une matrice, pas nécessairement carrée. Attention : Matlab effectuant des opérations numériques, avec une précision donnée, il est capable d’« inverser » des matrices non inversibles mais donnera un message d’erreur dans ce cas là. De là, il est possible de résoudre des systèmes linéaires à l’aide de ces instructions pré-définies
| Exercice 3.6 Générer trois matrices et deux vecteurs en utilisant la séquence d’instructions suivante : >> A = [1,2,3;4,5,6;7,8,10]; >> B = A(1 :2, :); >> C = A( :,1 :2); >> D = A(1 :2,1 :2); >> x = [1;2;3]; >> y = x(1 :2); et interpréter ces instructions au vu des résultats obtenus. Parmi les systèmes matriciels suivants, quels sont ceux qui admettent des solutions? une solution unique? x = Az ; x = Cz ; y = Bz ; y = Dz Résoudre ces systèmes en utilisant l’instruction inv ou mldivide. Utiliser l’aide en ligne pour interpréter le résultat obtenu. |
La fonction eig, qu’on a déjà vue, permet de calculer à la fois les vecteurs propres et valeurs propres d’une matrice carrée donnée. On rappelle que la syntaxe est la suivante :
>> [VectPropres, ValPropres] = eig(A);
et que ValPropres est une matrice diagonale, VectPropres étant une matrice dont les colonnes sont les vecteurs propres correspondant aux valeurs propres trouvées, dans le même ordre.
| Exercice 3.7 1. partant de la matrice A ci-dessous résoudre numériquement le système différentiel en utilisant la diagonalisation de A, et le fait que la solution s’exprime de fa¸con explicite : x(t) = PetDP?1x(0) . On pourra par exemple tracer le résultat obtenu. 2. Résoudre le problème précédent de fa¸con systématique : écrire une fonction masolexp.m prenant comme variables d’entrée la matrice A, la condition initiale, et un ensemble de valeurs de t, et retournant comme variables de sortie les coordonnées de la solution x(t). 3. Dans le cas bidimensionnel, compléter cette fonction en trac¸ant les valeurs ainsi obtenues comme des points dans le plan. On pourra essayer les matrices suivantes , et il est conseillé de prendre des valeurs du temps sur une échelle logarithmique (par exemple avec une instruction du type T=log(linspace(0.01,1,100)), qui prend les logarithmes de 100 valeurs régulièrement espacées entre 0.01 et 1). Un exemple de trajectoire avec la matrice A1 (et trois conditions initiales différentes) se trouve en figure 3.1. |
Un résultat classique d’algèbre linéaire montre qu’étant donnés n vecteurs orthogonaux deux à deux et normalisés dans un espace Euclidien (ou de Hilbert) de dimension N > n, il est toujours possible de trouver N ? n vecteurs normalisés de ce même espace tels que la famille complète des N vecteurs forme une base orthonormée de l’espace considéré.
La procédure de Gram-Schmidt fournit un procédé constructif permettant d’obtenir ces vecteurs. Notons v1, vn ? E les n premiers vecteurs. Soit w ? E un vecteur n’appartenant pas au sous-espace de E engendré par {v1, vn}. En posant
et
il est facile de vérifier que la famille {v1, vn+1} est une famille orthonormée dans E.
L’orthogonalisation de Gram-Schmidt utilise ce principe autant de fois qu’il le faut pour compléter la famille initiale et engendrer ainsi une base orthonormée de E.
La version matricielle de cette construction est la suivante. En mettant en colonne les composantes des n vecteurs initiaux dans la base canonique, on obtient une matrice M à N lignes et n colonnes. Cette matrice
3.2. PSEUDO-INVERSE
est isométrique, c’est à dire telle que M?M = In ; elle n’est par contre pas unitaire, car MM? 6= IN. La procédure de Gram-Schmidt permet alors de concaténer itérativement N ? n colonnes à M, de sorte à obtenir une matrice U unitaire, c’est à dire telle que U?U = UU? = IN.
Exercice 3.8 Ecrire une fonction GramSchmidtComplete.m prenant comme variable d’entrée une matrice isométrique N × n (avec n < N), et retournant une matrice unitaire complétée à la Gram-Schmidt. La fonction vérifiera que la matrice d’entrée est bien isométrique.
Pour utiliser la procédure de Gram-Schmidt, il est nécessaire de savoir générer à chaque itération un vecteur n’appartenant pas au sous espace engendré par les vecteurs initiaux. L’instruction randn permettra de générer un vecteur qui vérifiera presque suˆrement cette propriété.
Théorie De fa¸con générale, la théorie de la diagonalisation s’applique aux endomorphismes (et donc aux matrices carrées). Dans le cas de morphismes entre espaces vectoriels différents, il reste possible d’effectuer des opérations similaires, même lorsque les espaces sont de dimensions différentes.
Theor´ eme 3.1` Soit M ? MK(m,n) une matrice m × n à coefficients dans le corps K, avec K = R ou K = C. Alors il existe une factorisation de la forme :
M = U?V?,
ou` U est une matrice unitaire m × m sur K, ? est une matrice m × n dont les coefficients diagonaux sont des réels positifs ou nuls et tous les autres sont nuls (c’est donc une matrice diagonale dont on impose que les coefficients soient positifs ou nuls), et V?est la matrice adjointe de V , matrice unitaire n×n sur K. On appelle cette factorisation la décomposition en valeurs singulières de M.
– La matrice V contient un ensemble de vecteurs de base orthonormés pour M, dits vecteurs d’entrée; – La matrice U contient un ensemble de vecteurs de base orthonormés pour M, dits vecteurs de sortie; – La matrice ? contient les valeurs singulières de la matrice M.
3.2.3 Construction de la SVD :
La construction de cette matrice peut se faire de la fa¸con suivante : on considère la matrice
A = M?M ,
qui est par construction semi-définie positive, donc Hermitienne. Elle est donc diagonalisable, avec des valeurs propres réelles positives ou nulles, et il existe une matrice D ? M(n,n) et une matrice de passage V ? M(n,n) telles que
A = V DV?1.
V peut être choisie unitaire, de sorte que V?1 = V?. La matrice diagonale D prend quant à elle la forme
,
D1 ? M(r,r) étant une matrice diagonale dont les termes diagonaux sont strictement positifs (r est le rang de A et D). Il est usuel d’ordonner les valeurs diagonales de D par ordre décroissant.
On note V1 ? M(n,r) la matrice dont les colonnes sont les vecteurs propres de A de valeur propre non nulle (les valeurs propres étant rangées en ordre décroissant), et V2 ? M(n,n ? r) la matrice formée des autres vecteurs propres. Notons que V1 n’est pas unitaire : on a mais.
Posons maintenant (D1 étant inversible)
On a alors, mais. Il est par contre possible de montrer qu’en complétant U1 par m ? r colonnes orthogonales deux à deux, et orthogonales aux colonnes de U1, on obtient une matrice unitaire U. On montre alors
U? = MV ,
ou` ? s’obtient en complétant D1 par des zéros, d’ou` le résultat.
Remarque : cette démonstration de l’existence de la décomposition est constructive, elle fournit un algorithme de calcul de la svd qu’on va programmer plus loin. Ceci dit, cet algorithme est connu pour être numériquement instable (en présence de petites valeurs singulières), c’est pourquoi les bibliothèques lui préfèrent des algorithmes plus élaborés.
| Exercice 3.9 Ecrire une fonction masvd.m, prenant en entrée une matrice M, et calculant sa décomposition en valeurs singulières. – Diagonaliser la matrice A = M?M. Réorganiser le résultat (les matrices D et V ) de sorte que les valeurs propres soient classées par ordre décroissant (utiliser la fonction sort). Tracer les valeurs propres classées (en utilisant l’instruction bar plutôt que plot). – Extraire la matrice V1 et la matrice D1, et en déduire la matrice U1. – Compléter U1 en une matrice unitaire U : on pourra pour cela utiliser la procédure d’orthogonalisation de Gram-Schmidt. Compléter ?. – La fonction retournera les matrices U,V , et ?. |
Remarque : L’analyse en composantes principales présente de grosses similarités avec la SVD. En effet,
étant donné un tableau de données X, on effectue une diagonalisation de la matrice X?X (ou XX?), après avoir centré les lignes ou les colonnes de X. L’opération finale de l’ACP est une projection des données sur les sous-espaces engendrés par un certain nombre de vecteurs propres.
Soient E et F deux espaces vectoriels, et soit f ? L(E,F) une application linéaire de E dans F. Lorsque f n’est pas inversible, la notion de pseudo-inverse fournit un substitut intéressant.
Théorie
Definition 3.1´ Une application linéaire g ? L(F,E) est appelée pseudoinverse de f lorsqu’elle satisfait les deux conditions
f ? g ? f = f g ? f ? g = g
On montre dans ce cas les propriétés suivantes
– E est la somme directe du noyau de f et de l’image de g
– F est la somme directe du noyau de g et de l’image de f
– f induit un isomorphisme de l’image de g sur l’image de f
– g induit un isomorphisme de l’image de f sur l’image de g, qui est l’inverse de l’isomorphisme précédent. Il s’agit bien d’une extension de la notion d’inverse dans la mesure ou` si f admet effectivement un inverse, celui-ci est aussi un, et même l’unique, pseudo-inverse de f.
Interprétation La pseudo-inverse d’une matrice A ? M(m,n) peut être interprétée de fac¸on assez intuitive en termes de systèmes linéaires. En considérant le système linéaire
AX = Y ,
A†Y est, parmi les vecteurs X qui minimise l’erreur en moyenne quadratique kY ? AXk, celui dont la norme kXk est minimale.
– Si le rang de A est égal à m, alors AA? est inversible, et A† = A?(AA?)?1
– Si le rang de A est égal à n, alors A?A est inversible, et A† = (A?A)?1A? Si le rang de A est égal à m = n, alors A est inversible, et A† = A?1.
3.3. DEVOIR
Construction des pseudo-inverses La description géométrique qui vient d’être donnée pour les pseudo-inverses est une propriété caractéristique de ceux-ci. En effet, si on introduit des supplémentaires K du noyau de f (dans E) et I de l’image de f (dans F), il est possible de leur associer un unique pseudo-inverse g de f. Les restrictions de g aux espaces I et image de f sont parfaitement définies : application nulle sur I, et inverse de l’isomorphisme induit par f sur K, sur Im(f). Les deux propriétés de pseudo-inverses sont alors facilement vérifiées en séparant selon les couples d’espaces supplémentaires. Cette construction montre qu’en général il n’y a pas unicité du pseudo-inverse associé à une application linéaire.
| Proposition 3.1 Supposons maintenant que E et F soient deux espaces Hermitiens (Euclidiens dans le cas réel). Soit A la matrice de f par rapport à deux bases données de E et F, soit A†la matrice du pseudo-inverse g. Alors A†est l’unique matrice qui satisfait les conditions : 1. AA†A = A 2. A†AA† = A† 3. AA†est une matrice Hermitienne : (AA†)? = AA†. 4. A†A est une matrice Hermitienne : (A†A)? = A†A. ou` on a noté M?la conjuguée Hermitienne de la matrice M. |
Calcul de la pseudo-inverse Notons k le rang de A ? M(m,n), une matrice m × n. A peut donc
être décomposée sous la forme
A = BC ,
ou` B ? M(m,k) et C ? M(k,n). On a alors
Proposition 3.2 Avec les notations ci-dessus, la pseudo-inverse de A s’écrit
A† = C?(CC?)?1(B?B)?1B?
Si k = m, B peut être la matrice identité. De même, si k = n, C peut être l’identité. Ainsi, la formule se simplifie.
En pratique, la pseudo-inverse se calcule à partir de la décomposition en valeurs singulières :
A = U?V?,
ou` ? est diagonale. On a alors
A† = V ?†U?,
ou` ?† est obtenue à partir de ? par transposition, et remplacement des éléments non nuls par leur inverse.
Exercice 3.10 En utilisant soit votre fonction masvd, soit la fonction svd, construire une fonction mapseudoinverse.m.
Les devoirs sont à rendre par binôme. Finaliser les fonctions masolexp.m, GramSchmidtComplete.m, masvd.m et mapseudoinverse.m demandées dans les exercices des sections précédentes.
Les devoirs de cette section sont à transmettre sous forme d’archive, par exemple sous forme de fichier noms , nom ou noms . L’archive devra contenir les fichiers des programmes, mais aussi un petit compte-rendu décrivant le travail qui a été fait, le cas échéant ce qui n’a pas pu être fait
ou terminé, ainsi que quelques explications sur l’utilisation des fonctions. On pourra aussi ajouter des fichiers de graphiques (format jpg par exemple).
26 3. ALGEBRE LIN` EAIRE´
| 4 |
L’interpolation est une opération mathématique permettant de construire une courbe à partir de la donnée d’un nombre fini de points, ou une fonction à partir de la donnée d’un nombre fini de valeurs. La solution du problème d’interpolation passe par les points prescrits, et, suivant le type d’interpolation, il lui est demandé de vérifier des propriétés supplémentaires.
Le type le plus simple d’interpolation est l’interpolation linéaire, qui consiste à ”joindre les points” donnés. Il s’agit d’une procédure locale, à l’inverse de l’interpolation globale telle que l’interpolation polynômiale.
L’interpolation doit être distinguée de l’approximation de fonction, aussi appelée régression, qui consiste à chercher la fonction la plus proche possible, selon certains critères, d’une fonction donnée. Dans le cas de la régression, il n’est en général plus imposé de passer exactement par les points donnés initialement. Ceci permet de mieux prendre en compte le cas des erreurs de mesure, et c’est ainsi que l’exploitation de données expérimentales pour la recherche de lois empiriques relève plus souvent de la régression linéaire, ou plus généralement de la méthode des moindres carrés.
L’interpolation polynomiale consiste à utiliser un polynôme unique (et non des segments comme on le verra plus loin), de degré aussi grand que nécessaire, pour estimer localement l’équation représentant la courbe afin de déterminer la valeur entre les échantillons. Comme on le verra, l’interpolation polynômiale a tendance à devenir instable lorsque le degré du polynôme croˆ?t.
Le prototype de l’interpolation polynômiale est l’interpolation de Lagrange, décrite et analysée ci-dessous.
| Theor´ eme 4.1` Soient {(xk,yk),k = 0 n} n + 1 points donnés, tels que les xk soient tous deux à deux différents. Il existe un et un seul polynoˆme de degré n, noté pn, tel que pn(xk) = yk pour tout k = 0, n. Ce polynôme est de la forme , (4.1) ou` les `k sont des polynômes de degré n donnés par . (4.2) |
Ce résultat peut se démontrer de différentes fa¸cons. La plus simple est sans doute d’effectuer le calcul explicite : le polynôme p recherché doit satisfaire p(xk) = yk pour tout k = 0, n, ce qui conduit au système
| qui s’écrit sous forme matricielle | ||||
|
?1 1 ? ?1 ? ? ? ? 1 | x0x1x2 xn | x20x21x22x2n | xn0??a0? ?y0? xn1??a1? ?y1? xn2????a2?? = ??y2?? ?? .. ?? ?? ??? ?? . ? ? ?? xnn an yn 27 |
Ce système admet une solution unique si et seulement si le déterminant de cette matrice, appelé déterminant de Vandermonde est non-nul. On peut alors utiliser le résultat suivant
Comme nous avons supposé que les xk sont deux à deux disjoints, nous sommes assurés de l’existence et unicité de la solution. Finalement, il suffit de vérifier que la solution proposée satisfait bien les conditions exigées, et l’unicité fait le reste. ?
Matlab n’implémente pas de fonction effectuant explicitement l’interpolation de Lagrange. Ceci dit, il est possible d’utiliser une alternative, basée sur l’approximation par moindres carrés, qu’on verra plus loin. On se contente à ce point du résultat suivant, dont la preuve est immédiate.
| Proposition 4.1 Soient {(xk,yk),k = 0 n} n + 1 points donnés, tels que les xk soient tous deux à deux différents. Soit qn le polynôme de degré n qui minimise l’erreur d’approximation en moyenne quadratique n qn = argminX|yk ? q(xk)|2. q k=0 Alors qn = pn, le polynôme d’interpolation de Lagrange. |
La fonction polyfit de Matlab effectue l’interpolation polynômiale par moindres carrés. On peut donc l’utiliser pour résoudre le problème d’interpolation de Lagrange, et illustrer ses caractéristiques. On utilise ensuite la fonction polyval pour évaluer le polynôme ainsi obtenu en des points bien choisis.
| Exercice 4.1 Se renseigner sur les fonctions polyfit et polyval. Etudier l’exemple suivant, qui illustre l’utilisation de ces fonctions. >> n = 4; >> x = rand(n+1,1); Génération de 5 abscisses tirées au hasard >> x = sort(x); Réarrangement dans l’ordre croissant >> y = randn(n+1,1); Génération d’ordonnées aléatoires >> P = polyfit(x,y,n); Calcul du polynôme interpolateur >> xmin = min(x); xmax = max(x); >> xx = linspace(xmin,xmax,50); >> yy = polyval(P,xx); >> plot(xx,yy); >> hold on; >> plot(x,y,’.r’); Effectuer cet exercice plusieurs fois (en changeant aussi la valeur de n). Sauvegarder quelques figures significatives, et commenter ce que l’on observe dans un fichier compte rendu. |
Un exemple se trouve en Figure 4.1. On peut observer que le polynôme interpolant peut s’éloigner très significativement des données. Ce phénomène est à rapprocher du phénomène de Runge, que l’on rencontre lorsque l’on s’intéresse aux propriétés de convergence de l’approximation d’une fonction par ses polynômes de Lagrange.
Une question naturelle pour le mathématicien est d’analyser la précision de l’interpolation dans le cas d’une fonction connue : en d’autres termes, étant donnée une fonction f, quelques valeurs ponctuelles f(xk), et le polynôme pn de degré n obtenu par interpolation de lagrange à partir de celles-ci, quelle est la différence entre f et pn ? On déduit facilement du théorème de Rolle le résultat suivant.
4.1. INTERPOLATION
Fig. 4.1 – Interpolation polynômiale : polynôme de degré 5; les données initiales sont représentées par des étoiles.
| Proposition 4.2 Soit f ? Cn+1une fonction sur un intervalle I = [a,b], et soit pn le polynôme de Lagrange associé aux données (x0,y0), (xn,yn)}. Alors, ?x, il existe ?x tel que ou` ?n+1est le polynôme de degré n + 1 défini par. |
On en déduit facilement la borne
,
de sorte que la précision de l’approximation dépend, comme on peut s’y attendre, de la régularité de f, mais aussi de la position des points xk. Sans hypothèse supplémentaire, on a en se plac¸ant dans un intervalle [a,b]
|?n+1(x)| ? (b ? a)n+1,
alors qu’en se pla¸cant dans le cas de points xk régulièrement espacés dans l’intervalle [a,b] on peut montrer que
.
Dans un cas comme dans l’autre, cette borne n’est pas très précise.
Si f est une fonction analytique (c’est à dire si sa série entière converge dans un domaine de rayon R), alors la croissance des dérivées successives f(n+1)(?x) en fonction de n peut être estimée, et il est possible de montrer que pn converge vers f dans un intervalle dont la largeur dépend de R. L’exercice qui suit est un exemple classique montrant la (non) convergence de l’approximation par polynôme de Lagrange.
| Exercice 4.2 On considère la fonction u? : [?1,1] ? R définie par , ou` ? ? R+ est un réel positif fixé, et on s’intéresse à ses approximations polynômiales par interpolation de Lagrange. Ecrire une fonction erreurinterpol.m, prenant comme variable d’entrée la valeur de ? et un ensemble de points d’échantillonnage x, et évaluant la précision de l’approximation par polynôme interpolateur de Lagrange : 1. Calcul des valeurs de u? sur les points x. 2. Evaluation du polynôme interpolateur 3. Génération d’un vecteur de points régulièrement espacés xx (utiliser la fonction linspace), et évaluation de u? en ces points. 4. Evaluation de l’approximation de la fonction u? en ces points 5. Calcul de l’erreur relative : norme de l’erreur divisé par la norme de la référence ku?k. La fonction retournera l’erreur relative en variable de sortie. On effectuera plusieurs essais, avec des points x tirés aléatoirement entre -1 et 1, puis régulièrement espacés; on fera quelques courbes significatives, avec diverses valeurs de ? et n, qu’on commentera dans le compte rendu. |
Interpolation linéaire par morceaux Dans le cas d’une interpolation linéaire, on constitue une courbe d’interpolation qui est une succession de segments. Entre deux points p1 = (x1,y1) et p2 = (x2,y2), l’interpolation est donnée par la formule suivante
y = p · (x ? x1) + y1
ou` la pente p s’exprime comme
Exercice 4.3 Programmer une fonction Matlab interpollin.m, prenant comme variables d’entrée les coordonnées {(x0,y0), (xn,yn)}, ainsi qu’une liste de valeurs uk, calculant une interpolation linéaire par morceaux, la tra¸cant dans une figure, et retournant comme variable de sortie la liste de valeurs interpolées.
On sauvegardera un (ou plusieurs) exemples significatifs de figure, qu’on commentera dans le compte rendu.
Interpolation spline L’interpolation “spline” généralise l’interpolation linéaire par morceaux en rempla¸cant les morceaux affines par des morceaux polynômiaux. La plus “classique” est l’interpolation par fonctions spline “cubiques”, qui sont donc des polynômes de degré 3.
4.1. INTERPOLATION
Fig. 4.2 – Interpolation spline : linéaire (rouge) et cubique (bleu); les données initiales sont représentées par des étoiles.
| Definition 4.1´ Etant donnés n + 1 points {(x0,y0), (xn,yn)} (appelés “noeuds”), une fonction interpolante spline cubique passant par ces points est une fonction polynômiale de degré 3 par morceaux ?S0(x), x ? [x0,x1] ??? S1(x), x ? [x1,x2] S(x) = ··· ??? Sn?1(x), x ? [xn?1,xn] telle que 1. S(xk) = f(xk) pour k = 0, n 2. S est continue sur les noeuds, Sk?1(xk) = Sk(xk), pour k = 1, ,n ? 1 3. S est deux fois continuˆment différentiable ux noeuds : et , pour k = 1, ,n ? 1. |
Pour déterminer une telle fonction d’interpolation, on doit déterminer 4n paramètres (chaque polynôme de degré 3 est caractérisé par 4 paramètres). Les conditions d’interpolation donnent n+1 équations linéaires, et les conditions de continuité en donnent 3n ? 3, donc au total 4n ? 2 conditions. Les deux conditions restantes sont au choix de l’utilisateur; par exemple, les splines cubiques “naturelles” supposent
S00(x0) = S00(xn) = 0 .
| Exercice 4.4 Se documenter sur la fonction interp1, et tester les commandes suivantes. >> x = rand(50,1); >> x = sort(x); >> y = randn(size(x)); >> plot(x,y,’.’); >> xmin = min(x); xmax = max(x); >> x2 = linspace(xmin,xmax,100); >> y2 = interp1(x,y,x2,’linear’); >> hold on; plot(x2,y2,’r’); Comparer le résultat de l’interpolation linéaire effectuée dans l’exercice 4.3 avec celui de l’interpolation effectuée avec l’option ’linear’ de interp1, puis l’option ’spline’. Tracer un exemple significatif en comparant l’interpolation linéaire par morceaux et l’interpolation cubique par morceaux. Sauvegarder la figure, et la commenter dans le compte rendu. Un exemple se trouve en Figure 4.2 |
Fig. 4.3 – Régression : linéaire (rouge/trait plein) et polynôme d’ordre 5 (vert/tirets); les données initiales sont représentées par des points.
Dans le problème de régression, il s’agit encore une fois de décrire des données par une fonction s’en approchant, mais on relâche cette fois la contrainte de “passer par les points exactement” par une contrainte de type “moindres carrés”.
On peut là encore faire de la régression linéaire ou polynômiale, ou de la régression polynômiale “par morceaux”.
Commen¸cons par le cas unidimensionnel. Supposons données des observations {(x0,y0), (xn,yn)}, que l’on suppose liées par une relation linéaire
,
ou` a et b sont deux réels, et ou` les k représentent un “bruit”, c’est à dire une erreur.
Comme on l’a déjà vu, les coefficients a et b peuvent être estimés par moindres carrés, en résolvant le problème minX|yk ? axk ? b|2 (4.4) a,b k
En égalant à zéro les dérivées de cette expression par rapport à a et b, on obtient
Un exemple de régression linéaire se trouve dans la figure 4.3 (avec un exemple correspondant de régression polynômiale, voir plus bas).
Dans une problématique de statistique descriptive, il arrive fréquemment que l’on désire savoir si deux variables sont liées linéairement ou pas. La décision peut se faire via un test sur le coefficient de corrélation de Pearson.
,
ou` sx et sy sont les écarts-types de x et y
.
4.2. REGRESSION´
Lemme 4.2 Supposons données des observations {(x1,y1), (xn,yn)}. Sous l’hypothèse d’indépendance de x et y, le coefficient de Pearson suit une loi de Student-t à n ? 1 degrés de liberté.
Ainsi, après avoir effectué une régression linéaire, il est possible de tester la significativité de la relation linéaire en comparant la valeur obtenue pour ? à sa distribution sous hypothèse de décorrélation.
| Exercice 4.5 Ecrire une fonction regrlin.m prenant en entrée deux vecteurs de même dimension, et retournant les nombres a et b estimés, ainsi que l’écarttype de l’erreur estimée ˆ et le coefficient de corrélation de Pearson. La fonction représentera aussi graphiquement les données et la droite de régression sur le même graphique. Elle représentera aussi l’histogramme des erreurs ˆk dans un autre graphique. |
Exercice 4.6 Télécharger le fichier se trouvant sur http ?
et l’importer dans Matlab en utilisant soit l’instruction importdata (se renseigner), soit l’utilitaire d’importation se trouvant dans l’onglet File du menu de Matlab . Ce fichier contient des données relatives au poids corporel et poids du cerveau dans un certain nombre d’espèces animales.
En trac¸ant l’une contre l’autre, puis le logarithme de l’une contre le logarithme de l’autre, estimer les paramètres de ce qui vous semble être un bon modèle reliant ces deux quantités.
Remarque : La fonction importdata génère des structures, c’est à dire des variables possédant plusieurs champs. Par exemple,
>> tmp = importdata(’’); retourne une variable tmp possédant en général trois champs : qui contient les données elles mêmes, mais aussi tmp.textdata et tmp.colheaders. On pourra examiner les contenus de ces champs pour se familiariser avec cette fonction.
On considère maintenant le problème multidimensionnel : on suppose données des observations scalaires, mais des variables vectorielles (Xk,yk), et on cherche à modéliser ces données sous la forme
Ici, chaque Xk est un vecteur de dimension N, que l’on va écrire en ligne, et ? est un vecteur (colonne)
à N lignes. A partir de là on peut former une matrice X ? M(K,N), et écrire le problème sous la forme , ou même
en ayant posé
De là, la solution des moindres carrés au problème de régression multiple s’obtient en résolvant
,
problème dont on a vu que la solution est donnée par la pseudo-inverse Z† de Z :
?ˆ = Z†Y .
Exercice 4.7 Programmer une fonction regrlinmult.m prenant en entrée une matrice de données X et un vecteur d’observations Y , vérifiant leurs tailles respectives, et retournant le vecteur ? et le scalaire ? de la régression, ainsi que l’écart-type de l’erreur de régression. La fonction tracera les valeurs prédites par le modèle en fonction des vraies valeurs, l’erreur de prédiction en fonction des valeurs, et affichera l’histogramme des erreurs de régression.
| Exercice 4.8 Télécharger le fichier se trouvant sur http ? et l’importer dans Matlab en utilisant soit l’instruction importdata (se renseigner), soit l’utilitaire d’importation se trouvant dans l’onglet File du menu de Matlab . Ce fichier contient des données relatives à la hauteur de neige en montagne en fonction de divers paramètres environnementaux (pente, altitude, rugosité, ). Après avoir effectué une régression linéaire de la hauteur en fonction de quelques unes des variables, effectuer une régression linéaire multiple. Dans tous les cas, on conservera les valeurs de l’écart-type de l’erreur de régression, et on comparera les résultats obtenus. |
La fonction polyfit que nous avons déjà rencontrée plus haut effectue la régression polynômiale : étant donnée une famille d’observations {(xk,yk),k = 0, K}, et un ordre fixé N, polyfit produit le polynôme P de degré N x ? P(x) solution du problème de régression
!
On a déjà vu que lorsque N = K, la régression par moindres carrés co¨?ncide avec l’interpolation polynômiale, qui est malheureusement souvent instable. Dans le cadre de problèmes de régression, on se place souvent dans le cas.
Exercice 4.9 Programmer une fonction regrpolyn.m prenant en entrée deux vecteurs x et y, un degré de polynôme N, faisant les vérifications nécessaires, et effectuant la régression polynômiale (en utilisant polyval). En outre, la fonction affichera dans une fenêtre graphique les données ainsi que le polynôme obtenu, et dans une autre fenêtre l’histogramme des erreurs de régression. La fonction affichera à l’écran l’écart-type de l’erreur.
Exercice 4.10 En utilisant les données , utiliser la fonction regrpolyn pour effectuer une régression polynômiale de différents degrés. On tracera l’évolution de l’erreur en fonction du degré du polynôme, et on interprétera les résultats obtenus.
Finaliser les exercices de ce chapitre. Les devoirs sont à rendre (par binôme) sous forme d’une archive contenant
– Les fichier Matlab des fonctions développées; celles ci devront être suffisamment commentées pour qu’un utilisateur puisse les utiliser sans effort.
– Un bref compte-rendu, décrivant dans chaque cas l’approche suivie pour résoudre le problème posé, la syntaxe de la fonction programmée si nécessaire; le compte-rendu pourra également inclure des illustrations (sauvegardes de figures au format JPEG (extention .jpg) par exemple) si vous le jugez utile.
– En cas de difficultés pour inclure des figures dans le texte, on pourra plus simplement ajouter à l’archive les fichiers .jpg et préciser dans le compte-rendu quelle figure on doit regarder.
| 5 |
Un générateur de nombres pseudo-aléatoires est un algorithme qui génère une séquence de nombres présentant certaines propriétés du hasard. Par exemple, les nombres sont supposés être approximativement indépendants, et il est potentiellement difficile de repérer des groupes de nombres qui suivent une certaine règle (comportements de groupe).
Cependant, les sorties d’un tel générateur ne sont évidemment pas aléatoires; elles s’approchent seulement des propriétés idéales des sources complètement aléatoires. De vrais nombres aléatoires peuvent être produits avec du matériel qui tire parti de certaines propriétés physiques stochastiques (bruit d’une résistance par exemple).
La raison pour laquelle on se contente d’un rendu pseudo-aléatoire est :
– d’une part, il est difficile d’obtenir de “vrais” nombres aléatoires et, dans certaines situations, il est possible d’utiliser des nombres pseudo-aléatoires, en lieu et place de vrais nombres aléatoires.
– d’autre part, ce sont des générateurs particulièrement adaptés à une implémentation informatique, donc plus facilement et plus efficacement utilisables.
Les méthodes pseudo-aléatoires sont souvent employées sur des ordinateurs, dans diverses tâches comme la méthode de Monte-Carlo, la simulation ou les applications cryptographiques. Une analyse mathématique rigoureuse est nécessaire pour déterminer le degré d’aléa d’un générateur pseudo-aléatoire.
La plupart des algorithmes pseudo-aléatoires essaient de produire des sorties qui sont uniformément distribuées. Une classe très répandue de générateurs utilise une congruence linéaire. D’autres s’inspirent de la suite de Fibonacci en additionnant deux valeurs précédentes ou font appel à des registres à décalage dans lesquels le résultat précédent est injecté après une transformation intermédiaire.
Rappelons tout d’abord la définition probabiliste de l’histogramme.
Definition 5.1´ Soient (X1, XN) des variables aléatoires i.i.d.
(indépendantes identiquement distribuées) à valeurs dans un ensemble E. Soit E = U1 ? U2 ? ··· ? UKune partition de E en K sous-ensembles. L’histogramme correspondant est la collection de variables aléatoires (H1, HK) définies par
N
Hk= XIUk(Xn) .
n=1
Il existe des résultats mathématiques précisant les propriétés de l’histogramme. Par exemple, dans la limite des grandes dimensions, on a le résultat de convergence suivant :
| Lemme 5.1 Avec les notations précédentes, on a lim Hk = µ(Uk) . N?? |
En statistiques, un histogramme est un graphe permettant de représenter la répartition d’une variable continue. L’histogramme est un moyen simple et rapide pour représenter la distribution de cette variable. Nous utiliserons l’histogramme pour visualiser la distribution d’une collection de nombres (pseudo) aléatoires (voir par exemple la Figure 5.1. En pratique, il est important de savoir régler le nombre
35
Fig. 5.1 – Histogramme
K de classes à utiliser. K ne doit être ni trop grand (auquel cas l’histogramme n’est pas représentatif, ls classes ayant souvent des effectifs trop faibles), ni trop petit (auquel cas l’histogramme ne représente pas la distribution étudiée de fa¸con assez fine). Il existe des règles empiriques reliant le nombre de classes à la taille de l’échantillon, par exemple
?
K ? N , ou.
Ceci dit, l’histogramme étant avant tout un outil de visualisation, il est conseillé de l’utiliser pour plusieurs valeurs de N.
Comme on l’a déjà vu, la commande Matlab hist permet de tracer l’histogramme d’une population de nombres. Elle permet également de retourner les effectifs des classes créées, ainsi que les valeurs correspondantes de la variable. Par exemple, la suite d’instructions
>> X = randn(5000,1);
>> [NN,XX] = hist(X,20);
effectue une partition du domaines de valeurs des nombres créés en 20 intervalles, et compte les effectifs de ces 20 classes, mais ne trace pas l’histogramme. Celui-ci peut être obtenu par >> bar(XX,NN);
Les générateurs congruentiels génèrent des nombres entiers positifs inférieurs à une certaine valeur maximale N, en utilisant une congruence modulo N. Les suites de nombres ainsi obtenues sont évidemment périodiques avec une période inférieure ou égale à N. Les plus simples sont les générateurs de Lehmer, introduits en 1948
Definition 5.2´ Un générateur congruentiel linéaire est défini par une initialisation X0et une relation de congruence de la forme
Xn+1 = (aXn+ b)[modN] ,
ou` a et b sont deux entiers positifs, et N est un entier positif égal à la plus grande valeur souhaitée.
La période est relativement courte, de l’ordre de la base N choisie. Les valeurs par défaut sont a = 16807, et N = 231 ? 1.
Il existe bien d’autres exemples, notamment
– Les générateurs de Fibonacci, définis par la récurrence
Xn = (Xn?1 + Xn?2)[modN] ,
avec une initialisation adaptée, c’est à dire la donnée de X1 et X2. Il faut que ceux ci soient suffisamment grands pour que les suites de nombres ainsi générées soient de bonne qualité.
– Les générateurs de Fibonacci généralisés, définis par
Xn+1 = (Xn?k + Xn?`+c)[modN] ,
avec c fixé, et une initialisation adaptée, c’est à dire la donnée de X1,X2, Xk (si on suppose k > `).
– Générateurs non-linéaires : par exemple
Xn+1 = Xn(Xn + 1)[mod2k] ,
pour un k assez grand, et une initialisation bien choisie X2 = 2[mod4]
| Exercice 5.1 1. Ecrire une fonction fibogen, qui génère des nombres pseudo-aléatoires uniformément distribués entre 0 et 1 avec l’algorithme de Fibonacci. La fonction prendra en entrée la base N, les valeurs initiales X1 et X2, ainsi que le nombre de valeurs à générer. 2. Ecrire une fonction lincongen, qui génère des nombres pseudo-aléatoires uniformément distribués entre 0 et 1 avec l’algorithme congruentiel linéaire ci-dessus. La fonction prendra en entrée la base N, les entiers a et b, la valeurs initiale X1, ainsi que le nombre de valeurs à générer. 3. Dans les deux cas, étudier la qualité du générateur obtenu. On pourra par exemple calculer un histogramme, et calculer l’erreur de celui-ci par rapport à la distribution attendue, et étudier l’évolution de cette erreur en fonction du nombre N de nombres aléatoires générés. |
Dans cet exercice, il sera utile d’effectuer des tests avec les générateurs ainsi construits, pour différents choix de paramètres, et d’interpréter les résultats obtenus.
Comme on l’a vu, l’instruction rand permet de générer des séquences de nombres pseudo-aléatoires aussi “indépendants” que possible, et distribués selon une loi proche de la loi U([0,1]). Les paramètres sont optimisés pour fournir des nombres d’aussi bonne qualité que possible.
rand propose trois méthodes, correspondant à trois algorithmes différents, que l’on peut choisir en utilisant l’instruction
>> rand(methode,etat)
La variable etat est un entier positif caractérisant l’état courant du système (en gros, l’initialisation du générateur), et la variable methode est l’un des trois choix suivants
– ’state’ : le générateur par défaut, dont la période peut théoriquement atteindre 21492.
– ’seed’ : l’ancien générateur par défaut, basé sur un gégérateur congruentiel multiplicatif, et dont la période vaut 231 ? 2.
– ’twister’ : utilise le twister de Mersenne, l’un des meilleurs générateurs actuels (plus complexe, et donc aussi plus lent que les autres), dont la période vaut (219937 ? 1)/2. L’instruction
>> s = rand(methode);
retourne une chaˆ?ne de caractères contenant l’état du système pour la méthode choisie.
Une fois la méthode et l’état initial choisis, on utilise le générateur comme on l’a déjà vu :
>> X = rand(m,n);
génère une matrice à m lignes et n colonnes contenant des nombres pseudo-aléatoires avec une distribution aussi proche que possible de celle d’un vecteur i.i.d. U([0,1]).
randn génére des matrices pseudo-aléatoires approchant une loi i.i.d. N(0,1). L’instruction
>> X = s*randn(M,N) + m;
ou` m et s sont des réels, génère M × N nombres pseudo-aléatoires i.i.d. suivant approximativement une loi N(m,s2). Pour randn, seuls deux méthodes sont disponibles : ’seed’ et ’state’.
On pourra se documenter aussi sur randperm (permutations aléatoires).
Les méthodes ci-dessus permettent de simuler des distributions uniformes d’entiers entre 0 et N ?1 (ou` N est très grand), et donc par division par N de réels (plutôt des rationnels) dans [0,1]. On va maintenant voir comment utiliser ces derniers pour simuler d’autres distributions.
Pourquoi se fatiguer quand quelqu’un d’autre a déjà fait le travail? La Stixbox, disponible sur http
contient des implémentations de la plupart des densités, des fonctions de répartition, des fonctions quantiles, et des générateurs de nombres aléatoires des lois classiques. La plupart sont basés sur l’inversion de la fonction de répartition, sauf quelques uns qui utilisent le rejet.
Prenons l’exemple de la distribution ?2. Une variable aléatoire X suit une loi ?2 à k degrés de liberté si la densité de X notée fX est :
,
ou` ? est la fonction gamma d’Euler.
Stixbox implémente quatre fonctions reliées à la distribution ?2 :
– rchisq(n,k) : génére n nombres pseudo-aléatoires distribués suivant une loi ?2 à k degrés de liberté.
– dchisq(x,k) : évalue la densité de probabilités d’une variable ?2 à k degrés de liberté au point x. – pchisq(x,k) : évalue la fonction de répartition d’une variable ?2 à k degrés de liberté au point x.
– qchisq(x,k) : évalue la réciproque de la fonction de répartition d’une variable ?2 à k degrés de liberté au point x.
La suite d’instructions suivante permet de générer la figure 5.2.
>> y = rchisq(10000,k);
>> subplot(2,1,1);
>> hist(y,50);
>> x = linspace(0,20,1000);
>> subplot(2,1,2);
>> plot(x,dchisq(x,k));
Fig. 5.2 – Densité de la distribution ?2 à 3 degrés de liberté (bas), et histogramme correspondant (sur 10 000 réalisations, en haut)
On peut également tracer la fonction de répartition et sa réciproque (voir figure 5.3) via les instructions
>> x = linspace(0,20,1000);
>> subplot(2,1,1);
>> plot(x,pchisq(x,k));
>> subplot(2,1,2);
>> x = linspace(0,1,100);
>> plot(x,qchisq(x,k));
La syntaxe est similaires pour un certain nombre d’autres distributions classiques (Fisher, Gamma, Gumbel, ).
Fig. 5.3 – Fonction de répartition d’une loi de ?2 à 3 degrés de liberté (haut), et sa réciproque (bas)
La premi‘ere méthode est basée sur la fonction de répartition de la variable aléatoire que l’on désire simuler, et exploite le résultat suivant :
Lemme 5.2 Soit F la fonction de répartition d’une variable aléatoire sur R, et soit U ? U([0,1]) une variable aléatoire uniforme sur [0,1]. Alors la variable aléatoire V = F?1(U) admet F pour fonction de répartition.
Cette méthode est facile à mettre en oeuvre lorsque la réciproque de la fonction de répartition est connue. Si tel n’est pas le cas, on peut toutefois se baser sur une approximation de celle-ci.
Remarque : Il est possible de réinterpréter la fonction VAdiscrete précédemment développée en termes de cette approche. Supposons en effet que nous voulions simuler une variable aléatoire discrète prenant les valeurs x1 ? x2 ? ··· ? xN, avec probabilités pk,k = 1, N. La fonction de répartition correspondante est comme indiqué en figure 5.4, ou` on pose
.
Fig. 5.4 – Fonction de répartition d’une variable aléatoire discrète
La fonction de répartition est constante par morceaux, et n’est pas inversible à proprement parler. On convient donc d’associer à tout f ? [0,1] la valeur X(f) = xk telle que xk?1 ? f < xk.
Il est facile de vérifier qu’on obtient bien de la sorte la distributiuon désirée.
Exercice 5.2 Ecrire une fonction monranexp.m générant des réalisations d’une variable aléatoire exponentielle de paramètre ?, de densité
??(x) = ?e??x, x ? R+.
On explicitera tout d’abord la fonction de répartition, puis on utilisera celle-ci. Les variables d’entrée de la fonction seront ? et le nombre de nombres aléatoires
à générer.
| Exercice 5.3 Ecrire une fonction générant des réalisations d’une variable aléatoire de Cauchy, de densité . On explicitera tout d’abord la fonction de répartition, puis on utilisera celleci. La variable d’entrée de la fonction sera le nombre de nombres aléatoires à générer. |
Exercice 5.4 Ecrire une fonction monrandn1.m générant un nombre donné de réalisations d’une variable aléatoire Gaussienne de moyenne a et écart-type s donnés, en utilisant la fonction de répartition. On se documentera sur les fonctions erf et erfinv de Matlab .
Dans le cas de variables Gaussiennes (pour lesquelles la réciproque de la fonction de répartition n’a pas d’expression explicite), il existe une autre méthode, basée sur le résultat suivant :
| Lemme 5.3 Soient U,V ? U([0,1]) deux variables aléatoires indépendantes. Soient X et Y définies par Alors X et Y sont deux variables aléatoires normales centrées réduites (c’est à dire X,Y ? N(0,1)) indépendantes. |
Cette transformation vient du fait que, dans un système Cartésien à deux dimensions ou` les coordonnées X et Y sont deux variables aléatoires indépendantes normales, les variables aléatoires R2 et ? sont elles aussi indépendantes et peuvent s’écrire
R2 = ?2lnU ? = 2?V .
Exercice 5.5 Ecrire une fonction monrandn2.m générant un nombre donné de réalisations d’une variable aléatoire Gaussienne de moyenne a et écart-type s donnés, en utilisant l’algorithme de Box-Mu¨ller.
Le principe des méthodes de rejet est d’utiliser une variable aléatoire que l’on sait simuler pour simuler une autre variable aléatoire plus complexe, en sélectionnant les “bonnes réalisations” de la première. Plus précisément, on se base sur le résultat suivant
Proposition 5.1 Soient ? et ?0deux densités de probabilités telles qu’il existe une constante positive K vérifiant
?(x) ? K?0(x) ?x .
Soient X une variable aléatoire de densité ?0, et soit U ? U([0,1]) une variable aléatoire uniformément distribuée entre 0 et 1, indépendante de X.
Soit E l’évènement
E = {KU?0(X) < ?(X)} (])
Alors, la loi conditionnelle de X sachant E a pour densité ?.
De là on tire un algorithme simple pour simuler des réalisations d’une variable aléatoire de densité ? :
– Tirer au hasard X et U suivant la densité ?0 et la loi U(0,1) respectivement.
– Calculer KU?0(X) et ?0(X).
– Si la condition (]) est remplie, conserver la valeur obtenue.
Prenons l’exemple d’une distribution de Cauchy vue précédemment
,
et d’une Gaussienne centrée réduite de densité
| Exercice 5.6 Montrer que . Utiliser cette remarque pour construire un générateur de nombres pseudoaléatoires Gaussiens (centrés réduits) à partir d’un générateur de distribution de Cauchy (le votre, ou celui de la Stixbox), et de la fonction Matlab rand. |
Notons que les méthodes que nous avons utilisées au premier chapitre pour simuler des distributions aléatoires dans un disque on une sphère sont des versions simplifiées de la méthode du rejet. On va maintenant en voir des versions plus efficaces.
Exercice 5.7 Utiliser la méthode du rejet pour simuler un vecteur aléatoire uniformément distribué dans un disque de rayon R0 et de centre (a,b) du plan, à partir de la fonction rand. On pourra commencer par le cas a = b = 0 et R0 = 1, et on utilisera le fait que pour un tel vecteur aléatoire, en passant en coordonnées polaires, les variables radiale R et angulaire ? sont indépendantes, et distribuées respectivement avec une densité ?R(r) = 2r,r ? [0,1] et selon une loi uniforme U(0,2?).
Passons maintenant au cas multivarié. On rappelle qu’étant donnée une matrice C ? M(N,N) symétrique définie positive, et un vecteur b ? RN, la densité de probabilités Gaussienne multivariée associée est la fonction de N variables définie par
,
ou` on a noté h·,·i le produit scalaire dans RN. b est la moyenne, et C est la matrice de variance-covariance : si X = (X1, XN) est un vecteur aléatoire associé à cette distribution, on a
bk = E{Xk} , Ck` = E{(Xk ? bk)(X` ? b`)} , k,` = 1, N .
Pour simuler de tels vecteurs aléatoires Gaussiens, on se base sur la décomposition de Cholesky de la matrice de Variance-Covariance.
Lemme 5.4 (Factorisation de Cholesky) Soit A une matrice symétrique définie positive, il existe au moins une matrice réelle triangulaire inférieure L telle que :
A = LLT
Remarque : On peut également imposer que les éléments diagonaux de la matrice L soient tous positifs, et la factorisation correspondante est alors unique.
Supposons maintenant que W = {W1, WN} soit un vecteur aléatoire Gaussien centré réduit : ses composantes Wk sont des variables aléatoires centrées réduites dindépendantes. Soit C = LLT une décomposition de Cholesky de la matrice définie positive C ? M(N,N), et posons X = LW. Alors
E{XkX`} = XLkmL`nE{WmWn} = XLkmL`m = LLTk` = Ck` .
m,n m
Des combinaisons linéaires de variables aléatoires Gaussiennes étant toujours Gaussiennes, on a donc montré
Lemme 5.5 Soit W = {W1, WN} soit un vecteur aléatoire Gaussien centré réduit, soit B ? RN un vecteur, soit C ? M(N,N) une matrice symétrique définie positive, soit C = LLT une décomposition de Cholesky de C. Alors le vecteur aléatoire défini par
X = LW + B
est un vecteur aléatoire Gaussien centré, de matrice de covariance C.
Fig. 5.5 – 2000 réalisations d’un vecteur aléatoire Gaussien dans le plan
| Exercice 5.8 Ecrire une fonction qui, partant d’une matrice de variancecovariance donnée C ? M(N,N) (symétrique et définie positive), engendre M réalisations d’un vecteur aléatoire Gaussien centré, de matrice de variancecovariance C. Attention : se document sur chol! Pour visualiser les résultats obtenus, on pourra par exemple se placer dans le cas bidimensionnel, et tester avec une matrice diagonale, ou la matrice et tracer dans le plan les points obtenus. On pourra faire de même en trois dimensions. |
Un exemple de visualisation se trouve en Figure 5.5; on peut y remarquer la forme “ellipso¨?dale” du nuage de points ainsi généré.
5.4. UNE APPLICATION : INTEGRATION DE MONTE-CARLO´
Il peut arriver que l’on ait à évaluer numériquement des intégrales difficilement calculables par des méthodes d’intégration déterministes (par exemple, en grande dimension). On peut alors avoir recours à des méthodes d’intégration utilisant des nombres pseudo-aléatoires.
Par exemple, supposons que l’on ait à calculer une intégrale d-dimensionnelle de la forme
,
D ? Rd est un domaine, et f est telle RDf(x)dx = 1. f peut alors être interprétée comme une densité de probabilités, et l’intégrale comme une espérance mathématique. Si on sait générer des vecteurs pseudoaléatoires suivant la loi de f, on peut alors approximer l’intégrale sous la forme
,
ou` les Xn sont N vecteurs pseudo-aléatoires indépendants suivant f. La loi des grands nombres comme 1/ N quand N augmente. Par contre l’erreur est proportionnelle à var{g(Xn)} et il faut simuler Nd variables aléatoires pour obtenir N vecteurs aléatoires de dimension d.
adaptées du site web du cours Matlab /Octave de l’EPFL : http matlab 2Développé par la société The MathWorks
Alternative : utiliser le moteur de recherche de Matlab disponible dans le help.
Inutile car Matlab possède déjà une fonction mean et une fonction var.
à ne pas confondre avec les instructions ldivide et rdivide, qui effectuent des divisions “point par point” de matrices de même taille.
[5] John von Neumann insista sur ce fait avec la remarque suivante : “Quiconque considère des méthodes arithmétiques pour produire des nombres aléatoires est, bien suˆr, en train de commettre un péché”.
[6] Robert R. Coveyou, du Oak Ridge National Laboratory écrivit dans un article que “la génération de nombres aléatoires est trop importante pour être confiée au hasard”.
