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%
vues
depuis le 2012-01-01:
54'250
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).
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.
(les quatre opérations fondamentales ont le même coût)
H2. Un accès mémoire
une opération arithmétique
H3. Un contrôle de
comparaison
une opération arithmétique
H4. Un seul processeur
Nous
notons l'ensemble
des données de taille n et
T(n) le
coût (en temps) de l'algorithme sur la donnée ou le
jeu de données de taille n.
-
La "complexité au meilleur" est donnée par la
fonction:
(57.1)
C'est
le plus petit temps qu'aura
à exécuter l'algorithme sur un jeu de données
(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):
(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é 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:
(57.3)
se résumera à avoir une complexité d'ordre .
-
La "complexité moyenne":
(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
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 .
Une mauvaise complexité est du type
ou .
Exemples:
E1. Évaluation d'un polynôme:
(57.5)
L'évaluation directe
de la valeur de P(x)
conduit à une complexité
(57.6)
Nous devons à Horner
un algorithme plus performant qui utilise une factorisation du polynôme
sous la forme:
(57.7)
Nous pouvons facilement montrer
que cette factorisation maintient le même nombre d'additions
mais réduit le nombre de multiplications à (n).
La complexité qui en découle est donc O(n).
Le gain est incontestablement important. De plus, cette factorisation
évite tout calcul de puissance.
E2. Un autre exemple connu de complexité algorithmique
est la recherche d'une information dans une colonne triée.
Un algorithme simple appelé "recherche dichotomique" consiste à prendre
la cellule 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, à l'étape
2, et plus généralement de à l'étape k.
Au pire, la recherche se termine quand il n'y a plus qu'une seule
cellule de la colonne à explorer, autrement dit quand k est
tel que .
Nous en déduisons le nombre maximal d'étapes: c'est
le plus petit k tel
que ,
soit ,
soit:
(57.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 .
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:
- Addition:
telle que
- Multiplication:
telle que
- Déterminant (par
la méthode directe de Cramer). Nous renvoyons le lecteur
au chapitre d'Algèbre Linéaire pour le détail
des méthodes de calcul du déterminant d'une matrice.
Nous pouvons alors montrer que la complexité d'un déterminant
d'ordre n est n multiplications, n-1
additions plus n
fois la complexité d'un déterminant d'ordre n-1.
Par cumul, nous arrivons à:
(57.9)
En faisant l'hypothèse
que l'ordinateur utilisé effectue une opération élémentaire
en
secondes (ce qui est déjà une bonne machine), nous
obtenons alors les temps de calculs suivants pour plusieurs valeurs
de n:
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 .
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 ,
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 il
est trivial de vérifier que :
ce problème est donc dans NP. Il semblerait (nous
n'avons pas vérifié ce résultat et
ni la possibilité
de le faire) que la complexité du meilleur algorithme de
factorisation en nombres premiers soit du type:
(57.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:
(57.11)
Nous notons aussi {x} pour désigner la partie
fractionnaire de x. Nous avons ainsi:
(57.12)
Soit .
Alors nous avons les propriétés suivantes:
P1. ,
,
où
P2. ,
lorsque
P3. ,
si
P4.
P5. si
,
si
P6. si
P7. Si ,
alors [x / a] représente le nombre d'entiers
positifs inférieurs ou égaux à x qui sont divisibles
par a.
Démonstrations:
La première partie de P1 est simplement la définition
de [x] sous forme algébrique. Les deux autres parties
sont des réarrangements de la première partie. Dans ce cas, nous
pouvons écrire
(57.13)
où
.
Pour P2, la somme est vide (ne contient aucun terme autrement
dit) pour et,
dans ce cas, on adopte la convention selon laquelle la somme vaut
0. Alors, pour ,
la somme compte le nombre d'entiers positifs n qui sont
plus petits ou égaux à x. Ce nombre est évidemment
[x].
La démonstration de P3 sera
supposée évidente.
Pour prouver P4, nous écrivons:
,
(57.14)
où n et m sont
des entiers et où et .
Alors:
(57.15)
En écrivant ,
où ,
nous avons
(57.16)
où .
Il s'ensuit que:
(57.17)
si et:
si .
Et nous obtenons du même coup la démonstration
P5.
Pour démontrer P6, nous
écrivons:
(57.18)
où ,
et:
(57.19)
où .
Nous obtenons ainsi:
(57.20)
puisque .
Par ailleurs:
(57.21)
et nous avons ainsi le résultat.
Pour la dernière partie,
nous observons que, si sont
tous les entiers positifs qui
sont divisibles par a,
il suffit de prouver que .
Puisque ,
alors:
(57.22)
c'est-à-dire:
(57.23)
soit le résultat attendu.
C.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:
(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:
(57.25)
Et l'astuce consiste à poser que:
(57.26)
où il faut prendre comme valeur initiale (ce
qui évidemment est en totale contradiction avec les développements
mais l'on voit alors mieux pourquoi c'est une pseudo-démonstration...).
C.Q.F.D.
Exemple:
Soit à calculer:
(57.27)
Nous prendrons donc 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:
(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 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 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:
(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:
Figure: 57.1 - Principe illustré de l'algorithme d'Archimède
Le
principe de l'algorithme d'Archimède est le suivant: Soit le
périmètre d'un polygone régulier de n
côtés inscrit
dans un cercle de rayon 1/2. Archimède arrive à montrer
par induction que:
(57.30)
Nous avons pour le périmètre
d'un
n-polygone:
et
(57.31)
Avec:
(57.32)
Donc:
(57.33)
Il suffit d'un ordinateur ensuite et de plusieurs itérations
pour
évaluer avec une bonne précision la valeur de .
Évidemment, on utilise l'algorithme d'Héron pour
calculer la racine...
C.Q.F.D.
Remarque: Il existe un très grand nombre d'algorithmes
pour calculer .
Celui présenté ci-dessus, sans être le plus
esthétique, serait historiquement le premier.
CALCUL
DU NOMBRE D'EULER
Hors
la constante ,
il existe d'autres constantes mathématiques importantes
qu'il faut pouvoir générer à l'ordinateur
(de nos jours les valeurs constantes sont stockées telles
quelles et ne sont plus recalculées
systématiquement).
Parmi celles-ci, se trouve le "nombre d'Euler" noté e (cf.
chapitre d'Analyse Fonctionnelle). Voyons comment calculer
ce nombre:
Soit
la série de Taylor (cf. chapitre
sur les Suites Et Séries), pour une fonction indéfiniment
dérivable f donnée
par:
(57.34)
Comme
(cf. chapitre de Calcul Différentiel Et Intégral):
(57.35)
nous
avons:
(57.36)
Donc
en résumé:
(57.37)
Cette
relation donne un algorithme pour calculer l'exponentielle à
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:
(57.38)
Et d'après les propriétés
des logarithmes:
(57.39)
Si n est
très grand (mais alors très...) alors:
(57.40)
Lorsque ,
la limite inférieure est négligeable et alors (approximation
qui nous est très utile dans le chapitre de Mécanique
Statistique):
(57.41)
Après une petite simplification
élémentaire, nous obtenons:
(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:
(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:
(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:
(57.45)
sont
les coefficients du système d'équations, et
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:
(57.46)
Ensuite,
nous résolvons l'équation ,
ce qui donne:
(57.47)
Nous
en déduisons:
(57.48)
La
transformation entre les deux systèmes:
(57.49)
s'effectue
simplement en multipliant chaque coefficient de la première égalité
par et
en soustrayant les résultats obtenus des coefficients correspondants
de la seconde égalité. Dans ce cas, l'élément est
appelé "pivot".
TROIS ÉQUATIONS À TROIS INCONNUES
Examinons
encore le cas des systèmes de trois équations à trois inconnues:
(57.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:
(57.51)
et
ensuite résoudre l'équation:
(57.52)
puis
la deuxième:
(57.53)
et
enfin:
(57.54)
Revenons
à la transformation des systèmes. Elle s'effectue en deux étapes:
1. Dans la première,
nous choisissons comme
pivot et nous éliminons les coefficients et de
la manière suivante: il faut multiplier chaque coefficient de
la première égalité par et
soustraire les résultats obtenus de la deuxième égalité, ainsi devient
nul. De même, en multipliant les coefficients de la première équation
par et
en soustrayant les résultats obtenus de la troisième égalité, disparaît.
Le système d'équation s'écrivant alors:
(57.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 comme
pivot. Cette méthode de résolution peut paraître
compliquée,
mais elle a l'avantage de pouvoir être généralisée
et être appliquée
à la résolution de systèmes de n équations
à
n inconnues.
N ÉQUATIONS
À N INCONNUES
Pour
simplifier l'écriture, les coefficients seront toujours notés et
non pas ,
etc. lors de chaque étape du calcul.
Soit
le système linéaire (nous pourrions très bien
le représenter sous
la forme d'une matrice augmentée afin d'alléger les écritures):
(57.56)
Il
faut choisir comme
pivot pour éliminer .
Ensuite, l'élimination de s'effectue
en prenant comme
pivot. Le dernier pivot à considérer est bien évidemment ,
il permet d'éliminer .
Le système prend alors la forme:
(57.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
tel
que:
(57.58)
Il
est amené à par
échange des première et j-ème lignes. L'élimination du reste
de la première colonne peut alors être effectuée. Ensuite, on recommence
avec les n-1 équations
qui restent.
POLYNÔMES
L'ensemble
des polynômes à 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:
(57.59)
alors:
(57.60)
avec:
(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:
(57.62)
mais en imposant cette fois-ci
.
La division s'écrira
donc nous le savons:
(57.63)
avec:
(57.64)
ou sinon .
Il est normalement connu d'avance (car démontré dans le chapitre
de Calcul Algébrique) que nous avons:
(57.65)
et:
(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:
(57.67)
Pour démontrer l'expression des différents ,
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.
Exemple:
Soient:
(57.68)
donc de ce que nous avons dit précédemment, il vient
(point de départ):
(57.69)
En utilisant le fait que (pour rappel):
(57.70)
nous avons donc immédiatement:
(57.71)
Ensuite (toujours en procédant de façon identique):
(57.72)
et enfin:
(57.73)
Donc de manière générale:
(57.74)
comme:
(57.75)
le premier reste est donc:
(57.76)
Ensuite:
(57.77)
Le deuxième reste
est alors:
(57.78)
etc.
Nous continuons jusqu'à
ce que .
RÉGRESSIONS
ET INTERPOLATIONS
Les régressions (ou "interpolations")
sont des outils très utiles aux statisticiens, ingénieurs,
informaticiens souhaitant établir une loi de corrélation
entre deux (ou plus) variables dans un contexte d'études
et d'analyse ou d'extrapolation.
Il existe un grand nombre
de méthodes d'interpolation: de la simple résolution
d'équations
du premier degré (lorsque uniquement deux points d'une
mesure sont connus) aux équations permettant d'obtenir à partir
d'un grand nombre de points des informations essentielles à l'établissement
d'une loi (ou fonction) de régression linéaire,
polynomiale, 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:
(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:
(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:
(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):
(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:
(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:
(57.84)
la relation suivante:
(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):
(57.86)
Soit sous la forme plus explicite que nous utiliserons
plus loin (cf.
chapitre de Statistiques):
(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:
(57.88)
ce qui donne b sous
la forme:
(57.89)
MÉTHODE
DES MOINDRES CARRÉS
Du
fait de l'erreur sur y,
les points expérimentaux, de coordonnées ,
ne se situent pas exactement sur la droite théorique. Il
faut donc trouver l'équation de la droite expérimentale
qui passe le plus près possible de ces points.
La "méthode
des moindres carrés" (DMC) 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 et
les valeurs calculées théoriques de .
Nous parlons alors de "méthode des moindres
carrés des écarts d'ordonnées":
(57.90)
où
n est le
nombre de points et:
(57.91)
d'où autrement écrit:
(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:
(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
est la grandeur scalaire qui fait office de multiplicateur de
Lagrange.
Soit après simplification:
(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:
(57.95)
Le système devient:
(57.96)
De la deuxième équation, nous tirons:
(57.97)
En remplaçant dans la première, nous obtenons:
(57.98)
De là nous avons:
(57.99)
Ainsi, les expressions
de la pente et de l'ordonnée à l'origine de l'équation
recherchée
sont:
(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:
(57.101)
soit sous forme discrète:
(57.102)
ainsi que par construction de la méthode la relation suivante:
(57.103)
Maintenant, nous faisons l'hypothèse que chaque valeur mesurée
est entachée d'un résidu tel que:
(57.104)
Soit en soustrayant les deux dernières relations:
(57.105)
Maintenant, passons par un résultat intermédiaire. Rappelons
que nous avons obtenu plus haut:
(57.106)
En remplaçant b par sa valeur:
(57.107)
nous avons alors:
(57.108)
Multipliant la deuxième relation ci-dessus par et
retranchant de la première, nous obtenons:
(57.109)
Soit après réarrangement:
(57.110)
Revenons maintenant à:
(57.111)
Si nous mettons le tout au carré et en sommant pour toutes les
observations, nous avons:
(57.112)
soit:
(57.113)
Or, nous venons de montrer avant que le double produit était
nul. Donc:
(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:
(57.115)
où SCT est la "somme des carrés
totale", SCE la "somme
des carrés expliquée" et SCR la "somme
des carrés résiduelle".
Cette dernière relation se trouve également souvent sous la forme
suivante dans la littérature:
(57.116)
Notons maintenant les sans
erreurs d'une manière différente et appelons cela le "modèle
linéaire a priori":
(57.117)
D'où l'égalité que nous réutiliserons plusieurs
fois:
(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:
(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:
(57.120)
Cette dernière relation peut
se représenter (en gros...)
graphiquement sous la forme suivante:
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:
(57.121)
qui est simplement une autre manière d'écrire la
décomposition de la variance (implicite):
(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):
(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:
(57.124)
Ou sinon puisque nous avons démontré plus
haut que (se rappeler que la variance indiquée est la variance
biaisée...
):
(57.125)
nous pouvons aussi écrire le coefficient de corrélation
sous la forme:
(57.126)
Donc nous en déduisons):
(57.127)
Montrons que ceci est égal à (notation souvent utilisée dans
la littérature spécialisée):
(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:
(57.129)
et puisque nous avons montré que:
(57.130)
Alors:
(57.131)
C.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:
Figure: 57.3 - Quelques valeurs du coefficient de corrélation linéaire
(source: Wikipédia)
Enfin, indiquons que nous retrouvons aussi très souvent
le coefficient de corrélation linéaire sous la forme suivante dans
les logiciels et la littérature:
(57.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):
(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, est
connu sans erreur, et que est
une réalisation d'une variable aléatoire que nous
noterons dorénavant et
la droite théorique des moindres carrés s'écrira
maintenant:
(57.134)
où 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:
et
(57.135)
donc autrement dit:
(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):
(57.137)
et puisque par hypothèse , il
vient immédiatement que:
(57.138)
Raison pour laquelle le modèle s'appelle "modèle
linaéire gaussien"....
Explicitement, cela donne donc:
(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:
(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
!!!!!!!!):
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:
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):
(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 et 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 . |
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 |
|
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. |
|
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):
(57.142)
et b comme une réalisation de la variable aléatoire donnée
par:
(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:
(57.144)
nous pouvons mettre A sous la forme:
(57.145)
et B:
(57.146)
Nous en déduisons les espérances pour A:
(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:
(57.148)
avec la même conclusion.
C.Q.F.D.
Donc l'espérance de A et B sont bien des
estimateurs sans biais (donc de variance minimale pour rappel)
de .
Comme ce sont des estimateurs, dans la littérature spécialisée, sont
souvent notés et
dès lors il vient la notation courante alternative:
(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:
(57.150)
Comme par hypothèse nous avons tous les qui
sont égaux, et qu'il n'y a pas d'autocorrélation,
nous pouvons alors écrire:
(57.151)
Soit si n est assez grand nous écrirons:
(57.152)
Avant de déterminer la variance de B, rappelons
que par hypothèse:
(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):
(57.154)
Dès lors:
(57.155)
En rappelant la relation de Huyghens (cf.
chapitre de Statistiques):
(57.156)
Nous avons finalement:
(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 . É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:
(57.158)
puisque la loi Normale est centrée pour les résidus
donc ...
et que le résidu est une variable aléatoire implicitement
dépendante
de la somme de deux variables aléatoires que sont A et B d'où la
minoration de deux fois l'erreur standard.
Indiquons aussi que dans la pratique nous notons fréquemment
ce dernier résultat en mélangeant les notations de
l'aspect aléatoire
et déterministe (donc en notant tout avec des minuscules):
(57.159)
où SEE signifie en anglais "Standard
Error of Estimate" (Erreur
Standard de l'Estimation). Il s'agit en français de "l'erreur
standard de la régression" 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:
(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 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:
(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:
(57.162)
est en fait l'erreur standard de la moyenne A,
alors par analogies:
(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):
(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:
(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:
(57.166)
Nous avons alors in extenso le fameux intervalle
de confiance du coefficient de corrélation:
(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...):
(57.168)
Et que nous avons vu juste plus haut que:
(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):
(57.170)
et donc:
(57.171)
Le test d'hypothèse que nous voulons faire est alors:
(57.172)
et donc équivalent à:
(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:
(57.174)
ainsi que:
(57.175)
ce qui nous amène en mixant les deux à avoir:
(57.176)
d'où:
(57.177)
Or, nous avons aussi montré que si n est assez grand:
(57.178)
Mais si n est assez petit, nous écrirons donc:
(57.179)
Donc:
(57.180)
et avec l'hypothèse nulle ,
il vient:
(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.
Exemples:
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):
(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:
(57.183)
alors:
(57.184)
Le problème étant contourné, nous avons maintenant en utilisant
les propriétés de la variance:
(57.185)
d'où:
(57.186)
Donc nous avons:
(57.187)
Maintenant, rappelons (cf. chapitre de
Statistiques) que:
(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:
(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:
(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:
(57.191)
Donc:
(57.192)
Or comme:
(57.193)
Alors:
(57.194)
et donc:
(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:
(57.196)
Donc:
(57.197)
Soit:
(57.198)
C'est une raison pour laquelle beaucoup de statisticiens notent
directement et sans détours:
(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:
(57.200)
et il en découle en bilatéral, un intervalle de
confiance à un
niveau donné et
pour un x fixé par:
(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 choisi.
Mais pour que cela ait un sens, il faut que toutes les hypothèses
du modèle construit soient satisfaites.
Raison pour laquelle de nombreux logiciels de statistiques donnent
le graphique suivant lors de régressions linéaires
(nous y voyons bien l'intervalle de confiance):
Figure: 57.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:
(57.202)
avec la variance:
(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:
(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:
(57.205)
Ce qui nous donne avec un logiciel comme Minitab
15:
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:
(57.206)
où n est le nombre de points. Mais cette fois-ci, nous écrivons:
(57.207)
d'où autrement écrit:
(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:
(57.209)
Soit après simplification:
(57.210)
Enfin:
(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:
(57.212)
Dans un échantillon de n individus, nous mesurons pour :
Observations |
|
|
... |
|
1 |
|
|
... |
|
2 |
|
|
... |
|
|
|
|
|
|
n |
|
|
... |
|
Au fait, pour estimer les paramètres (valeurs
estimées que nous noterons afin
de respecter cette fois-ci la tradition) l'approche est très
simple car elle consiste juste en une généralisation
de la méthode des
moindres carrés que nous avons 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:
(57.213)
avec le modèle théorique estimé:
(57.214)
Donc il nous faut minimiser:
(57.215)
Les parenthèses ci-dessus peuvent être réécrites
sous la forme:
(57.216)
Nous avons le vecteur des résidus qui vaut donc:
(57.217)
Comme nous le savons, la méthode des moindres carrés consiste à trouver
le vecteur qui
minimise:
(57.218)
Soit explicitement:
(57.219)
à remarquer que nous avons:
(57.220)
et:
(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:
(57.222)
Dérivons cette dernière "fonction objet F" à l'ordre
du vecteur (il
s'agit d'une dérivée intérieure composante par composante). Ce
que nous écrivons:
(57.223)
Maintenant, passons d'une écriture vectorielle à une écriture
matricielle pure:
(57.224)
Nous cherchons donc qui
annule cette dérivée. Donc nous devons résoudre l'équation suivante:
(57.225)
Soit:
(57.226)
Rappelons avant d'aller plus loin que:
(57.227)
Donc:
(57.228)
Puisque l'algèbre linéaire est associative, écrivons sans les
parenthèses:
(57.229)
Nous ne pouvons évidemment pas simplifier à droite
et gauche par car
comme il ne s'agit pas d'une matrice carrée, ce terme est
obligatoirement non-inversible. La seule chose que nous puissions
faire c'est identifier
les éléments tels que:
(57.230)
impose forcément que:
(57.231)
Nous trouvons alors que si la matrice carrée est
inversible alors:
(57.232)
La matrice 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 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:
(57.233)
Donc en ne supposant que 2 observations, nous avons alors:
(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 en :
(57.235)
Nous avons dans le cas présent:
(57.236)
Nous avons donc:
(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 (cf.
chapitre Algèbre Linéaire)
(57.238)
Donc:
(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...):
(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:
(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:
(57.242)
Ce que nous pouvons considérer sous la forme:
(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:
(57.244)
Mais sous cette convention d'écriture, nous avons
alors:
(57.245)
et écrivons cela sous forme vectorielle:
(57.246)
Maintenant, utilisons la technique du maximum de
vraisemblance (cf. chapitre de Statistique)
et cherchons la valeur des coefficients qui maximise donc:
(57.247)
Ce que nous pouvons écrire sous forme matricielle:
(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:
(57.249)
Maintenant cherchons l'expression de qui
maximise la log-vraisemblance. Il vient alors:
(57.250)
Utilisons la propriété de la transposée
démontrée dans le chapitre d'Algèbre Linéaire:
(57.251) Il vient:
(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):
(57.253)
Soit:
(57.254)
C'est-à-dire après réarrangement:
(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:
(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:
(57.257)
Nous avons quelque chose d'intéressantà observer:
(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é:
(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:
(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:
(57.261)
où 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 qui
minimisent l'erreur:
(57.262)
et nous réattaquons avec des dérivées partielles pour chaque coefficient:
(57.263)
Soit après simplification et un petit réarrangement:
(57.264)
De même:
(57.265)
Soit après simplification et un petit réarrangement:
(57.266)
Et enfin:
(57.267)
Soit après simplification et un petit réarrangement:
(57.268)
Donc en utilisant la notation de l'algèbre linéaire, nous avons
finalement le système suivant à résoudre:
(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é ,
qui nous donne à l'opposé la probabilité d'appartenir à la
classe ,
que nous nommerons "régression
logistique" ou "régression
logit" ou encore "régression
binomiale" (très
souvent utilisée
dans le cadre des réseaux de
neurones formels que nous verrons plus loin). Ensuite, dans une
deuxième étape,
nous définissons
pour un cas binaire une valeur "cutoff".
Par exemple, si nous prenons un cutoff de 0.5 alors les cas pour
lesquels appartiendront à la
classe 1 (et inversement dans le cas contraire).
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.
Exemple:
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.
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 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!):
(57.270)
pour une variable de prédiction (explicative) x où P 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:
(57.271)
S'il y a plusieurs variables prédictives nous avons alors:
(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:
(57.273)
il en découle donc la fonction de distribution:
(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:
(57.275)
Indiquons que si a est posé comme unitaire, et b comme
nul, nous avons alors la "loi logistique
standard" donnée par:
(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:
(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:
(57.278)
il vient alors:
(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!):
(57.280)
Ainsi, nous voyons que si nous posons:
(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:
(57.282)
peut par ailleurs être transformée de la
façon suivante:
(57.283)
d'où:
(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...:
(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):
(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:
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:
(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:
et
(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:
(57.289)
Exemple:
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 .
Quand le crédit est inférieur à 27'500, la proportion de bons débiteurs
est de .
Quand le montant des crédits est inférieur à 28'000 la proportion
de bons débiteurs est de .
Quand le montant des crédits est inférieur à 28'500, la proportion
de bons débiteurs est de ,
et pour des montants inférieurs à 30'000 la proportion est de .
Nous poserons pour cette régression logistique que est
un bon risque de crédit et est
un mauvais risque. Ensuite, nous créons le tableau suivant qui
est un récapitulatif des données de toutes les succursales:
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:
Figure: 57.9 - Pourcentage cumulé des bons débiteurs en fonction du crédit
Une fois ceci fait, nous utilisons la transformation en logit:
(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:
Figure: 57.10 - Logit des bons débiteurs en fonction du montant du crédit
avec pour équation:
(57.291)
La fonction logistique avec sa représentation vient alors
immédiatement (les unités de x sont en
milliers de francs!):
(57.292)
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):
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):
(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):
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:
(57.294)
mais que nous avons affiné avwec un logiciel de
statistiques spécialisé afin d'obtenir:
(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):
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:
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:
Figure: 57.15 - Formules du tableur pour construire la matrice de confusion
ce qui donne explicitement:
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:
Figure: 57.17 - Matrice de confusion de Minitab 16.1.1
Maintenant, observons la courbe ROC donnée par Minitab:
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:
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é:
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:
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):
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é):
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:
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:
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):
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:
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:
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 est:
(57.296)
Ce qui est juste puisque lorsque nous
sommes en A et lorsque nous
sommes en B. Donc et
le point M parcourt tout le segment [AB]. Par construction,
si A et B étaient des masses physiques égales à l'unité, 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:
Figure: 57.29 - Courbe de Bézier avec un point supplémentaire
1ère étape:
- Soit le
barycentre de (A,1-t) et (B, t) et
où décrit
[AB].
- Soit le
barycentre de (B,1-t) et (C, t) et où décrit
[BC].
2ème étape:
- Soit M(t) le
barycentre de .
Par construction, M(t) se
situe donc à la même proportion du segment que par
rapport au segment [AB] ou par
rapport au segment [BC].
La courbe obtenue est alors l'enveloppe des segments :
en tout point M, la tangente à la courbe est donc le segment .
M(t) décrit
alors une Courbe de Bézier de degré 2, qui, par construction
commence en A et se finit en C, et a pour tangentes [AB]
en A et [BC] en C.
C'est en fait un arc de parabole (que nous pourrions noter très
logiquement [ABC] ):
Figure: 57.30 - Arc de parabole [ABC]
Par le même schéma, nous pouvons définir une courbe de Bézier
de n points soit .
C'est ce que nous appelons "l'algorithme
de Casteljau".
Ainsi, soit:
(57.297)
Nous avons:
(57.298)
La récurrence se terminant pour:
(57.299)
Ainsi, pour nous
avons:
(57.300)
Soit:
(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:
(57.302)
Nous avons par la relation de récurrence:
(57.303)
où nous avons éliminé les termes contenant des
points non déterminés.
Nous avons donc:
(57.304)
Il vient alors:
(57.305)
Et donc:
(57.306)
soit sous forme vectorielle (plus conforme à la
notation mathématique d'usage):
(57.307)
et sous forme matricielle:
(57.308)
Par le même raisonnement, nous avons pour une courbe de Bézier
d'ordre 4:
(57.309)
soit sous forme vectorielle:
(57.310)
Ce qui correspond de manière générique à:
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:
(57.311)
Nous remarquons d'abord aisément la proportionnalité suivante:
(57.312)
et si nous regardons de plus près les coefficients, nous
remarquons que nous avons aussi:
(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:
(57.314)
Ainsi, les polynômes de Bernstein sont donnés par:
(57.315)
et finalement les courbes de Bernstein d'ordre n par:
(57.316)
Au fait, si nous avions noté plus haut la somme sous la forme
suivante:
(57.317)
Nous aurions alors les polynômes de Bernstein qui seraient
donnés
(ce qui est plus respectueux des traditions...) par:
(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:
(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 est
défini après coup à la souris en utilisant
des poignées
appelées "torseurs" dans
le langage de spécialistes Adobe... Voici un exemple pris
d'un de ces logiciels fait avec un tracé à la plume
de 5 points (soit 4 splines):
Figure: 57.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.
Exemple:
Un cercle, dessiné par un logiciel graphique est en pratique
composé de 4 arcs de Bézier. Pour observer cette particularité, il suffit de dessiner un cercle
avec Illustrator par exemple, puis de le sélectionner pour voir
apparaître les points de contrôle des arcs de Bézier qui le forment.
Nous allons nous intéresser à la meilleure façon de choisir les
points de contrôle de ces arcs de sorte qu'ils ressemblent à des
quarts de cercle, puis nous observerons la différence entre le
dessin produit et un vrai cercle:
Figure: 57.33 - Exemple de construction d'un cercle avec des courbes de Bézier
Prenons le quart de cercle de rayon 1 centré à l'origine:
Il est approché par un arc de Bézier dont les points de contrôle
sont .
Les extrémités de l'arc de Bézier étant et ,
il est naturel de choisir et .
L'intuition nous amène à choisir et et
il reste à trouver une valeur positive de k de sorte que
l'arc de Bézier ressemble à un arc de cercle.
Nous obtenons ainsi l'équation paramétrique de l'arc de Bézier
(57.320)
Soit:
(57.321)
Nous pouvons par exemple chercher la valeur de k pour
laquelle l'arc passe par le point:
(57.322)
en .
Il devient alors très simple à partir de l'équation paramétrique
de déterminer k. Il s'agit simplement pour x (ou
pour y) d'une simple équation à une inconnue.
MÉTHODE
D'EULER
Il s'agit là de la
méthode numérique la plus simple (elle est triviale
et dans l'idée assez proche de la méthode de Newton
même si leur objectif n'est pas le même). En fait
elle ne fournit 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
(h étant le "pas"). Nous notons
la valeur exacte et ,
la valeur approchée.
Il y a plusieurs méthodes
pour procéder (comme souvent):
1. Graphiquement:
Nous nous déplaçons
d'un pas h en x en suivant le vecteur de pente f(x,y)
Par construction, nous savons
(cf. chapitre de Calcul
Différentiel Et Intégral) que (nous adoptons
une notation un peu particulière dans ce contexte):
(57.323)
qui correspond donc bien
à la pente (non instantanée bien sûr!) en de
la "courbe intégrale" passant par ce point. D'où:
(57.324)
2. Analytiquement:
Nous remplaçons dans
la dernière relation
par .
Nous obtenons alors:
(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
une fonction connue sous forme explicite ou sous forme tabulée,
et supposons qu'un certain nombre de valeurs:
(57.326)
en sont données. Les
points
sont appelés les "points d'appui".
"Interpoler" f signifie
estimer les valeurs de f pour des abscisses situées
entre
et ,
c'est-à-dire dans l'intervalle d'interpolation, par une
fonction approximante ,
qui vérifie les "conditions de collocations" (rien à
voir avec votre colocataire!):
Figure: 57.34 - Illustration du concept d'interpolation
La fonction p s'appelle "fonction
de collocation" sur les .
Lorsque p est un polynôme, nous parlons de "polynôme de collocation"
ou de "polynôme d'interpolation".
"Extrapoler" f signifie
approcher f(x) par p(x) pour des abscisses situées "hors" de l'intervalle d'interpolation.
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é:
(57.327)
et nous avons connaissance
des points suivants (dont vous remarquerez l'ingéniosité des
points choisis par les auteurs de ces lignes...):
(57.328)
Nous en déduisons donc le
système d'équations:
(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:
(57.330)
Voyons le cas général:
Théorème:
Soient
des points d'appui, avec
si .
Alors il existe un polynôme
de degré inférieur ou égal à n,
et un seul, tel que
pour .
Démonstration:
Posons:
(57.331)
Les conditions de collocation:
(57.332)
s'écrivent donc:
(57.333)
Il s'agit d'un système
de n+1 équations à n+1 inconnues.
Son déterminant s'écrit:
(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 ,
nous considérons le déterminant:
(57.335)
qui correspond donc au système
(pour rappel):
(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):
(57.337)
Ce dernier polynôme
peut s'écrire:
(57.338)
Ce qui s'écrit:
(57.339)
Comme les
sont dans l'énoncé de notre problème tous
différents
tels que
alors le système a une solution unique! Ce qui prouve qu'il
existe toujours alors un polynôme d'interpolation.
C.Q.F.D.
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:
(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:
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 nous
associons un polynôme
de degré n et non nul en mais
nul en tous les autres points mesurés et
que tous les autres polynômes associés
aux autres points 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:
(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 (vous
pouvez le voir sur le graphique). Les coefficients sont
des constantes à déterminer.
Maintenant, rien nous nous empêche de sommer
puisque pour chaque nous
aurons toujours la bonne valeur (puisque
les autres polynômes
y sont nuls). Soit en généralisant
l'écriture, la somme devient:
(57.342)
En injectant respectivement dans
la relation précédente,
il vient:
(57.343)
Nous en déduisons alors immédiatement:
(57.344)
En substituant les valeurs des constantes dans l'expression
initiale de la somme, nous obtenons:
(57.345)
Ce qui peut s'écrire sous forme condensée:
(57.346)
Où le terme:
(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 ne
possède pas de solution algébrique (sauf accident...)
si est
un polynôme de degré supérieur à 4.
Il
existe un grand nombre d'algorithmes permettant de calculer les
racines de l'équation avec
une précision théorique arbitraire. Nous n'en verrons que les principaux.
Attention, la mise en oeuvre de tels algorithmes nécessite toujours
une connaissance approximative de la valeur cherchée et celle du
comportement de la fonction près de la racine. Un tableau des valeurs
de la fonction et sa représentation graphique permettent souvent
d'acquérir ces connaissances préliminaires.
Si
l'équation à résoudre est mise sous la forme ,
nous traçons les courbes représentant g et
h.
Les racines de l'équation étant
données par les abscisses des points d'intersection des deux courbes.
MÉTHODE
DES PARTIES PROPORTIONNELLES
La
mise en oeuvre, sur calculatrice, de cette méthode est particulièrement
simple. Les conditions à vérifier étant seulement:
-
f est
continue
-
f est
monotone dans un voisinage de la racine
Dans
un petit intervalle, nous pouvons remplacer une courbe par un segment
de droite. Il y a plusieurs situations possibles mais en voici
une
particulière généralisable facilement à n'importe
quoi:
Figure: 57.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):
(57.348)
d'où:
(57.349)
Si
,
nous pouvons négliger f(a) au
dénominateur et il vient:
(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 .
Si est
assez petit, nous arrêtons le calcul et affichons x1 et
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
-
nous retournons au point (2)
MÉTHODE
DE LA BISSECTION
La
condition préalable à satisfaire pour cette méthode est de trouver
un intervalle tel
que:
1. f(x) est
continue sur [a,b]
2.
Il
faut encore fixer qui
est défini comme la borne supérieure de l'erreur
admissible.
La
méthode consiste à appliquer successivement les 4 étapes suivantes:
1.
Calcul de
2.
Évaluation de f(x)
3.
Si alors
le travail est terminé, il faut afficher x et
f(x)
4.
Sinon on procède comme suit:
-
on remplace a par
x si
-
on remplace b par
x si
ou
-
on retourne en (1)
L'étape
(3) impose la condition pour
l'arrêt des calculs. Il est parfois préférable de choisir un autre
critère de fin de calcul. Celui-ci impose à la solution calculée
d'être confinée dans un intervalle de longueur contenant
.
Ce test s'énonce comme suit:
3'.
Si ,
le travail est terminé et est
affiché. Il est bien sûr évident que
MÉTHODE
DE LA SÉCANTE (REGULA FALSI)
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:
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.
Si
est
le point de coordonnées ,
alors les points sont
alignés sur la sécante. La proportion suivante (Thalès) est donc
vraie:
(57.351)
nous
en déduisons:
(57.352)
La
méthode consiste à appliquer successivement les étapes suivantes:
1.
Calcul de
2.
Évaluation de
3.
Si ,
le travail est terminé. Il faut afficher
4.
Sinon nous procédons comme suit:
-
nous remplaçons a par
si
-
nous remplaçons b par
si
ou
-
nous retournons en (1)
La
condition (3) peut être remplacée par la condition:
3'.
Si ,
alors le travail est terminé et nous affichons
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:
Figure: 57.38 - Illustration de la méthode de Newton dans le plan
Si est
une approximation de la racine ,
nous remarquons que en
est une meilleure. est
l'intersection de la tangente à la courbe en et
de l'axe des abscisses. est
encore une meilleure approximation de ,
est
obtenu de la même manière que mais
à partir de .
La méthode de Newton consiste
en la formalisation de cette constatation géométrique.
Pour utiliser cette technique,
rappelons que si nous prenons une fonction f qui
est dérivable en ,
alors nous pouvons la réécrire sous la forme (cf.
chapitre de Calcul Différentiel Et Intégral):
(57.353)
où est
la dérivée de f en
et
est
une fonction qui tend vers 0 comme pour
lorsque
x tend
vers (c'est
un terme correctif qui sous-tend la suite des termes du développement
de Taylor).
En appliquant ce résultat
à la résolution de ,
nous obtenons:
(57.354)
La fonction empêche
la résolution de cette équation par rapport à .
En négligeant le terme ,
l'équation se réécrit:
(57.355)
et se résout aisément par
rapport à .
Pour voir cela, commençons par :
(57.356)
Mais
ne
satisfait pas, en général, l'égalité .
Mais comme nous l'avons déjà souligné, est
plus petit que si
la fonction f satisfait
à certaines conditions.
La
méthode de Newton consiste à remplacer l'équation:
(57.357)
par:
(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
et
il
faut que:
1.
La fonction soit deux fois dérivable
2. La dérivée f ' ne s'annule pas (monotonie)
3.
La deuxième dérivée soit continue et ne s'annule pas (pas de point
d'inflexion)
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 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 et
.
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
par
le quotient différentiel:
(57.359)
où
h doit
être choisi suffisamment petit pour que la différence:
(57.360)
soit
elle aussi suffisamment petite.
L'itération
s'écrit alors:
(57.361)
Si
la méthode de résolution est convergente, l'écart
entre et diminue
à chaque itération. Ceci est assuré, par exemple,
si l'intervalle [a,b] contenant ,
voit sa longueur diminuer à chaque étape. La méthode
de Newton est intéressante car la convergence est quadratique:
(57.362)
alors
que la convergence des autres méthodes est linéaire:
(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 est
réduit de moitié à chaque étape du calcul:
(57.364)
Pour
démontrer la convergence quadratique de la méthode de Newton, il
faut utiliser les développements limités de f et
f ' au
voisinage de :
(57.365)
Mais:
(57.366)
donc:
(57.367)
En
soustrayant à
gauche et à droite de l'égalité et en mettant les
deux termes du second membre au même dénominateur, il vient:
(57.368)
et
dès que est
assez petit, le dénominateur peut être simplifié.
(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);
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:
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):
(57.370)
Nous avons alors sur la base de ce principe le développement
au deuxième ordre:
(57.371)
Il alors en négligeant les termes d'ordre
supérieur et en soustrayant et simplifiant les deux séries
ci-dessus:
(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:
(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:
(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):
(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:
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 et
.
Nous supposons dans ce cas que la fonction f est
à valeurs positives:
(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 en
n sous-intervalles
dont les bornes sont .
Les longueurs de ces sous-intervalles sont .
Nous construisons les rectangles dont les côtés sont et
.
Figure: 57.43 - Approche de l'aire sous une courbe par des
rectangles
inférieurs gauche
L'aire
de ces rectangles vaut:
(57.377)
Si
les sont
suffisamment petits, est
une bonne approximation de l'aire cherchée approchée par
la gauche.
Nous pouvons recommencer cet exercice en choisissant et
comme
côtés des rectangles (donc l'approche par la droite).
Nous obtenons alors:
(57.378)
La
figure correspondante est la suivante:
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:
(57.379)
Si
nous avons n rectangles,
h vaut
alors .
Les aires et deviennent:
(57.380)
MÉTHODE
DES TRAPÈZES
Afin
d'augmenter la précision des calculs, il est possible de
calculer:
(57.381)
Dans
le cas où tous les intervalles sont de longueur égale, vaut:
(57.382)
Ce que nous retrouvons souvent dans la littérature scolaire, sous
la forme:
(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:
(57.384)
et:
(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:
ou encore depuis la version 2010 du même logiciel (dont l'interface
a complétement changé):
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:
(57.386)
où les
sont
des variables qui influent sur la valeur de Z, et les
poids respectifs de ces variables modélisant l'importance relative
de chacune de ces variables sur la valeur de la fonction économique.
Les contraintes
relatives aux variables s'expriment par le système linéaire suivant:
(57.387)
Sous forme
générale et matricielle ce genre de problème s'écrit:
(57.388)
Exemple:
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):
(57.389)
La fonction
économique à maximiser étant:
(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":
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
maximisant la fonction économique.
Or, l'équation Z est représentée
par une droite de pente constante (-1.2) dont tous les points fournissent
la même valeur Z pour la fonction économique.
En particulier, la droite
passe
par l'origine et 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 plan)
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.
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":
(57.391)
avec:
(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":
(57.393)
Donc pour 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:
(57.394)
La contrainte A1 devient alors:
(57.395)
et la contrainte A2 devient respectivement:
(57.396)
Dès lors, le problème se résume à maximiser Z avec
les contraintes:
(57.397)
Partons donc d'une solution réalisable évidente
qui au vu des contraintes est trivialement:
(57.398)
Dès lors, avec le système:
(57.399)
nous trouvons immédiatement:
(57.400)
Les paramètres dans l'état actuel peuvent se résumer à:
(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:
(57.402)
c'est-à-dire (nous
pensons que c'est ainsi que Z croîtra le plus
vite). Nous parlons alors de comme
"direction du pivot".
Nous gardons donc et
nous faisons croître avec
le système qui se réduit alors à:
(57.403)
Avec donc et pour
commencer,
nous avons:
(57.404)
et nous voyons que les contraintes sont
toujours respectées, il en va de même si vaut
2, 3, 4, 5, ... et ce jusqu'à 31, car dès lors:
(57.405)
et l'une des variables d'écart étant devenue négative,
les contraintes
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 tout
en maintenant la condition lorsque
?
Et la réponse est assez simple:
(57.406)
et donc il s'agit de:
(57.407)
Nous parlons alors parfois du "pas
du pivot". Nous avons alors la solution actuelle:
(57.408)
Ce qui donne:
(57.409)
Graphiquement parlant, voilà à quoi correspond
ce que nous venons de faire:
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:
(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
en fonction de puisque
nous avions pour rappel:
(57.411)
Pour la suite, il nous faut exprimer ainsi
que Z en fonction de puisque
nous venons d'obtenir pour rappel:
(57.412)
Avant d'obtenir le nouveau système fonction de ,
faisons quelques manipulations algébriques:
(57.413)
ce qui nous donne après simplification:
(57.414)
et dès lors, il vient:
(57.415)
ce qui nous donne après simplification:
(57.416)
et nous avons de même:
(57.417)
Donc au final, le système est:
(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 et
donc cela nous amène à poser .
La valeur la plus contraignante de qui
permet de respecter toujours les contraintes est
dès lors:
(57.419)
Et pour cette valeur, nous avons:
(57.420)
La fonction économique initiale prend alors la valeur:
(57.421)
ce qui correspond donc graphiquement à:
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:
(57.422)
Et nous allons donc réécrire le système:
(57.423)
avec cette fois-ci ainsi
que Z en fonction de .
Nous avons alors d'abord:
(57.424)
et nous avons:
(57.425)
Et donc pour Z il vient (les coefficients y sont
tous négatifs donc nous devinons la suite...):
(57.426)
Nous avons donc le nouveau système:
(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:
(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:
(57.429)
C'est-à-dire
aussi, en utilisant des notations matricielles:
(57.430)
où
les matrices correspondent,
respectivement, aux coefficients des niveaux d'activité dans
la fonction objectif, aux coefficients techniques des activités
et aux seconds membres des contraintes.
Nous allons voir maintenant 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:
(57.431)
peut
être remplacée par l'équation:
(57.432)
impliquant une variable supplémentaire, ,
appelée donc "variable d'écart", et soumise à la contrainte
de non-négativité, .
Bien
évidemment, dans un cas contraire tel où le système est du type:
(57.433)
Nous
écrirons:
(57.434)
impliquant
donc également une variable supplémentaire et soumise à la contrainte
de non-négativité, .
Ce
travail de mise en forme standard nous permet donc de retrouver
un système d'équations linéaires à résoudre
(nous avons vu précédemment
sur le site 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 ,
une contenant les variables initiales D et
l'autre comportant les variables d'écart B tel
que:
(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:
(57.436)
d'ajouter
une variable d'écart tel que:
(57.437)
où
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 .
Les variables associées aux colonnes de la matrice D seront
appelées les "variables hors-base"; il s'agit des variables
.
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 (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 ".
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:
(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:
(57.439)
ou
bien aussi:
(57.440)
Si
la matrice B est
une matrice de base, elle est non singulière et admet donc une
matrice inverse .
En multipliant cette équation, à gauche et à droite, par nous
obtenons:
(57.441)
Le
système d'équations aura alors été mis sous une forme résolue en
.
Pour
obtenir une forme canonique utilisable associée à la base B,
correspondant aux variables de base, il ne reste plus qu'à éliminer
les variables de base de l'expression de la fonction économique.
Écrivons
cette fonction en séparant les vecteurs et
,
nous obtenons:
(57.442)
Nous
pouvons alors facilement exprimer en
fonction de .
En utilisant le système d'équations mis sous forme résolue en ,
nous avons dans un premier temps:
(57.443)
que
nous substituons dans la fonction économique, pour obtenir:
(57.444)
Nous
regroupons les termes en et
nous avons:
(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 de
la ligne .
Soit
un problème de programmation linéaire sous forme standard:
(57.446)
La
matrice A à
m lignes
(autant qu'il y a de contraintes) et n colonnes,
avec .
Si nous sélectionnons m variables
de base et si nous annulons les variables
hors base, la matrice A:
(57.447)
et
le système se réduit à:
(57.448)
La
matrice B est
de dimension .
Si elle définit une base, elle admet une matrice inverse .
Une solution du système est donc:
(57.449)
Si
l'expression est
non négative ,
nous avons une solution admissible qui vérifie les contraintes
et que nous appellerons un "programme
de base":
(57.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":
(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 alors
le programme de base est optimal.
Proposition
2: La solution d'un problème de P.L., si elle existe, peut toujours
être trouvée en un programme de base.
La
démonstration précise de ce résultat est assez délicate. Nous pouvons
cependant en avoir une intuition en considérant, une fois de plus,
la notion de coût réduit.
En
effet, pour un programme de base donné, considérons la forme canonique
utilisable associée à la base. Sur cette forme nous pouvons vérifier
que, ou bien le programme de base est optimal (tous les coûts
réduits
sont ),
ou bien que le programme de base peut être amélioré et remplacé
par un nouveau programme de base donnant à z une
valeur plus petite (un coût réduit est négatif et la variable
hors-base associée peut être augmentée jusqu'à ce qu'une ancienne
variable de base s'annule). Comme il y a un nombre fini de programmes
de
base (au plus égal au nombre ),
la solution de P.L. se trouve nécessairement en un programme de
base.
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:
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:
(57.452)
avec les contraintes:
(57.453)
ce que nous pouvons en effet vérifier visuellement.
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:
(57.454)
où pour rappel h et k sont
des variables et
sont
fixés et où nous avons la matrice hessienne:
(57.455)
que les américains spécialistes du
domaine ont pour habitude (malheureuse à mon avis...) de
noter:
(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:
(57.457)
Ce qui nous donne une forme plus condensée
et technique de la série de Taylor autour de
:
(57.458) En changeant un peu la notation:
(57.459)
Nous retrouvons donc l'expression d'usage d'une fonction de évaluée
en série de Taylor centrée en .
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:
(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 .
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:
(57.461)
Et rappelons que tous les termes sont
des constantes car il s'agit soit de la fonction évaluée sur
le point particulier , 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:
(57.462)
et en revenant aux notations d'usage du domaine des méthodes
numériques, nous avons alors:
(57.463)
Et donc le gradient devant être nul, nous avons:
(57.464)
et après un premier réarrangement:
(57.465)
et un deuxième réarrangement:
(57.466)
ce qui se note souvent:
(57.467)
et par les américains...:
(57.468)
Enfin avant de passer à un exemple concret il est important
que le lecteur se souvienne de la relation vu juste plus haut:
(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:
(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.
Exemple:
Nous cherchons un extrema local de la fonction "baleine à bosse" représentée
plus haut:
(57.471)
avec le point de départ (arbitraire):
(57.472)
Pour effectuer la recherche, nous calculons d'abord le gradient:
(57.473)
et la matrice hessienne:
(57.474)
Nous avons alors:
(57.475)
et donc:
(57.476)
et nous recommençons (bon on va se passer de tous les
détails maintenant car faut pas exagérer non plus...):
(57.477)
et donc:
(57.478)
et nous recommençons (avec encore moins de détails):
(57.479)
et donc:
(57.480)
et encore une fois (avec encore moins de détails):
(57.481)
et donc:
(57.482)
et encore une fois (avec encore moins de détails):
(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:
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:
|
|
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!):
(57.484)
Nous cherchons donc qui
minimisent la somme des carrés
des écarts entre les valeurs expérimentales et théoriques
tels que:
(57.485)
avec donc:
(57.486)
Notons pour la suite comme le veut la tradition dans le domaine:
(57.487)
Nous avons alors la notation courante:
(57.488)
Maintenant, imaginons que nous ayons trouvé une valeur
du couple qui
donne ce minimum et notons le 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:
(57.489)
Dès lors, il en découle immédiatement
que:
(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):
(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:
(57.492)
et la solution compatible nous amène donc à bien évidemment à:
(57.493)
De la même manière, nous avons:
(57.494)
Nous avons donc au final les deux relations suivantes:
(57.495)
Étant donné que pour la solution compatible nous
avons:
(57.496)
Il s'ensuit que dans ce cas que la deuxème relation devient:
(57.497)
où H est la matrice Hessienne (cf.
chapitre de Suites Et Séries) ce que les américains notent simplement:
(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:
(57.499)
Maintenant, rappelons la relation de base de la méthode
de Newton-Raphson obtenu plus haut:
(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!):
(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:
(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:
(57.503)
Donc en multipliant simplement par un facteur ½ (ce
qui ne change rien au résultat). Nous avons alors:
(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:
(57.505)
Nous commençons avec un couple qui nous semble proche
de la solution cherchée:
(57.506)
Nous avons alors:
(57.507)
et donc nous avons:
(57.508)
Pour la suite, il vient alors:
(57.509)
Ce que nous pouvons donc récrire sous la forme:
(57.510)
Nous avons aussi in extenso:
(57.511)
Ensuite, nous appliquons la relation démontrée plus haut:
(57.512)
Soit:
(57.513)
et après une petite simplification mineure:
(57.514)
Soit:
(57.515)
et donc le prochain couple pour l'itération sera:
(57.516)
Ce qui correspond bien évidemment aux valeurs pour la
première itération:
(57.517)
Nous n'allons pas refaire aussi explicitement les autres itérations.
Donc voilà ce que cela donne au final:
i |
|
0 |
|
1
|
|
2
|
|
3
|
|
4
|
|
Tableau: 57.10 - Itérations Gauss-Newton
avec donc pour solution
locale à la 4ème itération:
(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:
(57.519)
et nous écrirons la fonction de la façon suivant
pour la méthode de Newton-Raphson:
(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):
(57.521)
et:
(57.522)
Nous avons donc:
(57.523)
qui devient:
(57.524)
Soit:
(57.525)
et donc le prochain couple pour l'itération sera:
(57.526)
Ce qui correspond bien évidemment aux valeurs pour la
première itération:
(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 |
|
0 |
|
1
|
|
2
|
|
3
|
|
4
|
|
5 |
|
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();
>restart;rand();
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();
>restart;readlib(randomize):randomize():rand();
Épreuve
a priori réussie (au fait, il nous faudrait faire un beaucoup
plus grand nombre d'essais afin de bien vérifier que le
générateur ne
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):
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:
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]:
(57.528)
Soit:
(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 :
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:
(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 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:
(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...):
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):
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):
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:
- 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é.
- 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é.
- 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 d'une
variable aléatoire. Pour cet exemple,
l'intervalle de confiance de l'espérance
d'une variable aléatoire
est parfaitement déterminé à partir
de la moyenne et de la variance calculées sur l'échantillon
(cf.
chapitre de Statistiques).
Nous considérons un échantillon de la variable aléatoire
composé de estimations:
(57.532)
La moyenne arithmétique de l'échantillon est:
(57.533)
et son écart-type (estimateur de maximum de vraisemblance non
biaisé):
(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:
(57.535)
où S est une autre notation traditionnelle dans certains
domaines de la statistique pour la notation de l'écart-type empirique
(cf. chapitre de Statistiques). Nous
avons alors pour l'intervalle de confiance à 95% de l'espérance:
(57.536)
Soit:
(57.537)
Ce qui donne:
(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:
(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):
(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.
Exemple:
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:
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:
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):
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:
(57.541)
et le fait que la racine se trouve dans l'intervalle:
(57.542)
Nous calculons:
(57.543)
Si alors
la racine est dans l'intervalle sinon
elle est dans l'intervalle .
Nous avons donc ramené le problème à une taille inférieure.
Nous arrêterons
l'algorithme quand la précision sera suffisante.
L'algorithme
Maple 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!!
|