News :: Erreur :: Statistiques :: ClearType :: Imprimer ::

MÉTHODES NUMÉRIQUES | FRACTALES | SYSTÈMES LOGIQUES | CODES CORRECTEURS CRYPTOGRAPHIE | AUTOMATES | INFORMATIQUE QUANTIQUE

L'informatique théorique est la branche des sciences qui s'occupe du développement d'algorithmes et d'outils théoriques applicables à l'informatique pour la résolution de problèmes formels lié à la simulation de phénomènes physiques ou au traitement et l'échange d'informations. (Larousse)

Dernière mise-à-jour de ce chapitre: 05.01.2012 13:11
Version: 2.5 Revision 6 | Rédacteur: Vincent Isoz  | Avancement: ~50%
vues depuis le 01.01.2012: 776

LISTE DES SUJETS TRAITÉS SUR CETTE PAGE | DISCUTER DE CETTE PAGE

L'analyse numérique est une discipline des mathématiques. Elle s'intéresse tant aux fondements théoriques qu'à la mise en pratique des méthodes permettant de résoudre, par des calculs purement numériques, des problèmes d'analyse mathématique.

Définition: "L'analyse numérique" est l'étude des algorithmes permettant de résoudre les problèmes de mathématiques continues (distinguées des mathématiques discrètes).

Cela signifie qu'elle s'occupe principalement de répondre numériquement à des questions à variable réelle ou complexe comme l'algèbre linéaire numérique sur les champs réels ou complexes, la recherche de solutions numériques d'équations différentielles et d'autres problèmes liés survenant dans les sciences physiques et l'ingénierie.

Certains problèmes de mathématique continue peuvent être résolus de façon exacte par un algorithme. Ces algorithmes sont appelés alors "méthodes directes". Des exemples sont l'élimination de Gauss-Jordan pour la résolution d'un système d'équations linéaires et l'algorithme du simplexe en programmation linéaire (voir plus loin).

Cependant, aucune méthode directe n'est connue pour certains problèmes (et il est même démontré que pour une classe de problèmes dits "NP complets", il n'existe aucun algorithme fini de calcul direct en temps polynomial). Dans de tels cas, il est parfois possible d'utiliser une méthode itérative pour tenter de déterminer une approximation de la solution. Une telle méthode démarre depuis une valeur devinée ou estimée grossièrement et trouve des approximations successives qui devraient converger vers la solution sous certaines conditions. Même quand une méthode directe existe cependant, une méthode itérative peut être préférable car elle est souvent plus efficace et même souvent plus stable (notamment elle permet le plus souvent de corriger des erreurs mineures dans les calculs intermédiaires).

L'utilisation de l'analyse numérique est grandement facilitée par les ordinateurs. L'accroissement de la disponibilité et de la puissance des ordinateurs depuis la seconde moitié du 20ème siècle a permis l'application de l'analyse numérique dans de nombreux domaines scientifiques, techniques et économiques, avec souvent des effets révolutionnaires.

Lors de simulations numériques de systèmes physiques, les conditions initiales sont très importantes dans la résolution d'équations différentielles (voir les différents chapitres du site ou apparaissent des effets chaotiques). Le fait que nous ne puissions les connaître avec exactitude fait que les résultats des calculs ne peuvent jamais être parfaitement exacts (nous le savons très bien pour la météo qui en est l'exemple connu le plus flagrant). Cet effet est une conséquence des résultats de la physique fondamentale (basée sur des mathématiques pures) qui démontre que l'on ne peut connaître parfaitement un système en y effectuant des mesures puisqu'elles perturbent directement ce dernier (principe d'incertitude de Heisenberg) et elles font l'objet de la théorie du Chaos (classique ou quantique).

Avec les nouveaux outils informatiques à disposition en ce début du 21ème siècle, il est devenu pratique et passionnant de connaître les méthodes numériques afin de s'amuser dans certains logiciels (OpenGL, 3D Studio Max, Maple, Matlab, Mathematica, Comsol, etc.) à simuler sous forme graphique 2D ou 3D des systèmes physiques.

Définition: Un "algorithme" est une suite finie de règles à appliquer dans un ordre déterminé à un nombre fini de données pour arriver, en un nombre fini d'étapes (dont la quantité, ou réciproquement le temps d'exécution est définie par le terme "coût"), à un certain résultat, et cela indépendamment du type de données

Les algorithmes sont intégrés dans des calculateurs par l'intermédiaire de "programmes".

Définition: Un "programme" est la réalisation (l'implémentation) d'un algorithme au moyen d'un langage donné (sur une architecture donnée). Il s'agit de la mise en oeuvre du principe.

Axiomes de la programmation (anecdotique):

A1. Plus nous écrivons de code, plus nous produirons d'erreurs

A2. Il n'existe pas de programmes sans de possibles erreurs (dû au programme lui-même, à l'électronique sous-jacente ou le plus souvent à l'utilisateur même)

Basiquement voici la démarche minimale à suivre lors du développement d'un produit informatique (au niveau du code):

M1. Auditer les besoins présents et anticiper les besoins futurs

M2. Définir les objectifs

M3. Calculer la faisabilité, le risque d'erreur, la durée nécessaire au développement

M4. Créer les algorithmes en langage formel (cela comprend la gestion des erreurs!)

M5. Optimiser la complexité et contrôler les algorithmes

M6. Choix d'une stratégie de développement (modulable ou autre)

M7. Traduire les algorithmes dans la technologie choisie (ce choix est un sujet assez sensible...)

M8. Tests (maintenance préventive)

Lors du développement d'un programme scientifique, il peut être intéressant, voire même rigoureux de s'attarder la notion de "complexité" de ce dernier. Sans aller trop loin, voyons un peu de quoi il s'agit:

COMPLEXITÉ

La complexité d'un algorithme est la mesure du nombre d'opérations fondamentales qu'il effectue sur un jeu de données. La complexité est donc exprimée comme une fonction de la taille du jeu de données.

Les hypothèses permettant un calcul de cette complexité sont:

H1. (les quatre opérations fondamentales ont le même coût)

H2. Un accès mémoire une opération arithmétique

H3. Un contrôle de comparaison une opération arithmétique

H4. Un seul processeur

Nous notons  l'ensemble des données de taille n et T(n) le coût (en temps) de l'algorithme sur la donnée ou le jeu de données de taille n.

- La "complexité au meilleur" est donnée par la fonction:

  (57.1)

C'est le plus petit temps qu'aura exécuter l'algorithme sur un jeu de données (de lignes de code) de taille fixée, ici à n dont le coût (la durée) d'exécution est C(d). C'est une borne inférieure de la complexité de l'algorithme sur un jeu de données de taille n.

- La "complexité au pire":

  (57.2)

C'est le plus grand temps qu'aura à exécuté l'algorithme sur un jeu de données de taille fixée, ici à n. Il s'agit donc d'un maximum, et l'algorithme finira toujours avant d'avoir effectué  opérations. Cependant, cette complexité peut ne pas refléter le comportement usuel de l'algorithme, le pire cas pouvant ne se produire que très rarement, mais il n'est pas rare que le cas moyen soit aussi mauvais que le pire cas.

- La "complexité moyenne":

  (57.3)

Il s'agit de la moyenne des complexités de l'algorithme sur des jeux de données de taille n (en toute rigueur, il faut bien évidemment tenir compte de la probabilité d'apparition de chacun des jeux de données). Cette moyenne reflète le comportement général de l'algorithme si les cas extrêmes sont rares ou si la complexité varie peu en fonction des données. Cependant, la complexité en pratique sur un jeu de données particulier peut être nettement plus importante que la complexité moyenne, dans ce cas la complexité moyenne ne donnera pas une bonne indication du comportement de l'algorithme.

En pratique, nous ne nous intéressons qu'à la complexité au pire, et à la complexité moyenne

Définition: Un algorithme est dit "algorithme optimal" si sa complexité est la complexité minimale parmi les algorithmes de sa classe.

Comme nous l'avons fait sous-entendre précédemment, nous nous intéressons quasi exclusivement à la complexité en temps des algorithmes. Il est parfois intéressant de s'intéresser à d'autres de leurs caractéristiques, comme la complexité en espace (taille de l'espace mémoire utilisé), la largeur de bande passante requise, etc.

Pour que le résultat de l'analyse d'un algorithme soit pertinent, il faut avoir un modèle de la machine sur laquelle l'algorithme sera implémenté (sous forme de programme). Nous prenons habituellement comme référence, la "machine à accès aléatoire (RAM)" et à processeur unique, où les instructions sont exécutées l'une après l'autre, sans opérations simultanées et sans processus stochastiques (contrairement aux possibles machines quantiques futures).

Les algorithmes usuels peuvent être classés en un certain nombre de grandes classes de complexité dont l'ordre O varie d'une certaine manière:

- Les algorithmes sublinéaires dont la complexité est en général de l'ordre O(log(n))

- Les algorithmes linéaires en complexité O(n) et ceux en complexité en O(n log(n)) sont considérés comme rapides.

- Les algorithmes polynomiaux en O(n^k) pour  sont considérés comme lents, sans parler des algorithmes exponentiels (dont la complexité est supérieure à tout en n) que l'on s'accorde à dire impraticables dès que la taille des données est supérieure à quelques dizaines d'unités.

Remarque: Une bonne complexité est du type O(n^k) pour . Une mauvaise complexité est du type ou .

Exemples:

E1. Évaluation d'un polynôme:

  (57.4)

L'évaluation directe de la valeur de P(x) conduit à une complexité

  (57.5)

Nous devons à Horner un algorithme plus performant qui utilise une factorisation du polynôme sous la forme:

  (57.6)

Nous pouvons facilement montrer que cette factorisation maintient le même nombre d'additions mais réduit le nombre de multiplications à (n). La complexité qui en découle est donc O(n). Le gain est incontestablement important. De plus, cette factorisation évite tout calcul de puissance.

E2. Un autre exemple connu de complexité algorithmique est la recherche d'une information dans une colonne triée. Un algorithme simple appelé "recherche dichotomique" consiste à prendre la cellule au milieu de la colonne et de voir si nous y trouvons la valeur recherchée. Sinon, la recherche doit continuer dans la partie supérieure ou inférieure du tableau (cela dépend de l'ordre lexicographique).

L'algorithme est récursif et permet de diminuer par deux, à chaque étape, la taille de l'espace de recherche. Si cette taille, initialement, est de n cellules dans la colonne, elle est de n/2 à l'étape 1,  à l'étape 2, et plus générale à  à l'étape k.

Au pire, la recherche se termine quand il n'y a plus qu'une seule cellule de la colonne à explorer, autrement dit quand k est tel que . Nous en déduisons le nombre maximal d'étapes: c'est le plus petit k tel que , soit , soit:

  (57.7)

à comparer avec une recherche séquentielle dans une colonne de 25'000 données dont la complexité linéaire est O(n) soit 25'000 alors qu'avec la méthode dichotomique, nous avons une complexité sublinéaire . Le gain est donc considérable!

E3. Soit A et B deux matrices carrées de dimension n. Les principales opérations sur ces matrices ont les complexités suivantes (nous laisserons le soin au lecteur de vérifier car c'est vraiment trivial):

- Lecture (itérations): O(n²)

- Calcul de la trace:

- Addition: tel que

- Multiplication: tel que

- Déterminant (par la méthode directe de Cramer). Nous renvoyons le lecteur au chapitre d'Algèbre Linéaire pour le détail des méthodes de calcul du déterminant d'une matrice. Nous pouvons alors montrer que la complexité d'un déterminant d'ordre n est n multiplications, n-1 additions plus n fois la complexité d'un déterminant d'ordre n-1. Par cumul, nous arrivons à:

  (57.8)

En faisant l'hypothèse que l'ordinateur utilisé effectue une opération élémentaire en secondes (ce qui est déjà une bonne machine). Nous obtenons alors les temps de calculs suivants pour plusieurs valeurs de n:

n	5	10	15	20	50
t

     Tableau: 57.1  - Temps de calcul d'un déterminant

d'où la nécessité de faire un calcul de complexité avant de mettre l'algorithme en route (à moins que vous ne travailliez exclusivement pour les générations futures, à condition qu'il y en ait encore...).

Finalement, pour résumer un peu, nous distinguons quelques types de complexités classiques: O(1) indépendant de la taille de la donnée, O(log(n)), complexité logarithmique, O(n) complexité linéaire, O(n log(n)), complexité quasi-linéaire, O(n²), complexité quadratique, O(n³), complexité cubique, O(kⁿ) complexité exponentielle, O(n!), complexité factorielle.

NP-COMPLÉTUDE

Nous allons introduire pour la culture générale le concept de NP-complétude, c'est-à-dire que nous allons tenter de définir sans trop de formalisme (comme à l'habitude) la classe des problèmes dit "NP-complets". Ces problèmes sont ceux pour lesquels personne à l'heure actuelle ne connaît d'algorithme efficace (i.e. seulement des algorithmes exponentiels). Ce sont des problèmes vraiment difficiles par opposition à ceux pour lesquels nous connaissons des algorithmes de complexité polynomiale.

Définitions:

D1. Les problèmes de "classe P" sont de bons problèmes dans le sens où le calcul de leur solution est faisable dans un temps raisonnable avec un algorithme à complexité polynomiale. Plus formellement, ce sont les problèmes pour lesquels nous pouvons construire une machine déterministe (voir la remarque après les définitions) dont le temps d'exécution est de complexité polynomiale (le sigle "P" signifiant "Polynomial Time").

D2. Les problèmes de "classe E" sont des problèmes dans le sens où le calcul de leur solution est faisable dans un temps exponentiel par nature du type . Plus formellement, ce sont les problèmes pour lesquels nous pouvons construire une machine déterministes dont le temps d'exécution est de complexité exponentielle (le sigle "E" signifiant "Exponential Time").

D3. Les problèmes de la "classe NP" sont ceux pour lesquels nous pouvons construire une machine de Turing non déterministe (voir la remarque après les définitions) dont le temps d'exécution est de complexité polynomiale (le sigle "NP" provient de "Nondeterministic Polynomial time" et non de "Non Polynomial").

Il importe de remarquer à ce stade de la discussion que la classe P est incluse dans la classe NP, nous écrivons , car si nous pouvons construire une machine déterministe pour résoudre efficacement (en un temps au pire polynomial) un problème, nous pouvons certainement (du moins dans notre esprit) en construire une non déterministe qui résout aussi efficacement le même problème. Par ailleurs, il ne faut pas croire que ce concept de divination optimale qu'est la machine déterministe permet de tout résoudre puisqu'il existe en informatique théorique d'autres types de problèmes n'appartenant pas à la classe NP qui sont les problèmes indécidables.

Pour savoir si un problème donné appartient ou non à la classe NP, il s'agit simplement de l'exprimer sous la forme dont la réponse est soit OUI, soit NON. Le problème appartient alors à la classe NP si par définition, nous arrivons à démontrer à l'aide d'un algorithme de complexité polynomiale que n'importe quelle instance OUI du problème est bel et bien correcte. Nous n'avons pas à nous préoccuper des instances NON du problème puisque la machine non déterministe, par définition, prend toujours la bonne décision (lorsque celle-ci existe).

Par exemple, le problème consistant à trouver un cycle hamiltonien (cycle qui passe une et une seule fois par tous les sommets du graphe - voir chapitre de Théorie Des Graphes) dans un graphe appartient à NP puisque, étant donné un cycle, il est trivial de vérifier en temps linéaire qu'il contient bien une et une seule fois chaque sommet.

Un autre exemple de problème difficile des mathématiques est la factorisation d'un entier en produit de facteurs premiers. Nous ne savons pas à ce jour s'il existe un algorithme polynomial qui réussisse cette opération. Autrement dit, nous ne savons pas si ce problème est dans P. En revanche, étant donnés des nombres premiers il est trivial de vérifier que : ce problème est donc dans NP. Il semblerait (nous n'avons pas vérifié ce résultat et ni la possibilité de le faire) que la complexité du meilleur algorithme de factorisation en nombres premiers soit du type:

  (57.9)

il resterait donc du travail à faire (si un internaute pouvait nous fournir les détails qui ont amené à ce résultat, nous sommes preneurs).

Parmi l'ensemble des problèmes NP, il en existe un sous-ensemble qui contient les problèmes les plus difficiles: nous les appelons les problèmes "NP-complet" N.P.C. Ainsi, un problème NP-complet possède la propriété que tout problème dans NP peut être transformé en celui-ci en temps polynomial.

La raison essentielle pour laquelle les problèmes NPC sont les plus difficiles parmi les problèmes de NP est que ces premiers peuvent toujours servir comme des sous-routines pour solutionner ces derniers. Cette réduction à un ou des sous-routines assez facilement faisable puisque réalisable, si elle existe, en temps polynomial. Un problème NPC est donc complet en ce sens qu'il contient l'essentiel de la complexité des problèmes appartenant à NP, et qu'une solution polynomiale à ce problème implique une solution polynomiale à tous les problèmes NP.

Autrement formulé, les problèmes NPC ont une complexité exponentielle et ils ont tous la même classe de complexité (modulo les polynômes).

Finalement, ce qu'il importe de bien comprendre et de retenir de toute cette théorie, son idée maîtresse, est que si nous trouvons un jour un algorithme de complexité polynomiale pour un seul de ces problèmes vraiment difficiles que sont les problèmes NPC, alors d'un seul coup NP devient égal à P et tous les problèmes difficiles deviennent faciles !

Pour résumer, être dans P, c'est trouver une solution en un temps polynomial, tandis qu'être dans NP, c'est prouver en un temps polynomial qu'une proposition de réponse est une solution du problème. Ainsi, tout problème qui est dans P se trouve dans NP. Un champ de recherche majeur des mathématiques actuelles est l'investigation de la réciproque: a-t-on P=NP? Autrement dit, peut-on trouver en un temps polynomial ce que l'on peut prouver en temps polynomial?

Passons maintenant à l'étude de quelques applications types des méthodes numériques dont il est très souvent fait usage dans l'industrie. Nous irons du plus simple au plus compliqué et sans oublier que beaucoup de méthodes ne se trouvant pas dans ce chapitre peuvent parfois être trouvées dans d'autres sections du site!

PARTIE ENTIÈRE

Le plus grand entier inférieur ou égal à un nombre réel x est par [x], qui se lit "partie entière de x".

Ainsi, le nombre M est entier si et seulement si [M]=M. De même, le naturel A est divisible dans l'ensemble des naturels par le naturel B si et seulement si:

  (57.10)

Nous notons aussi {x} pour désigner la partie fraction de x; on a ainsi:

  (57.11)

Soit . Alors nous avons les propriétés suivantes:

P1. , ,  où

P2. , lorsque

P3. , si

P4.

P5.  si ,  si

P6. si

P7. Si , alors [x / a] représente le nombre d'entiers positifs inférieurs ou égaux à x qui sont divisibles par a.

Démonstrations:

La première partie de P1 est simplement la définition de [x] sous forme algébrique. Les deux autres parties sont des réarrangements de la première partie. Dans ce cas, nous pouvons écrire

  (57.12)

où .

Pour P2, la somme est vide pour  et, dans ce cas, on adopte la convention selon laquelle la somme vaut 0. Alors, pour , la somme compte le nombre d'entiers positifs n qui sont plus petits ou égaux à x. Ce nombre est évidemment [x].

La démonstration de P3 sera supposée évidente.

Pour prouver P4, nous écrivons:

,   (57.13)

où n et m sont des entiers et où  et . Alors:

  (57.14)

En écrivant , où , nous avons

  (57.15)

où .

Il s'ensuit que:

0 si -1 si   (57.16)

et on obtient la démonstration P5.

Pour démontrer P6, nous écrivons:

  (57.17)

où , et:

  (57.18)

où . Nous obtenons ainsi:

  (57.19)

puisque . Par ailleurs:

  (57.20)

et nous avons ainsi le résultat.

Pour la dernière partie, nous observons que, si  sont tous les entiers positifs  qui sont divisibles par a, il suffit de prouver que . Puisque  , alors:

  (57.21)

c'est-à-dire:

  (57.22)

soit le résultat attendu.

C.Q.F.D.

ALGORITHME D'HÉRON

Soit à calculer la racine carrée:

  (57.23)

Il existe un algorithme dit "algorithme d'Héron" qui permet de calculer la valeur de cette racine carrée.

Démonstration:

  (57.24)

Nous obtenons alors la relation:

  (57.25)

C.Q.F.D.

Exemple:

Soit à calculer:

  (57.26)

Nous prenons :

Dans le cas de la racine cubique, la démonstration est semblable et nous obtenons:

  (57.27)

Signalons encore que le lecteur pourra trouver dans le chapitre de Théorie des Nombres la méthode utilisée pendant l'antiquité (du moins une analogie) et utilisant les fractions continues.

ALGORITHME D'ARCHIMÈDE

Le calcul de la constante universelle "pi" notée  est très certainement le plus grand intérêt de l'algorithmique puisque l'on retrouve cette constante un peu partout en physique et mathématique (nous pouvons vous conseiller un très bon ouvrage sur le sujet).

Nous rappelons que nous n'en avons pas donné la valeur ni en géométrie, ni dans les autres sections de ce site jusqu'à maintenant. Nous allons donc nous attacher à cette tâche.

Nous définissons en géométrie le nombre   dit "pi", quelque soit le système métrique utilisé (ce qui fait son universalité), comme le rapport de la moitié du périmètre d'un cercle par son rayon tel que:

  (57.28)

Nous devons le premier algorithme du calcul de cette constante à Archimède (287-212 av. J.-C.) dont voici la démonstration.

Démonstration:

Soit un n-polygone inscrit dans un cercle:

Figure: 57.1 - Principe illustré de l'algorithme d'Archimède

Le principe de l'algorithme d'Archimède est le suivant: Soit le périmètre d'un polygone régulier de n côtés inscrit dans un cercle de rayon 1/2. Archimède arrive par induction que:

  (57.29)

Nous avons pour périmètre d'un n-polygone:

 et   (57.30)

Avec:

  (57.31)

Donc:

  (57.32)

Il suffit d'un ordinateur ensuite et de plusieurs itérations pour évaluer avec une bonne précision la valeur de .  Évidemment, on utilise l'algorithme d'Héron pour calculer la racine...

C.Q.F.D.

Remarque: Il existe un très grand nombre d'algorithmes pour calculer . Celle présentée ci-dessus, sans être la plus esthétique, est historiquement la première.

CALCUL DU NOMBRE D'EULER

Parmi la constante , il existe d'autres constantes mathématiques importantes qu'il faut pouvoir générer à l'ordinateur (de nos jours les valeurs constantes sont stockées telles quelles et ne sont plus recalculées systématiquement). Parmi celles-ci, se trouve le "nombre d'Euler" noté e (cf. chapitre d'Analyse Fonctionnelle). Voyons comment calculer ce nombre:

Soit la série de Taylor, pour une fonction indéfiniment dérivable f donnée par (cf chapitre sur les Suites Et Séries):

  (57.33)

Comme (cf. chapitre de Calcul Différentiel Et Intégral):

  (57.34)

nous avons:

  (57.35)

Donc en résumé:

  (57.36)

Cette relation donne un algorithme pour calculer l'exponentielle  à un ordre n donnée de précision.

CALCUL DE LA FACTORIELLE (FORMULE DE STIRLING)

Évidemment, la factorielle pourrait être calculée avec une simple itération. Cependant, ce genre de méthode génère un algorithme à complexité exponentielle ce qui n'est pas le mieux. Il existe alors une autre méthode:

Soit, la définition de la factorielle:

  (57.37)

Et d'après les propriétés des logarithmes:

  (57.38)

Si n est très grand (mais alors très...) alors:

  (57.39)

Lorsque , la limite inférieure est négligeable et alors (approximation qui nous est très utile dans le chapitre de Mécanique Statistique):

  (57.40)

Après une petite simplification élémentaire, nous obtenons:

  (57.41)

Cette dernière relation est utile si l'on suppose bien évidemment que la constante d'Euler est une valeur stockée dans la machine...

SYSTÈMES D'ÉQUATIONS LINÉAIRES

Il existe de nombreuses méthodes de résolution de systèmes d'équations linéaires. La plupart d'entre elles ont été mises au point pour traiter des systèmes particuliers. Nous en étudierons une, appelée la "méthode de réduction de Gauss", qui est bien adaptée à la résolution des petits systèmes d'équations (jusqu'à 50 inconnues).

UNE ÉQUATION À UNE INCONNUE

Considérons le cas le plus simple: celui d'une équation à une inconnue:

  (57.42)

a et b sont les coefficients de l'équation et x en est l'inconnue. Résoudre cette équation, c'est déterminer x en fonction de a et b. Si a est différent de 0 alors:

  (57.43)

est la solution de l'équation. Si a est nul et si b est différent de 0 alors l'équation n'admet pas de solutions. Si a et b sont nuls, alors l'équation possède une infinité de solutions.

DEUX ÉQUATIONS À DEUX INCONNUES

Un système (linéaire) de deux équations à deux inconnues s'écrit:

  (57.44)

 sont les coefficients du système d'équations,  et  en sont les inconnues.

Pour résoudre le système, nous procédons comme suit:

À l'aide de manipulations algébriques (addition ou soustraction des différentes égalités entre elles - manipulations autorisées par l'indépendance linéaire des vecteurs-lignes) nous transformons le système en un autre donné par:

  (57.45)

Ensuite, nous résolvons l'équation , ce qui donne:

  (57.46)

Nous en déduisons:

  (57.47)

La transformation entre les deux systèmes:

  (57.48)

s'effectue simplement en multipliant chaque coefficient de la première égalité par  et en soustrayant les résultats obtenus des coefficients correspondants de la seconde égalité. Dans ce cas, l'élément  est appelé "pivot".

TROIS ÉQUATIONS À TROIS INCONNUES

Examinons encore le cas des systèmes de trois équations à trois inconnues:

  (57.49)

Nous pouvons par la suite des opérations élémentaires (cf. chapitre d'Algèbre Linéaire) réduire le système linéaire précédent en le système suivant:

  (57.50)

et ensuite résoudre l'équation:

  (57.51)

puis la deuxième:

  (57.52)

et enfin:

  (57.53)

Revenons à la transformation des systèmes. Elle s'effectue en deux étapes:

1. Dans la première, nous choisissons  comme pivot et nous éliminons les coefficients  et  de la manière suivante: il faut multiplier chaque coefficient de la première égalité par  et soustraire les résultats obtenus de la deuxième égalité, ainsi  devient nul. De même, en multipliant les coefficients de la première équation par  et en soustrayant les résultats obtenus de la troisième égalité,  disparaît. Le système d'équation s'écrivant alors:

  (57.54)

2. La deuxième étape consiste à traiter le système de deux équations à deux inconnues formé des deuxième et troisième égalités du système précédent, et ce, en choisissant comme pivot. Cette méthode de résolution peut paraître compliquée, mais elle a l'avantage de pouvoir être généralisée et être appliquée à la résolution de systèmes de n équations n inconnues.

N ÉQUATIONS À N INCONNUES

Pour simplifier l'écriture, les coefficients seront toujours notés  et non pas , etc. lors de chaque étape du calcul.

Soit le système linéaire (nous pourrions très bien le représenter sous la forme d'une matrice augmentée afin d'alléger les écritures):

  (57.55)

Il faut choisir  comme pivot pour éliminer . Ensuite, l'élimination de  s'effectue en prenant  comme pivot. Le dernier pivot à considérer est bien évidemment , il permet d'éliminer . Le système prend alors la forme:

  (57.56)

En résolvant d'abord la dernière équation, puis l'avant-dernière et ainsi de suite jusqu'à la première.

Cette méthode, appelée "méthode de résolution de Gauss" ou encore "méthode du pivot" doit cependant être affinée pour éviter les pivots de valeur 0. L'astuce consiste à permuter l'ordre dans lequel sont écrites les équations pour choisir comme pivot le coefficient dont la valeur absolue est la plus grande. Ainsi, dans la première colonne, le meilleur pivot est l'élément  tel que:

  (57.57)

Il est amené à  par échange des première et j-ème lignes. L'élimination du reste de la première colonne peut alors être effectuée. Ensuite, on recommence avec les n-1 équations qui restent.

POLYNÔMES

L'ensemble des polynômes et à coefficients réels a été étudié dans le chapitre d'Analyse Fonctionnelle en détails. Nous allons ici traiter de l'aspect numérique de quelques problèmes liés aux polynômes.

Mis à part l'addition et la soustraction de polynômes que nous supposerons comme triviaux (optimisation de la complexité mis à part comme le schéma de Horner), nous allons voir comment multiplier et diviser deux polynômes.

Voyons d'abord comment multiplier deux polynômes:

Soit:

  (57.58)

alors:

  (57.59)

avec:

  (57.60)

C'était simple...

Un tout petit peu plus difficile maintenant: la division euclidienne des polynômes (cf. chapitre de Calcul Algébrique).

Reprenons:

  (57.61)

mais en imposant cette fois-ci .

La division s'écrira donc nous le savons:

  (57.62)

avec ou sinon .

Il est connu d'avance que nous avons bien évidemment:

et deg(r(x))<m.

Nous avons donc par définition q(x) qui est le quotient de la division et r(x) le reste de la division euclidienne de f(x) par g(x).

Rien ne nous interdit donc de poser de la manière la plus générale qui soit:

  (57.63)

Exemple:

Soit:

  (57.64)

donc:

  (57.65)

Nous avons donc:

  (57.66)

Ensuite:

  (57.67)

et enfin:

  (57.68)

Donc de manière générale:

  (57.69)

comme:

  (57.70)

le premier reste est donc:

  (57.71)

Ensuite:

  (57.72)

Le deuxième reste est alors:

  (57.73)

etc.

Nous continuons jusqu'à ce que .

RÉGRESSIONs ET INTERPOLATIONS

Les régressions (ou "interpolations") sont des outils très utiles aux statisticiens, ingénieurs, informaticiens souhaitant établir une loi de corrélation entre deux (ou plus) variables dans un contexte d'études et d'analyse ou d'extrapolation.

Il existe un grand nombre de méthodes d'interpolation: de la simple résolution d'équations du premier degré (lorsque uniquement deux points d'une mesure sont connus) aux équations permettant d'obtenir à partir d'un grand nombre de points des informations essentielles à l'établissement d'une loi (ou fonction) de régression linéaire, polynomiale ou encore logistique.

RÉGRESSION LINÉAIRE À UNE VARIABLE EXPLICATIVE

Nous présentons ici plusieurs algorithmes (méthodes) utiles et connus dans les sciences expérimentales (nous en avons déjà parlé lors de notre étude des statistiques). L'objectif ici est de chercher à exprimer la relation linéaire entre deux variables x et y indépendantes par un "modèle linéaire" (ML) le plus simple possible (sinon quoi il faudrait des centaines de pages pour présenter le sujet!).

- x est la variable indépendante ou "explicative" appelée aussi "covariable" ou "variable prédictive". Les valeurs de x sont fixées par l'expérimentateur et sont supposées connues sans erreur

- y est la variable dépendante ou "expliquée" (exemple: réponse de l'analyseur). Les valeurs de y sont entachées d'une erreur de mesure. L'un des buts de la régression sera précisément d'estimer cette erreur.

Nous cherchons une relation de la forme:

  (57.74)

C'est l'équation d'une droite (fonction affine), d'où le terme de "régression linéaire" avec a qui est appelé dans ce cadre d'étude: "coefficient de régression".

Dans la vie réelle, les relations linéaires constituent vraiement une exception. Objectivement, on ne s'y intéresse que dans les salles de cours et les livres de cours parce qu'elles sont plus faciles à comprendre. La majorité des relations sont non linéaires dans la réalité et même non continues dans certaines situations... de plus, ce n'est pas parce qu'elles sont linéaires dans l'intervalle de mesures effectuées qu'elles le sont à plus petite échelle ou à plus grande échelle!

DROITE DE RÉGRESSION

Il existe aussi une autre manière commune de faire une régression linéaire du type:

  (57.75)

qui consiste à se baser sur les propriétés de la covariance et de l'espérance (cf. chapitre de Statistiques) et très utilisée entre autres en finance (mais aussi dans n'importe quel domaine où l'on fait un peu de statistiques).

Soit x, y deux variables dont l'une dépend de l'autre (souvent c'est y qui dépend de x) nous avons selon la propriété de linéarité de la covariance qui est, rappelons-le:

  (57.76)

la relation suivante:

  (57.77)

Il vient donc pour la pente (nous réutiliserons cette relation lors de l'étude du rendement d'un portefeuille selon le modèle de Sharpe dans le chapitre d'Économétrie):

  (57.78)

Soit sous la forme plus explicite que nous utiliserons plus loin (cf. chapitre de Statistiques):

  (57.79)

Pour déterminer l'ordonnée à l'origine nous utilisons les propriétés de l'espérance démontrées dans le chapitre de Statistiques:

  (57.80)

ce qui donne b sous la forme:

  (57.81)

MÉTHODE DES MOINDRES CARRÉS

Du fait de l'erreur sur y, les points expérimentaux, de coordonnées , ne se situent pas exactement sur la droite théorique. Il faut donc trouver l'équation de la droite expérimentale qui passe le plus près possible de ces points.

La "méthode des moindres carrés" (DMC) consiste à chercher les valeurs des paramètres a, b qui rendent minimale la somme des carrés des écarts e_i résiduels (SS_r: Sum of Squared Residuals) entre les valeurs observées  et les valeurs calculées théoriques de :

  (57.82)

où n est le nombre de points et:

  (57.83)

d'où autrement écrit:

  (57.84)

Cette relation fait apparaître la somme des carrés des écarts comme une fonction des paramètres a,b. Lorsque cette fonction est minimale (extrémale), les dérivées par rapport à ses paramètres s'annulent:

  (57.85)

Remarque: Cette méthode de recherche de minimum (optimisation) est nommée "méthode des multiplicateurs de Lagrange" dans le monde de l'économétrie. Dans notre exemple est la grandeur scalaire qui fait office de multiplicateur de Lagrange.

Soit après simplification:

  (57.86)

Le système ci-dessus est dit appelé "système des équations normales".

C'est un système linéaire de deux équations à deux inconnues. Notons pour simplifier:

  (57.87)

Le système devient :

  (57.88)

De la deuxième équation, nous tirons :

  (57.89)

En remplaçant dans la première, nous obtenons :

  (57.90)

De là nous avons :

  (57.91)

Ainsi, l'expression de la pente et de l'ordonnée à l'origine de l'équation recherchée est:

  (57.92)

Nous avons donc maintenant pour la droite des moindres carrés:

  (57.93)

soit sous forme discrète:

  (57.94)

ainsi que par construction de la méthode la relation suivante:

  (57.95)

Maintenant, nous faisons l'hypothèse que chaque valeur mesurée est entachée d'un résidu tel que:

  (57.96)

Soit en soustrayant les deux dernières relations:

  (57.97)

Maintenant, passons par un résultat intermédiaire. Rappelons que nous avons obtenu plus haut:

  (57.98)

En remplaçant b par sa valeur:

  (57.99)

nous avons alors:

  (57.100)

Multipliant la deuxième relation ci-dessus par  et retranchant de la première, nous obtenons:

  (57.101)

Soit après réarrangement:

  (57.102)

Revenons maintenant à:

  (57.103)

Si nous mettons le tout au carré et en sommant pour toutes les observations, nous avons:

  (57.104)

soit:

  (57.105)

Or, nous venons de montrer avant que le double produit était nul. Donc:

  (57.106)

Cette dernière relation est appelée "équation d'analyse de la variance". En fait, il s'agit de sommes de carrés. Il faudrait diviser par n pour obtenir des variances biaisées.

Cette dernière relation s'écrit aussi souvent:

  (57.107)

où SCT est la "somme des carrés totale", SCE la "somme des carrés expliquée" et SCR la "somme des carrés résiduelle".

Cette dernière relation se trouve également souvent sous la forme suivante dans la littérature:

  (57.108)

Notons maintenant les  sans erreurs d'une manière différente et appelons cela le "modèle linéaire a priori":

  (57.109)

Il est effectivement important dans la pratique de différencier le modèle a priori qui ne prend pas en compte les erreurs du modèle réel entaché d'erreurs!

Nous avons alors:

  (57.110)

ce qui se représente graphiquement par:

Figure: 57.2 - Représentation graphique de SCT, SCE, SCR

La dernière relation est souvent noté dans la littérature sous la forme plus pédagogique suivante:

  (57.111)

qui est simplement une autre manière d'écrire:

  (57.112)

et il vient alors immédiatement la relation utilisée parfois dans la pratique pour calculer les résidus (connaissant les valeurs calculées et les valeurs mesurées):

  (57.113)

Il est important de se rappeler (ou de constater) que les relations ci-dessus entre SCT, SCE et SCR ne sont valables que dans le cas d'un modèle linéaire!

Rappelons maintenant que dans le chapitre de Statistiques, nous avions démontré que le coefficient de corrélation s'exprimait par:

  (57.114)

ou sinon puisque comme démontré plus (se rappeler que la variance indiquée est la variance biaisée... il faut faire attention aux notations):

  (57.115)

nous pouvons aussi écrire le coefficient de corrélation sous la forme:

  (57.116)

Donc nous en déduisons:

  (57.117)

Montrons que ceci est égal à (notation souvent utilisée dans la littérature spécialisée):

  (57.118)

Démonstration:

Nous partons donc de:

  (57.119)

et puisque nous avons montré que:

  (57.120)

Donc:

  (57.121)

C.Q.F.D.

et dans le cadre des modèles de régression voici quelques cas typiques de la valeur de ce coefficient:

Figure: 57.3 - Quelques valeurs du coefficient de corrélation linéaire (source: Wikipédia)

Enfin, indiquons que nous retrouvons aussi très souvent le coefficient de corrélation linéaire sous la forme suivante dans les logiciels et la littérature:

  (57.122)

Cette dernière forme met mieux en évidence que si la somme des carrés des résidus SCR sont nuls, le modèle est parfaitement modélisable par une relation linéaire dans l'intervalle d'étude considéré.

MODÈLE LINÉAIRE GAUSSIEN

Nous admettrons que, pour un individu i prélevé au hasard dans la population,  est connu sans erreur, et que  est une réalisation d'une variable aléatoire que nous noterons dorénavant  et la droite théorique des moindres carrés s'écrira maintenant:

  (57.123)

où  est par hypothèse un résidu identiquement distribué et indépendant (pas d'auto-corrélation) pour chaque point i selon une loi Normale centrée (de moyenne nulle et d'écart-type égal pour tout k) tel que:

 et   (57.124)

donc:

  (57.125)

où nous avons le résidu qui est donc donné par la différence entre l'ordonnée théorique et l'ordonnée mesurée:

  (57.126)

et puisque par hypothèse , il vient immédiatement que:

  (57.127)

Les hypothèses précédentes concernant les moments des résidus sont appelées "hypothèses de Gauss-Markov" et l'hypothèse particulière d'égalité des variances s'appelle comme nous l'avons vu dans le chapitre de Statistique "l'homoscédasticité" (ou le fait que les variances ne soient pas égales s'appelle pour rappel "l'hétéroscédasticité").

Nous avons de par la propriété de l'espérance (cf. chapitre de Statistiques):

  (57.128)

Alors sous les hypothèses ci-dessus, nous allons montrer que a et b sont des estimateurs sans biais (cf. chapitre de Statistiques) de  et  et qu'il est possible d'estimer l'écart-type à partir de SCR ce qui est un résultat important!

Conformément au modèle adopté, a est à considérer maintenant comme une réalisation de la variable aléatoire donnée par (démontré plus haute comme étant le rapport de la covariance et de la variance):

  (57.129)

et b comme une réalisation de la variable aléatoire donnée par:

  (57.130)

Tenant compte de ce que:

  (57.131)

nous pouvons mettre A sous la forme:

  (57.132)

et B:

  (57.133)

Nous en déduisons les espérances pour A:

  (57.134)

ce qui montre que notre hypothèse initiale pour l'expression de A n'est pas trop mauvais.

Respectivement, pour B il vient:

  (57.135)

Donc A et B sont bien des estimateurs sans biais de . Comme ce sont des estimateurs, dans la littérature spécialisée, A et B sont souvent notés et dès lors il vient la notation courante alternative:

  (57.136)

Nous devons enfin calculer les variances de A et B en utilisant ces propriétés (cf. chapitre de Statistiques) et les hypothèses sur les résidus, nous avons:

  (57.137)

Comme par hypothèse nous avons tous les   qui sont égaux, et qu'il n'y a pas d'auto-corrélation, nous pouvons alors écrire:

  (57.138)

Soit:

  (57.139)

Avant de déterminer la variance de B, rappelons que par hypothèse:

  (57.140)

Donc il en découle immédiatement (cf. chapitre de Statistiques):

  (57.141)

Dès lors:

  (57.142)

En rappelant la relation de Huyghens (cf. chapitre de Statistiques):

  (57.143)

Nous avons finalement:

  (57.144)

Le problème réside maintenant dans la détermination de . Évidemment pour ce faire nous allons être obligés de passer par un estimateur statistique.

Nous savons que nous pouvons écrire selon ce qui a été vu dans le chapitre de Statistique en ce qui concerne les estimateurs:

  (57.145)

puisque la loi Normale est centrée pour les résidus donc ... et que le résidu est une variable aléatoire implicitement dépendante de la somme de deux variables aléatoires que sont A et B d'où la minoration de deux fois l'erreur-standard.

Indiquons aussi que dans la pratique nous notons fréquemment ce dernier résultat en mélangeant les notations de l'aspect aléatoire et déterministe:

  (57.146)

où SEE signifie en anglais "Standard Error of Estimate" (Erreur Standard de l'Estimation). Il s'agit en français de "l'erreur standard de la régression".

Nous avons donc les estimateurs non biaisés des variances de A et de B:

  (57.147)

relations qui ne sont donc valable que pour une régression linéaire à une variable explicative (et sous les hypothèses de construction du modèle linéaire Gaussien).

Ce qui est sympa connaissant ces variances, c'est que nous pouvons aussi estimer la variance de la variable expliquée de notre régression facilement (en utilisant les propriétés de la variance vues dans le chapitre de Statistiques).

Il serait intéressant aussi de faire de l'inférence statistique sur l'espérance des paramètres A et B (donc la pente et l'ordonnée à l'origine) étant donné leur espérance empirique connue. Dès lors, nous avons démontré dans le chapitre de Statistique l'intervalle de confiance suivant:

  (57.148)

Il s'ensuit en faisant un parallègle digne de l'ingénieur... que comme A est un estimateur non biaisé de la moyenne de la pente a et que:

  (57.149)

est en fait l'erreur standard de la moyenne A, alors par analogies:

  (57.150)

et alors (c'est à prendre avec des pincettes et il vaut mieux utiliser les développements qui vont suivra par la suite):

  (57.151)

ce qui donne donc l'intervalle de confiance de la pente d'un modèle linéaire Gaussien à une variable explicative (c'est ce que donne MS Excel pour chaque coefficient). La démarche est la même pour l'ordonnée à l'origine et dans le cas d'une régression linéaire ayant plusieurs variables explicatives, assimilées alors au concept de "degrés de liberté" (d.d.l.) ou en anglais de "degrees of freedom" (df), nous avons:

  (57.152)

et là encore, c'est l'intervalle de confiance qu'utilisent beaucoup de logiciels de statistiques (dont la tableur MS Excel) pour chaque coefficient d'une régression linéaire ayant df variables explicatives.

INTERVALLE DE CONFIANCE DES VALEURS PRÉDICTIVES

Nous souhaiterions pour chaque valeur mesurée de la variable expliquée, connaître l'intervalle de confiance. En d'autres termes, nous souhaiterions connaître la variance de l'estimateur statistique de Y :

  (57.153)

Malheureusement, nous allons nous casser les dents, car la covariance est difficile à calculer (a et b ne sont pas indépendants comme le montrent les expressions que nous avons obtenues plus haut).

Par contre, en étant bon observateur, nous voyons que si nous utilisons le résultat vu plus haut:

  (57.154)

alors:

  (57.155)

Le problème étant contourné, nous avons maintenant en utilisant les propriétés de la variance:

  (57.156)

d'où:

  (57.157)

Donc nous avons:

  (57.158)

Maintenant, rappelons (cf. chapitre de Statistiques) que:

  (57.159)

Dès lors:

  (57.160)

Il est important de remarquer que dans la relation précédente, y est considéré comme la moyenne de la variable expliquée! Dans la pratique il faut aussi vérifier que ce rapport suit effectivement une loi Normale selon un modèle Gaussien pour pouvoir faire les intervalles de confiance et tests statistiques qui s'en suivente.

Rappelons maintenant que nous avons démontré dans le chapitre de Statistique que:

    (57.161)

suit donc une loi de Student de degré de liberté k et U une loi du khi-deux de degré de liberté k.

Maintenant revenons à l'expression du Z précédemment obtenu et rappelons que:

  (57.162)

Donc:

  (57.163)

Or comme:

  (57.164)

Alors:

  (57.165)

et donc:

  (57.166)

Correspond à une somme de carrés lois Normales centrées réduites. Et donc conformément à ce que nous avions démontré dans le chapitre de Statistique il vient que:

  (57.167)

Donc:

  (57.168)

Soit:

  (57.169)

C'est une raison pour laquelle beaucoup de statisticiens notent directement et sans détours:

  (57.170)

Ce qui n'est pas forcément évident au premier coup d'oeil. Raison pour laquelle, suite à la demande d'un lecteur, nous avons détaillé un peu de manière exagérée, le mécanisme qui se cache derrière cette implication.

Bon ceci étant, fait, nous avons donc:

  (57.171)

et il en découle en bilatéral, un intervalle de confiance à un niveau  donné par:

  (57.172)

Remarque: Le même type de développement peut être fait pour la pente et l'ordonnée à l'origine. Raison pour laquelle des logiciels comme MS Excel, SPSS, Minitabl, Statistica, etc. donnent la valeur de la distribution de Student ainsi que l'intervalle de confiance à un niveau de choisi. Mais pour que cela ait un sens, il faut que toutes les hypothèses du modèle construit soient satisfaites.

Raison pour laquelle de nombreux logiciels de statistiques donnent le graphique suivant lors de régressions linéaires (nous y voyons bien l'intervalle de confiance).:

Figure: 57.4 - Capture du logicie Minitab 15 pour l'intervalle de confiance

Le lecteur aura remarqué que:

1. Il est fort pénible sans logiciels d'obtenir l'intervalle de confiance pour la droite des moindres carrés tracé puisqu'il faudrait le calculer pour chaque point...

2. L'intervalle de confiance est courbe ce qui est parfois considéré comme du bon sens, du moins dans la version temporelle de la régression : plus l'échéance de la prévision est lointaine, moins elle est sûre.

La valeur vraie de Y est donc donnée par:

  (57.173)

avec la variance:

  (57.174)

qui est indépendant de l'estimateur Y. Dès lors, la variance de la différence entre Y et y (donc entre estimateur et réel) a pour variance:

  (57.175)

Dès lors il est d'usage de considérer que l'intervalle de prédiction (à ne pas confondre avec l'intervalle de confiance de l'estimateur) est pris comme étant:

  (57.176)

Ce qui nous donne avec un logiciel comme Minitab 15:

Figure: 57.5 - Capture du logiciel Minitab 15 pour l'intervalle de confiance et de prédiction

où nous distinguons bien l'intervelle de confiance (en rouge) de l'intervalle de prévision (en vert).

RÉGRESSION LINÉAIRE À UNE VARIABLE EXPLICATIVE FORCÉE PAR L'ORIGINE

Un cas très fréquent et demandé dans les laboratoires (et globalement dans d'autres départements) des entreprises, est de forcer la régression linéaire à passer par l'origine.

Nous allons voir ici que l'approche est seulement une variante simplifiée de la méthode des moindres carrés.

Nous utilisons comme précédemment:

  (57.177)

où n est le nombre de points. Mais cette fois-ci, nous écrivons:

  (57.178)

d'où autrement écrit:

  (57.179)

Cette relation fait apparaître la somme des carrés des écarts comme une fonction du paramètre a. Lorsque cette fonction est minimale (extrémale), les dérivées par rapport à ses paramètres s'annulent:

  (57.180)

Soit après simplification:

  (57.181)

Enfin:

  (57.182)

Vous pouvez vérifier aussi avec un tableur quelconque (MS Excel par exemple) que les calculs correspondent bien.

RÉGRESSION LINÉAIRE MULTIPLE

Bien évidemment, dans certaines situations la régression linéaire est trop simpliste ou simplement pas adaptée. Le cas le plus typique qui nous concerne dans ce qui va suivre étant les situations où nous avons plusieurs variables explicatives.

Le principe de la régression linéaire multiple est relativement simple. Nous voulons déterminer la variable expliquée y à partir de p-1 variables explicatives reliées par une loi linéaire de la forme générale:

  (57.183)

Dans un échantillon de n individus, nous mesurons  pour :

Au fait, pour estimer les paramètres  (valeurs estimées que nous noterons  afin de respecter cette fois-ci la tradition) l'approche est très simple car elle consiste juste en une généralisation de la méthode des moindres carrés que nous avons vu plus haut pour la régression linéaire simple.

Donc en fin de compte nous récrivons la relation du Sum of Squared Residuals vu plus haut en modifiant un tout petit peu puisque maintenant nous avons du multilinéaire:

  (57.184)

avec:

  (57.185)

Donc il nous faut minimiser:

  (57.186)

Le système peut se réécrire:

  (57.187)

Nous avons le vecteur des résidus qui vaut donc:

  (57.188)

Comme nous le savons, la méthode des moindres carrés consiste à trouver le vecteur  qui minimise:

  (57.189)

Soit explicitement:

  (57.190)

a remarquer que nous avons:

  (57.191)

et:

  (57.192)

puisque chacun des éléments de la multiplication est un simple vecteur.

Nous avons alors (attention! ne pas oublier que certaines multiplications dans la relation à suivre sont des produits scalaires!!!) la fonction quadratique multivariée de coefficients de vecteur:

  (57.193)

Dérivons cette dernière "fonction objet F" à l'ordre du vecteur  (il s'agit d'une dérivée intérieure composante par composante). Ce que nous écrivons:

  (57.194)

Maintenant, passons d'une écriture vectorielle à une écriture matricielle pure:

  (57.195)

Nous cherchons donc  qui annule cette dérivée. Donc nous devons résoudre l'équation suivante:

  (57.196)

Soit:

  (57.197)

Rappelons avant d'aller plus loin que:

  (57.198)

Donc:

  (57.199)

Puisque l'algèbre linéaire est associative, écrivons sans les parenthèses:

  (57.200)

Nous ne pouvons évidemment pas simplifier à droite et gauche par  car comme il ne s'agit pas d'une matrice carrée, ce terme est obligatoirement pas inversible. La seule chose que nous puissions faire c'est identifier les éléments tels que:

  (57.201)

impose forcément que:

  (57.202)

Nous trouvons alors que si la matrice carrée  est inversible alors:

  (57.203)

Pour montrer que cela semble juste, retrouvons les résultats de la régression linéaire simple:

  (57.204)

donc avec . Nous avons alors:

  (57.205)

En utilisant la relation démontrée dans le chapitre d'Algèbre Linéaire pour calculer en toute généralité l'inverse d'une matrice en :

  (57.206)

Nous avons dans le cas présent:

  (57.207)

Nous avons donc:

  (57.208)

vu que nous avons une matrice carrée de dimension 2 seulement, le calcul des quatre déterminants se réduit au fait à sélectionner les composantes de :

  (57.209)

Donc:

  (57.210)

Donc à un changement de notation près pour les indices et les mesures expérimentales, nous retrouvons bien le résultat que nous avions obtenu lors de notre étude de la régression linéaire simple qui était (pour rappel...):

  (57.211)

Maintenant, il nous faut un indicateur de qualité en ce qui concerne notre régression multilinéaire. Rappelons que dans le cadre de notre étude de la régression linéaire à une variable explicative, nous avons démontré que le coefficient de corrélation linéaire pouvait être écrit sous la forme:

  (57.212)

et au fait celui-ci est applicable aussi directement à la régression linéaire multiple puisque qu'il ne présuppose par le nombre de variables explicatives!!

RÉGRESSION LOGISTIQUE (LOGIT)

Bien souvent, les données statistiques disponibles sont relatives à des caractères qualitatifs. Or, comme nous allons le voir, les méthodes d'inférence traditionnelles ne permettent pas de modéliser et d'étudier ce type de caractères. Des méthodes spécifiques doivent être utilisées tenant compte par exemple de l'absence de continuité des variables traitées ou de l'absence d'ordre naturel entre les modalités que peut prendre le caractère qualitatif. Ce sont ces méthodes spécifiques les plus usuelles qui seront l'objet du texte qui suit.

Comme nous l'avons vu plus haut, la régression linéaire simple a donc pour but de modéliser la relation entre une variable dépendante quantitative et une variable explicative quantitative.

Lorsque la "variable de classe" Y à expliquer est binaire (oui-non, présence-absence,0-1,etc.) nous approchons dans un premier temps celle-ci par une fonction de probabilité , qui nous donne à l'opposé la probabilité d'appartenir à la classe , que nous nommerons "régression logistique" ou "régression logit" ou encore "régression binomiale" (très souvent utilisée dans le cadre des réseaux de neurones formels que nous verrons plus loin) . Ensuite, dans une deuxième étape, nous définissons pour un cas binaire une valeur "cutoff". Par exemple, si nous prenons un cutoff de 0.5 alors les cas pour lesquels  appartiendront à la classe 1 (et inversement dans le cas contraire).

Considérons par exemple la variable dichotomique: fin des études. Celle-ci prend deux modalités (en cours, a fini). L'âge est une variable explicative de cette variable et nous cherchons à modéliser la probabilité d'avoir terminé ses études en fonction de l'âge.

Exemple:

Pour construire le graphique suivant, nous avons calculé et représenté en ordonnées, pour des jeunes d'âge différent, le pourcentage de ceux qui ont arrêté leurs études.

Figure: 57.6 - Partie de la population en études en fonction de l'äge

Mais comment obtenons-nous pareil graphique avec une variable dichotomique??? Au fait c'est relativement simple... Imaginez un échantillon de 100 individus. Pour ces 100 individus supposez pour un âge donné que 70% "a fini" et 30% "en cours". Eh bien la courbe représente simplement la proportion des deux classes pour l'âge donné. Il est même parfois indiqué les grosseurs des classes avec des cercles sur toute la longueur des asymptotes horizontales pour bien signifier qu'il s'agit d'une variable dichotomique.

Les points sont distribués selon une courbe en S (une sigmoïde): il y a deux asymptotes horizontales car la proportion est comprise entre 0 et 1. Nous voyons immédiatement qu'un modèle linéaire serait manifestement inadapté.

Cette courbe évoquera pour certains, à juste titre, une courbe cumulative représentant une fonction de répartition (d'une loi Normale par exemple, mais d'autres lois continues ont la même allure). Ainsi, pour ajuster une courbe à cette représentative, nous pourrions nous orienter vers la fonction de répartition d'une loi Normale, et au lieu d'estimer les paramètres a et b de la régression linéaire, nous pourrions estimer les paramètres  de la loi Normale (qui est très similaire à la loi logistique comme nous le démontrerons plus loin). Nous parlons alors d'un "modèle Probit".

La loi qui nous intéresse cependant est donc la loi logistique. Contrairement à la loi Normale, nous savons calculer l'expression de sa fonction de répartition dichotomique (probabilité cumulée) qui est du type (c'est son premier avantage!):

  (57.213)

pour une variable de prédiction x où P est donc la probabilité d'avoir un 1. Nous voyons immédiatement que cette dernière relation étant la primitive de la fonction de distribution (voir la démo un peu plus bas), que prenant x de moins l'infini à plus l'infini que nous avons bien 1. Il s'agit donc bien d'une fonction de répartition!

S'il y a plusieurs variables prédictives nous avons alors:

  (57.214)

Lorsque nous optons pour cette fonction de répartition de la loi logistique, nous obtenons le modèle de régression logistique ou "modèle Logit" pour le choix de la "fonction de lien" et c'est là son deuxième avantage le plus important: nous pouvons faire des statistiques sur une variable binaire comme si nous faisions une simple régression linéaire multiple!

Ainsi, nous estimerons la probabilité cumulée d'avoir fini ses études pour un individu d'âge x par (il existe plusieurs manières d'écrire cette loi suivant les habitudes et le contexte):

  (57.215)

il en découle la fonction de distribution:

  (57.216)

Nous pouvons calculer aussi l'espérance de la fonction de distribution en appliquant ce qui a déjà été vu au chapitre de Statistiques mais une partie de cette intégrale ne peut être résolue que numériquement par contre... si nous posons:

  (57.217)

comme étant la variable aléatoire alors nous pouvons calculer formellement:

  (57.218)

qui après une intégration numérique donne 0. Nous obtenons alors aussi le résultat suivant:

  (57.219)

Ainsi, nous voyons que si nous posons:

  (57.220)

Nous retombons sur une fonction de répartition ayant parfaitement les mêmes caractéristiques qu'une loi Normale centrée réduite (moyenne nulle et variance unitaire).

Exemple:

La fonction sigmoïde (de répartition) est présentée ci-dessous pour  :

Figure: 57.7 - Illustration de la fonction sigmoïde

Les paramètres a, b sont ajustés selon le principe du maximum de vraisemblance (cf. chapitre de Statistiques). De plus, ces paramètres doivent généralement être ajustés de manière itérative, à l'aide d'un programme auquel nous fournissons des valeurs initiales, et qui optimise ces valeurs de manière récurrente (nous n'entrerons pas dans ces détails qui dépassent le cadre théorique de ce site à ce jour).

La dernière relation:

  (57.221)

peut par ailleurs être transformée de la façon suivante :

  (57.222)

d'où:

  (57.223)

Ce que certains écrivent aussi...:

  (57.224)

Le résultat de cette dernière transformation est appelé "logit". Il est égal au logarithme de "l'odds" P/1-P.

Donc lorsque les coefficients a et b ont été déterminés, l'expression précédente permet de déterminer P connaissant x facilement (il s'agit de résoudre une équation linéaire) et inversement! Par ailleurs, puisque x est une variable dichotomique les coefficients sont très facilement interprétables.

Revenons un peu sur l'odds car il est possible d'introduire la notion de fonction logistique en faisant la démarche inverse de celle présentée ci-dessus (soit de commencer par la définition de l'odds pour aller jusqu'au logit) et ceci peut parfois même s'avérer plus pédagogique.

Supposons que nous connaissons la taille (hauteur) d'une personne pour prédire si la personne est un homme ou une femme. Nous pouvons donc parler de probabilité d'être un homme ou une femme. Imaginons que la probabilité d'être un homme pour une hauteur donnée est 90%. Alors l'odds d'être un homme est:

  (57.225)

Dans notre exemple, l'odds sera de 0.90/0.10 soit 9. Maintenant, la probabilité d'être une femme sera donc de 0.10/0.90 soit 0.11. Cette asymétrie des valeurs est peu parlante parce que l'odds d'être un homme devrait être l'opposé de l'odds d'être une femme idéalement. Nous résolvons justement cette asymétrie à l'aide du logarithme naturel. Ainsi, nous avons:

 et    (57.226)

Ainsi, le logit (logarithme de l'odds) est exactement l'opposé de celui d'être une femme de par la propriété du logarithme:

  (57.227)

Exemple:

Imaginons qu'une banque souhaite faire un scoring de ses débiteurs. Comme elle a plusieurs succursales (la banque) elle construit les tables de données (fictives...) suivantes pour chacune (toutes les succursales ne sont donc pas représentées):

- 1^ère succursale:

- 2^ème succursale:

- 3^ème succursale:

Nous pouvons voir que la proportion totale des bons débiteurs dans les trois succursales est de .

Quand le crédit est inférieur à 27'500, la proportion de bons débiteurs est de . Quand le montant des crédits est inférieur à 28'000 la proportion de bons débiteurs est de .

Quand le montant des crédits est inférieur à 28'500, la proportion de bons débiteurs est de , et pour des montants inférieurs à 30'000 la proportion est de .

Nous poserons pour cette régression logistique que  est un bon risque de crédit et  est un mauvais risque. Ensuite, nous créons le tableau suivant qui est un récapitulatif des données de toutes les succursales:

Ce qui donne graphiquement en Kilo-francs:

Figure: 57.8 - Pourcentage cumulé des bons débiteurs en fonction du crédit

Une fois ceci fait, nous utilisons la transformation en logit:

  (57.228)

Ce qui donne:

Une régression linéaire par la méthode des moindres carrés donne:

Figure: 57.9 - Logit des bons débiteurs en fonction du crédit

avec pour équation:

  (57.229)

La fonction logistique avec sa représentation vient alors immédiatement (les unités de x sont en millier de francs!):

  (57.230)

Ainsi, il est possible de dire dans cet exemple, qu'elle est la proportion P de bons ou mauvais payeurs en fonction d'une valeur de crédit X plus petite ou égale à une certaine valeur donnée. Puisque 0 est un mauvais risque de crédit, nous voyons que plus les crédits sont élevés, moins le risque est gros (dans ce cas fictif...). Par ailleurs, avec un logiciel comme Minitab, la différence entre les calculs (grossiers) effectués ici à la main et ceux effectués avec l'outil de régression logistique binaire du logiciel sont de l'ordre de 10% (car évidemment... Minitab utilise le concept des estimateurs de maximum de vraisemblance vus dans le chapitre de Statistiques pour déterminer les coefficients et la constante)

COEFFICIENT DE CORRÉLATION (DÉTERMINATION) GÉNÉRALISÉ

Nous avons vu plus haut ainsi que dans le chapitre de Statistiques plusieurs formulations du coefficient de corrélation linéaire, qui mis au carré était alors nommé coefficient de détermination de Pearson.

Nous avons également vu dans le chapitre de Statistiques le coefficient de corrélation non-paramétrique de Spearman mais celui-ci aussi est uniquement valable dans le cadre d'une relation linéaire (affine).

Afin de généraliser cela, une définition empirique plus générale a été donnée (je n'en ai jamais vu la démonstration ni le nom de la personne qui est l'origine de cette définition...):

  (57.231)

Ainsi, pour n'importe quelle dépendance, on cherche les paramètres des équations souhaitées avec les techniques de recherche opérationnelle (voir plus bas dans ce chapitre) les valeurs des coefficients de la fonction qui minimisent le coefficient de corrélation.

Ainsi dans le cas particulier de la régression linéaire puisque nous avons:

  (57.232)

Nous retrouvons alors:

  (57.233)

INTERPOLATION POLYNOMIALE

Il existe de nombreuses techniques d'interpolation de polynômes plus ou moins complexes et élaborées. Nous nous proposons ici d'en présenter quelques-unes dans l'ordre croissant de difficulté et de puissance d'application.

COURBES DE BÉZIER (B-SPLINE)

L'ingénieur russe Pierre Bézier (Peugeot), aux débuts de la Conception Assistée par Ordinateur (C.A.O), dans les années 60, a donné un moyen de définir des courbes et des surfaces à partir de points. Ceci permet la manipulation directe, géométrique, des courbes sans avoir à donner d'équation à la machine!!

Le thème des Courbes de Bézier est une notion à multiples facettes, vraiment très riche, au croisement de nombreux domaines mathématiques très divers: Analyse, Cinématique, Géométrie Différentielle, Géométrie Affine, Géométrie Projective, Géométrie Fractale, Probabilités, ...

Les Courbes de Bézier sont par ailleurs devenues incontournables dans leurs applications concrètes dans l'industrie, l'infographie, ...

Voilà l'approche mathématique de cette technique:

D'abord, nous savons que l'équation d'une droite que nous noterons dans le domaine M (par respect de tradition) joignant deux points  est:

  (57.234)

Ce qui est juste puisque lorsque  nous sommes en A et lorsque  nous sommes en B. Donc  et le point M parcoure tout le segment [AB]. Par construction, si A et B étaient des masses physiques égales à l'unité,  représente le barycentre (centre de gravité) du système pour un t donné.

Par définition, le segment [AB] est la "courbe de Bézier de degré 1" avec points de contrôle A et B et les Polynômes 1-t et t sont les "polynômes de Bernstein de degré 1".

Construisons maintenant une courbe paramétrée en rajoutant une 2ème étape à ce qui précède:

Figure: 57.10 - Courbe de Bézier avec un point supplémentaire

1ère étape:

- Soit  le barycentre de (A,1-t) et (B, t) et où  décrit [AB].
- Soit  le barycentre de (B,1-t) (C, t) et où  décrit [BC].

2ème étape:

- Soit M(t) le barycentre de .

Par construction, M(t) se situe donc à la même proportion du segment  que  par rapport au segment [AB] ou par rapport au segment [BC].

La courbe obtenue est alors l'enveloppe des segments : en tout point M, la tangente à la courbe est donc le segment .

M(t) décrit alors une Courbe de Bézier de degré 2, qui, par construction commence en A et se finit en C, et a pour tangentes [AB] en A et [BC] en C.

C'est en fait un arc de parabole (que nous pourrions noter très logiquement [ABC] ):

Figure: 57.11 - Arc de parabole [ABC]

Par le même schéma, nous pouvons définir une courbe de Bézier de n points  soit . C'est ce que nous appelons "l'algorithme de Casteljau". Ainsi, soit:

  (57.235)

Nous avons:

  (57.236)

La récurrence se terminant pour:

  (57.237)

Ainsi, pour  nous avons trivialement:

  (57.238)

Soit:

  (57.239)

Ainsi, nous avons forcément avec deux points l'équation d'une droite.

Considérons maintenant M(t) la courbe de Bézier d'ordre 3, nous avons donc les points définis par:

  (57.240)

Nous avons par la relation de récurrence:

  (57.241)

où nous avons éliminé les termes contenant des points non déterminés.

Nous avons donc:

  (57.242)

Il vient alors:

  (57.243)

Et donc:

  (57.244)

soit sous forme vectiorelle (plus conforme à la notatin mathématique d'usage):

  (57.245)

et sous forme matricielle:

  (57.246)

Par le même raisonnement, nous avons pour une courbe de Bézier d'ordre 4:

  (57.247)

soit sous forme vectorielle:

  (57.248)

Ce qui correspond de manière générique à:

Figure: 57.12 - Courbe de Bézier d'ordre 4

Et là nous pouvons creuser un peu les coefficients en comparant les coefficients de Bézier des courbes d'ordre 2, 3 et 4.

D'abord, reprenons la courbe de Bézier précédente:

  (57.249)

Nous remarquons d'abord aisément la proportionnalité suivante:

  (57.250)

et si nous regardons de plus près les coefficients, nous remarquons que nous avons aussi:

  (57.251)

Il ne s'agit ni plus ni moins que du triangle de Pascal!! Et donc les constantes sont simplement les coefficients binomiaux (cf. chapitre de Calcul Algébrique) donnés par pour l'ordre n dans notre cas par:

  (57.252)

Ainsi, les polynômes de Bernstein sont donnés par:

  (57.253)

et finalement les courbes de Bernstein par d'ordre n:

  (57.254)

Au fait, si nous avions noté plus haut la somme sous la forme suivante:

  (57.255)

Nous aurions alors les polynômes de Bernstein qui seraient donnés (ce qui est plus respectueux des traditions...):

  (57.256)

C'est une relation très pratique car elle permet de calculer facilement et rapidement le polynôme correspondant à une courbe de Bézier d'ordre n.

Nous avons alors:

  (57.257)

Un exemple très connu des courbes de Bézier d'ordre 3 est l'outil Plume des produits phares Adobe Photoshop ou Adobe Illustrator. Effectivement, ces deux outils créent une succession de courbes de Bézier d'ordre 3 dont le point  est défini après coup à la souris en utilisant des poignées appelées "torseurs" dans le langage de spécialistes Adobe... Voici un exemple pris d'un de ces logiciels fait avec un tracé à la plume de 5 points (soit 4 splines):

Figure: 57.13 - Exemples de splines avec un logiciel de dessin

Tant que l'utilisateur ne bouge pas les poignées de points alors tous les points sont alignés sur la droite. Nous avons alors l'impression d'avoir une spline d'ordre 2.

Exemple:

Un cercle, dessiné par un logiciel graphique est en pratique composé de 4 arcs de Bézier. Pour

observer cette particularité, il suffit de dessiner un cercle avec Illustrator par exemple, puis de le sélectionner pour voir apparaître les points de contrôle des arcs de Bézier qui le forment.

Nous allons nous intéresser à la meilleure façon de choisir les points de contrôle de ces arcs de sorte qu'ils ressemblent à des quarts de cercle, puis nous observerons la différence entre le dessin produit et un vrai cercle:

Figure: 57.14 - Exemple de construction d'un cercle avec des courbes de Bézier

Prenons le quart de cercle de rayon 1 centré à l'origine:

Il est approché par un arc de Bézier dont les points de contrôle sont . Les extrémités de l'arc de Bézier étant  et , il est naturel de choisir  et .

L'intuition nous amène à choisir   et  et il reste à trouver une valeur positive de k de sorte que l'arc de Bézier ressemble à un arc de cercle.

Nous obtenons ainsi l'équation paramétrique de l'arc de Bézier

  (57.258)

Soit:

  (57.259)

Nous pouvons par exemple chercher la valeur de k pour laquelle l'arc passe par le point:

  (57.260)

en . Il devient alors très simple à partir de l'équation paramétrique de déterminer k. Il s'agit simplement pour x (ou pour y) d'une simple équation à une inconnue.

MÉTHODE D'EULER

Il s'agit là de la méthode numérique la plus simple (elle est triviale et dans l'idée assez proche de la méthode de Newton même si leur objectif n'est pas le même). En fait elle ne fournit une approximation (au sens très large du terme) d'une fonction f(x) dont la dérivée est connue.

Ici les points choisis sont équidistants, c'est-à-dire (h étant le "pas"). Nous notons la valeur exacte et , la valeur approchée.

Il y a plusieurs méthodes pour procéder (comme souvent):

1. Graphiquement:

Nous nous déplaçons d'un pas h en x en suivant le vecteur de pente f(x,y)

Par construction, nous savons (cf. chapitre de Calcul Différentiel Et Intégral) que (nous adoptons une notation un peu particulière dans ce contexte):

  (57.261)

qui correspond donc bien à la pente (non instantané bien sûr!) en de la "courbe intégrale" passant par ce point. D'où:

  (57.262)

2. Analytiquement:

Nous remplaçons dans la dernière relation par . Nous obtenons alors:

  (57.263)

appelée "équation aux différences pour la méthode d'Euler".

L'application en est triviale et ne nécessite pas d'exemples particuliers.

POLYNÔME DE COLLOCATION

Soit une fonction connue sous forme explicite ou sous forme tabulée, et supposons qu'un certain nombre de valeurs:

  (57.264)

en sont données. Les points sont appelés les "points d'appui".

"Interpoler" f signifie estimer les valeurs de f pour des abscisses situées entre et , c'est-à-dire dans l'intervalle d'interpolation, par une fonction approximante , qui vérifie les "conditions de collocations" (rien à voir avec votre colocataire !):

Figure: 57.15 - Illustration du concept d'interpolation

La fonction p s'appelle "fonction de collocation" sur les . Lorsque p est un polynôme, nous parlons de "polynôme de collocation" ou de "polynôme d'interpolation".

"Extrapoler" f signifie approcher f(x) par p(x) pour des abscisses situées "hors" de l'intervalle d'interpolation.

Quand nous connaissons un polynôme de degré n en n+1 points, nous pouvons donc connaître par une méthode simple (mais pas très rapide - mais il existe plusieurs méthodes) complètement ce polynôme.

Pour déterminer le polynôme, nous allons utiliser les résultats exposés précédemment lors de notre étude des systèmes d'équations linéaires. Le désavantage de la méthode présentée ici est qu'il faut deviner à quel type de polynôme nous avons affaire et savoir quels sont les bons points qu'il faut choisir...

Un exemple particulier devrait suffire à la compréhension de cette méthode, la généralisation en étant assez simple (voir plus loin).

Soit un polynôme du second degré:

  (57.265)

et nous avons connaissance des points suivants (dont vous remarquerez l'ingéniosité des points choisis par les auteurs de ces lignes...):

  (57.266)

Nous en déduisons donc le système d'équations:

  (57.267)

Système qui une fois résolu dans les règles de l'art nous donne:

  (57.268)

Voyons le cas général:

Théorème: Soient des points d'appui, avec si . Alors il existe un polynôme de degré inférieur ou égal à n, et un seul, tel que pour .

Démonstration:

Posons:

  (57.269)

Les conditions de collocation:

  (57.270)

s'écrivent donc:

  (57.271)

Il s'agit d'un système de n+1 équations à n+1 inconnues.

Son déterminant s'écrit:

  (57.272)

relation que nous appelons "déterminant de Vandermonde". Nous savons que si le système à une solution, le déterminant du système doit être non nul (cf. chapitre d'Algèbre linéaire).

Montrons par l'exemple (en reprenant un polynôme du même degré que celui que nous avons utilisé plus haut) que le déterminant se calcule selon la relation suivante précédente (le lecteur généralisera par récurrence):

Donc dans le cas , nous considérons le déterminant:

  (57.273)

qui correspond dont au système (pour rappel):

  (57.274)

Calculons ce déterminant suivant la colonne 1 (en faisant usage des cofacteurs comme démontré dans le chapitre d'Algèbre linéaire):

  (57.275)

Ce dernier polynôme peut s'écrire:

  (57.276)

Ce qui s'écrit:

  (57.277)

Comme les sont dans l'énoncé de notre problème tous différents tels que alors le système a une solution unique! Ce qui prouve qu'il existe toujours alors un polynôme d'interpolation.

C.Q.F.D.

RECHERCHE DES RACINES

Bien des équations rencontrées en pratique ou en théorie ne peuvent pas être résolues exactement par des méthodes formelles ou analytiques. En conséquence, seule une solution numérique approchée peut être obtenue en un nombre fini d'opérations.

Évariste Galois a démontré, en particulier, que l'équation  ne possède pas de solution algébrique (sauf accident...) si  est un polynôme de degré supérieur à 4.

Il existe un grand nombre d'algorithmes permettant de calculer les racines de l'équation  avec une précision théorique arbitraire. Nous n'en verrons que les principaux. Attention, la mise en oeuvre de tels algorithmes nécessite toujours une connaissance approximative de la valeur cherchée et celle du comportement de la fonction près de la racine. Un tableau des valeurs de la fonction et sa représentation graphique permettent souvent d'acquérir ces connaissances préliminaires.

Si l'équation à résoudre est mise sous la forme , nous traçons les courbes représentant g et h. Les racines de l'équation  étant données par les abscisses des points d'intersection des deux courbes.

MÉTHODE DES PARTIES PROPORTIONNELLES

La mise en oeuvre, sur calculatrice, de cette méthode est particulièrement simple. Les conditions à vérifier étant seulement:

- f est continue

- f est monotone dans un voisinage de la racine

Dans un petit intervalle, nous pouvons remplacer une courbe par un segment de droite. Il y a plusieurs situations possibles mais en voici une particulière généralisable facilement à n'importe quoi:

Figure: 57.16 - Approximation d'une courbe par un segment de droite

Sur cette figure, nous tirons à l'aide des théorèmes de Thalès (cf. chapitre de Géométrie Euclidienne):

  (57.278)

d'où:

  (57.279)

Si , nous pouvons négliger f(a) au dénominateur et il vient:

  (57.280)

L'algorithme consiste donc à réaliser les étapes suivantes:

1. Choisir a et b, calculer f(a) et f(b) 

2. Déterminer . Si  est assez petit, nous arrêtons le calcul et affichons x et f(x)

3. Sinon nous procédons comme suit:

- Nous remplaçons b par a et f(b) par f(a) 

- nous remplaçons a par x et f(a)  par f(x)

- nous retournons au point (2)

MÉTHODE DE LA BISSECTION

La condition préalable à satisfaire pour cette méthode est de trouver un intervalle  tel que:

1.  f(x) est continue sur [a,b]

2.

Il faut encore fixer  qui est défini comme la borne supérieure de l'erreur admissible.

La méthode consiste à appliquer successivement les 4 étapes suivantes:

1. Calcul de

2. Evaluation de  f(x)

3. Si  alors le travail est terminé, il faut afficher x et  f(x)

4. Sinon on procède comme suit:

- on remplace a par x si

- on remplace b par x si  ou

- on retourne en (1)

L'étape (3) impose la condition  pour l'arrêt des calculs. Il est parfois préférable de choisir un autre critère de fin de calcul. Celui-ci impose à la solution calculée d'être confinée dans un intervalle de longueur  contenant . Ce test s'énonce comme suit:

3'. Si , le travail est terminé et  est affiché. Il est bien sûr évident que

MÉTHODE DE LA SÉCANTE (REGULA FALSI)

Figure: 57.17 - Illustration de la méthode de la sécante

Les conditions préalables sont les suivantes:

Il faut déterminer un intervalle [a,b] tel que:

1.  f(x) est continue sur [a,b]

2.

Si  est le point de coordonnées , alors les points  sont alignés sur la sécante. La proportion suivante (Thalès) est donc vraie:

  (57.281)

nous en déduisons:

  (57.282)

La méthode consiste à appliquer successivement les étapes suivantes:

1. Calcul de

2. Évaluation de

3. Si , le travail est terminé. Il faut afficher

4. Sinon nous procédons comme suit:

- nous remplaçons a par  si

- nous remplaçons b par  si  ou

- nous retournons en (1)

La condition (3) peut être remplacée par la condition:

3'. Si , alors le travail est terminé et nous affichons

MÉTHODE DE NEWTON

Considérons la figure suivante:

Figure: 57.18 - Illustration de la méthode de Newton

Si  est une approximation de la racine , nous remarquons que  en est une meilleure.  est l'intersection de la tangente à la courbe en  et de l'axe des abscisses.  est encore une meilleure approximation de ,  est obtenu de la même manière que  mais à partir de .

La méthode de Newton consiste en la formalisation de cette constatation géométrique.

Pour utiliser cette technique, rappelons que si nous prenons une fonction f qui est dérivable en , alors nous pouvons la réécrire sous la forme:

  (57.283)

où  est la dérivée de f en  et  est une fonction qui tend vers 0 comme  pour  lorsque x tend vers  (c'est un terme correctif qui sous-tend la suite des termes du développement de Taylor).

En appliquant ce résultat à la résolution de , nous obtenons:

  (57.284)

La fonction  empêche la résolution de cette équation par rapport à . En négligeant le terme , l'équation se réécrit:

  (57.285)

et se résout aisément par rapport à :

  (57.286)

Mais  ne satisfait pas, en général, l'égalité . Mais comme nous l'avons déjà souligné,  est plus petit que  si la fonction f satisfait à certaines conditions.

La méthode de Newton consiste à remplacer l'équation:

  (57.287)

par:

  (57.288)

et à résoudre itérativement cette équation.

Les conditions suivantes sont suffisantes pour assurer la convergence de la méthode:

Dans un intervalle [a,b] comprenant  et  il faut que:

1. La fonction soit deux fois dérivable

2. La dérivée f ' ne s'annule pas (monotonie)

3. La deuxième dérivée soit continue et ne s'annule pas (pas de point d'inflexion)

La condition (2) est évidente, en effet si  alors l'itération peut conduire à une erreur de calcul (singularité).

La condition (3) est moins évidente, mais le graphique suivant illustre un cas de non-convergence. Dans ce cas, le processus à une boucle calculant alternativement  et .

Figure: 57.19 - Exemple de cas de non-covergence avec la méthode de Newton

Si la fonction f est donnée analytiquement, sa dérivée peut-être déterminée analytiquement. Mais dans bien des cas, il est utile, voire indispensable de remplacer  par le quotient différentiel:

  (57.289)

où h doit être choisi suffisamment petit pour que la différence:

  (57.290)

soit elle aussi suffisamment petite.

L'itération s'écrit alors:

  (57.291)

Convergence de la méthode:

Si la méthode de résolution est convergente, l'écart entre  et  diminue à chaque itération. Ceci est assuré, par exemple, si l'intervalle [a,b] contenant , voit sa longueur diminuer à chaque étape. La méthode de Newton est intéressante car la convergence est quadratique:

  (57.292)

alors que la convergence des autres méthodes est linéaire:

  (57.293)

Considérons, par exemple, la méthode de la bissection vue précédemment. À chaque itération la longueur de l'intervalle [a,b] diminue de moitié. Ceci nous assure que l'écart  est réduit de moitié à chaque étape du calcul:

  (57.294)

Pour démontrer la convergence quadratique de la méthode de Newton, il faut utiliser les développements limités de f et f ' au voisinage de :

  (57.295)

Mais:

  (57.296)

donc:

  (57.297)

En soustrayant  à gauche et à droite de l'égalité et en mettant les deux termes du second membre au même dénominateur, il vient:

  (57.298)

et dès que  est assez petit, le dénominateur peut être simplifié.

  (57.299)

ce qui montre bien que la convergence est quadratique.

AIRES ET SOMMES DE RIEMANN

Considérons la figure suivante:

Figure: 57.20 - Illustration d'un intervalle sous une courbe

Nous désirons calculer l'aire comprise entre l'axe x, la courbe de f et les droites d'équations  et . Nous supposons dans ce cas que la fonction f est à valeurs positives:

  (57.300)

Ce problème, dans sa généralité, est difficile voire impossible à résoudre analytiquement. Voici donc quelques méthodes numériques permettant le calcul approché de cette aire (ces méthodes sont utilisées parfois dans les entreprises par les employés qui n'ont que des tableurs du type MS Excel ou OpenOffice Calc pour calculer des intégrales).

MÉTHODE DES RECTANGLES

Nous subdivisons l'intervalle  en n sous-intervalles dont les bornes sont . Les longueurs de ces sous-intervalles sont . Nous construisons les rectangles dont les côtés sont  et .

Figure: 57.21 - Approche de l'intégrale de Riemann par des rectangles inférieurs

L'aire de ces rectangles vaut:

  (57.301)

Si les  sont suffisamment petits,  est une bonne approximation de l'aire cherchée. Nous pouvons recommencer cet exercice en choisissant  et  comme côtés des rectangles. Nous obtenons alors:

  (57.302)

La figure correspondante est la suivante:

Figure: 57.22 - Approche de l'intégrale de Riemann par des rectangles supérieurs

Encore une fois, l'aire de ces rectangles approche l'aire cherchée. Afin de simplifier la programmation, il est utile de choisir des intervalles de longueur identique:

  (57.303)

Si nous avons n rectangles, h vaut alors . Les aires  et  deviennent:

  (57.304)

MÉTHODE DES TRAPÈZES

Afin d'augmenter la précision des calculs, il est possible de calculer:

  (57.305)

Dans le cas où tous les intervalles sont de longueur égale,  vaut:

  (57.306)

Il existe une foule d'autres méthodes permettant la résolution de ce problème (dont la méthode de Monte-Carlo que nous verrons plus loin).

Dans le cas où la fonction f n'est pas à valeurs positives, nous ne parlons plus d'aire mais de "somme de Riemann". Les sommes à calculer sont alors:

  (57.307)

et:

  (57.308)

Tous les calculs doivent être conduits de la même manière, mais les résultats peuvent être positifs, négatifs ou nuls.

PROGRAMMATION LINÉAIRE

L'objectif de la programmation linéaire (P.L.) est de trouver la valeur optimale d'une fonction linéaire sous un système d'équations d'inégalités de contraintes linéaires. La fonction à optimiser est baptisée "fonction économique" (utilisée en économie dans le cadre d'optimisations) et on la résout en utilisant une méthode dite "méthode simplexe" (voir plus loin) dont la représentation graphique consiste en un "polygone des contraintes". 

La recherche opérationnelle a pour domaine l'étude de l'optimisation de processus quels qu'ils soient. Il existe de nombreux algorithmes s'inspirant des problèmes du type exposés lors de notre étude de la programmation linéaire. Nous nous attarderons en particulier sur l'algorithme le plus utilisé qui est "l'algorithme du simplexe".

Lorsque l'on peut modéliser un problème sous forme d'une fonction économique à maximiser dans le respect de certaines contraintes, alors on est typiquement dans le cadre de la programmation linéaire.

Soit une fonction économique Z telle que:

  (57.309)

où les  sont des variables qui influent sur la valeur de Z, et les  les poids respectifs de ces variables modélisant l'importance relative de chacune de ces variables sur la valeur de la fonction économique.

Les contraintes relatives aux variables s'expriment par le système linéaire suivant:

  (57.310)

Sous forme générale et matricielle ce genre de problème s'écrit:

  (57.311)

Voyons un exemple qui consiste à résoudre le problème simple suivant:

Une usine fabrique 2 pièces P1 et P2 usinées dans deux ateliers A1 et A2. Les temps d'usinage sont pour P1 de 3 heures dans l'atelier A1 et de 6 heures dans l'atelier A2 et pour P2 de 4 heures dans l'atelier A1 et de 3 heures dans l'atelier A2.

Le temps de disponibilité hebdomadaire de l'atelier A1 est de 160 heures et celui de l'atelier A2 de 180 heures.

La marge bénéficiaire est de 1'200.- pour une pièce P1 et 1'000.- pour une pièce P2.

Quelle production de chaque type doit-on fabriquer pour maximiser la marge hebdomadaire?

Le problème peut se formaliser de la façon suivante:

  (57.312)

La fonction économique étant:

  (57.313)

à maximiser.

Résolution graphique du problème (ou méthode du "polygone des contraintes"): les contraintes économiques et de signe sont représentées graphiquement par des demi-plans. Les solutions, si elles existent, appartiennent donc à cet ensemble appelé "région des solutions admissibles":

Figure: 57.23 - Illustration d'un problème de recherche opérationnelle simple avec région des admissibles

Pour trouver les coordonnées des sommets, on peut utiliser le graphique si les points sont faciles à déterminer.

Il s'agit donc de chercher à l'intérieur de ce domaine (connexe), le couple maximisant la fonction économique.

Or, l'équation Z est représentée par une droite de pente constante (-1.2) dont tous les points fournissent la même valeur Z pour la fonction économique.

En particulier, la droite passe par l'origine et donne une valeur nulle à la fonction économique. Pour augmenter la valeur de Z et donc la fonction économique, il suffit d'éloigner de l'origine (dans le quart de plan) la droite de pente -1.2.

Pour respecter les contraintes, cette droite sera déplacée, jusqu'à l'extrême limite où il n'y aura plus qu'un point d'intersection (éventuellement un segment) avec la région des solutions admissibles.

Figure: 57.24 - Recherch de solutions graphiquement avec la droite de la fonction économique

La solution optimale se trouve donc nécessairement sur le pourtour de la région des solutions admissibles et les parallèles formées par la translation de la fonction économique s'appellent les "droites isoquantes" ou "droites isocoûts"...

Voyons maintenant comment résoudre ce problème de manière analytique avant de passer à la partie théorique.

Nous avons donc le "système canonique":

  (57.314)

avec:

  (57.315)

Nous introduisons d'abord la "variable d'écart" afin de transformer les 3 inégalités par des égalités. Le système d'équations devient alors une "forme standard":

  (57.316)

La situation peut se résumer dans le tableau suivant (nous omettons la représentation de la variable d'écart dans le tableau-matrice qui ne sert qu'à égaliser les équations):

Nous déterminons maintenant le pivot (voir plus loin la méthode du pivot), pour cela nous choisissons la colonne où le coefficient économique est le plus grand. Ici c'est la colonne 1.

Ensuite, nous effectuons les procédures suivantes:

1. Le pivot est remplacé par son inverse

2. On divise les éléments de la ligne du pivot (pivot exclu) par le pivot

3. On divise les éléments de la colonne du pivot (pivot exclu) par le pivot mais on change leur signe ensuite

4. Pour les autres éléments de la première ligne: élément de la ligne 1 diminué de l'élément correspondant sur la ligne de pivot multiplié par 3/6 (rapport des valeurs dans la colonne de pivot)

Nous obtenons dès lors:

Ce qui donne:

Nous n'atteignons la solution optimale que lorsque tous les éléments de la marge sont négatifs ou nuls. Il faut donc continuer (car il reste 400 dans la colonne ) ... ici, on atteint déjà l'optimum au troisième tableau, mais ce n'est pas une généralité (le pivot est 2.5 cette fois). On recommence donc les opérations:

Ce qui donne:

Le processus est terminé car tous les termes de la fonction économique sont négatifs. Le programme optimum est donc de et pour un résultat de:

  (57.317)

ALGORITHME DU SIMPLEXE

Pour mettre en oeuvre cet algorithme, nous devons poser le problème sous une forme "standard" et introduire la notion de "programme de base" qui est l'expression algébrique correspondant à la notion de "point extrême du polyèdre des programmes admissibles" étudiée lors de la programmation linéaire (noté ci-après P.L.). En effet, nous verrons que la solution d'un problème du type P.L. si elle existe, peut toujours être obtenue en un programme de base. La méthode du simplexe va donc consister à trouver un premier programme de base puis à construire une suite de programmes de base améliorant constamment la fonction économique et donc conduisant à l'optimum.

Un problème de P.L. est donc mis sous sa "forme standard" s'il implique la recherche du minimum de la fonction objectif sous des contraintes ayant la forme d'équations linéaires et de conditions de non négativité des variables, c'est-à-dire s'il se pose sous la forme que nous avons vue lors de notre étude de la programmation linéaire:

  (57.318)

C'est-à-dire aussi, en utilisant des notations matricielles:

  (57.319)

où les matrices  correspondent, respectivement, aux coefficients des niveaux d'activité dans la fonction objectif, aux coefficients techniques des activités et aux seconds membres des contraintes.

Nous allons voir maintenant comme un problème général de P.L. peut toujours être ramené à une forme standard. La notion de "variable d'écart" est essentielle pour effectuer cette "réduction".

Chercher le maximum d'une fonction f(x) revient à chercher le minimum de la fonction de signe opposé -f(x) . D'autre part, une contrainte qui se présente comme une inéquation:

  (57.320)

peut être remplacée par l'équation:

  (57.321)

impliquant une variable supplémentaire, , appelée donc "variable d'écart", et soumise à la contrainte de non-négativité, .

Bien évidemment, dans un cas contraire tel où le système est du type:

  (57.322)

Nous écrirons:

  (57.323)

impliquant donc également une variable supplémentaire et soumise à la contrainte de non-négativité, .

Ce travail de mise en forme standard nous permet donc de retrouver un système d'équations linéaires à résoudre (nous avons vu précédemment sur le site comme résoudre ce genre de système avec l'algorithme du pivot).

La matrice A qui représente les composantes du système d'équations peut s'exprimer dans différentes variantes en fonction de la base vectorielle choisie (voir le chapitre d'Analyse vectorielle dans la section d'algèbre). Nous allons introduire la notion de "forme canonique utilisable" associée au choix d'une base et nous montrerons que cette reformulation du système de contraintes va nous permettre de progresser vers l'optimum.

La matrice A peut, après introduction des variables d'écart se décompenser en deux sous-matrices , une contenant les variables initiales D et l'autre comportant les variables d'écart B tel que:

  (57.324)

Les colonnes de la matrice B sont bien évidemment, par définition de la méthode, des colonnes unités, linéairement indépendantes. Ces colonnes forment une base de l'espace vectoriel des colonnes à m éléments (ou dimensions) - le nombre de lignes du système. Nous appelons B la "matrice de la base".

Les variables associées aux composantes colonnes de la matrice B seront dès maintenant appelées les "variables de bases". Dans ce cas, les variables de base sont donc essentiellement les variables d'écart . Les variables associées aux colonnes de la matrice D seront appelées les "variables hors-base"; il s'agit des variables .

En résumé, tout P.L. une fois mis sous forme standard est tel que:

- il existe une sous-matrice carrée de la matrice A des coefficients techniques, qui est appelée matrice de base et qui est égale à la matrice carrée unité I de dimension  (effectivement il y autant de variables d'écart que de lignes dans le système d'équations original - au nombre de m - et autant de colonnes puisque chaque variable d'écart à un indice différent).

- les variables de base associées n'apparaissent pas dans l'expression de la fonction économique.

- le second membre des contraintes est constitué d'éléments non négatifs.

Nous disons que le problème est mis sous "forme canonique utilisable associée à la base B, correspondant aux variables de base ".

Il est maintenant commode d'introduire les notations suivantes:

  (57.327)

qui sont respectivement les vecteurs des variables de base et le vecteur des variables hors-base.

Ainsi, le système d'équations décrivant les contraintes peut s'écrire indifféremment:

  (57.328)

ou bien aussi:

  (57.329)

Si la matrice B est une matrice de base, elle est non singulière et admet donc une matrice inverse . En multipliant cette équation, à gauche et à droite, par  nous obtenons:

  (57.330)

Le système d'équations aura alors été mis sous une forme résolue en .

Pour obtenir une forme canonique utilisable associée à la base B, correspondant aux variables de base, il ne reste plus qu'à éliminer les variables de base de l'expression de la fonction économique.

Écrivons cette fonction en séparant les vecteurs  et , nous obtenons:

  (57.331)

Nous pouvons alors facilement exprimer  en fonction de . En utilisant le système d'équations mis sous forme résolue en , nous avons dans un premier temps:

  (57.332)

que nous substituons dans la fonction économique, pour obtenir:

  (57.333)

Nous regroupons les termes en  et nous avons:

  (57.334)

Nous avons alors toujours système d'équations mais ne comportant plus d'inégalités mais au contraire des égalités !!! Reste plus qu'à démontrer que la solution de ce système dit "programme base" par la méthode du pivot est optimale.

Définition: nous appelons coût réduit de l'activité hors base j, le coefficient correspondant  de la ligne .

Soit un problème de programmation linéaire sous forme standard:

  (57.335)

La matrice A à m lignes (autant qu'il y a de contraintes) et n colonnes, avec . Si nous sélectionnons m variables de base et si nous annulons les  variables hors base, la matrice A:

  (57.336)

et le système se réduit à:

  (57.337)

La matrice B est de dimension . Si elle définit une base, elle admet une matrice inverse . Une solution du système est donc:

  (57.338)

Si l'expression  est non négative , nous avons une solution admissible qui vérifie les contraintes et que nous appellerons un programme de base:

  (57.339)

Le problème de programmation linéaire, s'écrit aussi sous la forme suivante, que nous appelons "forme canonique utilisable associée au programme de base":

  (57.340)

Définition: Nous appelons "coût réduit" de l'activité hors base j, le coefficient correspondant  de chaque ligne j de l'expression

À partir des développements effectués précédemment nous pouvons énoncer le résultat suivant:

Proposition 1: si dans la forme canonique utilisable associée à un programme de base, tous les coûts réduits sont  alors le programme de base est optimal.

Proposition 2: La solution d'un problème de P.L., si elle existe, peut toujours être trouvée en un programme de base.

La démonstration précise de ce résultat est assez délicate. Nous pouvons cependant en avoir une intuition en considérant, une fois de plus, la notion de coût réduit.

En effet, pour un programme de base donné, considérons la forme canonique utilisable associée à la base. Sur cette forme nous pouvons vérifier que, ou bien le programme de base est optimal (tous les coûts réduits sont ), ou bien que le programme de base peut être amélioré et remplacé par un nouveau programme de base donnant à z une valeur plus petite (un coût réduit est négatif et la variable hors-base associée peut être augmentée jusqu'à ce qu'une ancienne variable de base s'annule). Comme il y a un nombre fini de programmes de base (au plus égal au nombre ), la solution de P.L. se trouve nécessairement en un programme de base.

MÉTHODES DE MONTE-CARLO

La méthode de Monte-Carlo est un moyen très efficace de contourner les problèmes mathématiques et physiques les plus complexes. Elle trouve ses applications dans des domaines variés dont voici quelques exemples:

- Problèmes de neutronique liés à la bombe atomique

- Calculs d'intégrales ou de paramètres divers de variables aléatoires (finance, risque)

- Résolution d'équations elliptiques ou paraboliques

- Résolution de systèmes linéaires

- Résolution de problèmes d'optimisation (recherche opérationnelle, gestion de projets)

Il existe donc deux types de problèmes qui peuvent être traités par la méthode de Monte-Carlo: les problèmes probabilistes, qui ont un comportement aléatoire, et les problèmes déterministes, qui n'en ont pas.

Pour ce qui est du cas probabiliste, il consiste à observer le comportement d'une série de nombres aléatoires qui simule le fonctionnement du problème réel et à en tirer les solutions.

Pour le cas déterministe, le système étudié est complètement défini et on peut en principe prévoir son évolution, mais certains paramètres du problème peuvent être traités comme s'il s'agissait de variables aléatoires. Le problème déterministe devient alors probabiliste et ré solvable de façon numérique. On parle alors d'estimation de "Monte-Carlo" ou d'une approche de "MC élaborée".

La dénomination de méthode de "Monte-Carlo" date des alentours de 1944. Des chercheurs isolés ont cependant utilisé bien avant des méthodes statistiques du même genre: par exemple, Hall pour la détermination expérimentale de la vitesse de la lumière (1873), ou Kelvin dans une discussion de l'équation de Boltzmann (1901), mais la véritable utilisation des méthodes de Monte-Carlo commença avec les recherches sur la bombe atomique.

Au cours de l'immédiat après-guerre, Von Neumann, Fermi et Ulam avertirent le public scientifique des possibilités d'application de la méthode de Monte-Carlo (par exemple, pour l'approximation des valeurs propres de l'équation de Schrödinger). L'étude systématique en fut faite par Harris et Hermann Khan en 1948. Après une éclipse due à une utilisation trop intensive pendant les années 1950, la méthode de Monte-Carlo est revenue en faveur pour de nombreux problèmes: en sciences physiques, en sciences économiques, pour des prévisions électorales, etc., bref, partout où il est fructueux d'employer des procédés de simulation.

GÉNÉRATION DES VARIABLES ALÉATOIRES

Le mieux pour comprendre la méthode de Monte-Carlo c'est de faire des exemples. Mais pour cela, il faut d'abord d'avoir un très bon générateur de nombres aléatoires (ce qui est très difficile). C'est un domaine très délicat et sensible donc des normes internationales ont été édictées (ISO 28640:2010).

Prenons comme exemple le générateur de Maple:

>rand();

>restart;rand();

Nous voyons donc que la fonction par défaut de générateur de nombres aléatoires de Maple est à utiliser avec la plus grande prudence puisqu'une réinitialisation du système suffit à retrouver des valeurs aléatoires... égales. Il s'agit donc d'un "générateur pseudo-aléatoire" permettant de faire des simulations appelées parfois "pseudo Monte-Carlo".

Cependant il existe des libraires spécialisées dans Maple tel que:

>restart;readlib(randomize):randomize():rand();

>restart;readlib(randomize):randomize():rand();

Épreuve a priori réussie (au fait, il nous faudrait faire un beaucoup plus grand nombre d'essais afin de bien vérifier que le générateur ne suive pas une loi de distribution connue... ce qui n'est malheureusement jamais le cas).

Les fonctions ALEA( ) et ALEA.ENTRE.BORNES( ) de MS Excel sont aussi des générateurs pseudo-aléatoires dont voici un échantillon de 100 simulations (évidemment dans MS Excel le graphique ci-dessous changera à chaque fois que vous activerez la touche F9 du clavier):

Figure: 57.25 - Illustration d'une séquence de nombres pseudo-aléatoires avec MS Excel

Il peut malheureusement arriver avec les nombres pseudo-aléatoires que les nombres générés se présentent en grappes, c'est-à-dire en séries de nombres rapprochés les uns des autres, ce qui nuit à l'efficacité de la simulation de Monte-Carlo.

Une technique empirique consiste à faire appel à des séquences de nombres générés sur la base d'algorithmes qui balaient à coup sur l'intervalle [0,1]. Nous parlons alors de "nombres quasi-aléatoires" permettant de faire des simulation appelées parfois "quasi Monte-Carlo". Avec MS Excel, il est possible de créer une fonction qui remplacera les générateurs pseudo-aléatoires que sont les fonctions ALEA( ) ou ALEA.ENTRE.BORNES( ).

Voici donc un exemple de fonction en V.B.A. (Visual Basic for Application) qui génère des nombres quasi-aléatoires appelés "séquence de Fauré":

Function SequenceFaure(n) As Double
Dim f As Double, sb As Double
Dim i As Integer, n1 As Integer, n2 As Integer
n1 = n
sb = 1 / 2
Do While n1 > 0
n2 = Int(n1 / 2)
i = n1 - n2 * 2
f = f + sb * i
sb = sb / 2
n1 = n2
Loop
SequenceFaure = f
End Function

Ce qui donnera la séquence suivante pour un échantillon de 100 simulations:

Figure: 57.26 - Illustration d'une séquence de nombres quasi-aléatoires avec MS Excel

où nous voyons bien que la séquence couvre bien la surface comprise 0 et 1 (nous disons alors qu'elle couvre plus rapidement la surface d'intégration). Cette technique est parfois appréciée car elle a pour avantage de conserver les valeurs de la simulation à chaque fois que l'on relance la simulation (donc dans MS Excel le graphique ci-dessous ne changera pas quand vous activerez la touche F9 du clavier).

Par contre les générateurs de séquence ont une grande faiblesse: ils ne sont applicables (à ma connaissance du moins) que pour des problèmes de simulations avec une seule et unique variable aléatoire (typiquement du pricing d'options selon Black & Scholes). Effectivement si nous avons plusieurs variables aléatoires (et c'est le cas le plus courant!), alors les variables sont artificiellement corrélées (coefficient de corrélation à 1) car elles parcourent toutes la surface comprise entre 0 et 1 de la même manière. Donc une bonne simulation avec plusieurs variables est une simulation dont les variables traitées ont un coefficient de corrélation qui tend vers zéro.

De plus, les générateurs de séquence nécessitent des algorithmes qui sont très gourmands lorsqu'il y a de nombreuses variables par rapport à un générateur pseudo-aléatoire, raison pour laquelle dans la majorité des situations, on préférera cette bonne vieille méthode.

Une fois le générateur créé et testé, nous pouvons voir quelques résultats de la méthode de Monte-Carlo. Ainsi, dans le calcul des intégrales, celle-ci s'avère très utile et très rapide en termes de vitesse de convergence.

CALCUL D'UNE INTÉGRALE

Soit à calculer l'intégrale d'une fonction f définie et positive sur l'intervalle [a,b]:

  (57.341)

Soit:

  (57.342)

Nous considérons le rectangle englobant de la fonction sur [a,b] défini par les points  . Nous tirons un grand nombre N de points au hasard dans ce rectangle. Pour chaque point, nous testons s'il est au-dessous de la courbe. Soit F la proportion de points situés au-dessus, nous avons:

  (57.343)

L'algorithme Maple est donné par:

intmonte:=proc(f,a,b,N)
local i,al,bl,m,F,aleaabs,aleaord,estaudessous;
m:=round(max(a,b)*10^4);
al:=round(a*10^4);
bl:=round(b*10^4);
aleaabs:=rand(al..bl);
aleaord:=rand(0..m);
F:=0;
for i from 1 to N do
     estaudessous:=(f(aleaabs()/10^4)-aleaord()/10^4)>=0;
     if estaudessous then
          F:=F+1;
     fi
od:
RETURN((b-a)*max(a,b)*F/N)
end:

CALCUL DE PI

Pour le calcul de  le principe est le même et constitue donc à utiliser la proportion du nombre de points dans un quartier de cercle (cela permet de simplifier l'algorithme en se restreignant aux coordonnées strictement positives) inscrit dans un carré relativement au nombre de points totaux (pour tester si un point est à l'extérieur du cercle, nous utilisons bien évidemment le théorème de Pythagore) tel que:

  (57.344)

L'algorithme Maple est donné par:

estalinterieur:=proc(x,y) x^2+y^2<1 end:
calculepi:=proc(N)
local i,F,abs,ord,alea,erreur,result;
alea:=rand(-10^4..10^4);
F:=0;
for i from 1 to N do
     abs:=alea()/10^4;ord:=alea()/10^4;
       if estalinterieur(abs,ord) then
            F:=F+1;
       fi
od;
RETURN(4*F/N)
end:

MODÉLISATION

L'application la plus courante de la méthode par Monte-Carlo est certainement l'étude de variables aléatoires. Par ailleurs, cette méthode fait partie intégrante de la norme ISO 31010 de gestion du risque sous le nom "analyse de Monte-Carlo" tellement elle est courante et utile. De nombreuses entreprises font de la modélisation de Monte-Carlo avec un tableur comme MS Excel (même les multinationales!) et dans une moindre mesure avec des logiciels comme @Risk, CrystalBall, TreeAge ou encore Isograph.

Les avantages de méthode dans la modélisation de variables aléatoires sont les suivants:

- On peut intégrer n'importe quelle distribution dans une variable d'entrée, y compris des empiriques!

- Les modèles sont très simples à développer et peuvent être étendus à mesure des besoins

- Toutes les influences ou relations se produisant dans la réalité peuvent être représentées

- L'analyse de la sensibilité (cf. chapitre Techniques De Gestion) peut être appliquée

- Les modèles sont aisément compréhensibles et fournissent une mesure de l'exactitude d'un résultat

- De nombreux logiciels sont disponible et peu onéreux

Considérons un cas simple mais concret (très utilisé dans les entreprises) d'un petit projet de deux tâches notées A et B qui se succèdent sans marge libre. La durée de chacune des tâches a été estimée conformément à la recommandation du Project Management Institute avec une loi beta (cf. chapitre de Statistiques) comme l'apprennent tous les responsables de projets lors de leur cursus de formation (cf. chapitre Techniques De Gestion).

Pour cet exemple, la tâche A a une durée optimiste de 5 jours et pessimiste 8 jours. La tâche B une durée optimiste de 1 jour et pessimiste de 4 jours. Nous souhaiterions dans MS Excel à l'aide d'une simulation de pseudo Monte-Carlo (basé donc obligatoirement sur une variable pseudo-aléatoire) indiquer les trois informations traditionnelles minimales suivantes:

- Un tableau avec les 3 colonnes (durée de A, de B et somme des deux) de 10'000 simulations
- La fonction de distribution de la somme des deux variables aléatoires sous forme graphique
- La convergence du 95ème centile des 100 premières simulations (utile pour le sujet d'après).

Nous construisons alors le tableau suivant sur 10'000 lignes (la capture d'écran ne prend que les premières...):

Figure: 57.27 - Mise en place d'une petite simulaiton de Monte-Carlo avec MS Excel

où toutes les cellules de la colonne A contiennent la fonction suivante (MS Excel 2010):

=BETA.INVERSE.N(ALEA.ENTRE.BORNES(1;9999)/10000;3+RACINE(2);3-RACINE(2);5;8)

et toutes les cellules de la colonne B contiennent la fonction suivante (MS Excel 2010):

=BETA.INVERSE.N(ALEA.ENTRE.BORNES(1;9999)/10000;3+RACINE(2);3-RACINE(2);1;4)

et enfin la colonne C contient la fonction suivante:

=A2+B2

évidemment les valeurs évidemment dans MS Excel changeront à chaque fois que vous activerez la touche F9 du clavier.

Cela nous donne alors pour l'histogramme (dont je ne vais pas détailler la construction car il s'agit d'un sujet élémentaire en bureautique):

Figure: 57.28 - Distribution obtenue de la somme des variables aléatoires

et la convergence du 95ème centile sur les 100 premières simulations (car le problème étant simple, le système converge suffisamment rapidement pour ne pas avoir besoin d'en prendre plus de 100 pour l'exemple):

Figure: 57.29 - Illustration de la convergence du 95ème centile

Évidemment dans MS Excel les graphiques ci-dessus changeront à chaque fois que vous activerez la touche F9 du clavier.

Il existe d'autres techniques de réduction de la variance (ou in extenso: de l'écart-type) que la méthode de quasi Monte-Carlo permettant de réduire la variance des estimateurs de Monte-Carlo dans certaines conditions particulières:

1. Une de ces techniques est l'usage des "variables antithétiques" qui consiste très simplement (la programmation de cette technique est du niveau du lycée) à décoreller les simulations pour rendre la covariance entre les variables négatives et ainsi de réduire la variance (puisque comme nous l'avons vu dans le chapitre de Statistiques, la variance de la somme de deux variables aléatoires fait apparaître un terme de covariance). Malheureusement, cette technique fonctionne de manière satisfaisante qu'avec des distributions symétriques ce qui fait qu'à ma connaissance elle n'est pas implémentée dans les logiciels de simulation disponibles sur le marché.
   
   
2. Il existe aussi la technique "d'échantillonnage stratifié" qui consiste à découper l'espace des pré-images de la variables aléatoires en intervalles réguliers (la programmation de cette technique est aussi du niveau du lycée). Cette technique fonctionne très bien lorsque le nombre de simulations doit être faible mais seulement dans le cas d'une une unique variable. Raison pour laquelle est n'est pas implémentée à ma connaissance dans les logiciels de simulation disponibles sur le marché.
   
   
3. Il existe une généralisation de l'échantillonage stratifié (la programmation de cette technique est aussi du niveau du lycée) pour les simulations comportant plusieurs variables et qui se nomme "Latin Hypercube" (abrégé "LHS" pour Latin Hypercube Stratification). Cette technique assure donc que chaque n-uplet de variables aléatoires (correspondant à un espace à n -dimensions) utilise une pré-image unique à chaque itération, d'où le nom de la technique (Latin: faire référence aux carrés magiques où chaque valeur apparaît de manière unique, Hypercube: car il s'agit d'une généralisation à n dimensions d'un carré magique). Certains logiciels de simulation disponibles sur le marchés implémentent cette technique (@Risk, CrystalBall).

Pour résumer, que ce soit la technique des générateurs de séquence de Fauré, des variables antithétiques, des variables de contrôle, de l'échantillonage stratifié ou de Latin Hypercube même si ces techniques sont toutes faciles à programmer, la métode utilisant les variables pseudo aléatoires est privilégiée car est la plus adaptée à la majorité des situations courantes de l'économie mondiale.

BOOTSTRAPPING

En Statistiques, les techniques de bootstrap sont des méthodes d'Inférence statistique requérant des calculs informatiques relativement intensifs. L'objectif est de connaître certaines indications sur une statistique : son estimation bien sûr, mais aussi la dispersion (variance, écart-type), des intervalles de confiance voire un test d'hypothèse. Cette méthode est basée sur des simulations, comme les méthodes de Monte-Carlo, à la différence près que le bootstrap ne nécessite pas d'information supplémentaire que celle disponible dans l'échantillon. En général, il est basé sur de  nouveaux échantillons obtenus par tirage avec remise à partir de l'échantillon initial (on parle de rééchantillonnage).

Nous allons illustrer le principe du bootstrap (dit aussi "bootstrapping") sur l'exemple de l'intervalle de confiance de l'espérance  d'une variable aléatoire. Pour cet exemple, l'intervalle de confiance de l'espérance d'une variable aléatoire est parfaitement déterminé à partir de la moyenne et de la variance calculées sur l'échantillon (cf. chapitre de Statistiques).

Nous considérons un échantillon de la variable aléatoire composé de  estimations:

  (57.345)

La moyenne arithmétique de l'échantillon est:

  (57.346)

et son écart-type (estimateur de maximum de vraisemblance non biaisé):

  (57.347)

Comme nous sommes dans la situation d'une moyenne empirique connue et d'une variance empirique connue, pour faire le calcul d'un intervalle de confiance, nous avons alors démontré dans le chapitre de Statistiques qu'il fallait utiliser:

  (57.348)

où S est une autre notation traditionnelle dans certains domaines de la statistique pour la notation de l'écart-type empirique (cf. chapitre de Statistiques). Nous avons alors pour l'intervalle de confiance à 95% de l'espérance:

  (57.349)

Soit:

  (57.350)

Ce qui donne:

  (57.351)

L'intervalle de confiance peut être également calculé par bootstrap. Il est alors obtenu par l'algorithme suivant:

À partir de l'échantillon initial, nous simulons de nouveaux échantillons, appelés "répliques", de taille n, par tirages aléatoires avec remise. Par exemple avec la série précédente, nous pourrions obtenir la réplique suivante:

  (57.352)

dans laquelle certaines valeurs de l'échantillon initial ne figurent pas, et où d'autres apparaissent plusieurs fois. Plusieurs échantillons sont ainsi simulés. Nous pouvons ainsi former un nombre de répliques (arrangements) égal à (cf. chapitre de Probabilités):

  (57.353)

Pour chaque échantillon simulé, une moyenne est calculée (plusieurs milliers de moyennes!). L'intervalle de confiance à 95% est défini sur cet ensemble de moyennes typiquement à l'aide du calcul des centiles (via les fonctions d'un tableur ou d'un langage de programmation).

Évidemment pour chaque ensemble de plusieurs milliers de valeurs, les centiles ne seront pas les mêmes donc il est même possible de créer un intervalle de confiance pour les centiles eux-mêmes!

Il est très facile (au même titre que la méthode de Monte-Carlo) de créer des répliques avec des tableurs (de type MS Excel) sans faire de la programmation informatique! En plus la technique du bootstrap est très puissante car elle ne fait appel à aucune hypothèse sur la distribution statistique sous-jacente. Le domaine le plus courant et simple d'application du bootstrapping est la gestion de projets où lors de réunions avec une dizaine de ressources chacune estime la durée d'une tâche ou d'une phase.

Le bootstrap peut donc être appliqué à tout estimateur autre que la moyenne, tel que la médiane, le coefficient de corrélation entre deux variables aléatoires ou la valeur propre principale d'une matrice de variance-covariance (pour l'analyse en composantes principales) et c'est là sa grande force!!! Effectivement, pour ces estimateurs, il n'existe pas de relation mathématique qui définisse l'erreur-standard ou l'intervalle de confiance. Les seules méthodes applicables sont des "méthodes de rééchantillonage" (resampling) comme en fait partie le bootstrapping.

Exemple:

Avec par exemple le tableau MS Excel 2007 et en s'interdisant de faire de la programmation VBA, nous construisons un petit tableau avec l'échantillon:

Figure: 57.30 - Échantillon de base

Nous souhaiterions pouvoir déterminer un intervalle de confiance de la médiane (nous faisons exprès de prendre un indicateur statistique pour lequel il n'existe pas d'intervalle de confiance analytique). Pour cela, nous calculons la médiane de plusieurs milliers de réplications où chaque réplication correspond à une ligne:

Figure: 57.31 - Médianes de répliques

avec la formule MS Excel suivante qu'il faut mettre dans F5 et tirer ensuite jusqu'à la fin de la feuille:

=MEDIANE(INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1);
INDEX($A$5:$A$14;ALEA.ENTRE.BORNES(1;10);1))

Nous pourrons donc en avoir 10 milliards pas plus... (comme MS Excel est limité à un million de lignes, cela coupe net toute discussion...).

Il suffit ensuite dans une cellule de votre choix de mettre:

=CENTILE(F5:F2003;0.025)

et dans une autre:

=CENTILE(F5:F2003;0.975)

ce qui avec 2'000 réplications donnera respectivement 7 et 29.5.

Avec des connaissances élémentaires du tableur, il est possible de montrer graphiquement la convergence de la moyenne de la médiane en fonction du nombre de réplications (ci-dessous avec les 100 premières réplications):

Figure: 57.32 - Convergence de la médiance en fonction du nombre de réplications

Évidemment, ce graphique aura un aspect différent à chaque fois que vous relancerez la simulation dans MS Excel en appuyant sur la touche F9 du clavier

DICHOTOMIE

La dichotomie consiste pour un objet de taille N à exécuter un algorithme de façon à réduire la recherche à un objet de taille N/2. On répète l'algorithme de réduction sur ce dernier objet. Ce type d'algorithme est souvent implémenté de manière récursive. Lorsque cette technique est utilisable, elle conduit à un algorithme très efficace et très lisible.

Un exemple simple est la recherche de la racine d'une fonction continue (nous avons déjà étudié différentes méthodes plus haut pour résoudre ce genre de problèmes). C'est-à-dire le point pour laquelle la fonction f s'annule.

Supposons qu'une fonction soit croissante et continue sur un intervalle [a,b] et tel la racine cherchée soit entre a et b. Nous avons donc par le fait que la fonction soit croissante dans l'intervalle:

  (57.354)

et le fait que la racine se trouve dans l'intervalle:

  (57.355)

Nous  calculons:

  (57.356)

Si  alors la racine est dans l'intervalle  sinon elle est dans l'intervalle . Nous avons donc ramené le problème à une taille inférieure. Nous arrêterons l'algorithme quand la précision sera suffisante.

L'algorithme Maple est donné par:

zero:=proc(f,a,b,pre)
local M;
M:=f((a+b)/2);
if abs(M)<pre then 
     RETURN((a+b)/2)
elif M>0 then
     zero(f,a,(a+b)/2,pre)
else zero(f,(a+b)/2,b,pre)
     fi
end:

et ce ne sont que quelques exemples auxquels la méthode est applicable!!

Noter cette page:

Noté par 14 visiteur(s).

12345

PageRater 0.5.0