f
loadingPage en cours de chargement
    ACCUEIL | TÉLÉCHARGER | QCM | DON | ANNONCES | CHAT | FORUM | LIVRE D'OR | PARTENAIRES | CONTACT | BLOG
 
  Rechercher
  separation
  Introduction
  Arithmétique
  Algèbre
  Analyse
  Géométrie
  Mécanique
  Électrodynamique
  Atomistique
  Cosmologie
  Chimie
  Informatique Théorique
  Maths. Sociales
  Ingénierie
  separation
  Biographies
  Références
  Liens
  separation
  Humour
  Serveur d'exercices
  separation
  Parrains
2 connectés
News News :: Erreur Erreur :: Statistiques Visiteurs :: ClearType ClearType :: Imprimer Imprimer :: Bookmark and Share

Informatique Théorique

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)

57. MÉTHODES (ANALYSES) NUMÉRIQUES (1/2)

Dernière mise à jour de ce chapitre: 2017-12-31 17:59:29 | {oUUID 1.791}
Version: 3.10 Révision 43 | Avancement: ~60%
viewsvues depuis le 2012-01-01: 54'250

Table des matières LISTE DES SUJETS TRAITÉS SUR 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, pour certains problèmes aucune méthode directe n'est connue (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 à partir d'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 implique 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 ces perturbations font l'objet de la théorie du Chaos (classique ou quantique).

Sens de la vie

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

Remarques:

R1. Beaucoup de méthodes numériques utilisées en informatique se basent sur des raisonnements qui ont déjà été étudiés dans d'autres sections de ce site, méthodes sur lesquelles nous ne reviendrons pas.

R2. Ce chapitre étant à la limite entre l'ingénierie et la mathématique appliquée, nous avons décidé de donner parfois des exemples d'applications des outils développés.

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 (anecdotiques):

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

A2. Il n'existe pas de programmes sans de possibles erreurs (dues 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. Choisir une stratégie de développement (modulable ou autre).

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

Remarque: Il faudrait dans l'étape (7) veiller à respecter les normes de nommage et de représentation du code ainsi que des conventions (traditions) de fréquence d'insertion des commentaires.

M8. Tests (maintenance préventive)

Remarque: Le débogage (la gestion des erreurs) d'un programme et les tests de fonctionnement doivent prendre autant, voire plus, de temps que le développement du programme lui-même.

Lors du développement d'un programme scientifique, il peut être intéressant, voire même rigoureux de s'attarder sur 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 dans le pire des cas sur un jeu de données. La complexité est donc exprimée comme une fonction de la taille du jeu de données.

Mesurer la complexité exacte est peu pertinent car souvent trop complexe vu la taille des programmes. Pour éviter de calculer en détails la complexité d'un algorithme, nous repèrons l'opération fondamentale. Ces opérations fondamentales peuvent être: une assignation, une comparaison entre deux variables, une opération arithmétique entre deux variables, etc.

Ainsi, les hypothèses d'usage permettant un calcul de la complexité sont:

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

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

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

H4. Un seul processeur

Nous notons equation 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:

equation   (57.1)

C'est le plus petit temps qu'aura à exécuter l'algorithme sur un jeu de données (des 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" (la plus intéressante pour le praticien):

equation   (57.2)

C'est le plus grand temps qu'aura à exécuter 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é equation opérations. Cependant, cette complexité peut ne pas refléter le comportement usuel de l'algorithme, le pire cas ne pouvant se produire que très rarement, mais il n'est pas rare que le cas moyen soit quasi aussi mauvais que le pire des cas.

Si l'expression de la complexité au pire d'un algorithme contient plusieurs termes (additions ou soustractions), on ne garde que le terme qui croît le plus vite. Ainsi, un algorithme ayant une complexité du type:

equation   (57.3)

se résumera à avoir une complexité d'ordre equation.

- La "complexité moyenne":

equation   (57.4)

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 de 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(nk) pour equation 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(nk) pour equation. Une mauvaise complexité est du type equation ou equation.

exempleExemples:

E1. Évaluation d'un polynôme:

equation   (57.5)

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

equation   (57.6)

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

equation   (57.7)

Nous pouvons facilement montrer que cette factorisation maintient le même nombre d'additions equation 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 située à mi-colonne et de voir si nous y trouvons la valeur recherchée. Sinon, la recherche doit continuer sur le même mode opératoire dans la partie supérieure ou inférieure du tableau (cela dépend de l'ordre lexicographique).

L'algorithme est récursif et permet à chaque étape, de diviser par deux, 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, equation à l'étape 2, et plus généralement de equation à 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 equation. Nous en déduisons le nombre maximal d'étapes: c'est le plus petit k tel que equation, soit equation, soit:

equation   (57.8)

à comparer avec une recherche séquentielle (utile lorsque le tri est trop coûteux en ressources). Par exemple, 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 equation. Le gain est donc considérable (à condition que les données soient triées)!

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

- Lecture (itérations): O(n2)

- Calcul de la trace: equation

- Addition:equation telle que equation

- Multiplication: equation telle que equation

- 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 à:

equation   (57.9)

En faisant l'hypothèse que l'ordinateur utilisé effectue une opération élémentaire en equation 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

equation

equation

equation

equation

equation

     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(n2) complexité quadratique, O(n3) complexité cubique, O(kn) 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 (c.-à-d. 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 equation. Plus formellement, ce sont les problèmes pour lesquels nous pouvons construire une machine déterministe 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").

Remarque: Contrairement aux machines déterministes qui exécutent une séquence d'instructions bien déterminée, les machines non déterministes ont la remarquable capacité de toujours choisir la meilleure séquence d'instructions qui mène à la bonne réponse lorsque celle-ci existe. Il va sans dire qu'une telle machine ne peut pas exister autrement que dans l'esprit d'un théoricien (à moins qu'avec les ordinateurs quantiques...)! Néanmoins, comme nous le verrons par la suite, ce concept abstrait n'est pas sans intérêt et constitue en fait la base de toute la théorie de la NP-complétude.

Il importe de remarquer à ce stade de la discussion que la classe P est incluse dans la classe NP, nous écrivons equation, 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 non 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 une 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 puisse réussir cette opération. Autrement dit, nous ne savons pas si ce problème est dans P. En revanche, étant donnés des nombres premiers equation il est trivial de vérifier que equation: 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:

equation   (57.10)

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

Remarque: Si un problème est dans NP, alors il existera un algorithme en temps exponentiel pour le résoudre mais le contraire n'est pas toujours vrai (il faut donc être prudent).

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 un 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 à une ou des sous-routines est assez facilement faisable puisque réalisable, si elle existe, en un 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: et donc finalement a-t-on P=NP? Autrement dit, peut-on trouver en un temps polynomial ce que l'on peut prouver en temps polynomial?

Remarque: Ce problème est si important en informatique qu'il fait partie (arbitrairement) des 7 problèmes du millénaire, dont la résolution est primée 1 million de dollars par le Clay Mathematic Institute.

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 exprimé par [x], qui se lit "partie entière de x" et selon la lnorme ISO 80000-2:2009 devrait s'écrire: int 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:

equation   (57.11)

Nous notons aussi {x} pour désigner la partie fractionnaire de x. Nous avons ainsi:

equation   (57.12)

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

P1. equation, equation, equation où equation

P2. equation, lorsque equation

P3. equation, si equation

P4. equation

P5. equation si equation, equation si equation

P6.equation si equation

P7. Si equation, 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

equation  (57.13)

equation.

Pour P2, la somme est vide (ne contient aucun terme autrement dit) pour equation et, dans ce cas, on adopte la convention selon laquelle la somme vaut 0. Alors, pour equation, 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:

equation, equation   (57.14)

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

equation   (57.15)

En écrivant equation, où equation, nous avons

equation   (57.16)

equation.

Il s'ensuit que:

equation   (57.17)

si equation et:

equation

si equation . Et nous obtenons du même coup la démonstration P5.

Pour démontrer P6, nous écrivons:

equation   (57.18)

equation, et:

equation   (57.19)

equation. Nous obtenons ainsi:

equation   (57.20)

puisque equation. Par ailleurs:

equation   (57.21)

et nous avons ainsi le résultat.

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

equation   (57.22)

c'est-à-dire:

equation   (57.23)

soit le résultat attendu.

equationC.Q.F.D.

Remarque: La méthode d'arrondi de valeurs réelles est donnée dans le chapitre d'Économie.

ALGORITHME D'HÉRON

Soit à calculer la racine carrée:

equation   (57.24)

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

Démonstration:

Voici une pseudo-démonstration car historiquement l'algorithme a été construit sur des considérations purement intuitives (car 100 ans avant J.C. l'algèbre n'existait pas...). Dans les classes le résultat est donné en tant que définition et on en étudie la convergence.

Donc pour la démonstration, nous procéderons ainsi:

equation   (57.25)

Et l'astuce consiste à poser que:

equation   (57.26)

où il faut prendre comme valeur initiale equation (ce qui évidemment est en totale contradiction avec les développements mais l'on voit alors mieux pourquoi c'est une pseudo-démonstration...).

equationC.Q.F.D.

exempleExemple:

Soit à calculer:

equation   (57.27)

Nous prendrons donc equation et cela nous donne le tableau d'itérations suivant:

Itération
xi /2
A/2xi
xi+1
Écart
1
5
0.5
5.50
~2.3
2
2.750
0.90909
3.659 090 909
~0.49
3
1.82954
1.3664
3.196 005 083
~0.033
4
1.59800
1.5644
3.162 455 624
~0.0002
5
1.58122
1.5810
3.162 277 665
~0.5 10-8 
Tableau: 57.2  - Itérations pour l'algorithme d'Héron

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

equation   (57.28)

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 equation 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 atteler à cette tâche.

Nous définissons en géométrie le nombre equation  dit "pi", quel que 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:

equation   (57.29)

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:

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

Le principe de l'algorithme d'Archimède est le suivant: Soit equation 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 à montrer par induction que:

equation   (57.30)

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

equation et equation   (57.31)

Avec:

equation   (57.32)

Donc:

equation   (57.33)

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

equationC.Q.F.D.

Remarque: Il existe un très grand nombre d'algorithmes pour calculer equation. Celui présenté ci-dessus, sans être le plus esthétique, serait historiquement le premier.

CALCUL DU NOMBRE D'EULER

Hors la constante equation, 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 (cf. chapitre sur les Suites Et Séries), pour une fonction indéfiniment dérivable f donnée par:

equation   (57.34)

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

equation   (57.35)

nous avons:

equation   (57.36)

Donc en résumé:

equation   (57.37)

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

Remarque: Pour diminuer la complexité de cet algorithme, la factorielle peut être calculée avec la formule exposée ci-après.

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:

equation   (57.38)

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

equation   (57.39)

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

equation   (57.40)

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

equation   (57.41)

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

equation   (57.42)

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" ou "pivot de Gauss", qui est bien adaptée à la résolution des petits systèmes d'équations (jusqu'à 50 inconnues).

Remarques:

R1. La validité de certaines des opérations que nous allons effectuer ici pour résoudre les systèmes linéaires se démontre dans le chapitre traitant de l'Algèbre Linéaire. Au fait, pour être bref, le tout fait appel à des espaces vectoriels dont les vecteurs-colonnes sont linéairement indépendants.

R2. Les systèmes admettent une solution si et seulement si (rappel) le rang de la matrice augmentée est inférieur ou égal au nombre d'équations (cf. chapitre d'Algèbre Linéaire).

UNE ÉQUATION À UNE INCONNUE

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

equation   (57.43)

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:

equation   (57.44)

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 solution. 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:

equation   (57.45)

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

Remarque: Les notations usitées ci-dessus n'ont rien à voir avec le calcul tensoriel.

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:

equation   (57.46)

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

equation   (57.47)

Nous en déduisons:

equation   (57.48)

La transformation entre les deux systèmes:

equation   (57.49)

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

TROIS ÉQUATIONS À TROIS INCONNUES

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

equation   (57.50)

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:

equation   (57.51)

et ensuite résoudre l'équation:

equation   (57.52)

puis la deuxième:

equation   (57.53)

et enfin:

equation   (57.54)

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

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

equation   (57.55)

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 equationcomme 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 equation et non pas equation, 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):

equation   (57.56)

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

equation   (57.57)

Et de résoudre 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 equation tel que:

equation   (57.58)

Il est amené à equation 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 à 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.

Mises à part l'addition et la soustraction de polynômes que nous supposerons comme triviales (optimisation de la complexité mise à 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:

Soient:

equation   (57.59)

alors:

equation   (57.60)

avec:

equation   (57.61)

C'était simple...

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

Reprenons:

equation   (57.62)

mais en imposant cette fois-ci equation.

La division s'écrira donc nous le savons:

equation   (57.63)

avec:

equation   (57.64)

ou sinon equation.

Il est normalement connu d'avance (car démontré dans le chapitre de Calcul Algébrique) que nous avons:

equation   (57.65)

et:

equation   (57.66)

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:

equation   (57.67)

Pour démontrer l'expression des différents equation, nous avons préféré pour des raisons pédagogiques de passser par un exemple particulier (voir ci-dessous) dont nous généraliserons le résultat.

exempleExemple:

Soient:

equation   (57.68)

donc de ce que nous avons dit précédemment, il vient (point de départ):

equation   (57.69)

En utilisant le fait que (pour rappel):

equation   (57.70)

nous avons donc immédiatement:

equation   (57.71)

Ensuite (toujours en procédant de façon identique):

equation   (57.72)

et enfin:

equation   (57.73)

Donc de manière générale:

equation   (57.74)

comme:

equation   (57.75)

le premier reste est donc:

equation   (57.76)

Ensuite:

equation   (57.77)

Le deuxième reste est alors:

equation   (57.78)

etc.

Nous continuons jusqu'à ce que equation.

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, logistique ou autre.

Listons les techniques les plus utilisées dans les entreprises et administrations (dont les modèles mathématiques ne sont pas tous présentés ici):

- Modèle de régression linéaire à une ou plusieurs variables explicatives basé sur la méthode des moindres carrés avec variables explicatives binaires ou continues avec réponse réelle. Présenté en détail dans le présent chapitre (implicitement ce modèle contient les interactions entre variables).

- Modèle de régression linéaire gaussien (approche statistique de la régression linéaire basée sur la méthode des moindres carrés) avec variables explicatives binaires ou continues avec réponse réelle. Présenté en détail dans le présent chapitre.

- Modèles de régression non-linéaires avec variables explicatives binaires ou continues avec réponse réelle. Présenté en détail dans le présent chapitre dans le cas où ils peuvent être ramenés à des cas linéaires ou non mais alors sans interactions des variables explicatives. Sinon basé sur des techniques de type quasi-Newton ou de Gauss-Newton présentées dans le présent chapitre.

- Modèle de régression polynomial par la méthode des B-spline et du polynôme de collocation. Présenté en détail dans le présent chapitre.

- Modèle de régression logistique (régression binomiale) avec variables explicatives binaires, nominales (catégorielles) ou ordinales ou continues avec réponse bornée entre 0 et 1. Présenté sommairement et naïvement dans le présent chapitre.

- Modèle de régression de comptage de Poisson (Poisson MLE, PMLE, GLM) ou binomial négatif (binomial MLE et QGPMLE) avec variables explicatives binaires, nominales (catégorielles) ou ordinales ou continues avec réponse entière positive.

- Modèle de régression linéaire orthogonal (ou régression de Deming) qui est utilisé comme complément au test-T pour données appariées pour vérifier la stabilité des instruments de mesure dans les laboratoires. Il s'agit d'un cas où la variable explicative et expliquée sont entachées d'une incertitude (je présenterai la démonstration quand le temps me le permettra).

- Modèle de régression quantile (très utilité dans le domaine médical et économique) basé sur la même idée que la régression par la méthode des moindres carrées mais où on ne minimise non pas la somme des carées des écarts à la moyenne mais la somme des écarts absolus à un quantile (la médiane ou autre).

...et ... parmi un certain nombre de ces approches nous différencions les modèles mathématiques prenant en compte les données censurées des données non censurées (ce qui fait au final un paquet de théories/modèles à étudier).

Notons enfin que dans la régression linéaire les variables explicatives forment une expression linéaire mais cela ne signifie pas qu'elles sont elles-mêmes linéaires. Ainsi, si nous considérons les deux expressions ci-dessous:

equation   (57.79)

la première est linéaire dans les paramètres mais la seconde ne l'est pas!

Enfin, touchons un mot sur une technique parfois utilisée pour une interprétation qualitative des influences des variables explicatives via leurs coefficients pour des régressions linéaires simples ou multivariées.

Quand les amplitudes de certaines variables explicatives (continues!) ont des ordres de grandeurs très différents cela pose alors un gros problème d'interpértation de l'influence de chacune des variables via la lecture de leur coefficient et pose aussi comme problème la précision des calculs dans les algorithmes à cause des différences de grandeur d'ordre et donc génère aussi des erreurs d'arrondis!

L'idée traditionnelle consiste alors à centrer-réduire les valeurs des variables explicative ce qui aide grandement à l'interprétation de l'influence de ces mêmes variables (mais il faut alors laisser de côté l'exploitation de la valeur numérique de la variable expliquée). Par contre attention au piège!!! Une fois le modèle théorique obtenu à partir des variables centrées-réduites, les nouvelles valeurs à expliquer devront être obtenues en ayant au préalable centré et réduit les nouvelles valeurs explicatives mais en déduisant par les anciennes moyennes et en réduisant par les anciens écart-types de respectivement chacune des variables explicatives injectées dans le modèle!

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 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" (et en économie "variable exogène"...). 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) appelée aussi en éconoome "variable endogène". 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:

equation   (57.80)

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 vraiment une exception car 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!

Cependant, dans la pratique nous nous arrangeons pour linéariser les fonctions soit par des transformations algébriques élémentaires comme celles qu'utilisent les tableurs (par exemple Microsoft Excel) pour la linéarisation d'une fonction logarithmique en faisant un simple changement de variables:

equation   (57.81)

ou encore pour les fonctions puissance et exponentielles en faisant aussi une petite manipulation algébrique avec les propriétés des logarithmes démontrées dans le chapitre d'Analyse Fonctionnelle (à condition que a soit strictement positif):

equation   (57.82)

soit en faisant des approximations en série de Taylor (cf. chapitre de Suites Et Séries).

Remarques: Si nous cherchons à déterminer la valeur de y pour un x non mesuré et se situant au-delà de l'étendue de mesure d'origine, nous parlons alors "d'intervalle de prédiction pour x".

DROITE DE RÉGRESSION

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

equation   (57.83)

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).

Soient 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:

equation   (57.84)

la relation suivante:

equation   (57.85)

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'Économie):

equation   (57.86)

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

equation   (57.87)

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:

equation   (57.88)

ce qui donne b sous la forme:

equation   (57.89)

MÉTHODE DES MOINDRES CARRÉS

Du fait de l'erreur sur y, les points expérimentaux, de coordonnées equation, 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) consistera dans le cadre d'étude particulier qui nous intéresse à chercher les valeurs des paramètres a, b qui rendent minimale la somme des carrés des écarts ei résiduels (SSr: Sum of Squared Residuals) entre les valeurs observées equation et les valeurs calculées théoriques de equation. Nous parlons alors de "méthode des moindres carrés des écarts d'ordonnées":

equation   (57.90)

n est le nombre de points et:

equation   (57.91)

d'où autrement écrit:

equation   (57.92)

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:

equation   (57.93)

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 (nous y reviendrons plus loin). Dans notre exemple equation est la grandeur scalaire qui fait office de multiplicateur de Lagrange.

Soit après simplification:

equation   (57.94)

Le système ci-dessus est appelé "système des équations normales". C'est un système linéaire de deux équations à deux inconnues. Notons pour simplifier:

equation   (57.95)

Le système devient:

equation   (57.96)

De la deuxième équation, nous tirons:

equation  (57.97)

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

equation   (57.98)

De là nous avons:

equation   (57.99)

Ainsi, les expressions de la pente et de l'ordonnée à l'origine de l'équation recherchée sont:

equation   (57.100)

Les deux dernières relations sont utilisées par la majorité des tableurs comme par exemple dans la version française de Microsoft Excel 11.8346 lors de l'utilisation de la fonction REGRESSION( ). Le terme b (l'ordonnée à l'origine) peut être obtenu directement avec la fonction ORDONNEE.ORIGINE( ) et a avec la fonction PENTE( ) et l'ensemble avec la fonction DROITEREG( ).

Voici si jamais un petit listing intéressant de quelques cas très pratiques avec ce tableur:

- Pour une droite:

a: =PENTE(y,x)
b: =ORDONNEE.ORIGINE(y,x)

- Pour une fonction logarithmique (nous retrouvons le changement de variable donné au début):

a: =INDEX(DROITEREG(y,LN(x)),1)
b: =INDEX(DROITEREG(y,LN(x)),1,2)

- Pour une fonction puissance (nous retrouvons le changement de variable donné au début):

a: =EXP(INDEX(DROITEREG(LN(y),LN(x),,),1,2))
b: =INDEX(DROITEREG(LN(y),LN(x),,),1)

- Pour une fonction exponentielle (nous retrouvons aussi le changement de variable donné au début):

a: =EXP(INDEX(DROITEREG(LN(y),x),1,2))
b: =INDEX(DROITEREG(LN(y),x),1)

Remarque: Il faut bien avoir en tête que la droite des moindres carrés, qui permet de résumer au mieux le nuage de points des observations, passe nécessairement par le centre de gravité du nuage, c'est-à-dire par un point moyen qui ne correspond que rarement à une observation (moyenne des abscisses et moyenne des ordonnées).

ANALYSE DE LA VARIANCE DE LA RÉGRESSION UNIVARIÉE

Avant de commencer il est important que le lecteur abandonne tout de suite le possible réflexe qui serait de tenter de ramener par analogies successives l'ANOVA de la régression que nous allons voir maintenant à l'ANOVA catégorielle que nous avons traité dans le chapitre de Statistiques!

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

equation   (57.101)

soit sous forme discrète:

equation   (57.102)

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

equation   (57.103)

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

equation   (57.104)

Soit en soustrayant les deux dernières relations:

equation   (57.105)

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

equation   (57.106)

En remplaçant b par sa valeur:

equation   (57.107)

nous avons alors:

equation   (57.108)

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

equation   (57.109)

Soit après réarrangement:

equation   (57.110)

Revenons maintenant à:

equation   (57.111)

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

equation   (57.112)

soit:

equation   (57.113)

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

equation   (57.114)

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:

equation   (57.115)

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:

equation   (57.116)

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

equation   (57.117)

D'où l'égalité que nous réutiliserons plusieurs fois:

equation   (57.118)

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:

equation   (57.119)

où la première somme après l'égalité est très souvent appelée "manque d'ajustement" ("lack of fit" en anglais). Plus explicitement:

equation   (57.120)

Cette dernière relation peut se représenter (en gros...) graphiquement sous la forme suivante:

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

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

equation   (57.121)

qui est simplement une autre manière d'écrire la décomposition de la variance (implicite):

equation   (57.122)

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):

equation   (57.123)

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 (était défini) par:

equation   (57.124)

Ou sinon puisque nous avons démontré plus haut que (se rappeler que la variance indiquée est la variance biaisée... ):

equation   (57.125)

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

equation   (57.126)

Donc nous en déduisons):

equation   (57.127)

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

equation   (57.128)

Remarque: Cette formulation du coefficient de corrélation est extrêmement utile car, contrairement à la formulation statistique, cette dernière se généralise immédiatement à la régression linéaire multiple que nous verrons un peu plus loin.

Démonstration:

Nous partons donc de:

equation   (57.129)

et puisque nous avons montré que:

equation   (57.130)

Alors:

equation   (57.131)

equationC.Q.F.D.

et dans le cadre des modèles de régression voici quelques cas typiques de la valeur de ce coefficient avec un corrélation linéaire pour les deux premières lignes et non linéaire pour la troisième ligne:

equation
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:

equation   (57.132)

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

Enfin, remarquons que l'ordonne à l'origine n'intervient pas dans la valeur du coefficient de corrélation puisque (propriété de bilinéarité de la covariance démontrée dans le chapitre de Statistiques):

equation   (57.133)

La suite viendra quand j'aurai trouvé une démonstration acceptable pédagogiquement des valeurs des degrés de liberté du test de Fisher de la régression linéaire....

MODÈLE LINÉAIRE GAUSSIEN

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

equation   (57.134)

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

equation et equation   (57.135)

donc autrement dit:

equation   (57.136)

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 (expérimentale):

equation   (57.137)

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

equation   (57.138)

Raison pour laquelle le modèle s'appelle "modèle linaéire gaussien".... Explicitement, cela donne donc:

equation   (57.139)

Nous allons opter pour le symbole ~ pour dire "suit la loi..." dans ce qui va suivre immédiatement afin d'éviter toute confusion:

equation   (57.140)

Ce qui équivaut graphiquement à avoir (donc normalement l'analyse statistique de la régression gaussienne ne se fait que si et seulement si nous avons prises plusieurs mesures de la variable dépendante pour des valeurs fixes et identiques des variables indépendante !!!!!!!!):

equation
Figure: 57.4 - Représentation graphique du principe de distribution Normale

Attention! Puisque le modèle est gaussien, la variable à expliquer a son domaine de définition qui est donc non borné (support de la loi Normale). Certaines entreprises (cabinets d'audits en particulier....) souhaitant parfois créer un simple modèle linéaire pour modéliser une probabilité (qui pour rappel est bornée [0,1]) typiquement de faillite/défaut en fonction de différents facteurs (variables explicatives) doivent donc au préalable transformer la probabilité en valeur Z de la loi Normale par simple inversion. Ce type d'approche est alors appelé "modèle linéaire Z-score".

Presque tous les logiciels d'analyse statistique proposent de faire un tracé (graphique) des résidus en fonction des valeurs de x. Ainsi, deux graphiques du type suivant améneraient à rejeter le modèle linéaire gaussien:

equation
Figure: 57.5 - Exemples de "plot" de résidus qui n'annoncent rien de bon...

Dans la figure ci-dessus, le graph en haut à gauche est ce dont on doit s'attendre pour pouvoir appliquer les tests statistiques du modèle linéaire gaussien. Le graphique en haut à droite montre que la variance des résidus n'est pas constante et qu'il y à violation de l'hypothèse d'homoscédasticité. Le graph en bas à gauche montre que la variance est constante mais cependant que notre modèle a des variables explicatives manquantes qui rajoutées expliqueraient ce glissement qui croit linéairement. Le graph en bas à droite indique une variance constante mais que le modèle serait plutôt non linéaire que linéaire.

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 Statistiques "l'homoscédasticité" (tandis que le fait que les variances ne soient pas égales s'appelle pour rappel "l'hétéroscédasticité").

Remarque: La majorité des logiciels (dont Microsoft Excel 11.8346) proposent un graphique qui montre les résidus en fonction des valeurs des ordonnées x. Bien évidemment, il vaut mieux que les points représentant les résidus ne soient pas trop divergents... sinon quoi l'hypothèse d'homoscédasticité est violée.

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

equation   (57.141)

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

Avant de voir la démonstration faisons un rappel et donnons quelques définitions des variables que nous avons déjà manipulées et des nouvelles que nous allons manipuler (si le vocabulaire semble technique au lecteur alors il devra au préalable lire ou relire le chapitre de Probabilités et de Statistiques):

Variable

Description

a

Pente (coefficient directeur) du modèle linéaire de la méhode des moindres carrés. Il s'agit d'une valeur ponctuelle (déterministe).

b

Ordonnée à l'origine du modèle linéaire de la méthode des moindres carrés. Il s'agit d'une valeur ponctuelle (déterministe).

A

Variable aléatoire de la pente (coefficient directeur) selon l'approche statistique du modèle gaussien et dont a est une réalisation. L'espérance de A étant un estimateur non biaisé de equation.

B

Variable aléatoire de l'ordonnée à l'origine selon l'approche statistique du modèle linéaire gaussien et dont b est une réalisation. L'espérance de B étant un estimateur non biaisé de equation

equation

Espérance non biaisée de la variable A représentant la pente (coefficient directeur) dans le cadre de l'approche statistique du modèle linéaire gaussien.

equation

Espérance non biaisée de la variable B représentant l'ordonnée à l'origine dans le cadre de l'approche statistique du modèle linéaire gaussien.

Tableau: 57.3 - Rappel des notations pour l'étude de la régression linéaire

Démonstration:

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ée plus haut comme étant le rapport de la covariance et de la variance):

equation   (57.142)

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

equation   (57.143)

Donc nous différencions les valeurs aléatoires et non aléatoires des coefficients par le passage de la notation minuscule à majuscule.

Tenant compte de ce que la variable expliquée théorique considérée comme la réalisation d'une variable aléatoire est donnée par:

equation   (57.144)

nous pouvons mettre A sous la forme:

equation   (57.145)

et B:

equation   (57.146)

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

equation   (57.147)

ce qui montre que notre hypothèse initiale pour l'expression de A n'est pas trop mauvaise... puisque l'espérance de A est visiblement un estimateur non biaisé de a.

Respectivement, pour B il vient:

equation   (57.148)

avec la même conclusion.

equationC.Q.F.D.

Donc l'espérance de A et B sont bien des estimateurs sans biais (donc de variance minimale pour rappel) de equation. Comme ce sont des estimateurs, dans la littérature spécialisée, equation sont souvent notés equation et dès lors il vient la notation courante alternative:

equation   (57.149)

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

equation   (57.150)

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

equation   (57.151)

Soit si n est assez grand nous écrirons:

equation   (57.152)

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

equation   (57.153)

donc par propriété de linéarité de la loi Normale, les variables aléatoires A et B suivent donc aussi une loi Normale.

Du rappel de cette hypothèse, il découle immédiatement (cf. chapitre de Statistiques):

equation   (57.154)

Dès lors:

equation   (57.155)

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

equation   (57.156)

Nous avons finalement:

equation   (57.157)

où évidemment la notation de la variance au dénominateur est très abusive (puisque x n'est pas une variable aléatoire) mais très pratique pour condenser les écritures.

Le problème réside maintenant dans la détermination de equation. É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 Statistiques en ce qui concerne les estimateurs:

equation   (57.158)

puisque la loi Normale est centrée pour les résidus donc equation... 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 (donc en notant tout avec des minuscules):

equation   (57.159)

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" qui s'obtient avec la version anglaise de Microsoft Excel via le carré fonction STEYX( ).

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

equation   (57.160)

relations qui ne sont donc valables que pour une régression linéaire à une variable explicative (et sous les hypothèses de construction du modèle linéaire Gaussien). Sachant par construction de l'hypothèse de départ que A et B suivent une loi Normale d'espérance respective equation et dont la variance est donnée juste ci-dessus, nous connaissons donc entièrement la distribution qui les caractérise.

Ce qui est sympathique 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 Statistiques l'intervalle de confiance suivant:

equation   (57.161)

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

equation   (57.162)

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

equation   (57.163)

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

equation   (57.164)

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 Microsoft Excel 11.8346 pour chaque coefficient). La démarche est la même pour l'ordonnée à l'origine.

Attention, si la variable explicative est aussi une variable aléatoire, nous utilisons alors:

equation   (57.165)

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", l'idée est la même mais les calculs sont plus longs (pas le courage de les faire...).

Enfin, rappelons que nous avons obtenu pour le coefficient de corrélation empirique:

equation   (57.166)

Nous avons alors in extenso le fameux intervalle de confiance du coefficient de corrélation:

equation   (57.167)

Il faut savoir qu'étant donné que calculer l'intervalle de confiance pour la pente ou pour le coefficient de corrélation est une démarche équivalente, nombreux sont les logiciels (Tanagra, Minitab, Microsoft Excel, etc.) qui donnent la valeur de la distribution de Student et la valeur critique de celle-ci uniquement pour la pente et supposent que le lecteur sait qu'il en est de même pour le coefficient de corrélation.

TEST DU COEFFICIENT DE CORRÉLATION DE PEARSON

Le  calcul obtenu plus haut pour l'intervalle de confiance du coefficient de corrélation est un peu pénible dans la pratique. C'est pour cette raison que de nombreux praticiens et logiciels de statistiques implémentent une alternative très simple communiquée uniquement sous la forme minimale qu'est la p-value.

Pour voir cette approche, rappelons que nous avons vu dans le chapitre de Statistiques que (cette fois-ci nous allons adopter la notation correcte...):

equation   (57.168)

Et que nous avons vu juste plus haut que:

equation   (57.169)

De la même manière, nous avons l'estimateur du coefficient de corrélation de Pearson qui est (en adaptant ici aussi les notations d'usage possibles que nous pouvons trouver dans la littérature):

equation   (57.170)

et donc:

equation   (57.171)

Le test d'hypothèse que nous voulons faire est alors:

equation   (57.172)

et donc équivalent à:

equation   (57.173)

l'hypothèse nulle étant évidemment que le coefficient de corrélation de Pearson est statistiquement significativement différent de zéro. Il s'agit donc d'un test bilatéral!

Pour trouver une forme simple du test, rappelons que nous avons obtenu:

equation   (57.174)

ainsi que:

equation   (57.175)

ce qui nous amène en mixant les deux à avoir:

equation   (57.176)

d'où:

equation   (57.177)

Or, nous avons aussi montré que si n est assez grand:

equation   (57.178)

Mais si n est assez petit, nous écrirons donc:

equation   (57.179)

Donc:

equation   (57.180)

et avec l'hypothèse nulle equation, il vient:

equation   (57.181)

Il faut faire attention à l'utilisation de ce test, appelé souvent et logiquement "test-t de Student de la pentre de régression univariée" (en anglais: "t-test for coefficient slope") , suivant si le coefficient de corrélation de Pearson est négatif ou positif et ne pas oublier qu'il est bilatéral.

exempleExemples:

E1. Nous avons calculé pour une série de données un coefficient de corrélation de Pearson positif R valant 0.298 et la variable explicative comporte 7 valeurs. Nous avons donc avec la versio anglaise de Microsoft Excel 11.8346 la p-value qui est (nous retrouvons exactement la même valeur qu'un logiciel comme Minitab 15.1.1):

=2*(1-T.DIST(0.298/SQRT((1-0.298^2)/(7-2));7-2;1))=2*(1-0.741869)=0.51626

Dans le cas présent nous ne rejetons donc pas l'hypothèse nulle comme quoi le coefficient de corrélation de Pearson est nul à un niveau de confiance de 5%.

E2. Nous avons calculé pour une série de données un coefficient de corrélation positif de Pearson R valant -0.084 et la variable explicative comporte 19 valeurs. Nous avons donc avec la version anglaise Microsoft Excel 11.8346 la p-value qui est (nous retrouvons exactement la même valeur qu'un logiciel comme Minitab 15.1.1):

=2* T.DIST((-0.084)/SQRT((1-(-0.084)^2)/(19-2));19-2;1)=2*(1-0.366)=0.74186

Dans le cas présent nous ne rejetons donc pas l'hypothèse nulle comme quoi le coefficient de corrélation de Pearson est nul à un niveau de confiance de 5%.

Ce petit piège fait que finalement on prend la valeur absolue du coefficient de corrélation de Pearson et on fait le calcul toujours d'une seule manière

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 (nous n'écrivons plus les indices pour gagner du temps):

equation   (57.182)

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:

equation   (57.183)

alors:

equation   (57.184)

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

equation   (57.185)

d'où:

equation   (57.186)

Donc nous avons:

equation   (57.187)

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

equation   (57.188)

et comme Y est distribué selon une loi Normale dont l'estimateur non biaisé de la moyenne ainsi que l'écart-type sont donnés par la relation antéprécédente il vient immédiatement:

equation   (57.189)

Dans la pratique il faut aussi vérifier que ce rapport suit effectivement une loi Normale pour pouvoir faire les intervalles de confiance et tests statistiques qui s'en suivent.

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

equation    (57.190)

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 obtenue et rappelons que:

equation   (57.191)

Donc:

equation   (57.192)

Or comme:

equation   (57.193)

Alors:

equation   (57.194)

et donc:

equation   (57.195)

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

equation   (57.196)

Donc:

equation   (57.197)

Soit:

equation   (57.198)

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

equation   (57.199)

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:

equation   (57.200)

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

equation   (57.201)

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 Microsoft Excel, IBM SPSS, Minitab, Statistica, etc. donnent la valeur de la distribution de Student ainsi que l'intervalle de confiance à un niveau de equation 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):

equation
Figure: 57.6 - Capture du logiciel Minitab 15 pour l'intervalle de confiance

Le lecteur aura remarqué que:

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

2. L'intervalle de confiance est courbe ce qui est parfois considéré comme relevant 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:

equation   (57.202)

avec la variance:

equation   (57.203)

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

equation   (57.204)

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:

equation   (57.205)

Ce qui nous donne avec un logiciel comme Minitab 15:

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

où nous distinguons bien l'intervalle 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:

equation   (57.206)

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

equation   (57.207)

d'où autrement écrit:

equation   (57.208)

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:

equation   (57.209)

Soit après simplification:

equation   (57.210)

Enfin:

equation   (57.211)

Vous pouvez vérifier aussi avec un tableur quelconque (Microsoft 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 indépendantes (donc en absence de "colinéarité"!) - et donc de p paramètres à déterminer - reliées par une loi linéaire de la forme générale:

equation   (57.212)

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

Observations

equation

equation

...

equation

1

equation

equation

...

equation

2

equation

equation

...

equation

equation

equation

equation

equation

equation

n

equation

equation

...

equation

Au fait, pour estimer les paramètres equation (valeurs estimées que nous noterons equation 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 vue plus haut pour la régression linéaire simple.

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

equation   (57.213)

avec le modèle théorique estimé:

equation   (57.214)

Donc il nous faut minimiser:

equation   (57.215)

Les parenthèses ci-dessus peuvent être réécrites sous la forme:

equation   (57.216)

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

equation   (57.217)

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

equation   (57.218)

Soit explicitement:

equation   (57.219)

à remarquer que nous avons:

equation   (57.220)

et:

equation   (57.221)

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 qui va suivre sont des produits scalaires!!!) la fonction quadratique multivariée de coefficients de vecteurs:

equation   (57.222)

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

equation   (57.223)

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

equation   (57.224)

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

equation   (57.225)

Soit:

equation   (57.226)

Rappelons avant d'aller plus loin que:

equation   (57.227)

Donc:

equation   (57.228)

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

equation   (57.229)

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

equation   (57.230)

impose forcément que:

equation   (57.231)

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

equation   (57.232)

La matrice equation que nous retrouverons souvent dans le domaine de la régression linéaire multiple et dans le domaine des plans d'expérience (cf. chapitre de Génie Industriel) est appelée "matrice d'information" et equation est appelée "matrice de dispersion" pour une raison qui paraîtra évidente un peu plus loin.

Remarque: Nous disons que la régression est "balancée et orthogonale" lorsque la matrice d'information est diagonale. Nous disons que la régression est "orthogonale" lorsque la sous-matrice de la matrice d'information excluant la première ligne et la première colonne est orthogonale. Nous disons que la régression est juste "balancée" lorsque toutes les valeurs de la première ligne et de la première colonne exceptée celle en étant à l'intersection sont nulles.

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

equation   (57.233)

Donc en ne supposant que 2 observations, nous avons alors:

equation   (57.234)

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 equation en equation:

equation   (57.235)

Nous avons dans le cas présent:

equation   (57.236)

Nous avons donc:

equation   (57.237)

et 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 equation (cf. chapitre Algèbre Linéaire)

equation   (57.238)

Donc:

equation   (57.239)

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...):

equation   (57.240)

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:

equation   (57.241)

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

Remarque: Le lecteur intéressé par la mise en pratique de ces résultats pourra, au même titre que pour la régression simple, se référer au serveur d'exercices - section Méthodes Numériques - où il y a des exemples pratiques avec Microsoft Excel.

Évidemment avec la régression linéaire multiple, nous pouvons maintenant, moyennant une astuce, faire de la régression linéaire de polynômes (nous verrons plus loin comment appliquer cependant directement la méthode des moindres carrés sur un polynôme). Effectivement, considérons un polynôme de la forme:

equation   (57.242)

Ce que nous pouvons considérer sous la forme:

equation   (57.243)

Donc dans des tableurs du type Microsoft Excel, il suffit d'utiliser l'utilitaire d'analyse de la régression avec en variable d'entrée la colonne x, une deuxième colonne que l'on aura pris soin de créer avec le carré de x et qui sera donc considérée comme la variable explicative w et enfin la troisième colonne que l'on aura aussi pris le soin de créer comme étant le cube de x et qui sera donc considérée comme la variable explicative z.

Nous pouvons aussi obtenir directement les coefficients de polynômes avec les fonctions déjà susmentionnées (mais vous n'aurez pas tous les résultats pertinents de l'Utilitaire d'Analyse). Par exemple pour un polynôme du deuxième degré:

a: =INDEX(REGRESSION(y,x^{1,2}),1)
b: =INDEX(REGRESSION(y,x^{1,2}),1,2)
c: =INDEX(REGRESSION(y,x^{1,2}),1,3)

et pour un polynôme du troisième degré:

a: =INDEX(REGRESSION(y,x^{1,2,3}),1)
b: =INDEX(REGRESSION(y,x^{1,2,3}),1,2)
c: =INDEX(REGRESSION(y,x^{1,2,3}),1,3)
d: =INDEX(REGRESSION(y,x^{1,2,3}),1,4)

C'est cette astuce qui permet de comprendre pourquoi et comment le coefficient de corrélation linéaire s'applique aussi aux polynômes dans la majorité des tableurs. Cependant nous pouvons avoir une approche plus directe qui ne nécessite pas de transformation mais qui dès lors est un peu plus longue.

Revenons maintenant sur le modèle linéaire gaussien avec une notation un peu plus rigoureuse et adaptée au cas multilinéaire et ce en particulier pour mettre en évidence la distinction entre estimateurs de la pente de la régression et valeurs exactes:

equation   (57.244)

Mais sous cette convention d'écriture, nous avons alors:

equation   (57.245)

et écrivons cela sous forme vectorielle:

equation   (57.246)

Maintenant, utilisons la technique du maximum de vraisemblance (cf. chapitre de Statistique) et cherchons la valeur des coefficients qui maximise donc:

equation   (57.247)

Ce que nous pouvons écrire sous forme matricielle:

equation   (57.248)

En prenant comme dans le chapitre de Statistique la log-vraisemblance pour faciliter la suite des calculs et en utilisant la propriété de la transposée des matrices démontrée dans le chapitre d'Algèbre Linéaire, il vient:

equation   (57.249)

Maintenant cherchons l'expression de equation qui maximise la log-vraisemblance. Il vient alors:

equation   (57.250)

Utilisons la propriété de la transposée démontrée dans le chapitre d'Algèbre Linéaire:

equation   (57.251)

Il vient:

equation   (57.252)

Pour des raisons qui paraîtront évident un peu plus loin nous choisisson la deuxième équivalence. Dès lors, nous avons (et en se rappelant que chercher un optimum revient à avoir la dérivée partielle nulle):

equation   (57.253)

Soit:

equation   (57.254)

C'est-à-dire après réarrangement:

equation   (57.255)

Soit exactement la même expressions que celle que nous avons obtenue juste un peu plus haut avec la méthode des moindres carrés dans le cas multilinéaire et qui était pour rappel:

equation   (57.256)

Ce qui montre que l'approche statistique par le modèle de maximum de vraisemblance du modèle linéaire permet de retomber dans le cas multilinéaire (et donc linéaire simple) sur les résultats de la méthode des moindres carrés. Cela peut sembler presque divin...

Enfin, indiquons que nous trouvons in extenso la fameuse "matrice chapeau" ou "matrice d'influence" H ("hat matrix" en anglais) qu relie à elle seule toute l'information entre les valeurs réelles et les valeurs théoriques (et ne dépendant que de X!) expliquées et donc in extenso l'erreur du modèle:

equation   (57.257)

Nous avons quelque chose d'intéressantà observer:

equation   (57.258)

Par définition, l'influence de l'observation i sur la régression est mesurée par la norme du vecteur et appelée "effet de levier" ("leverage score en anglais") et noté:

equation   (57.259)

C'est une méthode qualtitative (disponible dans de nombreux logiciels de statistiques) pour juger de l'influence de pointes qui pourraient être considérés comme abérrants. L'idée étant de comparer les valeurs des effets de levier entre eux. Un point admettent un levier dépassant deux fois sa moyenne (trois fois pour les petits échantillons) est suspect.

RÉGRESSION POLYNOMIALE

Nous allons voir maintenant comment déterminer par exemple le meilleur polynôme du deuxième degré qui passe par un nombre quelconque de points mais sans transformer la fonction contrairement à ce que nous venons de faire juste avant! Comme nous aimons bien la physique sur ce site, nous allons reprendre un cas classique de la cinématique afin de joindre l'utile à l'agréable...

Considérons donc que nous recherchons un polynôme du deuxième degré de la forme:

equation   (57.260)

sachant que la méthode est très facilement applicable à des polynômes de degré supérieur.

Relation qu'il est d'usage d'écrire dans le domaine de la régression polynomiale sous la forme suivante:

equation   (57.261)

i représente le nombre de points à notre disposition.

Pour la suite, nous allons nous baser encore une fois sur la méthode des moindres carrés. En d'autres termes, nous cherchons les coefficients equation qui minimisent l'erreur:

equation   (57.262)

et nous réattaquons avec des dérivées partielles pour chaque coefficient:

equation   (57.263)

Soit après simplification et un petit réarrangement:

equation   (57.264)

De même:

equation   (57.265)

Soit après simplification et un petit réarrangement:

equation   (57.266)

Et enfin:

equation   (57.267)

Soit après simplification et un petit réarrangement:

equation   (57.268)

Donc en utilisant la notation de l'algèbre linéaire, nous avons finalement le système suivant à résoudre:

equation   (57.269)

et donc il suffit de résoudre ce simple système linéaire soit à la main en utilisant les relations démontrées dans le chapitre d'Algèbre Linéaire, soit avec un simple tableur (du type Microsoft Excel).

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 feront 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 (non bornée) et une variable explicative quantitative.

Lorsque la "variable de classe" Y à expliquer est binaire (oui-non, présence-absence, 0-1, etc.), l'idée est de nous approcher dans un premier temps de celle-ci par une fonction de probabilité equation, qui nous donne à l'opposé la probabilité d'appartenir à la classe equation, 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 equation appartiendront à la classe 1 (et inversement dans le cas contraire).

Remarques:

R1. Au fait, la régression logistique n'est qu'une simple loi de distribution de probabilités dans le cas qui nous intéresse (nous verrons une autre régression logistique dans le chapitre d'Économie lors de notre étude des séries temporelles et encore une autre dans le chapitre de Dynamique des Populations).

R2. Il n'est évidemment pas possible d'appliquer systématiquement la régression logistique à n'importe quel type d'échantillon de données! Parfois il faut chercher ailleurs...

R3. Lorsque le nombre de modalités est égal à 2, nous parlons de "variable dichotomique" (oui-non) ou d'un "modèle dichotomique" ou encore de "variable indicatrice"; s'il est supérieur à 2, nous parlons de "variables polytomiques" (régression logistique polytomique). Le modèle logit est donc un "modèle dichotomique".

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.

exempleExemple:

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

equation
Figure: 57.8 - Partie de la population aux é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% sont "en cours". Eh bien la courbe représente simplement la proportion des deux classes pour l'âge donné. Il est même parfois indiqué la grosseur 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é (d'autant plus que la variable expliquée d'un modèle linéaire balaie tout l'ensemble des réels et ne se borne pas à l'intervalle [0, 1]).

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ésentation, 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 equation 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" (probability-unit).

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

equation   (57.270)

pour une variable de prédiction (explicative) xP est donc une probabilité bien évidemment comprise entre 0 et 1. Nous verrons un peu plus loin la raison historique de ce choix.

Nous voyons immédiatement que cette dernière relation étant la primitive de la fonction de distribution (voir la démonstration un peu plus bas) il s'agit donc bien d'une fonction de répartition puisque:

equation   (57.271)

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

equation   (57.272)

Lorsque nous optons pour la 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!

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) la fonction de répartition logistique suivante:

equation   (57.273)

il en découle donc la fonction de distribution:

equation   (57.274)

Évidemment, en fonction de la valeur de la probabilité nous associons à l'âge x le fait de ne pas avoir fini ses études (état associé à la valeur binaire: 0) ou de les avoir finies (état associé à la valeur binaire: 1).

Remarque: Après un petit changement de variable, nous retrouvons la loi logistique telle que définie sur Wikipédia:

equation   (57.275)

Indiquons que si a est posé comme unitaire, et b comme nul, nous avons alors la "loi logistique standard" donnée par:

equation   (57.276)

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:

equation   (57.277)

comme étant la variable aléatoire alors nous pouvons calculer formellement l'espérance de la loi logistique (le lecteur aura peut-être remarqué que c'est comme si nous posions a comme valant 1 et b comme valant 0). Effectivement, partant de:

equation   (57.278)

il vient alors:

equation   (57.279)

qui après une intégration numérique donne 0. Nous obtenons alors aussi le résultat suivant (si quelqu'un a le résultat analytique détaillé nous sommes preneurs!):

equation   (57.280)

Ainsi, nous voyons que si nous posons:

equation   (57.281)

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

La fonction de répartition:

equation   (57.282)

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

equation   (57.283)

d'où:

equation   (57.284)

Au fait c'est ici que réside l'astuce d'origine historique de la régression logistique. Nous transformons une variable P prenant ces valeurs dans 0 à 1 à l'aide du logarithme du rapport P/(1-P) en une variable prenant ses valeurs sur l'ensemble des réels et dès lors il est possible d'y associer une régression linéaire. Certes c'est empirique mais l'idée est bonne!

Ce que certains écrivent aussi...:

equation   (57.285)

Le résultat de cette dernière transformation est appelé "logit". Il est égal au logarithme de "l'odds" (sur lequel nous reviendrons très vite plus en détails):

equation   (57.286)

Il s'agit donc simplement du rapport d'un probabilité d'un événement sur la probabilité de l'événement complémentaire (ou contraire).

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.

Remarque: L'odds est également appelé "cote" par analogie à la cote des chevaux au Tiercé. Par exemple, si un cheval a 3 chances sur 4 d'être gagnant (donc in extenso 1 chance sur 4 d'être non-gagnant) sa cote est de 3 contre un 1, soit un odds égal à 3. On peut également introduire le concept de "odds ratio" (OR) pour le rapport des cotes qui est un indicateur très très utilisé en médecine. Ainsi, si la survenue d'un évènement dans un groupe A est p, et q dans le groupe B, le rapport des cotes est alors simplement donné par:

equation

L'odds ratio est toujours par construction supérieur ou égal à zéro. Si l'odds ratio est proche de 1, l'événement est indépendant du groupe, s'il est supérieur à 1, l'événement est plus fréquent dans le groupe A que dans le groupe B, s'il est inférieur à 1 l'événement est moins fréquent dans le groupe A que dans le groupe B et ainsi de suite...

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 partions de la taille (hauteur) d'une personne pour prédire si cette 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:

equation   (57.287)

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:

equation  et  equation   (57.288)

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

equation   (57.289)

exempleExemple:

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

- 1ère succursale:

Montant crédit

Payé

Non Payé

27'200

1

9

27'700

7

3

28'300

13

0

28'400

7

3

29'900

10

1

Tableau: 57.4  - Scoring débiteurs par montant de crédit succursale 1

- 2ème succursale:

Montant crédit

Payé

Non Payé

27'200

0

8

27'700

4

2

28'300

6

3

28'400

5

3

29'900

8

0

Tableau: 57.5  - Scoring débiteurs par montant de crédit succursale 2

- 3ème succursale:

Montant crédit

Payé

Non Payé

27'200

1

8

27'700

6

2

28'300

7

1

28'400

7

2

29'900

9

0

Tableau: 57.6  - Scoring débiteurs par montant de crédit succursale 3

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

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

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

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

Montant crédit

Proportion P

27'200

=2/27=0.0741

27'700

=17/24=0.7083

28'300

=26/30=0.8667

28'400

=19/27=0.7037

29'900

=27/28=0.9643

Tableau: 57.7  - Proportion des bons débiteurs

Ce qui donne graphiquement en Kilo-francs:

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

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

equation   (57.290)

Ce qui donne:

Montant crédit

Proportion P

Logit

27'200

0.0741

-2.5257

27'700

0.7083

0.8873

28'300

0.8667

1.8718

28'400

0.7037

0.8650

29'900

0.9643

3.2958

Tableau: 57.8  - Proportion des bons débiteurs et Logit

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

equation
Figure: 57.10 - Logit des bons débiteurs en fonction du montant du crédit

avec pour équation:

equation   (57.291)

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

equation   (57.292)

equation

Ainsi, il est possible de dire dans cet exemple, quelle 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).

Un logiciel come Minitab donnera une sortie sympathique qui est la "matrice de confusion". Elle permet de comparer le modèle à la réalité avec un cut-off traditionnel fixé à 50% (évidemment si le modèle correspondant parfaitement à l'échantillon, la matrice ci-dessous sera diagonale):

equation
Figure: 57.11 - Matrice de confusion de Minitab 16.1.1

où les bons débiteurs ont la valeur 1 (pas de risque de crédit) et les mauvais la valeur 0 (risque de crédit (nous détaillerons plus loi comment obtenir cette matrice avec un tableur). Enfin indiquons que dans des nombreux logiciels de Data Mining, il est d'usage de définir le "score" du modèle comme étant (dans le cas particulier mais facilement généralisable à une variable explicative):

equation   (57.293)

Remarque: Les logiciels de statistiques n'utilisent pas les méthode des moindres carrés pour déterminer les coefficients mais la méthode de la log-vraisemblance (cf. chapitre de Statistiques).

COURBEs ROC et lift

La "courbe ROC" (de l'anglais: "Receiver Operating Characteristic curve") est une mesure de la performance d'un classificateur binaire, c'est-à-dire d'un système qui a pour objectif de catégoriser des entités en deux groupes distincts sur la base d'une ou plusieurs de leurs caractéristiques. Graphiquement, on représente souvent la mesure ROC sous la forme d'une courbe qui donne le taux de vrais positifs (sensibilité : fraction des positifs qui sont détectés (correctement)) en fonction du taux de faux positifs (fraction des négatifs qui sont détectés (incorrectement)) pour ce même groupe. Les courbes ROC sont souvent utilisées en statistiques pour montrer les progrès réalisés grâce à un classificateur binaire lorsque le "seuil de discrimination" (cut-off) varie.

Remarque: Nous présentons ce sujet parce que de nombres logiciels de statistique renvoient une courbe ROC mais l'intérêt selon moi est objectivement très discutable. Nous montrerons cependant un outil plus utile juste après.

Voyons un exemple conrecte en reprenant notre exemple ci-dessus avec le tableur Microsoft Excel 14.0.7106. Dans un premier temps pour construire la courbe ROC il nous faut obtenir la matrice de confusion que nous avions présentée juste un peu plus haut avec Minitab. Pour obtenir celle-ci avec un tableur et sans code, voici une solution simple (mais ce n'est pas la solution la pluss condensée).

Le but va être d'abord de constuire le tableau suivant (ce n'est qu'une partie du tableau de l'exemple précédent puisqu'en réalité il y a 136 lignes de données) dont les colonnes B et C découlent des trois petits tableaux de l'exemple des crédits utilisé juste avant, la colonne A est juste la fréquence cumulée des individus 1/136... 2/136.... 3/136... et ainsi de suite):

equation
Figure: 57.12 - Liste à obtenir dans Microsoft Excel pour la courbe ROC

Voyons maintenant ce que contiennet les colonnes D, E et F qui sont directement liées au résultat obtenu plus haut et qui était pour rappel:

equation   (57.294)

mais que nous avons affiné avwec un logiciel de statistiques spécialisé afin d'obtenir:

equation   (57.295)

ce qui nous permet alors de constuire les trois colonnes précédemment mentionnées (ici nous n'avons que les premières 29 lignes car il suffit d'incrément les formules jusqu'en bas du tableau):

equation
Figure: 57.13 - Formules du tableur pour obtenir la probabilité du modèle et le score

La formule de la colonne F fait référence à la cellule M3 qui comme nous le verrons un peu plus loin contient la choix empirique de la valeur du cutoff que nous avions mentionnée lors de la présentation du modèle théorique de la régressin logistique (par défaut nous l'avons mis à 50%).

Ensuite, il nous faut construire les colonnes des vrais positifs, faux positifs, vrais négatifs, faux négatifs en utilisant des formules de tableurs élémentaires:

equation
Figure: 57.14 - Formules du tableur pour obtenir les vrais/faux positifs et négatifs

Une fois que nous avons ces données, nous allons pouvoir reconstuire le matrice de confusion qui nous avait été donnée par Minitab et qui nous sera indispensable pour élaborer la courbe ROC:

equation
Figure: 57.15 - Formules du tableur pour construire la matrice de confusion

ce qui donne explicitement:

equation
Figure: 57.16 - Valeurs correspondantes de la construction de la matrice de confusion

à comparer avec la matrice de confusion que nous avions donné et provenent de Minitab:

equation
Figure: 57.17 - Matrice de confusion de Minitab 16.1.1

Maintenant, observons la courbe ROC donnée par Minitab:

equation
Figure: 57.18 - Courbe ROC donnée par Minitab 16.1.1

Pourquoi l'abscisse représente-elle 1-% de faux positifs vous demanderez-vous alors qu'au niveau interprétation cela aurait été plus simple d'avoir simplement le taux de faux positfs? Eh bien pour deux raisons: la première c'est que les praticiens aiment bien les fonctions strictements croissantes... et la seconde raisons, la plus importante, c'est que si le classificateur binaire est parfaitement efficace, la surface est égale alors à 1 (c'est-à-dire à 100%). Ce qui il faut l'avouer est plus sympa que de dire qu'une surface nulle correspondant à une efficacité de 100%. Bon ceci étant dit... poursuivons.

Remarque: Certains logiciels (dont R) fournissent souvent que la surface sous la courbe. Ainsi plus celle-ci est proche de 1, meilleur est le classificateur!

Avant d'apprendre à interpréter ce graphique, comment obtenir cette même courbe dans un tableur? Eh bien il suffit simplement d'observer la colonne D de notre tableau Excel. Nous y avons logiquement 5 probabilités cumulées différentes (car il n'y avait que 5 tranches de crédits) qui sont respectivement:

0.235214057
0.481608302
0.777818743
0.813679766
0.991760974

Dès lors l'idée de la construction de la courbe ROC est de prendre chacune de ces probabilités cumulées comme valeurs de cutt-off. Ce qui donne respectivement:

equation
equation
equation
equation
equation
Figure: 57.19 - Valeurs de la matrice de confusion pour divers valeurs du cut off

Nous obtenons donc à chaque fois les coordonnées des points de la courbe ROC en fonction des valeurs de cut-off. Bon ceci étant fait, passons à l'interprétation en prenant le graphique ci-dessous complété:

equation
Figure: 57.20 - Courbes ROC (bleu: classificateur moyen, rouge: bon, vert: parfait, noir médiocre)

Alors comme le montre le graphique ci-dessous et l'intuition, un classificateur binaire parfait est celui dont le taux de vrais positifs est constant et toujours égal à 100%. Donc le classificateur binaire de notre exemple est moyennement bon.

Cependant de mon point de vue un bon classificateur binaire est celui qui pour une valeur donnée du cut-off maximise la somme du taux de vrais positifs et de vrais négatifs. Ainsi, avec un tableur comme celui de Microsoft Excel en mode "évolutionnaire" il est très simple de trouver cette valeur du cut-off qui maximise cet objectif. Mais cependant attention au piège: la valeur peut ne pas être unique, il peut s'agit d'un intervalle (ce qui est le cas dans notre exemple).

Enfin, si nous désignons par TP les vrais positifs (true positive), TN les vrais négatifs (true negative),FP les faux positifs (false positive) ou erreur de type I (dans une terminologie de théorie de la décision, ou de théorie des tests), et FN désigne les faux négatifs (false negative) ou erreur de type II.
On peut alors définir toute une batterie d'indicateurs permettant de juger de la qualité de notre prédicteur (ou plutôt de notre score):

  • TPR = TP / P = TP / (TP + FN) appelé "sensibilité", correspondant au taux de vrais positifs ("True Positive Rate" en anglais)
  • FPR = FP / N = FP / (FP + TN) correspondant au taux de faux positifs ("False Positive Rate" en anglais)
  • ACC = (TP + TN) / (P + N) appelé "précision" ("ACCuracy" en anglais)
  • SPC = TN / N = TN / (FP + TN) = 1 - FPR appelé "spécificité" ou taux de vrais négatifs
  • PPV = TP / (TP + FP) le taux de positifs prédits ("Positive Predictive Value" en anglais)
  • NPV = TN / (TN + FN) le taux de négatifs prédits ("Negative Predictive Value" en anglais)
  • FDR = FP / (FP + TP) correspondant au "False Discovery Rate" en anglais

Ceci étant fait, passons à la deuxième courbe.

Le principe de la "courbe lift" (que l'on pourrait traduire en français par "courbe levier") est très simple mais se base sur une hypothèse discutable (qui n'a heureusement cependant pas une grande importance) qui est que sans le modèle théorique (modèle logistique dans le cas présent) nous ne saurions absolument rien du comportement réel qu'auront les clients (débiteurs). L'hypothèse est de considérer que si nous prenions 50% des individus de notre échantillon, nous aurions 50% de vrais positifs (vrais mauvais débiteurs), si nous prenions 25% de notre échantillon, nous aurions 25% de vrais positifs (vrais mauvais débiteurs) et ainsi de suite. Cette hypothèse de départ (considérée par les praticiens comme le pire des cas) sera représentée graphiquement par la droite suivante:

equation
Figure: 57.21 - Courbe lift au pire (selon hypothèse)

Mais ça c'est juste pour avoir une référence empirique lorsque l'on travaille avec un seul modèle de classification statistique. En fait, face à la multiplicité des méthodes de modélisation statistiques, qui chacune possède ses propres indicateurs empiriques de qualité, les statisticiens ont recherché des critères généraux de performance d'un modèle (ceci dans l'idée évidemment de pouvoir en comparer plusieurs).

Ainsi, un choix empirique et assez intuitif consiste à dire qu'une méthode statistique est meilleure qu'un autre si pour un sous-ensemble donné de l'échantillon, son pouvoir prédictif (le nombre de vrais positifs cumulés par exemples) et meilleur ou non. C'est là le but de la courbe lift.

La courbe de lift est donc une variante de la courbe de ROC. La courbe lift classe les individus par score décroissant (à nouveau pour des raisons de facilitation d'interprétation de la surface sous la courbe et surtout pour avoir en premier le groupe cible de surveillance d'intérêt), en les regroupant par exemple par centiles, en déterminant le pourcentage d'événements d'intérêt dans chaque centile (normalement les vrais positifs), puis en traçant la courbe cumulative de ces pourcentages, de sorte qu'un point de ces coordonnées (n, m) sur la courbe signifie que les n% des individus ayant le plus fort score concentrent m% des événements. C'est la façon de construire cette courbe qui fait qu'on parle de performance prédictive du "ciblage marketing".

Voyons comment construire une telle courbe, toujours avec le même tableur et toujours les mêmes données. Le début de la construction consiste à reprendre certaines colonnes que nous avions utilisées pour la courbe ROC (31 premiers enregistrements sur les 136):

equation
Figure: 57.22 - Formules des premières colonnes à obtenir dans Microsoft Excel pour la courbe lift

où le lecteur aura peut-être remarqué que nous avons choisi empiriquement de fixer le cut-off à 50% (donc en réalité, et il ne faut jamais l'oublier, l'ensemble de la courbe lift change pour chaque valeur de cut-off!!!).

Ce qui donne explicitement (n'oubliez pas la colonne A qui comme pour la courbe ROC est triviale car représentant simplement l'effectif cumulé):

equation
Figure: 57.23 - Valeur explicites à obtenir dans Microsoft Excel pour la courbe lift

Une première différence maintenant par rapport à la courbe ROC, c'est que nous voulons les événements d'intérêt en premier (les mauvais débiteurs). Alors il faut trier la colonne score dans l'ordre décroissant. Ce qui donnera:

equation
Figure: 57.24 - Valeurs triées par ordre décroissant du score pour la courbe lift

Ensuite, nous ajoutons la colonne du % cumulé de vrais positifs:

equation
Figure: 57.25 - Ajout de la colonne de % cumulé de vrai positifs

Ce qui donne (une vue d'ensemble des quelques première lignes peut toujours être utile à la compréhension):

equation
Figure: 57.26 - Vue d'ensemble de toutes les colonnes construites

La dernière étape consiste maintenant simple à construire un graphique du nombre de % cumulé de vrai positifs (colonne G) par rapport au % cumulée de la taille de l'échantillon (colonne A), pour obtenir:

equation
Figure: 57.27 - Courbe lift (en bleu)

Alors évidemment nous ne pouvons pas dans le cas présent faire d'interprétation du pouvoir prédictif en terme de "lift" (ou "levier" en français) par rapport à un autre modèle statistique. Par contre nous pouvons comparer le pouvoir prédictif du modèle logistique utilisé dans le cas présent en termes de lift par rapport au cas considéré comme le pire (la diagonale pour rappel...). Donc dans le cas présent, si nous prenons le 20% des clients (débiteurs) ayant le score le plus grand, voilà ce que nous avons:

equation
Figure: 57.28 - Analyse lift

Donc dans cet exemple particulier, cela signifie que si nous concentrions notre analyse/étude/ciblage marketing ou autre que sur les 20% de l'échantillon étant le meilleur score (par exemple... pour des raisons de coût) , nous avons une surperformance de 37%/20%=1.85. Nous disons alors que le modèle a un lift de 1.85. L'idée ainsi lorsque nous comparons plusieurs modèles est de garder celui qui a le plus grand lift (levier).

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 un 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, Finance (courbe de Taux), ...

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

Voici 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 des traditions) joignant deux points equation est:

equation   (57.296)

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

Par définition, le segment [AB] est une "courbe de Bézier de degré 1" avec les 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étrique en rajoutant une 2ème étape à ce qui précède:

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

1ère étape:

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

2ème étape:

- Soit M(t) le barycentre de equation.

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

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

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] ):

equation
Figure: 57.30 - Arc de parabole [ABC]

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

equation   (57.297)

Nous avons:

equation   (57.298)

La récurrence se terminant pour:

equation   (57.299)

Ainsi, pour equation nous avons:

equation   (57.300)

Soit:

equation   (57.301)

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:

equation   (57.302)

Nous avons par la relation de récurrence:

equation   (57.303)

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

Nous avons donc:

equation   (57.304)

Il vient alors:

equation   (57.305)

Et donc:

equation   (57.306)

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

equation   (57.307)

et sous forme matricielle:

equation   (57.308)

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

equation   (57.309)

soit sous forme vectorielle:

equation   (57.310)

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

equation
Figure: 57.31 - 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:

equation   (57.311)

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

equation   (57.312)

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

equation   (57.313)

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 pour l'ordre n dans notre cas par:

equation   (57.314)

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

equation   (57.315)

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

equation   (57.316)

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

equation   (57.317)

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

equation   (57.318)

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:

equation   (57.319)

Remarque: Une courbe de Bézier est totalement modifiée dès que l'on déplace un point de contrôle. Nous disons alors que la méthode de Bézier est une "méthode globale".

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 equation 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):

equation
Figure: 57.32 - 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.

exempleExemple:

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:

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

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

equation

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

L'intuition nous amène à choisir  equation et equation 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

equation   (57.320)

Soit:

equation   (57.321)

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

equation   (57.322)

en equation. 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 qu'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 equation (h étant le "pas"). Nous notons equation la valeur exacte et equation, 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):

equation   (57.323)

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

equation   (57.324)

2. Analytiquement:

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

equation   (57.325)

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 equation une fonction connue sous forme explicite ou sous forme tabulée, et supposons qu'un certain nombre de valeurs:

equation   (57.326)

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

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

equation

equation
Figure: 57.34 - Illustration du concept d'interpolation

La fonction p s'appelle "fonction de collocation" sur les equation. 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.

Remarque: Il va sans dire que l'interpolation est un outil très important pour tous les chercheurs, statisticiens et autres.

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é:

equation   (57.327)

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

equation   (57.328)

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

equation   (57.329)

Système qui une fois résolu dans les règles de l'art (cf. chapitre d'Algèbre Linéaire) nous donne:

equation   (57.330)

Voyons le cas général:

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

Démonstration:

Posons:

equation   (57.331)

Les conditions de collocation:

equation   (57.332)

s'écrivent donc:

equation   (57.333)

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

Son déterminant s'écrit:

equation   (57.334)

relation que nous appelons "déterminant de Vandermonde". Nous savons que si le système a 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 equation, nous considérons le déterminant:

equation   (57.335)

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

equation   (57.336)

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

equation   (57.337)

Ce dernier polynôme peut s'écrire:

equation   (57.338)

Ce qui s'écrit:

equation   (57.339)

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

equationC.Q.F.D.

Il convient cependant de noter qu'il ne s'agit pas d'une méthode de régression polynomiale. Car avec une méthode de régression polynomiale, nous pourrions choisir un degré supérieur au nombre de points que nous avons!

MÉTHODE D'INTERPOLATION POLYNOMIALE DE LAGRANGE

La méthode d'interpolation polynomiale de Lagrange (très utilisée dans la pratique car l'algorithme est très simple et donc efficace) considère que nous avons initialement n + 1 points tels que:

equation   (57.340)

Et que nous cherchons donc un polynôme de collocation qui passe par tous les points. L'idée de l'interpolation polynomiale de Lagrange est alors simple et très astucieuse (comme toujours il fallait y penser...). Observez le graphique ci-dessous où nous avons 5 points par lesquels nous cherchons à faire passer un polynôme de colocation:

equation
Figure: 57.35 - Illustration du concept d'interpolation de Lagrange (source: Wikipédia)

Pour trouver le polynôme de colocation en rouge (supposez que nous ne le connaissons préalablement pas), l'idée est que pour chaque equation nous associons un polynôme de degré n et non nul en equation mais nul en tous les autres points mesurés equation et que tous les autres polynômes associés aux autres points equation y soient nul. Pour cela, comme nous le voyons dans le graphique ci-dessus, il faut donc que pour chaque point i nous ayons un polynôme associé qui ait n racines tel que:

equation   (57.341)

Donc nous voyons bien que nous avons 5 polynômes, chacun avec 4 racines (donc de degré 4) et qui sont respectivement nuls sur tous les points equation (vous pouvez le voir sur le graphique). Les coefficients equation sont des constantes à déterminer.

Maintenant, rien nous nous empêche de sommer puisque pour chaque equation nous aurons toujours la bonne valeur equation (puisque les autres polynômes y sont nuls). Soit en généralisant l'écriture, la somme devient:

equation   (57.342)

En injectant respectivement equation dans la relation précédente, il vient:

equation   (57.343)

Nous en déduisons alors immédiatement:

equation   (57.344)

En substituant les valeurs des constantes dans l'expression initiale de la somme, nous obtenons:

equation   (57.345)

Ce qui peut s'écrire sous forme condensée:

equation   (57.346)

Où le terme:

equation   (57.347)

est appelé "coefficients d'interpolations de polynômes de Lagrange".

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 equation ne possède pas de solution algébrique (sauf accident...) si equation est un polynôme de degré supérieur à 4.

Il existe un grand nombre d'algorithmes permettant de calculer les racines de l'équation equation 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 equation, nous traçons les courbes représentant g et h. Les racines de l'équation equation étant données par les abscisses des points d'intersection des deux courbes.

Remarque: Avant de résoudre numériquement l'équation equation, il faut vérifier que la fonction f choisie satisfasse certaines contraintes. Il faut par exemple, que la fonction soit strictement monotone au voisinage de la racine equation, lorsque la méthode de Newton est appliquée. Il est souvent utile, voire indispensable, de déterminer un intervalle [a,b] tel que:

- f est continue sur equation ou equation

- equation

- equation unique, equation

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 equation

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:

equation
Figure: 57.36 - 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):

equation   (57.348)

d'où:

equation   (57.349)

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

equation   (57.350)

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

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

2. Déterminer equation. Si equation est assez petit, nous arrêtons le calcul et affichons x1 et equation

3. Sinon nous procédons comme suit:

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

- nous remplaçons a par x1 et f(a)  par equation

- nous retournons au point (2)

MÉTHODE DE LA BISSECTION

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

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

2. equation

Il faut encore fixer equation 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 equation

2. Évaluation de  f(x)

3. Si equation 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 equation

- on remplace b par x si equation ou equation

- on retourne en (1)

L'étape (3) impose la condition equation 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 equation contenant equation. Ce test s'énonce comme suit:

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

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

La méthode de la sécante (ou "regula falsi" pour: régulièrement fausse) en méthodes numériques, est toujours un algorithme de recherche d'un zéro. Pour voir cela, considérons le schéma suivant:

equation
Figure: 57.37 - 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. equation

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

equation   (57.351)

nous en déduisons:

equation   (57.352)

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

1. Calcul de equation

2. Évaluation de equation

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

4. Sinon nous procédons comme suit:

- nous remplaçons a par equation si equation

- nous remplaçons b par equation si equation ou equation

- nous retournons en (1)

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

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

Remarque: Si l'intervalle [a,b] contient plusieurs racines, cette méthode converge vers l'une d'entre elles. Toutes les autres sont malheureusement perdues.

MÉTHODE DE NEWTON

Quelques années après sa découverte de théorie de la gravitation, Newton développa un technique extraordinaire qui permet de calculer les solutions d'une équation quelconque avec une rapidité (convergence phénomènale). Cette convergence surnaturellement rapide a ét utilisée pour démontrer certaines des résultats théoriques les plus marquants du 20ème siècle: le théorème de stabilité de Kolmogorov, le théorème de plongement isométrique de Nash... À elle seule, cette technique transcende la distinction entre mathématique pure et mathématique appliquée.

Pour étudier la méthode de Newton (appelée aussi "méthode de Newton-Raphson" ou encore "schéma d'approximation de Newton") dans le plan (donc à une variable explicative), considérons la figure suivante:

equation
Figure: 57.38 - Illustration de la méthode de Newton dans le plan

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

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 equation, alors nous pouvons la réécrire sous la forme (cf. chapitre de Calcul Différentiel Et Intégral):

equation   (57.353)

equation est la dérivée de f en equation et equation est une fonction qui tend vers 0 comme equation pour equation lorsque x tend vers equation (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 equation, nous obtenons:

equation   (57.354)

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

equation   (57.355)

et se résout aisément par rapport à equation. Pour voir cela, commençons par equation:

equation   (57.356)

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

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

equation   (57.357)

par:

equation   (57.358)

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 equation et equation 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)

Remarque: Il suffit souvent de vérifier les conditions (1) et (2) pour que le processus soit convergent.

La condition (2) est évidente, en effet si equation 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 a une boucle calculant alternativement equation et equation.

equation
Figure: 57.39 - Exemple de cas de non-convergence 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 equation par le quotient différentiel:

equation   (57.359)

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

equation   (57.360)

soit elle aussi suffisamment petite.

L'itération s'écrit alors:

equation   (57.361)

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

equation   (57.362)

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

equation   (57.363)

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 equation est réduit de moitié à chaque étape du calcul:

equation   (57.364)

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

equation   (57.365)

Mais:

equation   (57.366)

donc:

equation   (57.367)

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

equation   (57.368)

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

equation   (57.369)

ce qui montre bien que la convergence est quadratique.

Signalons enfin que nous étudierons la méthode de Newton à plusieurs variables dans le cadre de l'étude de la recherche (optimisation) non linéaire. C'est le choix pédagogique qui nous a semblé le plus judicieux.

Voici une application avec Maple 4.00b de cette méthode:

> with(plots): with(plottools):
Une fonction quelconque dont on cherche la racine
> f:=x->exp(x)*x^2-36;
> D(f)(x);
> Définition du point de départ
> x[0]:=3;
Définition du nombre d'étapes
> n:=7;
Équation de la tangente
> g:=x->f(x[i-1])+D(f)(x[i-1])*(x-x[i-1]);
> for i from 1 by 1 to n do;
> x[i]:=evalf(solve(g(x)=0,x));
> od;
> lines:={}:
> for i from 1 by 1 to n do;
> lines:=lines union {line([x[i-1],0],[x[i-1],f(x[i-1])],color=green), line([x[i-1],f(x[i-1])],[x[i],0],color=green)};
> od:
On peut jouer avec le x=.... pour mieux voir les itérations sur le graphique
> display({plot(f(x),x=2..3.01)} union lines);

equation
Figure: 57.40 - Application avec Maple 4.00b de la méthode de Newton

DÉRIVATIONS numériques

De nombreuses techniques de modélisation ou de résolution numériques que nous verrons plus loin utilisent les dérivées comme la recherche d'optimum (voir plus loin), les méthodes des éléments finis (voir plus loin). Par exemple, pour ne citer que le cas le plus connu, le solveur de Microsoft Office Excel des versions 2007 et antérieures propose quelques unes des dérivées numériques que nous allons étudier ici et réutiliser plus loin:

equation
Figure: 57.41 - Capture d'écran du solveur de Microsoft Excel 2003

Afin de permettre donc un traitement sur calculateur, les différentes dérivées présentes dansde nombreux algorithmes doivent être approchées numériquement. Pour ce faire, nous utilisons le principe des différences finies centrées qui s'appuie sur les développements en série de Taylor suivants (cf. chapitre Suites Et Séries):

equation   (57.370)

Nous avons alors sur la base de ce principe le développement au deuxième ordre:

equation   (57.371)

Il alors en négligeant les termes d'ordre supérieur et en soustrayant et simplifiant les deux séries ci-dessus:

equation   (57.372)

Relation que nous appelons "dérivée première centrée avec estimation tangente" (car nous négligeons tous les termes non linéaires). Nous retrouvons également souvent cette dernière relation sous la forme équivalente suivante:

equation   (57.373)

Maintenant, voyons ce que nous appelons les "dérivées première à droite" ou aussi "dérivées avant" ("forward derivate" en anglais) qui consistent simplement dans l'application de l'algorithme intuitif suivant:

equation   (57.374)

et accessoirement nous pouvons aussi définir les "dérivées première à gauche", appelées également "dérivées arrières" ("backward derivate" en anglais):

equation   (57.375)

Donc nous voyons que les dérivées centrales nécessitent plus de calculs mais sont aussi plus précises.

Nous pouvons également développer des relations plus élaborées en prenant des développements de Taylor à des ordres supérieurs, faire des moyennes entre différentes méthodes et donc c'est sans fin...

INTÉGRATIONS NUMÉRIQUES

Considérons la figure suivante:

equation
Figure: 57.42 - 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 equation et equation. Nous supposons dans ce cas que la fonction est à valeurs positives:

equation   (57.376)

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 Microsoft Excel ou OpenOffice Calc pour calculer des intégrales).

MÉTHODE DES RECTANGLES

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

equation
Figure: 57.43 - Approche de l'aire sous une courbe par des rectangles inférieurs gauche

L'aire de ces rectangles vaut:

equation   (57.377)

Si les equation sont suffisamment petits, equation est une bonne approximation de l'aire cherchée approchée par la gauche.

Nous pouvons recommencer cet exercice en choisissant equation et equation comme côtés des rectangles (donc l'approche par la droite). Nous obtenons alors:

equation   (57.378)

La figure correspondante est la suivante:

equation
Figure: 57.44 - Approche de l'aire sous une courbe par des rectangles supérieurs droite

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:

equation   (57.379)

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

equation   (57.380)

MÉTHODE DES TRAPÈZES

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

equation   (57.381)

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

equation   (57.382)

Ce que nous retrouvons souvent dans la littérature scolaire, sous la forme:

equation   (57.383)

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:

equation   (57.384)

et:

equation   (57.385)

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

PROGRAMMATION (OPTIMISATION) LINÉAIRE

L'objectif de la programmation linéaire (P.L.) est de trouver la valeur optimale d'une fonction linéaire soumise à un système d'équations d'inégalités constitué de contraintes elles aussi 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 du simplexe" (voir plus loin) dont la représentation graphique consiste en un "polygone des contraintes". 

Remarques:

R1. La programmation linéaire est beaucoup utilisée (pour ne citer que les cas les plus connus) dans la logistique (problème à flot maximal dit aussi "problème de transport"), la finance d'entreprise ou encore aussi en théorie de la décision lorsque nous devons résoudre un jeu à stratégie mixte (voir le chapitre de Théorie de la décision et des jeux pour un exemple pratique). C'est pour cette raison que Microsoft Excel 12.0 et antérieur intègre un outil appelé le "solveur" dans lequel il existe une option appelée "modèle supposé linéaire" qui alors impose l'utilisation du modèle du simplexe que nous allons voir ci-après:

equation

ou encore depuis la version 2010 du même logiciel (dont l'interface a complétement changé):

equation

R2. Dans le cadre de la résolution de problèmes où interviennent des produits de deux variables, nous parlons alors logiquement de "programmation quadratique" ou plus simplement de "programmation non linéaire". C'est typiquement le cas en économétrie dans la modélisation des portefeuilles (cf. chapitre d'Économie) ou dans le prévisionnel. Nous étudierons par ailleurs plus loin une version simplifiée et particulière des modèles correspondants qui sont: la méthode de Newton, de quasi-Newton, des gradients conjugués et la GRG non linéaire.

R3. Les programmations quadratique et linéaire sont réunies dans le cadre général de ce que nous appelons la "recherche opérationnelle".

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 de ceux 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:

equation   (57.386)

où les equation sont des variables qui influent sur la valeur de Z, et equation 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:

equation   (57.387)

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

equation   (57.388)

exempleExemple:

Une usine fabrique 2 types de 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 des ressources humaines (ouvriers) 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 (formulation canonique):

equation   (57.389)
equation

La fonction économique à maximiser étant:

equation   (57.390)

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":

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

Remarque: Dans le cas général, pour ceux qui aiment le vocabulaire des mathématiciens..., la donnée d'une contrainte linéaire correspond géométriquement à la donnée d'un demi-espace d'un espace à n dimensions (n étant le nombre de variables). Dans les cas élémentaires, l'ensemble des points de l'espace qui vérifient toutes les contraintes est un convexe limité par des portions d'hyperplan (voir le cas 2 variables, facile à illustrer), c'est pour cette raison que l'on parle alors "d'optimisation convexe". Si la fonction de coût est linéaire, l'extremum est à un sommet (facile à voir). L'algorithme du simplexe de base (voir plus loin) part d'un sommet et va au sommet d'à côté qui maximise localement le coût. Et recommence tant que c'est possible.

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 equation maximisant la fonction économique.

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

En particulier, la droite equation passe par l'origine et celle-ci fournit 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 planequation) la droite de pente -1.2. Évidemment on remarque alors très vite que la méthode du simplexe ne fonctionnera plus si le polygone des contraintes ne contient pas l'origine!

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

equation
Figure: 57.46 - Recherche 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 d'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":

equation   (57.391)
equation

avec:

equation   (57.392)

Nous introduisons d'abord les "variables d'écart" afin de transformer les 2 inégalités en des égalités. Le système d'équations prend alors une "forme standard":

equation   (57.393)

Donc pour equation fixés, les variables d'écarts dont les coefficients sont toujours unitaires, mesurent la distance à parcourir jusqu'aux sommets.

Il va sans dire que la technique des variables d'écart peut être utilisée pour les systèmes linéaires (ou non linéaires). Dès lors, un système d'optimisation contrainte avec inégalités, peut toujours être ramené à un système d'optimisation avec égalités.

Remarque: Il y a autant de variables d'écart que d'inéquations!

Pour la suite, nous avons remarqué, suite à une relecture de ce chapitre, que la technique des tableaux souvent présentée dans les livres et sur les sites Internet n'apportait finalement absolument rien à la compréhension profonde du mécanisme de résolution (même si pour programmer informatiquement la méthode cela est plus pratique). Comme le but de ce site est de démontrer toujours avec le maximum de détails le principe de fonctionnement des choses il va donc de soi que nous allons opter pour une approche purement algébrique. Voyons donc cela en revenant au système avec les variables d'écart et la fonction économique mais un peu réarrangée:

equation   (57.394)

La contrainte A1 devient alors:

equation   (57.395)

et la contrainte A2 devient respectivement:

equation   (57.396)

Dès lors, le problème se résume à maximiser Z avec les contraintes:

equation   (57.397)

Partons donc d'une solution réalisable évidente qui au vu des contraintes est trivialement:

equation   (57.398)

Dès lors, avec le système:

equation   (57.399)

nous trouvons immédiatement:

equation   (57.400)

Les paramètres dans l'état actuel peuvent se résumer à:

equation   (57.401)

Pour progresser, le but sera de faire croître Z et pour cela nous allons ne faire croître qu'une variable, en choisissant celle qui a le plus grand coefficient (poids) dans:

equation   (57.402)

c'est-à-dire equation (nous pensons que c'est ainsi que Z croîtra le plus vite). Nous parlons alors de equation comme "direction du pivot". Nous gardons donc equation et nous faisons croître equation avec le système qui se réduit alors à:

equation   (57.403)

Avec donc equation et pour commencerequation, nous avons:

equation   (57.404)

et nous voyons que les contraintes equation sont toujours respectées, il en va de même si equation vaut 2, 3, 4, 5, ... et ce jusqu'à 31, car dès lors:

equation   (57.405)

et l'une des variables d'écart étant devenue négative, les contraintes equation ne sont plus toutes satisfaites et donc cette solution n'est pas réalisable.

La question dans le cas général consiste à se demander jusqu'à combien (la valeur la plus contraignante, in extenso la plus petite) nous pouvons faire croître equation tout en maintenant la condition equationlorsque equation? Et la réponse est assez simple:

equation   (57.406)

et donc il s'agit de:

equation   (57.407)

Nous parlons alors parfois du "pas du pivot". Nous avons alors la solution actuelle:

equation   (57.408)

Ce qui donne:

equation   (57.409)

Graphiquement parlant, voilà à quoi correspond ce que nous venons de faire:

equation
Figure: 57.47 - Direction du pivot avec arrivée au point (30, 0)

Pour continuer à faire croître aussi simplement Z (en ne faisant croître qu'une variable), il nous faut
un nouveau système d'équations similaire au système initial:

equation   (57.410)

où nous avions exprimé les variables qui prennent une valeur non nulle en fonction des autres qui prennent une valeur nulle, c'est-à-dire equation en fonction de equation puisque nous avions pour rappel:

equation   (57.411)

Pour la suite, il nous faut exprimer equation ainsi que Z en fonction de equation puisque nous venons d'obtenir pour rappel:

equation   (57.412)

Avant d'obtenir le nouveau système fonction de equation, faisons quelques manipulations algébriques:

equation   (57.413)

ce qui nous donne après simplification:

equation   (57.414)

et dès lors, il vient:

equation   (57.415)

ce qui nous donne après simplification:

equation   (57.416)

et nous avons de même:

equation   (57.417)

Donc au final, le système est:

equation   (57.418)

à partir duquel nous itérons le processus (nous ne faisons croître qu'une variable dans Z en gardant les autres à 0). Quand nous ne pouvons plus faire augmenter Z car tous les coefficients sont négatifs, et bien c'est que nous sommes au maximum (merci la convexité). Voyons cela...

Dans Z le plus gros coefficient est maintenant equation et donc cela nous amène à poser equation. La valeur la plus contraignante de equation qui permet de respecter toujours les contraintes equation est dès lors:

equation   (57.419)

Et pour cette valeur, nous avons:

equation   (57.420)

La fonction économique initiale prend alors la valeur:

equation   (57.421)

ce qui correspond donc graphiquement à:

equation
Figure: 57.48 - Direction du pivot avec arrivée au point (16, 28) pour la deuxième itération

Donc nous voyons bien ci-dessus que nous sommes arrivés à la valeur optimale visible sur le graphique donné au début de l'exercice. Mais comment savoir que nous sommes arrivés au point final si nous n'avons pas de graphique ou que nous travaillons dans des dimensions supérieures?

Au fait, le processus est terminé soit quand tous les coefficients de la fonction économique sont négatifs ou que la valeur la plus contraignante qui respecte les contraintes est nulle!!! Voyons si c'est bien le cas! Nous avons donc dans le cas présent:

equation   (57.422)

Et nous allons donc réécrire le système:

equation   (57.423)

avec cette fois-ci equation ainsi que Z en fonction de equation. Nous avons alors d'abord:

equation   (57.424)

et nous avons:

equation   (57.425)

Et donc pour Z il vient (les coefficients y sont tous négatifs donc nous devinons la suite...):

equation   (57.426)

Nous avons donc le nouveau système:

equation   (57.427)

Comme tous les coefficient de Z sont négatifs, nous sommes bloqués car nous partirions dans la mauvaise direction. Il faut donc nous arrêter ici et nous adoptons au final la solution:

equation   (57.428)

La méthode des tableaux que l'on trouve souvent dans la littérature consiste à ne noter que les coefficients des variables du système dans un tableau, mais les transformations que l'on fait sont exactement celles que l'on vient de faire algébriquement (à part qu'elle masque le sens de la méthode).

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ée 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:

equation   (57.429)

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

equation   (57.430)

où les matrices equation 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 comment 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:

equation   (57.431)

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

equation   (57.432)

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

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

equation   (57.433)

Nous écrirons:

equation   (57.434)

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

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 comment 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écomposer en deux sous-matrices equation, une contenant les variables initiales D et l'autre comportant les variables d'écart B tel que:

equation   (57.435)

Remarque: Les variables d'écart sont des variables et non des constantes!! Il convient dans un système où les variables sont au nombre de n et les équations au nombre de m tel qu'une équation du système s'écrirait:

equation   (57.436)

d'ajouter une variable d'écart tel que:

equation   (57.437)

equation et sur chaque ligne m, la variable d'écart ajoutée devant être différente de celles déjà insérées dans le système. C'est la raison pour laquelle nous pouvons décomposer la matrice en deux-sous matrices.

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 equation. Les variables associées aux colonnes de la matrice D seront appelées les "variables hors-base"; il s'agit des variables equation.

Remarque: Rappelons que dans l'expression de la fonction économique, seules les variables hors-base apparaissent.

En résumé, toute P.L. une fois mise sous forme standard est telle 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 equation (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 equation".

Remarque: Nous pouvons intervertir les matrices (et donc changer de base canonique) B et D (ce qui revient à dire que la matrice des variables de base devient la matrice des variables hors-base et inversement).

Il est maintenant commode d'introduire les notations suivantes:

equation   (57.438)

qui sont respectivement le vecteur 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:

equation   (57.439)

ou bien aussi:

equation   (57.440)

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

equation   (57.441)

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

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 equation et equation, nous obtenons:

equation   (57.442)

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

equation   (57.443)

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

equation   (57.444)

Nous regroupons les termes en equation et nous avons:

equation   (57.445)

Nous avons alors toujours un 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 de 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 equation de la ligne equation.

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

equation   (57.446)

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

equation   (57.447)

et le système se réduit à:

equation   (57.448)

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

equation   (57.449)

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

equation   (57.450)

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":

equation   (57.451)

À 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 equation 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 equation), 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 equation), la solution de P.L. se trouve nécessairement en un programme de base.

PROGRAMMATION (OPTIMISATION) NON LINÉAIRE

Un programme d'optimisation non linéaire (N.L.O.) est une généralisation de la programmation linéaire (algorithme du simplexe) mais à des fonctions non linéaires et pouvant comporter aussi des contraintes et des fonctions économiques non linéaires.

Le but de ce qui va suivre est donc de comprendre dans les grandes lignes mais avec un niveau de rigueur acceptable les outils d'optimisation que proposent de nombreux tableurs comme les versions antérieurs à Microsoft Excel 2007:


equation
Figure: 57.49 - Exemple d'optimisation non linéaire du solveur de Microsoft Excel

Nous allons en particulier voir maintenant en quoi consiste la recherche Newton (sous-entendu: Gauss-Newton) avec les estimations Tangente et Quadratique. Après quoi nous étudierons la méthode de recherche des Gradients conjugués aussi avec les méthodes tangente et respectivement quadratique.

Remarque: Nous nous arrêterons à l'étude des modèles précités car il existe une quantité trop important de modèles empiriques pour citer par exemple que les plus connus des modèles (algorithmes): méthode de substitution, méthode des multiplicateurs de Lagrange, algorithme de Nelder-Mead, algorithme de Broyden–Fletcher–Goldfarb–Shanno (BFGS), algorithme d'annhiliation simultnaée (SANN ou SA), Méthodes de points intérieurs, et voir Wikipédia pour une liste plus compléte (il y en a plus d'une dizaine sans compter les variantes incluant des ajustements empiriques).

Nous le verrons plus loin, mais nous le devinons déjà, que le choix Tangente utilise une approximation linéaire de la tangente à la fonction à optimiser au point considéré alors que l'option Quadratique fera un estimation d'une fonction du deuxième degré au point considéré (typiquement une parabole). Si au point considéré, la fonction se modélise bien par une quadrique, alors l'option Quadratique peut faire économiser du temps en choisissant un meilleur point initial qui demandera moins de pas supplémentaires à chaque recherche. Si vous n'avez pas d'idée du comportement a priori de la fonction, la méthode Tangente est alors plus lente mais plus sure.

Un exemple connu dans la littérature pour introduire la recherche d'optimum de fonctions non linéaires, avant de passer à la partie intégrant les contraintes du système, est la fonction "baleine à bosse" qui consiste à trouver le minimum de:

equation   (57.452)

avec les contraintes:

equation   (57.453)

ce que nous pouvons en effet vérifier visuellement.

equation
Figure: 57.50 - Tracé de la fonction baleine à bosse avec minimum déjà visibles

Soit avec Maple 4.00b:

>plot3d(x^2*(4-2.1*x^2+1/3*x^4)+x*y+y^2*(-4+4*y^2),
x=-2..2,y=-1..1,contours=20,style=patchcontour,axes=boxed);

Comme nous pouvons le voir, cette fonction est également un excellent exemple de minimaux locaux multiples.

MÉTHODE DE NEWTON-RAPHSON (NEWTON QUADRATIQUE)

La "méthode de Newton-Raphson" est une technique permettant de chercher l'extremum d'une fonction ou également, comme nous le verrons plus loin lorsque nous comparerons dans un cas particulier d'application la méthode de Gauss-Newton à celle de Newton-Raphson, pour la régression non linéaire.

La méthode de Newton-Raphson, qui dans les versions antérieures à Microsoft Excel 2007 s'activait dans le solveur en sélectionnant l'option Newton et Quadratique, utilise les approximations de Taylor au deuxième ordre (donc avec les dérivées du deuxième ordre) pour avoir une fonction quadratique (parabole) qui converge si le point d'origine de la recherche est proche de l'optimum. Cette approximation est réitérée à chaque itération.

Pour commencer, rappelons que nous avons démontré dans le chapitre de Suites Et Séries qu'un développement de Taylor pour une fonction à deux variables pouvait s'écrire en approximation quadratique:

equation   (57.454)

où pour rappel h et k sont des variables et equation sont fixés et où nous avons la matrice hessienne:

equation   (57.455)

que les américains spécialistes du domaine ont pour habitude (malheureuse à mon avis...) de noter:

equation   (57.456)

la dernière expression étant la plus courante peut être très trompeuse avec la notation du laplacien.

Dans le domaine des méthodes numériques il est d'usage d'écrire la série de Taylor ci-dessus avec quelques changements de notation en posant d'abord:

equation   (57.457)

Ce qui nous donne une forme plus condensée et technique de la série de Taylor autour de equation :

equation   (57.458)

En changeant un peu la notation:

equation   (57.459)

Nous retrouvons donc l'expression d'usage d'une fonction de equation évaluée en série de Taylor centrée en equation.

Mais si nous cherchons un extrema local (appelé aussi parfois "point critique"), il faudra dans un premier temps que la dérivée de l'ensemble de la série de Taylor soit nulle. C'est-à-dire:

equation   (57.460)

et que le déterminant de la matrice hessienne soit positif (cf. chapitre Suites Et Séries). Et pour savoir si nous sommes sur un maximum local ou minimul local, il nous faudra regarder le signe de equation.

Récrivons la relation ci-dessus explicitement comme nous l'avions démontrée dans le chapitre de Suites Et Séries pour des raisons pédagogiques:

equation   (57.461)

Et rappelons que tous les termes equation sont des constantes car il s'agit soit de la fonction évaluée sur le point particulier equation, soit de la dérivée partielle évaluée en ce même point, soit de la dérivée partielle seconde toujours évaluée en ce même point, etc.

Donc le gradient donnera finalement:

equation   (57.462)

et en revenant aux notations d'usage du domaine des méthodes numériques, nous avons alors:

equation   (57.463)

Et donc le gradient devant être nul, nous avons:

equation   (57.464)

et après un premier réarrangement:

equation   (57.465)

et un deuxième réarrangement:

equation   (57.466)

ce qui se note souvent:

equation   (57.467)

et par les américains...:

equation   (57.468)

Enfin avant de passer à un exemple concret il est important que le lecteur se souvienne de la relation vu juste plus haut:

equation   (57.469)

Le fait que la gradient négatif apparaisse fait que la technique de Newton-Raphson (Newton Quadratique) appartient à la famille des techniques dites de "descente rapide".

Cette dernière égalité étant souvent noté chez les américains (... sans commentaires...) par:

equation   (57.470)

Évidemment un tableur comme Microsoft Excel ne pouvant pas déterminer les dérivées il va les calculer en utilisant la méthode numérique des dérivées à droite ou centrées comme nous l'avons présenté un peu plus haut.

exempleExemple:

Nous cherchons un extrema local de la fonction "baleine à bosse" représentée plus haut:

equation   (57.471)

avec le point de départ (arbitraire):

equation   (57.472)

Pour effectuer la recherche, nous calculons d'abord le gradient:

equation   (57.473)

et la matrice hessienne:

equation   (57.474)

Nous avons alors:

equation   (57.475)

et donc:

equation   (57.476)

et nous recommençons (bon on va se passer de tous les détails maintenant car faut pas exagérer non plus...):

equation   (57.477)

et donc:

equation   (57.478)

et nous recommençons (avec encore moins de détails):

equation   (57.479)

et donc:

equation   (57.480)

et encore une fois (avec encore moins de détails):

equation   (57.481)

et donc:

equation   (57.482)

et encore une fois (avec encore moins de détails):

equation   (57.483)

et les valeurs ne bougeront plus. Mais si nous regardons le graphique d'origine où nous avons mis en évidence le point de convergence par un point rouge:

equation
Figure: 57.51 - Mise en évidence du point de convergence dans la fonction baleine à bosse

nous constatons que ce système ne cherche pas un extremum global mais local comme nous l'avions déjà spécifié. En réalité, comme le lecteur pourra le tester lui-même, la convergence est très sensible au point de départ initial.

MÉTHODE DE GAUSS-NEWTON (NEWTON TANGENTE)

La méthode de Gauss-Newton est une approximation puissante sans les dérivées du deuxième ordre de la méthode de Newton-Raphson qui dans les versions antérieures à Microsoft Excel 2007 s'activait dans le solveur en sélectionnant l'option Newton et Tangente,

Pour aborder ce sujet, partons tout de suite accompagnés d'un exemple concret. Supposons que nous avons obtenu les données suivantes:

equation

equation

1

3.2939

2
4.2699
4
7.1749
5
9.3008
8
20.259
     Tableau: 57.9 - Données mesurées

et nous supposons que les données suivent le modèle théorique suivant (nous aurions pu faire n'importe quel autre choix!):

equation   (57.484)

Nous cherchons donc equation qui minimisent la somme des carrés des écarts entre les valeurs expérimentales et théoriques tels que:

equation   (57.485)

avec donc:

equation   (57.486)

Notons pour la suite comme le veut la tradition dans le domaine:

equation   (57.487)

Nous avons alors la notation courante:

equation   (57.488)

Maintenant, imaginons que nous ayons trouvé une valeur du couple equation qui donne ce minimum et notons le equation et en n'oubliant pas évidemment que cela ne sera qu'une solution locale! Considérons un cas particulier que nous appelons "solution compatible" et définie par le fait que le couple qui minimise la somme des carrées des erreurs est aussi tel que pour tout i nous ayons:

equation   (57.489)

Dès lors, il en découle immédiatement que:

equation   (57.490)

Avant d'aller plus loin, remarquons que par exemple que pour une composante j (ce qui correspond dans notre cas à chaque variable de la fonction théorique supposée a priori):

equation   (57.491)

où la dernière égalité condensée est souvent loi d'être triviale d'autant plus qu'elle fait usage du gradient d'un champ de vecteur (cf. chapitre de Calcul Vectoriel) que l'on retrouve rarement dans la pratique. Le lecteur déstabilisé pourra se reporter si besoin directement à l'exemple plus bas afin d'éclairer sa lanterne.

Nous en déduisons donc que:

equation   (57.492)

et la solution compatible nous amène donc à bien évidemment à:

equation   (57.493)

De la même manière, nous avons:

equation   (57.494)

Nous avons donc au final les deux relations suivantes:

equation   (57.495)

Étant donné que pour la solution compatible nous avons:

equation   (57.496)

Il s'ensuit que dans ce cas que la deuxème relation devient:

equation   (57.497)

H est la matrice Hessienne (cf. chapitre de Suites Et Séries) ce que les américains notent simplement:

equation   (57.498)

Donc nous pouvons approximer dans le cas de la solution compatible, la Hessienne qui contient des dérivées d'ordre deux par des dérivées premières.

Donc nous avons finalement dans ce cas particulier les deux relations qui sont à la base de méthode de Gauss-Newton:

equation   (57.499)

Maintenant, rappelons la relation de base de la méthode de Newton-Raphson obtenu plus haut:

equation   (57.500)

et pour information, toute technique mathématique (car elles sont nombreuses!), permettant de simplifier la matrice Hessien à droite de l'égalité fait alors partie de la famille des "méthodes de quasi-Newton".

Eh bien la méthode de Gauss-Newton qui nous intéresse ici et qui est donc une des techniques de la famille des méthodes de quasi-Newton consiste simplement dans un premier temps à se débarrasser des dérivées secondes de la Hessienne de la méthode de Newton-Raphson à droite de l'égalité à l'aide de la relation établie précédemment tel que (attention à se rappeler de l'abus d'écriture!):

equation   (57.501)

et dans un deuxième temps récrire le gradient à gauche de l'égalité à l'aide aussi de la relation précédemment établie. Ce qui nous donne:

equation   (57.502)

Le facteur 2 n'étant pas très esthétique, la quasi-totalité des ouvrages de référence optimisent le problème avec la relation de départ suivante:

equation   (57.503)

Donc en multipliant simplement par un facteur ½ (ce qui ne change rien au résultat). Nous avons alors:

equation   (57.504)

Rappelons encore une fois qu'un tableur comme Microsoft Excel ne pouvant pas déterminer les dérivées il va les calculer en utilisant la méthode numérique des dérivées à droite ou centrées comme nous l'avons présenté un peu plus haut.

Revenons maintenant à notre exemple du début! Nous avons donc:

equation   (57.505)

Nous commençons avec un couple qui nous semble proche de la solution cherchée:

equation   (57.506)

Nous avons alors:

equation   (57.507)

et donc nous avons:

equation   (57.508)

Pour la suite, il vient alors:

equation   (57.509)

Ce que nous pouvons donc récrire sous la forme:

equation   (57.510)

Nous avons aussi in extenso:

equation   (57.511)

Ensuite, nous appliquons la relation démontrée plus haut:

equation   (57.512)

Soit:

equation   (57.513)

et après une petite simplification mineure:

equation   (57.514)

Soit:

equation   (57.515)

et donc le prochain couple pour l'itération sera:

equation   (57.516)

Ce qui correspond bien évidemment aux valeurs pour la première itération:

equation   (57.517)

Nous n'allons pas refaire aussi explicitement les autres itérations. Donc voilà ce que cela donne au final:

i

equation

0

equation

1
equation
2
equation
3
equation
4
equation
     Tableau: 57.10 - Itérations Gauss-Newton

avec donc pour solution locale à la 4ème itération:

equation   (57.518)

Faisons pour clore une comparaison avec la méthode de Newton-Raphson pour la première itération en utilisant le même couple de départ. Rappelons encore une fois que pour cette méthode, les itérations sont basées sur la relation:

equation   (57.519)

et nous écrirons la fonction de la façon suivant pour la méthode de Newton-Raphson:

equation   (57.520)

où nous avons donc pour la première itération (nous retrouvons le même vecteur que pour la méthode de Gauss-Newton):

equation   (57.521)

et:

equation   (57.522)

Nous avons donc:

equation   (57.523)

qui devient:

equation   (57.524)

Soit:

equation   (57.525)

et donc le prochain couple pour l'itération sera:

equation   (57.526)

Ce qui correspond bien évidemment aux valeurs pour la première itération:

equation   (57.527)

Nous n'allons pas refaire aussi explicitement les autres itérations. Donc voilà ce que cela donne au final au niveau de la fonction à minimiser:

i

equation

0

equation

1
equation
2
equation
3
equation
4
equation
5
equation
     Tableau: 57.11 - Itérations Newton-Raphson

donc la méthode de Newton-Raphson converge dans ce cas particulier moins vite que celle de Gauss-Newton.

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 (ou tout autre problème de la même famille)

- 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)

- Création de tests statistiques (Anderson-Darling, Kolmogorov, Levene, Brown-Forsythe, etc.)

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ésoluble 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édiate 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 pour lequel des normes internationales ont été édictées (ISO 28640:2010).

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

>rand();

equation

>restart;rand();

 equation

Nous voyons donc que la fonction par défaut de générateur de nombres aléatoires de Maple 4.00b 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 4.00b telles que:

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

equation

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

equation

É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 suit pas une loi de distribution connue... ce qui n'est malheureusement jamais le cas).

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

equation
Figure: 57.52 - Illustration d'une séquence de nombres pseudo-aléatoires avec Microsoft Excel 14.0.6123

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 sûr l'intervalle [0,1]. Nous parlons alors de "nombres quasi-aléatoires" permettant de faire des simulations appelées parfois "quasi Monte-Carlo". Avec la version française de Microsoft Excel 11.8346, 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:

equation
Figure: 57.53 - Illustration d'une séquence de nombres quasi-aléatoires avec Microsoft Excel 11.8346

où nous voyons bien que la séquence couvre bien la surface comprise entre 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 Microsoft Excel 11.8346 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 égal à 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.

Remarque: Se référer à la norme internationale ISO 28640:2010 pour les ingénieurs ayant besoin d'implémenter des générateurs de nombres aléatoires dans leurs logiciels.

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]:

equation   (57.528)

Soit:

equation   (57.529)

la valeur maximale de la fonction entre les bornes [a,b].

Nous considérons le rectangle englobant la fonction sur [a,b] défini par les points equation:

equation
Figure: 57.54 - Principe de base du calcul de l'intégrale avec Monte-Carlo

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-dessous, nous avons:

equation   (57.530)

L'algorithme Maple 4.00b 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:

Remarque: Pour appeler cette procédure, il suffit d'écrire >intmonte(f,a,b,N) mais en remplaçant le premier argument passé en paramètre par l'expression d'une fonction et les autres arguments par des valeurs numériques bien évidemment.

CALCUL DE PI

Pour le calcul de equation le principe est le même et consiste 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:

equation   (57.531)

L'algorithme Maple 4.00b 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:
>evalf(calculepi(100));evalf(calculepi(1000));evalf(calculepi(10000));evalf(calculepi(100000));

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 Microsoft Excel (même les multinationales!) et dans une moindre mesure avec des logiciels comme @Risk, CrystalBall, TreeAge ou encore Isograph.

Les avantages de cette 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 disponibles 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 bêta (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 de 8 jours. La tâche B une durée optimiste de 1 jour et pessimiste de 4 jours. Nous souhaiterions dans Microsoft Excel à l'aide d'une simulation de pseudo Monte-Carlo (basée donc obligatoirement sur une variable pseudo-aléatoire) introduire 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...):

equation
Figure: 57.55 - Mise en place d'une petite simulation de Monte-Carlo avec Microsoft Excel 11.8346

où toutes les cellules de la colonne A contiennent la fonction suivante (version français de Microsoft Excel 14.0.6123):

=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 (version français de Microsoft Excel 14.0.6123):

=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 dans Microsoft Excel 14.0.6123 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):

equation
Figure: 57.56 - 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):

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

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

Remarque: Dans le cas de simulations de variables aléatoires, on peut dans les cas simples impliquant uniquement des sommes ou soustractions de variables aléatoires, comme c'est le cas pour l'exemple ci-dessus, déterminer l'espérance et l'écart-type du résultat analytiquement en utilisant la propriété de linéarité de l'espérance et de la variance (car normalement pour la variance de deux variables aléatoires indépendantes, la covariance est nulle). En analysant la différence entre la valeur analytique et celle obtenue par la simulation numérique, on peut corriger le décalage de certains autres indicateurs statistiques par simple ajout ou soustraction du différentiel. On parle alors de la technique des "variables de contrôle".

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écorréler 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 ne 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 variable aléatoire 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 unique variable. Raison pour laquelle elle 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'échantillonnage 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: fait 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é 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'échantillonnage stratifié ou de Latin Hypercube même si ces techniques sont toutes faciles à programmer, la méthode 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 distinguons en général, deux types de bootstrap:

- Les bootstraps qui ne font aucune hypothèse sur la loi de distribution des données analysées. Nous parlons alors de "bootstrap non-paramétrique". C'est le cas le plus courant et nous ferons un exemple uniquement pour celui-ci.

- Les bootstraps qui remplacent chacune des données mesurées par celles correspondantes à l'expression analytique de la loi distribution de probabilité supposée. Nous parlons alors de "bootstrap paramétrique". Une fois le remplacement de chacune de valeurs d'origine effectué, la démarche est exactement celle du bootstrap non-paramétrique.

Nous allons illustrer le principe du bootstrap (dit aussi "bootstrapping") sur l'exemple de l'intervalle de confiance de l'espérance equation 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 equation estimations:

equation   (57.532)

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

equation   (57.533)

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

equation   (57.534)

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:

equation   (57.535)

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:

equation   (57.536)

Soit:

equation   (57.537)

Ce qui donne:

equation   (57.538)

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 bootstrap", 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:

equation   (57.539)

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):

equation   (57.540)

Pour chaque échantillon simulé, une moyenne est calculée (donc nous aurons 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). Ce regroupement se nomme le "bagging" pour "bootstrap aggregating".

É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 Microsoft 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ééchantillonnage" (resampling) comme en fait partie le bootstrapping.

exempleExemple:

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

equation
Figure: 57.58 - É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:

equation
Figure: 57.59 - Médianes de répliques

avec la longue formule pour la version française de Microsoft Excel 14.0.6123 suivante qu'il faut mettre dans la cellule 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 Microsoft Excel 14.0.6123 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):

equation
Figure: 57.60 - Convergence de la médiane en fonction du nombre de réplications

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

Enfin indiquons qu'après avoir étudié des nombreux modèles de régression linéaire il n'empêche pas que dans de nombreux cas aucun modèle théorique ne soit adapté que ce soit pour interpoler ou pour in extenso extrapoler. Dès lors, si nous avons pour chaque valeur des abscisses un vingtaine de reproductions de la variable à expliquer (ce qui éléminie presque directement les applications dans le domaine de l'économie), nous pouvons faire du bootstrapping ce qui nous donnera des coefficients de régressions bootstrappés et aussi des valeurs interpolées ou extrapolées bootstrappées! C'est une technique extrêmement intéressante de régression non paramétrique dans la pratique! Bien que nous puissions faire cela dans des logiciels comme IBM SPSS (avec un module complémentaire) ou SAS avec un simple tableau en utilisant la technique ci-dessus mais en la généralisant nous pouvons faire des inférences très intéressantes.

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 lequel la fonction f s'annule.

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

equation   (57.541)

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

equation   (57.542)

Nous  calculons:

equation   (57.543)

Si equation alors la racine est dans l'intervalle equation sinon elle est dans l'intervalle equation. 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 4.00b 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!!


Haut de page

CHIMIE THERMIQUEMETHODES NUMERIQUES (2/2)


Like6   Dislike0
61.03% sur 100%
Notée par 39 visiteur(s)
12345
Commentaires: [0] 
 
 


W3C - HTMLW3C - CSS Firefox
Ce travail est dans le domaine public
2002-2024 Sciences.ch

Haut de page