# Programmer en R/Optimiser une fonction

* Introduction
** Faire ses premiers pas
** Installer R
** Choisir son environnement de travail
** Installer des librairies
** Se documenter
* Manipulation et mise en forme des données
** Manipuler les nombres
** Manipuler les vecteurs
** Manipuler les chaînes de caractères
** Manipuler les facteurs
** Manipuler les dates et les heures
** Manipuler un tableau de données
** Manipuler une liste
* Importation et exportation des données
** Importer un fichier CSV
** Résoudre les problèmes d'encodage des caractères
** Importer un fichier TSV
** Importer un fichier texte à largeur fixe
** Importer des données depuis un tableur
** Import un fichier JSON
** Extraire des données d'une page web
* Visualisation de données
** Comprendre la grammaire des graphiques
** Représenter une carte
** Visualiser une carte choroplèthe
* Modélisation statistique
** Estimer un modèle de régression linéaire
** Estimer un modèle de régression circulaire
* Communication
* Programmation avancée
** Programmer une boucle
** Programmer une fonction
* Outils mathématiques
** Optimiser une fonction
** Manipuler les matrices
** Les fonctions statistiques de base
** Les loi de probabilités, ajustement et test
* Introduction
** Faire ses premiers pas
** Installer R
** Choisir son environnement de travail
** Installer des librairies
** Se documenter
* Manipulation et mise en forme des données
** Manipuler les nombres
** Manipuler les vecteurs
** Manipuler les chaînes de caractères
** Manipuler les facteurs
** Manipuler les dates et les heures
** Manipuler un tableau de données
** Manipuler une liste
* Importation et exportation des données
** Importer un fichier CSV
** Résoudre les problèmes d'encodage des caractères
** Importer un fichier TSV
** Importer un fichier texte à largeur fixe
** Importer des données depuis un tableur
** Import un fichier JSON
** Extraire des données d'une page web
* Visualisation de données
** Comprendre la grammaire des graphiques
** Représenter une carte
** Visualiser une carte choroplèthe
* Modélisation statistique
** Estimer un modèle de régression linéaire
** Estimer un modèle de régression circulaire
* Communication
* Programmation avancée
** Programmer une boucle
** Programmer une fonction
* Outils mathématiques
** Optimiser une fonction
** Manipuler les matrices
** Les fonctions statistiques de base
** Les loi de probabilités, ajustement et test
[ Modifier le sommaire ]
Les méthodes d'optimisation sont des méthodes numériques de recherche du minimum ou d'un maximum d'une fonction.
Soit une fonction ƒ à valeurs réelles, définie sur un espace vectoriel de dimension n :
    y = ƒ(x$_1$, x$_2$, …, x$_n$)
On recherche le minimum de ƒ, c'est-à-dire un vecteur (x$_1$*, x$_2$*, …, x$_n$*) tel que la valeur de ƒ soit la plus petite possible.
Les méthodes de régression consistent à ajuster un modèle à des données. Elles consistent à minimiser une fonction d'écart, à l'optimiser.

## Méthode d'optimisation sans contrainte
### Fonction sans paramètre

Dans le cas d'une optimisation sans contrainte, R propose la commande optimise ou optimize (synonymes). Sa syntaxe globale est :
où
* fct est la fonction à minimiser, voir Programmation procédurale > Fonctions définies par l'utilisateur ;
* intervalle est l'intervalle sur lequel il effectue la recherche ;
* maximum = est un booléen : « vrai » si l'on cherche la maximum, « faux » (valeur par défaut) si l'on cherche le minimum.
Si par exemple on désire connaître le minimum de la fonction
    ƒ(x) = sin(x)
sur l'intervalle [-π ; π], on peut écrire
ce qui nous donne
Donc, le minimum est en −1,570,798, qui est une approximation de -π/2 ≃ −1,570,796, et la fonction y vaut -1. Le maximum est en 1,570,798, qui est une approximation de π/2, et la fonction y vaut 1.
On peut aussi utiliser la syntaxe
On peut lui indiquer la tolérance visée sur le résultat, avec l'option tol =. Par exemple, pour avoir une précision de 10$^{-6}$ :
Nous disposons également de la commande optim ; elle ne cherche que le minimum, elle ne permet pas par défaut de définir un intervalle, et il faut lui donner un point de départ ; par contre, elle permet de travailler avec des fonctions de plusieurs variables. Sa syntaxe globale est :
où
* x$_0$ est un vecteur de valeurs initiales ;
* f est la fonction à minimiser.
Si par exemple on désire connaître le minimum de la fonction
    ƒ(x) = 5x$^2$ + 3x - 4
on peut écrire (en utilisant la méthode de Horner pour calculer le polynôme)
Le solveur nous donne un message d'avertissement qui est sans conséquence.
L'abscisse de la solution est donnée par
qui est la solution exacte, et la valeur de la fonction en ce point est
Si l'on affiche la totalité de la solution, on obtient
qui nous indique :
* $count : le nombre de fois (26) qu'il a fallu appeler la fonction f ; c'est globalement le nombre d'itérations qui a été nécessaire pour trouver la solution ;
* $convergence : la qualité de la convergence ; « 0 » indique que la méthode a bien convergé ;
* $message : un message éventuel du solveur (aucun).
La commande permet donc de travailler avec une fonction de plusieurs variables, par exemple si l'on considère la fonction
    ƒ(x, y) = x$^2$ + y$^2$ - 10x - 8y + 41
on peut trouver son minimum avec
ce qui donne
La solution exacte est x = 5 et y = 4.

### Fonction avec paramètre

Si la fonction possède des paramètres, on doit les passer en utilisant leur nom exact. Par exemple, si la fonction est
    ƒ$_{a,b,c}$(x, y) = (x - a)$^2$ + (y - b)$^2$ + c
alors le script suivant permet de trouver le minimum :
ce qui donne

## Méthode d'optimisation avec contrainte linéaire

On recherche le minimum de ƒ, mais avec des contraintes linéaires, de type :
    inégalités : u$_i$×x$_i$ ≥ c$_i$
Le logiciel R propose la commande constrOptim, qui fait appel au solveur optim mais en ajoutant une fonction de pénalité pour gérer les contraintes. La syntaxe de base est :
où
* x$_0$ est le vecteur de départ, qui doit être dans la zone de faisabilité, c'est-à-dire vérifier les inégalités de contrainte ;
* f est la fonction ;
* NULL indique que l'on ne définit pas la fonction gradient ;
* u et c sont des vecteurs représentant les contraintes.
Par exemple, si l'on veut rechercher le minimum de la fonction y = sin(x) sur l'intervalle [-π ; π] :
* la fonction f est sin
* les contraintes linéaires sont 1×x ≥ -π, -1×x ≥ -π donc on a en notation matricielle u×x =c avec $\mathbf {u} ={\begin{pmatrix}1\\-1\end{pmatrix}}$ et $\mathbf {c} ={\begin{pmatrix}-\pi \\-\pi \end{pmatrix}}$ ;
* on part de x$_0$ = 0.
La résolution provoque des messages d'avertissement, qui sont sans conséquence ici. L'abscisse de la solution est donnée par
qui est une approximation de -π/2 ≃ -1,570,8 ; la valeur de la fonction en ce point est donnée par
Si l'on affiche la totalité de la solution, on obtient
qui nous indique, en plus des données d'optim :
* $outer.iterations : le nombre d'itération « extérieures », c'est-à-dire ici d'appels au solveur optim ;
* $barrier.value : c'est, lors de la dernière itération, la valeur du paramètre de pénalisation qui permet d'éviter de sortir des limites imposées.

## Les différents algorithmes d'optimisation

Jusqu'ici, nous avons utilisé les méthodes d'optimisation par défaut. Selon les problèmes que l'on aborde, on peut avoir intérêt à utiliser d'autres algorithmes.

### Fonction optimize

La fonction optimize utilise la méthode du nombre d'or associée à une interpolation parabolique.

### Fonctions optim et constrOptim
#### Méthode de Nelder-Mead

Par défaut, la fonction optim utilise la méthode de Nelder-Mead. Il consiste à prendre un simplexe — deux points sur une courbe, trois points sur une surface, … — et d'étendre ce simplexe vers les directions où la fonction est plus petite, puis à réduire ce simplexe lorsque le minimum est à l'intérieur.
Cet algorithme fonctionne lorsque la fonction n'est pas dérivable (mais elle doit être continue). Par contre, il n'est pas adapté aux problèmes à une dimension (fonction d'une seule variable), il vaut alors mieux utiliser la méthode de Brent ou bien la fonction optimize.

#### Méthode BFGS

Le deuxième algorithme proposé est celui de BFGS, avec l'option method = "BFGS". C'est un algorithme à direction de descente : à partir d'une point de départ, on se déplace sur le plan tangent à la surface, vers le bas. La direction de descente est l'opposé du gradient de la fonction.
Le pas de descente est fonction l'inverse de la courbure : si la courbure est importante, alors la surface s'écarte vite du plan tangent, il ne faut donc pas aller trop loin. L'algorithme calcule donc une approximation de l'inverse de cette dérivée seconde, ou plutôt, dans le cas d'une fonction de plusieurs variables, de la matrice hessienne (contenant les dérivées partielles d'ordre 2).
On peut fournir la fonction gradient, avec le troisième paramètre gr = … ; sinon, le solveur utilise une approximation par différences finies. Par exemple
ou bien
L'option hessian = TRUE retourne également la matrice hessienne calculée à la solution. On peut aussi calculer la matrice hessienne par
où solution$par est la solution trouvée par optimisation.
Le solveur propose une méthode BFGS modifiée, method = "L-BFGS-B", qui permet de définir des limites inférieures et supérieures aux variables.

#### Méthode du gradient conjugué

Le troisième algorithme proposé est celui du gradient conjugué, avec l'option method = "CG". Cette méthode ne calcule pas, et donc ne stocke pas, les matrices hessiennes ; elle permet donc de traiter des problèmes plus volumineux que la méthode BFGS, mais est moins robuste.
Le solveur accepte l'option hessian = TRUE.

#### Méthode du recuit simulé

Le quatrième algorithme proposé est celui du recuit simulé (simulated annealing), avec l'option method = "SANN". Le but de cette méthode est d'éviter de trouver un minimum local. Pour cela, une fois qu'une solution est trouvée, le solveur repart d'un autre point pour voir s'il converge vers une autre valeur. Cette méthode est utile lorsqu'il y a de nombreux minima locaux (surface « rugueuse »).

#### Méthode de Brent

Le quatrième algorithme proposé est la méthode de Brent, avec l'option method = "Brent". Ce solveur ne traite que des problèmes à une dimension, et est en fait similaire à la fonction optimize. En particulier, il faut définir l'intervalle de recherche avec les options lower = … et upper = …, mais il faut par contre toujours donner un point de départ.

#### Méthode de la barrière logarithmique

Pour la gestion des contrainte, la fonction constrOptim utilise la méthode de la barrière logarithmique. Cette méthode consiste à ajouter à la fonction ƒ une fonction de pénalité qui augmente lorsque l'on s'approche de la frontière. Si l'inégalité de contrainte s'écrit sous la forme
    h(x) ≥ 0
alors R prend pour fonction de pénalité
    μ×ln(h(x))
où μ est une constante positive. Cette constante vaut par défaut 10<spu>-4, mais peut être ajustée :
Une méthode consiste à faire des passes successives en diminuant au fur et à mesure la valeur de μ pour affiner le résultat, la fonction optimisée étant de plus en plus proche de la fonction initiale.

### Modèle linéaire généralisé

La régression des moindres carrés consiste à minimiser l'estimateur de la variance ; cette démarche est pertinente si le bruit suit une distribution normale ayant le même écart type. La commande glm permet d'utiliser d'autres distributions, avec le paramètre family.

## Régression non-linéaire

Notons tout d'abord qu'une régression linéaire est un modèle dans lequel on recherche des coefficients (β$_0$, β$_1$, …, β$_m$) ; c'est l'expression de ces coefficients qui doit être linéaire, du type :
    y = β$_0$ + β$_1$⋅ƒ$_1$(x) + β$_2$⋅ƒ$_2$(x) + … + β$_m$⋅ƒ$_m$(x)
où x est un vecteur (x$_1$, x$_2$, …, x$_n$). Les fonctions ƒ$_i$ ne sont pas nécessairement linéaires.
Un modèle non-linéaire est un modèle dans lequel les coefficients β$_i$ figurent eux-même de manière non-linéaire dans une formule. Par exemple, si l'on recherche un modèle gaussien :
    $y=\beta _{0}\cdot \exp \left(-\beta _{1}\cdot x^{2}\right)$
ce modèle n'est pas linéaire en β$_1$ — notons qu'on peut le rendre linéaire en considérant le logarithme de y.
La régression non-linéaire par la méthode des moindres carrés se fait avec la commande nls (nonlinear least squares). Par défaut, cette commande utilise la méthode de Gauss-Newton ; dans tous les cas, il s'agit d'algorithmes itératifs. La syntaxe de base est similaire à celle de lm, mais il faut introduire explicitement les paramètres à affiner dans la formule. Par exemple :
donne pour résultat :
La valeur de départ de l'algorithme itératif est choisie automatiquement. On peut lui indiquer des valeurs de départ avec le paramètre start, qui doit recevoir une liste nommée. par exemple, pour commencer avec β$_0$ = β$_1$ = 1 :

## Annexe : Formules symboliques

Le modèle de régression est une formule symbolique. Par exemple, si l'on écrit y ~ x, on utilise un modèle affine (droite avec ordonnée à l'origine) du type
    y = a + bx.
Pour une régression linéaire avec deux variables explicatives, on peut écrire y ~ x1 + x2 ; le signe plus n'indique pas ici l'addition, mais sépare les deux paramètres, le modèle sera donc
    y = a + b$_1$x$_1$ + b$_1$x$_2$.
Attention : les variables du modèle — x, x1, x2, y — doivent être les noms des variables utilisées dans la structure de données.
Si l'on veut un modèle linéaire sans ordonnée à l'origine, on modifie la formule en ajoutant un paramètre +0, par exemple y ~ x + 0.
Dans une formule symbolique, le signe plus indique que l'on combine linéairement deux variables. Si l'on veut appliquer une addition stricte, c'est-à-dire additionner deux termes sans les assortir d'un coefficient à définir par régression, il faut « inhiber l'interpréteur » avec la fonction I, par exemple :
indique un modèle
    x$^2$ + y$^2$ = a + b$_1$x + b$_2$y.
Formules symboliques dans R :
| Syntaxe | Modèle | Commentaires |
| Y ~ X | Y = β$_0$ + β$_1$X | |
| Y ~ 0 + X | Y = β$_1$X | |
| Y ~ X + I(X^2) | Y = β$_0$ + β$_1$X + β$_2$X$^2$ | |
| Y ~ X1 + X2 | Y = β$_0$ + β$_1$X$_1$ + β$_2$X$_2$ | modèle du premier ordre sans interaction croisée |
| Y ~ . | Y = β$_0$ + β$_1$X$_1$ + … | le point remplace toutes les colonnes des données (data.frame) non encore mentionnées dans la formule |
| Y ~ X1:X2 | Y = β$_0$ + β$_1$X$_1$X$_2$ | modèle ne contenant que des interactions croisée du premier ordre |
| Y ~ X1*X2 | Y = β$_0$ + β$_1$X$_1$ + β$_2$X$_2$ + β$_3$X$_1$X$_2$ | modèle du premier ordre complet équivalent de Y ~ X1 + X2 + X1:X2 |

## Notes et références

* ↑ le minimum est atteint à ƒ' = 0 soit 10x + 3 = 0 soit x = -3/10
* ↑ la fonction peut aussi s'écrire (x - 5)$^2$ + (y - 4)$^2$ - 25
* ↑ Kristel Van Steen, Using R for Linear Regression, lire en ligne